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


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

#include "real.hpp"
#include "sparkstack.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;
    }
}

enum BoolOp { Union, Intersection, Difference };

// 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,
                    BoolOp                 op)
{
    // 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
    auto                      t = std::chrono::system_clock::now();
    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 (op == Union) {
                if (bernstein::evalBernsteinPoly(phi1, x) < 0 || bernstein::evalBernsteinPoly(phi2, x) < 0)
                    volume += w * integrand(xmin + x * (xmax - xmin));
            } else if (op == Intersection) {
                if (bernstein::evalBernsteinPoly(phi1, x) < 0 && bernstein::evalBernsteinPoly(phi2, x) < 0)
                    volume += w * integrand(xmin + x * (xmax - xmin));
            } else if (op == Difference) {
                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);
    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;
}

template <int N, typename F>
void qConvMultiPloy2(real xmin, real xmax, const std::vector<xarray<real, N>>& phis, const F& integrand, int q)
{
    auto t = std::chrono::system_clock::now();

    ImplicitPolyQuadrature<N> ipquad(phis);

    // 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));
            bool in = true;
            for (int i = 0; i < phis.size(); ++i) {
                if (bernstein::evalBernsteinPoly(phis[i], x) >= 0) {
                    in = false;
                    break;
                }
            }
            if (in) 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 << "volume: " << volume << std::endl;
}

void fromCommon2Bernstein() {}

void testMultiPolys()
{
    const int                    N    = 3;
    const uvector<int, N>        P    = 3;
    real                         xmin = -1;
    real                         xmax = 1;
    std::vector<xarray<real, N>> phis; // 大量sphere
    // xarray<real, N>* phis = new xarray<real, N>[100];
    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.01 + rand() % 10 / 100.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 - r * r;
        };
        xarray<real, 3> phi(nullptr, P);
        algoim_spark_alloc(real, phi);
        bernstein::bernsteinInterpolate<N>([&](const uvector<real, N>& x) { return fphi(xmin + x * (xmax - xmin)); }, phi);
        phis.emplace_back(phi);
    }
    auto integrand = [](const uvector<real, 3>& x) { return 1.0; };
    // qConvMultiPloy<3>(phi0, phi1, -1.0, 1.0, 3, integrand, 3, "exampleC");
    std::cout << "testMultiPolys" << std::endl;
    qConvMultiPloy2<3>(-1.0, 1.0, phis, integrand, 10);
}

template <int N>
void power2BernsteinTensorProduct(const xarray<real, N>& phiPower, xarray<real, N>& phiBernstein)
{
    xarrayInit(phiBernstein);
    for (auto i = phiPower.loop(); ~i; ++i) {
        // phi.l(i) = powerFactors.l(i);
        real factorBase = phiPower.l(i);
        if (factorBase == 0) continue;
        auto                           traverseRange = phiPower.ext() - i();
        std::vector<std::vector<real>> decompFactors(N, std::vector<real>(max(traverseRange), 0.));
        for (int dim = 0; dim < N; ++dim) {
            // Sigma
            size_t      nDim      = phiPower.ext()(dim) - 1;
            const real* binomNDim = Binomial::row(nDim);
            for (int j = i(dim); j <= nDim; ++j) {
                const real* binomJ             = Binomial::row(j);
                decompFactors[dim][j - i(dim)] = binomJ[i(dim)] / binomNDim[i(dim)];
            }
        }
        xarray<real, N> subgrid(nullptr, traverseRange);
        algoim_spark_alloc(real, subgrid);

        for (auto ii = subgrid.loop(); ~ii; ++ii) {
            real factor = factorBase;
            for (int dim = 0; dim < N; ++dim) { factor *= decompFactors[dim][ii(dim)]; }
            phiBernstein.m(i() + ii()) += factor;
        }
    }
}

template <int N>
real powerEvaluation(const xarray<real, N>& phi, const uvector<real, N>& x)
{
    real res = 0;
    for (auto i = phi.loop(); ~i; ++i) {
        real item = phi.l(i);
        auto exps = i();
        for (int dim = 0; dim < N; ++dim) { item *= pow(x(dim), i(dim)); }
        res += item;
    }
    return res;
}

template <int N, typename F>
void qConvBernstein(const xarray<real, N>&                  phi,
                    real                                    xmin,
                    real                                    xmax,
                    const F&                                integrand,
                    int                                     q,
                    const std::vector<std::array<real, 4>>& halfFaces)
{
    auto t = std::chrono::system_clock::now();

    uvector<real, 3> testX(0., 0.75, 0.);
    real             testEvalBernstein = bernstein::evalBernsteinPoly(phi, testX);
    auto             vec1              = xarray2StdVector(phi);
    std::cout << "eval bernstein without interpolation:" << testEvalBernstein << std::endl;

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

    // Functional to evaluate volume and surface integrals of given integrand
    real volume;
    auto compute = [&](int q) {
        volume = 0.0;
        // compute volume integral over {phi < 0} using AutoMixed strategy
        if (!halfFaces.empty()) {
            ipquad.integrate(AutoMixed, q, [&](const uvector<real, N>& x, real w) {
                if (bernstein::evalBernsteinPoly(phi, x) >= 0) return;
                uvector<real, N> trueX = xmin + x * (xmax - xmin);
                // uvector<real, N> trueX = x;
                // std::cout << "trueX: " << trueX << std::endl;
                bool             in    = true;
                for (int i = 0; i < halfFaces.size(); ++i) {
                    uvector<real, 3> factors(halfFaces[i][0], halfFaces[i][1], halfFaces[i][2]);
                    if (sum(factors * trueX) + halfFaces[i][3] >= 0) {
                        in = false;
                        break;
                    }
                }
                if (in) volume += w * integrand(xmin + x * (xmax - xmin));
            });

        } else
            ipquad.integrate(AutoMixed, q, [&](const uvector<real, N>& x, real w) {
                if (bernstein::evalBernsteinPoly(phi, x) < 0) volume += w * integrand(xmin + x * (xmax - xmin));
            });
        // scale appropriately
        volume *= pow(xmax - xmin, N);
    };
    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;
}

// template <typename T, int N>
// void xarrayInit(xarray<T, N>& arr)
// {
//     for (auto i = arr.loop(); ~i; ++i) { arr.l(i) = 0; }
// }

void testSpherePowerDirectly()
{
    // a_x(x-c_x)^2 + a_y(y-c_y)^2 + a_z(x-c_z)^2 - r^2
    uvector<real, 3> xmin = -1, xmax = 1;
    uvector<real, 3> range = xmax - xmin;
    assert(all(range != 0));
    // uvector<real, 3> k    = 1 / range;
    // uvector<real, 3> bias = -k * xmin;
    uvector<real, 3> k    = range;
    uvector<real, 3> bias = xmin;
    real             a[3] = {1, 1, 1};
    real             c[3] = {0, 0, 0};
    real             r    = 1;
    uvector<int, 3>  ext  = 3;
    xarray<real, 3>  phiPower(nullptr, ext), phiBernstein(nullptr, ext);
    algoim_spark_alloc(real, phiPower, phiBernstein);
    xarrayInit(phiPower);
    auto v  = xarray2StdVector(phiPower);
    auto vv = xarray2StdVector(phiBernstein);
    for (int dim = 0; dim < 3; ++dim) {
        uvector<int, 3> idx  = 0;
        phiPower.m(idx)     += a[dim] * (bias(dim) - c[dim]) * (bias(dim) - c[dim]);
        idx(dim)             = 1;
        phiPower.m(idx)      = 2 * a[dim] * k(dim) * (bias(dim) - c[dim]);
        idx(dim)             = 2;
        phiPower.m(idx)      = a[dim] * k(dim) * k(dim);
    }
    phiPower.m(0)  -= r * r;
    auto integrand  = [](const uvector<real, 3>&) { return 1.0; };
    power2BernsteinTensorProduct(phiPower, phiBernstein);

    v = xarray2StdVector(phiPower);
    uvector<real, 3> testX(0., 0., 0.5);
    real             testEval = powerEvaluation(phiPower, testX);
    std::cout << "eval power:" << testEval << std::endl;

    real testEvalBernstein = bernstein::evalBernsteinPoly(phiBernstein, testX);
    auto vec               = xarray2StdVector(phiBernstein);
    std::cout << "eval bernstein:" << testEvalBernstein << std::endl;
    qConvBernstein(phiBernstein, -1, 1, integrand, 10, {});
}

void powerTransformation(const uvector<real, 3>& scale, const uvector<real, 3>& bias, xarray<real, 3>& phiPower)
{
    std::vector<std::vector<std::vector<real>>> dimOrderExpansion;
    const auto&                                 ext = phiPower.ext();
    for (int dim = 0; dim < 3; ++dim) {
        dimOrderExpansion.push_back(std::vector<std::vector<real>>(ext(dim)));
        for (int degree = 0; degree < ext(dim); ++degree) {
            const real* binomDegree = Binomial::row(degree);
            dimOrderExpansion[dim][degree].reserve(degree + 1);
            // 根据二项定理展开
            for (int i = 0; i <= degree; ++i) {
                dimOrderExpansion[dim][degree].push_back(binomDegree[i] * pow(scale(dim), i) * pow(bias(dim), degree - i));
            }
        }
    }
    for (auto i = phiPower.loop(); ~i; ++i) {
        // 迭代器必须按照坐标的升序进行访问，即，访问ijk时，(i-1)jk，i(j-1)k，ij(k-1)必须已经访问过
        real base     = phiPower.l(i);
        phiPower.l(i) = 0;
        auto exps     = i();
        for (MultiLoop<3> j(0, exps + 1); ~j; ++j) {
            real item = base;
            for (int dim = 0; dim < 3; ++dim) { item *= dimOrderExpansion[dim][exps(dim)][j(dim)]; }
            phiPower.m(j()) += item;
        }
    }
}

void testSpherePowerDirectlyByTransformation()
{
    // a_x(x-c_x)^2 + a_y(y-c_y)^2 + a_z(x-c_z)^2 - r^2
    uvector<real, 3> xmin = -1, xmax = 1;
    uvector<real, 3> range = xmax - xmin;
    assert(all(range != 0));
    // uvector<real, 3> k    = 1 / range;
    // uvector<real, 3> bias = -k * xmin;
    real            a[3] = {1, 1, 1};
    real            c[3] = {0, 0, 0};
    real            r    = 0.2;
    uvector<int, 3> ext  = 3;
    xarray<real, 3> phiPower(nullptr, ext), phiBernstein(nullptr, ext);
    algoim_spark_alloc(real, phiPower, phiBernstein);
    xarrayInit(phiPower);
    auto v  = xarray2StdVector(phiPower);
    auto vv = xarray2StdVector(phiBernstein);
    for (int dim = 0; dim < 3; ++dim) {
        uvector<int, 3> idx  = 0;
        phiPower.m(idx)     += a[dim] * (c[dim]) * (c[dim]);
        idx(dim)             = 1;
        phiPower.m(idx)      = 2 * a[dim] * (-c[dim]);
        idx(dim)             = 2;
        phiPower.m(idx)      = a[dim];
    }
    phiPower.m(0) -= r * r;
    powerTransformation(range, xmin, phiPower);

    auto integrand = [](const uvector<real, 3>& x) { return 1.0; };
    power2BernsteinTensorProduct(phiPower, phiBernstein);

    v = xarray2StdVector(phiPower);
    uvector<real, 3> testX(0., 0., 0.5);
    real             testEval = powerEvaluation(phiPower, testX);
    std::cout << "eval power:" << testEval << std::endl;

    real testEvalBernstein = bernstein::evalBernsteinPoly(phiBernstein, testX);
    auto vec               = xarray2StdVector(phiBernstein);
    std::cout << "eval bernstein:" << testEvalBernstein << std::endl;
    qConvBernstein(phiBernstein, -1, 1, integrand, 10, {});
}

void testCylinderPowerDirectly()
{
    // a_x(x-c_x)^2 + a_y(y-c_y)^2 + a_z(x-c_z)^2 - r^2
    uvector<real, 3> xmin = -1, xmax = 1;
    uvector<real, 3> range = xmax - xmin;
    assert(all(range != 0));
    uvector<real, 3> k    = range;
    uvector<real, 3> bias = xmin;
    real             c[3] = {0, 0, 0};
    real             r = 1, h = 0.8;
    uvector<int, 3>  ext = 3;
    xarray<real, 3>  phiPower(nullptr, ext), phiBernstein(nullptr, ext);
    algoim_spark_alloc(real, phiPower, phiBernstein);
    xarrayInit(phiPower);
    auto v   = xarray2StdVector(phiPower);
    auto vv  = xarray2StdVector(phiBernstein);
    real top = c[2] + h * 0.5, bottom = c[2] - h * 0.5;
    phiPower.m(uvector<int, 3>(2, 0, 2)) = -1;
    phiPower.m(uvector<int, 3>(0, 2, 2)) = -1;
    phiPower.m(uvector<int, 3>(2, 0, 1)) = top + bottom;
    phiPower.m(uvector<int, 3>(0, 2, 1)) = top + bottom;
    phiPower.m(uvector<int, 3>(2, 0, 0)) = -bottom * top;
    phiPower.m(uvector<int, 3>(0, 2, 0)) = -bottom * top;
    phiPower.m(uvector<int, 3>(1, 0, 2)) = 2 * c[0];
    phiPower.m(uvector<int, 3>(0, 1, 2)) = 2 * c[1];
    phiPower.m(uvector<int, 3>(1, 0, 1)) = -2 * c[0] * (top + bottom);
    phiPower.m(uvector<int, 3>(0, 1, 1)) = -2 * c[1] * (top + bottom);
    phiPower.m(uvector<int, 3>(1, 0, 0)) = 2 * c[0] * (top * bottom);
    phiPower.m(uvector<int, 3>(0, 1, 0)) = 2 * c[0] * (top * bottom);
    phiPower.m(uvector<int, 3>(0, 0, 2)) = -(c[0] * c[0] + c[1] * c[1] - r * r);
    phiPower.m(uvector<int, 3>(0, 0, 1)) = (bottom + top) * (c[0] * c[0] + c[1] * c[1] - r * r);
    phiPower.m(uvector<int, 3>(0, 0, 0)) = -bottom * top * (c[0] * c[0] + c[1] * c[1] - r * r);

    powerTransformation(range, xmin, phiPower);

    auto integrand = [](const uvector<real, 3>& x) { return 1.0; };
    power2BernsteinTensorProduct(phiPower, phiBernstein);

    v = xarray2StdVector(phiPower);
    uvector<real, 3> testX(0., 0., 0.75);
    real             testEval = powerEvaluation(phiPower, testX);
    std::cout << "eval power:" << testEval << std::endl;

    real testEvalBernstein = bernstein::evalBernsteinPoly(phiBernstein, testX);
    auto vec               = xarray2StdVector(phiBernstein);
    std::cout << "eval bernstein:" << testEvalBernstein << std::endl;
    std::cout << "testCylinderPowerDirectly" << std::endl;
    qConvBernstein(phiBernstein,
                   -1,
                   1,
                   integrand,
                   10,
                   {
                       {0, 0, 1,  -top  },
                       {0, 0, -1, bottom}
    });
}

void testCylinderAlgoim()
{
    const int                    N    = 3;
    const uvector<int, N>        P    = 3;
    real                         xmin = -1, xmax = 1;
    std::vector<xarray<real, N>> phis(3, xarray<real, N>(nullptr, P));
    auto                         fphiCircle = [](const uvector<real, 3>& xx) {
        real x = xx(0);
        real y = xx(1);
        real z = xx(2);
        return x * x + y * y - 1;
    };
    auto fphiUp = [](const uvector<real, 3>& xx) {
        real z = xx(2);
        return z - 0.5;
    };
    auto fphiFloor = [](const uvector<real, 3>& xx) {
        real z = xx(2);
        return -z - 0.5;
    };

    algoim_spark_alloc(real, phis[0]);
    algoim_spark_alloc(real, phis[1]);
    algoim_spark_alloc(real, phis[2]);
    bernstein::bernsteinInterpolate<N>([&](const uvector<real, N>& x) { return fphiCircle(xmin + x * (xmax - xmin)); },
                                       phis[0]);
    bernstein::bernsteinInterpolate<N>([&](const uvector<real, N>& x) { return fphiUp(xmin + x * (xmax - xmin)); }, phis[1]);
    bernstein::bernsteinInterpolate<N>([&](const uvector<real, N>& x) { return fphiFloor(xmin + x * (xmax - xmin)); }, phis[2]);
    auto integrand = [](const uvector<real, 3>& x) { return 1.0; };
    // // qConvMultiPloy<3>(phi0, phi1, -1.0, 1.0, 3, integrand, 3, "exampleC");
    std::cout << "testCylinderAlgoim" << std::endl;
    qConvMultiPloy2<3>(-1.0, 1.0, phis, integrand, 10);
}

void test6PlanesAlgoim()
{
    const int                    N    = 3;
    const uvector<int, N>        P    = 3;
    real                         xmin = -1, xmax = 1;
    std::vector<xarray<real, N>> phis(6, xarray<real, N>(nullptr, P));
    // auto                         fphiPlane1 = [](const uvector<real, 3>& xx) {
    //     real x = xx(0);
    //     real y = xx(1);
    //     real z = xx(2);
    //     return x + y + z - 1;
    // };
    // auto fphiPlane2 = [](const uvector<real, 3>& xx) {
    //     real x = xx(0);
    //     real y = xx(1);
    //     real z = xx(2);
    //     return x - y + z - 1;
    // };
    // auto fphiPlane3 = [](const uvector<real, 3>& xx) {
    //     real x = xx(0);
    //     real y = xx(1);
    //     real z = xx(2);
    //     return -x + y + z - 1;
    // };
    auto                         fphiPlane1 = [](const uvector<real, 3>& xx) {
        real x = xx(0);
        return x - 0.5;
    };
    auto fphiPlane2 = [](const uvector<real, 3>& xx) {
        real x = xx(0);
        return -x - 0.5;
    };
    auto fphiPlane3 = [](const uvector<real, 3>& xx) {
        real y = xx(1);
        return y - 0.5;
    };
    auto fphiPlane4 = [](const uvector<real, 3>& xx) {
        real y = xx(1);
        return -y - 0.5;
    };
    auto fphiPlane5 = [](const uvector<real, 3>& xx) {
        real z = xx(2);
        return z - 0.5;
    };
    auto fphiPlane6 = [](const uvector<real, 3>& xx) {
        real z = xx(2);
        return -z - 0.5;
    };
    algoim_spark_alloc(real, phis[0]);
    algoim_spark_alloc(real, phis[1]);
    algoim_spark_alloc(real, phis[2]);
    algoim_spark_alloc(real, phis[3]);
    algoim_spark_alloc(real, phis[4]);
    algoim_spark_alloc(real, phis[5]);
    bernstein::bernsteinInterpolate<N>([&](const uvector<real, N>& x) { return fphiPlane1(xmin + x * (xmax - xmin)); },
                                       phis[0]);
    bernstein::bernsteinInterpolate<N>([&](const uvector<real, N>& x) { return fphiPlane2(xmin + x * (xmax - xmin)); },
                                       phis[1]);
    bernstein::bernsteinInterpolate<N>([&](const uvector<real, N>& x) { return fphiPlane3(xmin + x * (xmax - xmin)); },
                                       phis[2]);
    bernstein::bernsteinInterpolate<N>([&](const uvector<real, N>& x) { return fphiPlane4(xmin + x * (xmax - xmin)); },
                                       phis[3]);
    bernstein::bernsteinInterpolate<N>([&](const uvector<real, N>& x) { return fphiPlane5(xmin + x * (xmax - xmin)); },
                                       phis[4]);
    bernstein::bernsteinInterpolate<N>([&](const uvector<real, N>& x) { return fphiPlane6(xmin + x * (xmax - xmin)); },
                                       phis[5]);
    auto integrand = [](const uvector<real, 3>& x) { return 1.0; };
    // // qConvMultiPloy<3>(phi0, phi1, -1.0, 1.0, 3, integrand, 3, "exampleC");
    std::cout << "testPlanesAlgoim" << std::endl;
    qConvMultiPloy2<3>(-1.0, 1.0, phis, integrand, 10);
}

void test8PlanesAlgoim()
{
    const int                    N    = 3;
    const uvector<int, N>        P    = 3;
    real                         xmin = -1, xmax = 1;
    std::vector<xarray<real, N>> phis(8, xarray<real, N>(nullptr, P));
    auto                         fphiPlane1 = [](const uvector<real, 3>& xx) {
        real x = xx(0);
        real y = xx(1);
        real z = xx(2);
        return x + y + z - 1;
    };
    auto fphiPlane2 = [](const uvector<real, 3>& xx) {
        real x = xx(0);
        real y = xx(1);
        real z = xx(2);
        return x - y + z - 1;
    };
    auto fphiPlane3 = [](const uvector<real, 3>& xx) {
        real x = xx(0);
        real y = xx(1);
        real z = xx(2);
        return -x + y + z - 1;
    };
    auto fphiPlane4 = [](const uvector<real, 3>& xx) {
        real x = xx(0);
        real y = xx(1);
        real z = xx(2);
        return -x - y + z - 1;
    };
    auto fphiPlane5 = [](const uvector<real, 3>& xx) {
        real x = xx(0);
        real y = xx(1);
        real z = xx(2);
        return x + y - z - 1;
    };
    auto fphiPlane6 = [](const uvector<real, 3>& xx) {
        real x = xx(0);
        real y = xx(1);
        real z = xx(2);
        return x - y - z - 1;
    };
    auto fphiPlane7 = [](const uvector<real, 3>& xx) {
        real x = xx(0);
        real y = xx(1);
        real z = xx(2);
        return -x + y - z - 1;
    };
    auto fphiPlane8 = [](const uvector<real, 3>& xx) {
        real x = xx(0);
        real y = xx(1);
        real z = xx(2);
        return -x - y - z - 1;
    };
    algoim_spark_alloc(real, phis[0]);
    algoim_spark_alloc(real, phis[1]);
    algoim_spark_alloc(real, phis[2]);
    algoim_spark_alloc(real, phis[3]);
    algoim_spark_alloc(real, phis[4]);
    algoim_spark_alloc(real, phis[5]);
    algoim_spark_alloc(real, phis[6]);
    algoim_spark_alloc(real, phis[7]);
    bernstein::bernsteinInterpolate<N>([&](const uvector<real, N>& x) { return fphiPlane1(xmin + x * (xmax - xmin)); },
                                       phis[0]);
    bernstein::bernsteinInterpolate<N>([&](const uvector<real, N>& x) { return fphiPlane2(xmin + x * (xmax - xmin)); },
                                       phis[1]);
    bernstein::bernsteinInterpolate<N>([&](const uvector<real, N>& x) { return fphiPlane3(xmin + x * (xmax - xmin)); },
                                       phis[2]);
    bernstein::bernsteinInterpolate<N>([&](const uvector<real, N>& x) { return fphiPlane4(xmin + x * (xmax - xmin)); },
                                       phis[3]);
    bernstein::bernsteinInterpolate<N>([&](const uvector<real, N>& x) { return fphiPlane5(xmin + x * (xmax - xmin)); },
                                       phis[4]);
    bernstein::bernsteinInterpolate<N>([&](const uvector<real, N>& x) { return fphiPlane6(xmin + x * (xmax - xmin)); },
                                       phis[5]);
    bernstein::bernsteinInterpolate<N>([&](const uvector<real, N>& x) { return fphiPlane7(xmin + x * (xmax - xmin)); },
                                       phis[6]);
    bernstein::bernsteinInterpolate<N>([&](const uvector<real, N>& x) { return fphiPlane8(xmin + x * (xmax - xmin)); },
                                       phis[7]);

    auto integrand = [](const uvector<real, 3>& x) { return 1.0; };
    // // qConvMultiPloy<3>(phi0, phi1, -1.0, 1.0, 3, integrand, 3, "exampleC");
    std::cout << "testPlanesAlgoim" << std::endl;
    qConvMultiPloy2<3>(-1.0, 1.0, phis, integrand, 10);
}

void compositePower(const std::vector<xarray<real, 3>>& powers,
                    int                                 powerIdx,
                    uvector<int, 3>                     powerSum,
                    real                                factor,
                    xarray<real, 3>&                    res)

{
    if (powerIdx == 0) {
        {
            uvector3 ext(1, 1, 1);
            for (auto& t : powers) { ext += t.ext() - 1; }
            assert(all(ext == res.ext()));
        }
        xarrayInit(res);
    }
    if (powerIdx == powers.size()) {
        res.m(powerSum) += factor;
        return;
    }
    auto& power = powers[powerIdx];
    for (auto i = power.loop(); ~i; ++i) {
        if (power.l(i) == 0) {
            // factor = 0;
            continue;
        }
        compositePower(powers, powerIdx + 1, powerSum + i(), factor * power.l(i), res);
    }
    // int a = 1;
}

void test8PlanesPowerDirectly()
{
    // a_x(x-c_x)^2 + a_y(y-c_y)^2 + a_z(x-c_z)^2 - r^2
    uvector<real, 3> xmin = -1, xmax = 1;
    uvector<real, 3> range = xmax - xmin;
    assert(all(range != 0));
    uvector<real, 3>             k    = range;
    uvector<real, 3>             bias = xmin;
    uvector<int, 3>              ext  = 2;
    std::vector<xarray<real, 3>> phiPowers(8, xarray<real, 3>(nullptr, ext));
    algoim_spark_alloc(real,
                       phiPowers[0],
                       phiPowers[1],
                       phiPowers[2],
                       phiPowers[3],
                       phiPowers[4],
                       phiPowers[5],
                       phiPowers[6],
                       phiPowers[7]);
    for (int i = 0; i < 8; i++) {
        xarrayInit(phiPowers[i]);
        if (i & 1) {
            phiPowers[i].m(uvector<int, 3>(0, 0, 1)) = 1;
        } else {
            phiPowers[i].m(uvector<int, 3>(0, 0, 1)) = -1;
        }
        if (i & 2) {
            phiPowers[i].m(uvector<int, 3>(0, 1, 0)) = 1;
        } else {
            phiPowers[i].m(uvector<int, 3>(0, 1, 0)) = -1;
        }
        if (i & 4) {
            phiPowers[i].m(uvector<int, 3>(1, 0, 0)) = 1;
        } else {
            phiPowers[i].m(uvector<int, 3>(1, 0, 0)) = -1;
        }
        phiPowers[i].m(uvector<int, 3>(0, 0, 0)) = -1;
    }
    uvector<int, 3> resExt = 1 + phiPowers.size();
    xarray<real, 3> phiPowerAll(nullptr, resExt), phiBernstein(nullptr, resExt);
    algoim_spark_alloc(real, phiPowerAll, phiBernstein);
    compositePower(phiPowers, 0, 0, 1, phiPowerAll);

    powerTransformation(range, xmin, phiPowerAll);

    auto integrand = [](const uvector<real, 3>& x) { return 1.0; };
    power2BernsteinTensorProduct(phiPowerAll, phiBernstein);

    uvector<real, 3> testX(0., 0., 0.75);
    real             testEval = powerEvaluation(phiPowerAll, testX);
    std::cout << "eval power:" << testEval << std::endl;

    real testEvalBernstein = bernstein::evalBernsteinPoly(phiBernstein, testX);
    auto vec               = xarray2StdVector(phiBernstein);
    std::cout << "eval bernstein:" << testEvalBernstein << std::endl;
    std::cout << "test8PlanesPowerDirectly" << std::endl;
    qConvBernstein(phiBernstein,
                   -1,
                   1,
                   integrand,
                   10,
                   {
                       {1,  1,  1,  -1},
                       {1,  1,  -1, -1},
                       {1,  -1, 1,  -1},
                       {1,  -1, -1, -1},
                       {-1, 1,  1,  -1},
                       {-1, 1,  -1, -1},
                       {-1, -1, 1,  -1},
                       {-1, -1, -1, -1}
    });
}

void test3PlanesAlgoim()
{
    const int                    N    = 3;
    const uvector<int, N>        P    = 3;
    real                         xmin = -1, xmax = 1;
    std::vector<xarray<real, N>> phis(3, xarray<real, N>(nullptr, P));
    auto                         fphiPlane1 = [](const uvector<real, 3>& xx) {
        real x = xx(0);
        real y = xx(1);
        return x + y;
    };
    auto fphiPlane2 = [](const uvector<real, 3>& xx) {
        real x = xx(0);
        real y = xx(1);
        return x - y;
    };
    auto fphiPlane3 = [](const uvector<real, 3>& xx) {
        real x = xx(0);
        return -x - 1;
    };


    algoim_spark_alloc(real, phis[0]);
    algoim_spark_alloc(real, phis[1]);
    algoim_spark_alloc(real, phis[2]);
    bernstein::bernsteinInterpolate<N>([&](const uvector<real, N>& x) { return fphiPlane1(xmin + x * (xmax - xmin)); },
                                       phis[0]);
    bernstein::bernsteinInterpolate<N>([&](const uvector<real, N>& x) { return fphiPlane2(xmin + x * (xmax - xmin)); },
                                       phis[1]);
    bernstein::bernsteinInterpolate<N>([&](const uvector<real, N>& x) { return fphiPlane3(xmin + x * (xmax - xmin)); },
                                       phis[2]);

    auto integrand = [](const uvector<real, 3>& x) { return 1.0; };
    // // qConvMultiPloy<3>(phi0, phi1, -1.0, 1.0, 3, integrand, 3, "exampleC");
    std::cout << "testPlanesAlgoim" << std::endl;
    qConvMultiPloy2<3>(-1.0, 1.0, phis, integrand, 10);
}

void test3PlanesPowerDirectly()
{
    // a_x(x-c_x)^2 + a_y(y-c_y)^2 + a_z(x-c_z)^2 - r^2
    uvector<real, 3> xmin = -1, xmax = 1;
    uvector<real, 3> range = xmax - xmin;
    assert(all(range != 0));
    uvector<real, 3>             k    = range;
    uvector<real, 3>             bias = xmin;
    uvector<int, 3>              ext  = 2;
    std::vector<xarray<real, 3>> phiPowers(3, xarray<real, 3>(nullptr, ext));
    algoim_spark_alloc(real, phiPowers[0]);
    algoim_spark_alloc(real, phiPowers[1]);
    algoim_spark_alloc(real, phiPowers[2]);
    for (int i = 0; i < 3; i++) { xarrayInit(phiPowers[i]); }
    phiPowers[0].m(uvector<int, 3>(1, 0, 0)) = 1;
    phiPowers[0].m(uvector<int, 3>(0, 1, 0)) = 1;
    phiPowers[1].m(uvector<int, 3>(1, 0, 0)) = 1;
    phiPowers[1].m(uvector<int, 3>(0, 1, 0)) = -1;
    phiPowers[2].m(uvector<int, 3>(1, 0, 0)) = -1;
    phiPowers[2].m(uvector<int, 3>(0, 0, 0)) = -1;
    uvector<int, 3> resExt                   = 1 + phiPowers.size();
    xarray<real, 3> phiPowerAll(nullptr, resExt), phiBernstein(nullptr, resExt);
    algoim_spark_alloc(real, phiPowerAll);
    algoim_spark_alloc(real, phiBernstein);
    compositePower(phiPowers, 0, 0, 1, phiPowerAll);
    auto v = xarray2StdVector(phiPowerAll);

    uvector<real, 3> testX(0., 0.75, 0);

    real testEval = powerEvaluation(phiPowerAll, testX);
    std::cout << "eval power before trans:" << testEval << std::endl;

    powerTransformation(range, xmin, phiPowerAll);

    auto vAfterTrans = xarray2StdVector(phiPowerAll);
    auto integrand   = [](const uvector<real, 3>& x) { return 1.0; };
    power2BernsteinTensorProduct(phiPowerAll, phiBernstein);
    // bernstein::normalise(phiBernstein);

    real testEval0 = powerEvaluation(phiPowers[0], testX);
    std::cout << "eval power0:" << testEval0 << std::endl;

    real testEval1 = powerEvaluation(phiPowers[1], testX);
    std::cout << "eval power1:" << testEval1 << std::endl;

    real testEval2 = powerEvaluation(phiPowers[2], testX);
    std::cout << "eval power2:" << testEval2 << std::endl;


    real testEvalBernstein = bernstein::evalBernsteinPoly(phiBernstein, testX);
    auto vec               = xarray2StdVector(phiBernstein);
    std::cout << "eval bernstein:" << testEvalBernstein << std::endl;
    std::cout << "test8PlanesPowerDirectly" << std::endl;
    qConvBernstein(phiBernstein,
                   -1,
                   1,
                   integrand,
                   10,
                   {
                       {1,  1,  0, 0 },
                       {1,  -1, 0, 0 },
                       {-1, 0,  0, -1},
    });
}

void testConicAlgoim()
{
    const int                    N    = 3;
    const uvector<int, N>        P    = 3;
    real                         xmin = -1, xmax = 1;
    std::vector<xarray<real, N>> phis(3, xarray<real, N>(nullptr, P));
    auto                         fphiPlane1 = [](const uvector<real, 3>& xx) {
        real x = xx(0);
        real y = xx(1);
        real z = xx(2);
        return x * x + y * y - z * z;
    };
    auto fphiPlane2 = [](const uvector<real, 3>& xx) {
        real x = xx(0);
        real y = xx(1);
        real z = xx(2);
        return x * z - y;
    };
    auto fphiPlane3 = [](const uvector<real, 3>& xx) {
        real x = xx(0);
        return -x - 1;
    };


    algoim_spark_alloc(real, phis[0]);
    algoim_spark_alloc(real, phis[1]);
    algoim_spark_alloc(real, phis[2]);
    bernstein::bernsteinInterpolate<N>([&](const uvector<real, N>& x) { return fphiPlane1(xmin + x * (xmax - xmin)); },
                                       phis[0]);
    bernstein::bernsteinInterpolate<N>([&](const uvector<real, N>& x) { return fphiPlane2(xmin + x * (xmax - xmin)); },
                                       phis[1]);
    bernstein::bernsteinInterpolate<N>([&](const uvector<real, N>& x) { return fphiPlane3(xmin + x * (xmax - xmin)); },
                                       phis[2]);

    auto integrand = [](const uvector<real, 3>& x) { return 1.0; };
    // // qConvMultiPloy<3>(phi0, phi1, -1.0, 1.0, 3, integrand, 3, "exampleC");
    std::cout << "testPlanesAlgoim" << std::endl;
    qConvMultiPloy2<3>(-1.0, 1.0, phis, integrand, 10);
}

void testQuarterSphere()
{
    uvector<real, 3> xmin = -1, xmax = 1;
    uvector<real, 3> range = xmax - xmin;
    assert(all(range != 0));
    uvector<real, 3>             k    = range;
    uvector<real, 3>             bias = xmin;
    std::vector<xarray<real, 3>> phiPowers;
    phiPowers.push_back(xarray<real, 3>(nullptr, 3));
    phiPowers.push_back(xarray<real, 3>(nullptr, 2));
    phiPowers.push_back(xarray<real, 3>(nullptr, 2));
    algoim_spark_alloc(real, phiPowers[0]);
    algoim_spark_alloc(real, phiPowers[1]);
    algoim_spark_alloc(real, phiPowers[2]);
    for (int i = 0; i < 3; i++) { xarrayInit(phiPowers[i]); }
    phiPowers[0].m(uvector<int, 3>(2, 0, 0)) = 1;
    phiPowers[0].m(uvector<int, 3>(0, 2, 0)) = 1;
    phiPowers[0].m(uvector<int, 3>(0, 0, 2)) = 1;
    phiPowers[0].m(uvector<int, 3>(0, 0, 0)) = -1;
    phiPowers[1].m(uvector<int, 3>(1, 0, 0)) = 1;
    phiPowers[1].m(uvector<int, 3>(0, 1, 0)) = 1;
    // phiPowers[1].m(uvector<int, 3>(0, 0, 1)) = 1;
    phiPowers[2].m(uvector<int, 3>(1, 0, 0)) = 1;
    phiPowers[2].m(uvector<int, 3>(0, 1, 0)) = -1;
    // phiPowers[2].m(uvector<int, 3>(0, 0, 1)) = -1;

    uvector<int, 3> resExt = 1 + 2 + 1 + 1;
    xarray<real, 3> phiPowerAll(nullptr, resExt), phiBernstein(nullptr, resExt);
    algoim_spark_alloc(real, phiPowerAll);
    algoim_spark_alloc(real, phiBernstein);
    compositePower(phiPowers, 0, 0, 1, phiPowerAll);
    auto v = xarray2StdVector(phiPowerAll);

    uvector<real, 3> testX(0., 0.75, 0);

    real testEval = powerEvaluation(phiPowerAll, testX);
    std::cout << "eval power before trans:" << testEval << std::endl;

    powerTransformation(range, xmin, phiPowerAll);

    auto vAfterTrans = xarray2StdVector(phiPowerAll);
    auto integrand   = [](const uvector<real, 3>& x) { return 1.0; };
    power2BernsteinTensorProduct(phiPowerAll, phiBernstein);
    // bernstein::normalise(phiBernstein);

    real testEval0 = powerEvaluation(phiPowers[0], testX);
    std::cout << "eval power0:" << testEval0 << std::endl;

    real testEval1 = powerEvaluation(phiPowers[1], testX);
    std::cout << "eval power1:" << testEval1 << std::endl;

    real testEval2 = powerEvaluation(phiPowers[2], testX);
    std::cout << "eval power2:" << testEval2 << std::endl;


    real testEvalBernstein = bernstein::evalBernsteinPoly(phiBernstein, testX);
    auto vec               = xarray2StdVector(phiBernstein);
    std::cout << "eval bernstein:" << testEvalBernstein << std::endl;
    std::cout << "testSphereQuarterDirectly" << std::endl;
    qConvBernstein(phiBernstein,
                   -1,
                   1,
                   integrand,
                   30,
                   {
                       {1, 1,  0, 0},
                       {1, -1, 0, 0},
    });
}

void test4PlanesAlgoim()
{
    const int                    N    = 3;
    const uvector<int, N>        P    = 3;
    real                         xmin = -1, xmax = 1;
    std::vector<xarray<real, N>> phis(4, xarray<real, N>(nullptr, P));
    auto                         fphiPlane1 = [](const uvector<real, 3>& xx) {
        real x = xx(0);
        real y = xx(1);
        return x + y;
    };
    auto fphiPlane2 = [](const uvector<real, 3>& xx) {
        real x = xx(0);
        real y = xx(1);
        return x - y;
    };
    auto fphiPlane3 = [](const uvector<real, 3>& xx) {
        real x = xx(0);
        return -x - 1;
    };
    auto fphiPlane4 = [](const uvector<real, 3>& xx) {
        real x = xx(0);
        real y = xx(1);
        real z = xx(2);
        // return x + y - z;
        return y - 1;
    };


    algoim_spark_alloc(real, phis[0]);
    algoim_spark_alloc(real, phis[1]);
    algoim_spark_alloc(real, phis[2]);
    algoim_spark_alloc(real, phis[3]);
    bernstein::bernsteinInterpolate<N>([&](const uvector<real, N>& x) { return fphiPlane1(xmin + x * (xmax - xmin)); },
                                       phis[0]);
    bernstein::bernsteinInterpolate<N>([&](const uvector<real, N>& x) { return fphiPlane2(xmin + x * (xmax - xmin)); },
                                       phis[1]);
    bernstein::bernsteinInterpolate<N>([&](const uvector<real, N>& x) { return fphiPlane3(xmin + x * (xmax - xmin)); },
                                       phis[2]);
    bernstein::bernsteinInterpolate<N>([&](const uvector<real, N>& x) { return fphiPlane4(xmin + x * (xmax - xmin)); },
                                       phis[3]);

    auto integrand = [](const uvector<real, 3>& x) { return 1.0; };
    // // qConvMultiPloy<3>(phi0, phi1, -1.0, 1.0, 3, integrand, 3, "exampleC");
    std::cout << "test4PlanesAlgoim" << std::endl;
    qConvMultiPloy2<3>(-1.0, 1.0, phis, integrand, 10);
}

void test4PlanesPowerDirectly()
{
    // a_x(x-c_x)^2 + a_y(y-c_y)^2 + a_z(x-c_z)^2 - r^2
    uvector<real, 3> xmin = -1, xmax = 1;
    uvector<real, 3> range = xmax - xmin;
    assert(all(range != 0));
    uvector<real, 3>             k    = range;
    uvector<real, 3>             bias = xmin;
    uvector<int, 3>              ext  = 2;
    std::vector<xarray<real, 3>> phiPowers(4, xarray<real, 3>(nullptr, ext));
    algoim_spark_alloc(real, phiPowers[0]);
    algoim_spark_alloc(real, phiPowers[1]);
    algoim_spark_alloc(real, phiPowers[2]);
    algoim_spark_alloc(real, phiPowers[3]);
    for (int i = 0; i < 4; i++) { xarrayInit(phiPowers[i]); }
    phiPowers[0].m(uvector<int, 3>(1, 0, 0)) = 1;
    phiPowers[0].m(uvector<int, 3>(0, 1, 0)) = 1;
    phiPowers[1].m(uvector<int, 3>(1, 0, 0)) = 1;
    phiPowers[1].m(uvector<int, 3>(0, 1, 0)) = -1;
    phiPowers[2].m(uvector<int, 3>(1, 0, 0)) = -1;
    phiPowers[2].m(uvector<int, 3>(0, 0, 0)) = -1;
    phiPowers[3].m(uvector<int, 3>(1, 0, 0)) = 1;
    phiPowers[3].m(uvector<int, 3>(0, 1, 0)) = 1;
    phiPowers[3].m(uvector<int, 3>(0, 0, 1)) = -1;
    uvector<int, 3> resExt                   = 1 + phiPowers.size();
    xarray<real, 3> phiPowerAll(nullptr, resExt), phiBernstein(nullptr, resExt);
    algoim_spark_alloc(real, phiPowerAll);
    algoim_spark_alloc(real, phiBernstein);
    compositePower(phiPowers, 0, 0, 1, phiPowerAll);
    auto v = xarray2StdVector(phiPowerAll);

    uvector<real, 3> testX(0., 0.75, 0.2);

    real testEval = powerEvaluation(phiPowerAll, testX);
    std::cout << "eval power before trans:" << testEval << std::endl;

    powerTransformation(range, xmin, phiPowerAll);

    real testEvalAfterTrans = powerEvaluation(phiPowerAll, testX);
    std::cout << "eval power after trans:" << testEvalAfterTrans << std::endl;

    auto vAfterTrans = xarray2StdVector(phiPowerAll);
    auto integrand   = [](const uvector<real, 3>& x) { return 1.0; };
    power2BernsteinTensorProduct(phiPowerAll, phiBernstein);
    // bernstein::normalise(phiBernstein);

    real testEval0 = powerEvaluation(phiPowers[0], testX);
    std::cout << "eval power0:" << testEval0 << std::endl;

    real testEval1 = powerEvaluation(phiPowers[1], testX);
    std::cout << "eval power1:" << testEval1 << std::endl;

    real testEval2 = powerEvaluation(phiPowers[2], testX);
    std::cout << "eval power2:" << testEval2 << std::endl;


    real testEvalBernstein = bernstein::evalBernsteinPoly(phiBernstein, testX);
    auto vec               = xarray2StdVector(phiBernstein);
    std::cout << "eval bernstein:" << testEvalBernstein << std::endl;
    std::cout << "test4PlanesPowerDirectly" << std::endl;
    qConvBernstein(phiBernstein,
                   -1,
                   1,
                   integrand,
                   10,
                   {
                       {1,  1,  0,  0 },
                       {1,  -1, 0,  0 },
                       {-1, 0,  0,  -1},
                       {1,  1,  -1, 0 },
    });
}

void testMultiScale()
{
    auto phi0 = [](const uvector<real, 3>& xx) {
        real x = xx(0) + 100000.;
        real y = xx(1);
        real z = xx(2) - 100000.;
        real r = 800000.;
        return x * x + y * y + z * z - r * r;
    };

    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) + 100000.;
        real y = xx(1) - 800000.;
        real z = xx(2) - 100000.;
        real r = 1000.;
        return x * x + y * y + z * z - r * r;
    };

    auto phi2 = [](const uvector<real, 3>& xx) {
        real z = xx(2);
        return -z - 100000.;
    };
    auto integrand = [](const uvector<real, 3>& x) { return 1.0; };
    qConvMultiPloy<3>(phi0, phi2, -1000000., 1000000., 3, integrand, 20, Difference);
}

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 testPower2Bernstein()
{
    uvector<real, 3> ext(3, 3, 3);
    xarray<real, 3>  power(nullptr, ext), bernstein(nullptr, ext);
    algoim_spark_alloc(real, power, bernstein);
    xarrayInit(power);
    power.m(uvector<int, 3>(1, 2, 0)) = 3;
    power.m(uvector<int, 3>(0, 1, 2)) = 4;
    power.m(uvector<int, 3>(1, 2, 2)) = 5;
    power.m(uvector<int, 3>(0, 0, 0)) = 3;
    uvector3 X(0.8, 0.7, 0.6);
    auto     evalPower = powerEvaluation(power, X);
    power2BernsteinTensorProduct(power, bernstein);
    auto evalBernstein = bernstein::evalBernsteinPoly(bernstein, X);
    std::cout << "eval power: " << evalPower << std::endl;
    std::cout << "eval bernstein: " << evalBernstein << std::endl;
}

void testMain()
{
    // testBooluarray();
    // testMultiPolys();
    // testMultiScale();
    // testSpherePowerDirectly();
    // testSpherePowerDirectlyByTransformation();
    // testCylinderPowerDirectly();
    // testCylinderAlgoim();

    // test8PlanesAlgoim();
    // test8PlanesPowerDirectly();


    // test4PlanesPowerDirectly();
    // test4PlanesAlgoim();

    // testQuarterSphere();

    // testPower2Bernstein();
}