#include <bitset>
#include <iostream>
#include <booluarray.hpp>


#include <cstddef>
#include <iostream>
#include <iomanip>
#include <fstream>
#include <vector>
#include "bernstein.hpp"
#include "quadrature_multipoly.hpp"

#include "real.hpp"
#include "uvector.hpp"
#include "vector"
#include "xarray.hpp"
#include<chrono>


using namespace algoim;

// Driver method which takes a functor phi defining a single polynomial in the reference
// rectangle [xmin, xmax]^N, of Bernstein degree P, along with an integrand function,
// and performances a q-refinement convergence study, comparing the computed integral
// with the given exact answers, for 1 <= q <= qMax.
template <int N, typename Phi, typename F>
void qConv1(const Phi&      phi,
            real            xmin,
            real            xmax,
            uvector<int, N> P,
            const F&        integrand,
            int             qMax,
            real            volume_exact,
            real            surf_exact)
{
    // Construct Bernstein polynomial by mapping [0,1] onto bounding box [xmin,xmax]
    xarray<real, N> phipoly(nullptr, P);
    algoim_spark_alloc(real, phipoly);
    bernstein::bernsteinInterpolate<N>([&](const uvector<real, N>& x) { return phi(xmin + x * (xmax - xmin)); }, phipoly);

    // Build quadrature hierarchy
    ImplicitPolyQuadrature<N> ipquad(phipoly);

    // Functional to evaluate volume and surface integrals of given integrand
    real volume, surf;
    auto compute = [&](int q) {
        volume = 0.0;
        surf   = 0.0;
        // compute volume integral over {phi < 0} using AutoMixed strategy
        ipquad.integrate(AutoMixed, q, [&](const uvector<real, N>& x, real w) {
            if (bernstein::evalBernsteinPoly(phipoly, x) < 0) volume += w * integrand(xmin + x * (xmax - xmin));
        });
        // compute surface integral over {phi == 0} using AutoMixed strategy
        ipquad.integrate_surf(AutoMixed, q, [&](const uvector<real, N>& x, real w, const uvector<real, N>& wn) {
            surf += w * integrand(xmin + x * (xmax - xmin));
        });
        // scale appropriately
        volume *= pow(xmax - xmin, N);
        surf   *= pow(xmax - xmin, N - 1);
    };

    // Compute results for all q and output in a convergence table
    for (int q = 1; q <= qMax; ++q) {
        compute(q);
        std::cout << q << ' ' << volume << ' ' << surf << ' ' << std::abs(volume - volume_exact) / volume_exact << ' '
                  << std::abs(surf - surf_exact) / surf_exact << std::endl;
    }
}

// Driver method which takes two phi functors defining two polynomials in the reference
// rectangle [xmin, xmax]^N, each of of Bernstein degree P, builds a quadrature scheme with the
// given q, and outputs it for visualisation in a set of VTP XML files
template <int N, typename F1, typename F2, typename F>
void qConvMultiPloy(const F1&              fphi1,
                    const F2&              fphi2,
                    real                   xmin,
                    real                   xmax,
                    const uvector<int, N>& P,
                    const F&               integrand,
                    int                    q,
                    std::string            qfile)
{
    // Construct phi by mapping [0,1] onto bounding box [xmin,xmax]
    xarray<real, N> phi1(nullptr, P), phi2(nullptr, P);
    algoim_spark_alloc(real, phi1, phi2);
    bernstein::bernsteinInterpolate<N>([&](const uvector<real, N>& x) { return fphi1(xmin + x * (xmax - xmin)); }, phi1);
    // bernstein::bernsteinInterpolate<N>([&](const uvector<real, N>& x) { return fphi2(xmin + x * (xmax - xmin)); }, phi2);

    // Build quadrature hierarchy
    ImplicitPolyQuadrature<N> ipquad(phi1, phi2);
    // ImplicitPolyQuadrature<N> ipquad(phi1);

    // Functional to evaluate volume and surface integrals of given integrand
    real volume, surf;
    auto compute = [&](int q) {
        volume = 0.0;
        surf   = 0.0;
        // compute volume integral over {phi < 0} using AutoMixed strategy
        ipquad.integrate(AutoMixed, q, [&](const uvector<real, N>& x, real w) {
            if (bernstein::evalBernsteinPoly(phi1, x) < 0 && bernstein::evalBernsteinPoly(phi2, x) < 0) volume += w * integrand(xmin + x * (xmax - xmin));
            // if (bernstein::evalBernsteinPoly(phi1, x) < 0) volume += w * integrand(xmin + x * (xmax - xmin));
        });
        // compute surface integral over {phi == 0} using AutoMixed strategy
        ipquad.integrate_surf(AutoMixed, q, [&](const uvector<real, N>& x, real w, const uvector<real, N>& wn) {
            surf += w * integrand(xmin + x * (xmax - xmin));
        });
        // scale appropriately
        volume *= pow(xmax - xmin, N);
        surf   *= pow(xmax - xmin, N - 1);
    };
    compute(q);
    std::cout << "q volume: " << volume << std::endl;
}

template <int N, typename F>
void qConvMultiPloy2(
                    real                   xmin,
                    real                   xmax,
                    const uvector<int, N>& P,
                    const F&               integrand,
                    int                    q,
                    std::string            qfile)
{
    // Construct phi by mapping [0,1] onto bounding box [xmin,xmax]

    // bernstein::bernsteinInterpolate<N>([&](const uvector<real, N>& x) { return fphi1(xmin + x * (xmax - xmin)); }, phi1);
    // bernstein::bernsteinInterpolate<N>([&](const uvector<real, N>& x) { return fphi2(xmin + x * (xmax - xmin)); }, phi2);

		std::vector<xarray<real, N>> phis; // 大量sphere
		for (int i = 0; i < 100; i++) {
			real centerX =  1.0 - rand() % 20 / 10.0, centerY = 1.0 - rand() % 20 / 10.0, centerZ = 1.0 - rand() % 20 / 10.0;	
			real r = 0.1 + rand() % 9 / 10.0;
			auto fphi = [&](const uvector<real, 3>& xx) {
				real x = xx(0) - centerX;
				real y = xx(1) - centerY;
				real z = xx(2) - centerZ;
				return x * x + y * y + z * z - 1;
    	};
			xarray<real, N> phi(nullptr, P);
			algoim_spark_alloc(real, phi);
			bernstein::bernsteinInterpolate<N>([&](const uvector<real, N>& x) { return fphi(xmin + x * (xmax - xmin)); }, phi);
		}
		auto t = std::chrono::system_clock::now();

		ImplicitPolyQuadrature<N> ipquad(phis);



    // Build quadrature hierarchy
    // ImplicitPolyQuadrature<N> ipquad(phi1, phi2);
    // ImplicitPolyQuadrature<N> ipquad(phi1);

    // Functional to evaluate volume and surface integrals of given integrand
    real volume, surf;
    auto compute = [&](int q) {
        volume = 0.0;
        surf   = 0.0;
        // compute volume integral over {phi < 0} using AutoMixed strategy
        ipquad.integrate(AutoMixed, q, [&](const uvector<real, N>& x, real w) {
            // if (bernstein::evalBernsteinPoly(phi1, x) < 0 && bernstein::evalBernsteinPoly(phi2, x) < 0) volume += w * integrand(xmin + x * (xmax - xmin));
            // if (bernstein::evalBernsteinPoly(phi1, x) < 0) volume += w * integrand(xmin + x * (xmax - xmin));
						volume += w * integrand(xmin + x * (xmax - xmin));
        });
        // compute surface integral over {phi == 0} using AutoMixed strategy
        // ipquad.integrate_surf(AutoMixed, q, [&](const uvector<real, N>& x, real w, const uvector<real, N>& wn) {
        //     surf += w * integrand(xmin + x * (xmax - xmin));
        // });
        // scale appropriately
        volume *= pow(xmax - xmin, N);
        surf   *= pow(xmax - xmin, N - 1);
    };
    compute(q);

		auto tAfter = std::chrono::system_clock::now();
		std::chrono::duration<double> elapsed_seconds = tAfter - t;
		std::cout << "elapsed time: " << elapsed_seconds.count() << "s\n";
    std::cout << "q volume: " << volume << std::endl;
}

void fromCommon2Bernstein() {

}

void testMultiPolys()
{
    // Visusalisation of a 2D implicitly-defined domain involving the intersection of two polynomials; this example
    // corresponds to the top-left example of Figure 15, https://doi.org/10.1016/j.jcp.2021.110720
    if (true) {
        auto phi0 = [](const uvector<real, 3>& xx) {
            real x = xx(0);
            real y = xx(1);
            real z = xx(2);
            return x * x + y * y + z * z - 1;
        };
        auto phi1 = [](const uvector<real, 3>& xx) {
            // real x = xx(0);
            // real y = xx(1);
            // real z = xx(2);
            // return x * x + y * y + z * z - 1;
            real x = xx(0);
            real y = xx(1);
            real z = xx(2);
            return x * y * z;
        };
        // auto phi0 = [](const uvector<real, 3>& xx) {
        //     real x = xx(0) * 2 - 1;
        //     real y = xx(1) * 2 - 1;
        //     return 0.014836540349115947 + 0.7022484024095262 * y + 0.09974561176434385 * y * y
        //            + x * (0.6863910464417281 + 0.03805619999999999 * y - 0.09440658332756446 * y * y)
        //            + x * x * (0.19266932968830816 - 0.2325190091204104 * y + 0.2957473125000001 * y * y);
        // };
        // auto phi1 = [](const uvector<real, 3>& xx) {
        //     real x = xx(0) * 2 - 1;
        //     real y = xx(1) * 2 - 1;
        //     return -0.18792528379702625 + 0.6713882473904913 * y + 0.3778666084723582 * y * y
        //            + x * x * (-0.14480813208127946 + 0.0897755603159206 * y - 0.141199875 * y * y)
        //            + x * (-0.6169311810674598 - 0.19449299999999994 * y - 0.005459163675646665 * y * y);
        // };
        auto integrand = [](const uvector<real, 3>& x) { return 1.0; };
        // qConvMultiPloy<3>(phi0, phi1, -1.0, 1.0, 3, integrand, 3, "exampleC");
        qConvMultiPloy2<3>(-1.0, 1.0, 3, integrand, 20, "exampleC");
        std::cout << "\n\nQuadrature visualisation of a 2D implicitly-defined domain involving the\n";
        std::cout << "intersection of two polynomials, corresponding to the top-left example of Figure 15,\n";
        std::cout << "https://doi.org/10.1016/j.jcp.2021.110720, written to exampleC-xxxx.vtp files\n";
        std::cout << "(XML VTP file format).\n";
    }
}

void testBitSet()
{
    std::bitset<10> set(128);
    set.set();
    std::cout << set << std::endl;
}

void testBooluarray()
{
    algoim::booluarray<2, 3> tmp;
    tmp(0) = true;
    tmp(1) = true;
    tmp(2) = true;
    std::cout << tmp.bits << std::endl;
}

void testMain()
{
    testBooluarray();
    testMultiPolys();
}