TasmanianSparseGrid.hpp

/*
 * Copyright (c) 2017, Miroslav Stoyanov
 *
 * This file is part of
 * Toolkit for Adaptive Stochastic Modeling And Non-Intrusive ApproximatioN: TASMANIAN
 *
 * Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
 *
 * 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
 *
 * 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions
 *    and the following disclaimer in the documentation and/or other materials provided with the distribution.
 *
 * 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse
 *    or promote products derived from this software without specific prior written permission.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,
 * INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
 * IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY,
 * OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA,
 * OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
 * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 *
 * UT-BATTELLE, LLC AND THE UNITED STATES GOVERNMENT MAKE NO REPRESENTATIONS AND DISCLAIM ALL WARRANTIES, BOTH EXPRESSED AND IMPLIED.
 * THERE ARE NO EXPRESS OR IMPLIED WARRANTIES OF MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE, OR THAT THE USE OF THE SOFTWARE WILL NOT INFRINGE ANY PATENT,
 * COPYRIGHT, TRADEMARK, OR OTHER PROPRIETARY RIGHTS, OR THAT THE SOFTWARE WILL ACCOMPLISH THE INTENDED RESULTS OR THAT THE SOFTWARE OR ITS USE WILL NOT RESULT IN INJURY OR DAMAGE.
 * THE USER ASSUMES RESPONSIBILITY FOR ALL LIABILITIES, PENALTIES, FINES, CLAIMS, CAUSES OF ACTION, AND COSTS AND EXPENSES, CAUSED BY, RESULTING FROM OR ARISING OUT OF,
 * IN WHOLE OR IN PART THE USE, STORAGE OR DISPOSAL OF THE SOFTWARE.
 */

 
#ifndef __TASMANIAN_SPARSE_GRID_HPP
#define __TASMANIAN_SPARSE_GRID_HPP
 
#include "tsgGridGlobal.hpp"
#include "tsgGridLocalPolynomial.hpp"
#include "tsgGridWavelet.hpp"
#include "tsgGridFourier.hpp"


namespace TasGrid{


class TasmanianSparseGrid{
public:
    TasmanianSparseGrid();
    TasmanianSparseGrid(const TasmanianSparseGrid &source);
    TasmanianSparseGrid(TasmanianSparseGrid &&source) = default;
    ~TasmanianSparseGrid() = default;

    TasmanianSparseGrid& operator=(TasmanianSparseGrid const &source);
    TasmanianSparseGrid& operator=(TasmanianSparseGrid &&source) = default;

    static const char* getVersion(); // human readable
    static int getVersionMajor();
    static int getVersionMinor();
    static const char* getLicense(); // human readable
    static const char* getGitCommitHash();
    static const char* getCmakeCxxFlags();
    static bool isOpenMPEnabled();
    static bool isCudaEnabled();
    static bool isHipEnabled();
    static bool isDpcppEnabled();

    void write(const char *filename, bool binary = mode_binary) const;
    void read(const char *filename); // auto-check if format is binary or ascii

    void write(std::ostream &ofs, bool binary = mode_binary) const;
    void read(std::istream &ifs, bool binary = mode_binary);

    void write(std::string const& fname, bool binary = mode_binary) const{ write(fname.c_str(), binary); }
    void read(std::string const& fname){ read(fname.c_str()); }

    void makeGlobalGrid(int dimensions, int outputs, int depth, TypeDepth type, TypeOneDRule rule,
                        std::vector<int> const &anisotropic_weights, double alpha = 0.0, double beta = 0.0,
                        const char* custom_filename = nullptr, std::vector<int> const &level_limits = std::vector<int>());
    void makeGlobalGrid(int dimensions, int outputs, int depth, TypeDepth type, TypeOneDRule rule,
                        const int *anisotropic_weights = nullptr, double alpha = 0.0, double beta = 0.0,
                        const char* custom_filename = nullptr, const int *level_limits = nullptr);
    void makeGlobalGrid(int dimensions, int outputs, int depth, TypeDepth type, CustomTabulated &&crule,
                        std::vector<int> const &anisotropic_weights, std::vector<int> const &level_limits = std::vector<int>());
    void makeGlobalGrid(int dimensions, int outputs, int depth, TypeDepth type, CustomTabulated &&crule,
                        const int *anisotropic_weights = nullptr, const int *level_limits = nullptr);

    void makeSequenceGrid(int dimensions, int outputs, int depth, TypeDepth type, TypeOneDRule rule,
                          std::vector<int> const &anisotropic_weights, std::vector<int> const &level_limits = std::vector<int>());
    void makeSequenceGrid(int dimensions, int outputs, int depth, TypeDepth type, TypeOneDRule rule,
                          const int *anisotropic_weights = nullptr, const int *level_limits = nullptr);

    void makeLocalPolynomialGrid(int dimensions, int outputs, int depth, int order, TypeOneDRule rule, std::vector<int> const &level_limits);
    void makeLocalPolynomialGrid(int dimensions, int outputs, int depth, int order = 1, TypeOneDRule rule = rule_localp, const int *level_limits = nullptr);

    void makeWaveletGrid(int dimensions, int outputs, int depth, int order, std::vector<int> const &level_limits);
    void makeWaveletGrid(int dimensions, int outputs, int depth, int order = 1, const int *level_limits = nullptr);

    void makeFourierGrid(int dimensions, int outputs, int depth, TypeDepth type,
                         std::vector<int> const &anisotropic_weights, std::vector<int> const &level_limits = std::vector<int>());
    void makeFourierGrid(int dimensions, int outputs, int depth, TypeDepth type, const int* anisotropic_weights = nullptr, const int* level_limits = nullptr);

    void copyGrid(const TasmanianSparseGrid *source, int outputs_begin = 0, int outputs_end = -1);
    void copyGrid(TasmanianSparseGrid const &source, int outputs_begin = 0, int outputs_end = -1){ copyGrid(&source, outputs_begin, outputs_end); }

    void updateGlobalGrid(int depth, TypeDepth type, std::vector<int> const &anisotropic_weights, std::vector<int> const &level_limits = std::vector<int>());
    void updateGlobalGrid(int depth, TypeDepth type, const int *anisotropic_weights = nullptr, const int *level_limits = nullptr);

    void updateSequenceGrid(int depth, TypeDepth type, std::vector<int> const &anisotropic_weights, std::vector<int> const &level_limits = std::vector<int>());
    void updateSequenceGrid(int depth, TypeDepth type, const int *anisotropic_weights = nullptr, const int *level_limits = nullptr);

    void updateFourierGrid(int depth, TypeDepth type, std::vector<int> const &anisotropic_weights, std::vector<int> const &level_limits = std::vector<int>());
    void updateFourierGrid(int depth, TypeDepth type, const int *anisotropic_weights = nullptr, const int *level_limits = nullptr);
    void updateGrid(int depth, TypeDepth type, std::vector<int> const &anisotropic_weights, std::vector<int> const &level_limits = std::vector<int>());
    void updateGrid(int depth, TypeDepth type, const int *anisotropic_weights = nullptr, const int *level_limits = nullptr);

    double getAlpha() const{ return (isGlobal()) ? get<GridGlobal>()->getAlpha() : 0.0; }
    double getBeta() const{ return (isGlobal()) ? get<GridGlobal>()->getBeta() : 0.0; }
    int getOrder() const{ return (isLocalPolynomial()) ? get<GridLocalPolynomial>()->getOrder() : ((isWavelet()) ? get<GridWavelet>()->getOrder() : -1); }
 

    int getNumDimensions() const{ return makeCall<BaseMethod::getNumDimensions>(); }
    int getNumOutputs() const{ return makeCall<BaseMethod::getNumOutputs>(); }
    TypeOneDRule getRule() const;
    const char* getCustomRuleDescription() const; // used only for Global Grids with rule_customtabulated

    int getNumLoaded() const{ return makeCall<BaseMethod::getNumLoaded>(); }
    int getNumNeeded() const{ return makeCall<BaseMethod::getNumNeeded>(); }
    int getNumPoints() const{ return makeCall<BaseMethod::getNumPoints>(); }

    std::vector<double> getLoadedPoints() const{ std::vector<double> x; getLoadedPoints(x); return x; }
    void getLoadedPoints(std::vector<double> &x) const{ x.resize(Utils::size_mult(getNumDimensions(), getNumLoaded())); getLoadedPoints(x.data()); }
    void getLoadedPoints(double x[]) const;

    std::vector<double> getNeededPoints() const{ std::vector<double> x; getNeededPoints(x); return x; }
    void getNeededPoints(std::vector<double> &x) const{ x.resize(Utils::size_mult(getNumDimensions(), getNumNeeded())); getNeededPoints(x.data()); }
    void getNeededPoints(double *x) const;

    std::vector<double> getPoints() const{ std::vector<double> x; getPoints(x); return x; }
    void getPoints(std::vector<double> &x) const{ x.resize(Utils::size_mult(getNumDimensions(), getNumPoints())); getPoints(x.data()); }
    void getPoints(double x[]) const; // returns the loaded points unless no points are loaded, then returns the needed points

    std::vector<double> getQuadratureWeights() const{ std::vector<double> w; getQuadratureWeights(w); return w; }
    void getQuadratureWeights(std::vector<double> &weights) const{ weights.resize((size_t) getNumPoints()); getQuadratureWeights(weights.data()); }
    void getQuadratureWeights(double weights[]) const; // static memory, assumes that weights has size getNumPoints()

    std::vector<double> getInterpolationWeights(std::vector<double> const &x) const;
    std::vector<double> getInterpolationWeights(double const x[]) const{ std::vector<double> w((size_t) getNumPoints()); getInterpolationWeights(x, w.data()); return w; }
    void getInterpolationWeights(const std::vector<double> &x, std::vector<double> &weights) const;
    void getInterpolationWeights(const double x[], double weights[]) const;

    std::vector<double> getDifferentiationWeights(std::vector<double> const &x) const;
    std::vector<double> getDifferentiationWeights(double const x[]) const {
        std::vector<double> w((size_t) getNumPoints() * (size_t) getNumDimensions());
        getDifferentiationWeights(x, w.data());
        return w;
    }

    void getDifferentiationWeights(const std::vector<double> &x, std::vector<double> &weights) const;
    void getDifferentiationWeights(const double x[], double weights[]) const;

    void loadNeededValues(std::vector<double> const &vals); // checks if vals has size num_outputs X getNumNeeded()
    void loadNeededValues(const double *vals);
    void loadNeededPoints(std::vector<double> const &vals) {loadNeededValues(vals);}
    void loadNeededPoints(const double *vals) {loadNeededValues(vals);}
    const double* getLoadedValues() const;

    void evaluate(std::vector<double> const &x, std::vector<double> &y) const;
    void evaluate(const double x[], double y[]) const;
    template<typename FloatType> void evaluateBatch(std::vector<FloatType> const &x, std::vector<FloatType> &y) const;
    void evaluateBatch(const double x[], int num_x, double y[]) const;
    void evaluateBatch(const float x[], int num_x, float y[]) const;
    template<typename FloatType> void evaluateBatchGPU(const FloatType gpu_x[], int cpu_num_x, FloatType gpu_y[]) const;
    template<typename FloatType> void evaluateFast(const FloatType x[], FloatType y[]) const{ evaluateBatch(x, 1, y); }
    template<typename FloatType> void evaluateFast(std::vector<FloatType> const &x, std::vector<FloatType> &y) const{ evaluateBatch(x, y); }
    void integrate(std::vector<double> &q) const;
    void integrate(double q[]) const;
    std::vector<double> integrate() const{
        std::vector<double> result;
        integrate(result);
        return result;
    }

    void differentiate(std::vector<double> const &x, std::vector<double> &jacobian) const;
    std::vector<double> differentiate(std::vector<double> const &x) const{
        std::vector<double> jacobian;
        differentiate(x, jacobian);
        return jacobian;
    }

    void differentiate(const double x[], double jacobian[]) const;

    bool isEmpty() const{ return std::holds_alternative<std::monostate>(base); }
    bool empty() const{ return std::holds_alternative<std::monostate>(base); }
    bool isGlobal() const{ return std::holds_alternative<GridGlobal>(base); }
    bool isSequence() const{ return std::holds_alternative<GridSequence>(base); }
    bool isFourier() const{ return std::holds_alternative<GridFourier>(base); }
    bool isLocalPolynomial() const{ return std::holds_alternative<GridLocalPolynomial>(base); }
    bool isWavelet() const{ return std::holds_alternative<GridWavelet>(base); }

    void setDomainTransform(std::vector<double> const &a, std::vector<double> const &b);
    void setDomainTransform(const double a[], const double b[]);
    bool isSetDomainTransfrom() const;
    void clearDomainTransform();
    void getDomainTransform(std::vector<double> &a, std::vector<double> &b) const;
    void getDomainTransform(double a[], double b[]) const;

    void setConformalTransformASIN(std::vector<int> const &truncation);
    bool isSetConformalTransformASIN() const;
    void clearConformalTransform();
    std::vector<int> getConformalTransformASIN() const;

    void clearLevelLimits(){ llimits.clear(); }
    std::vector<int> getLevelLimits() const{ return llimits; }

    void setAnisotropicRefinement(TypeDepth type, int min_growth, int output, const std::vector<int> &level_limits);
    void setAnisotropicRefinement(TypeDepth type, int min_growth, int output, const int *level_limits = nullptr);
    std::vector<double> getAnisotropicRefinement(TypeDepth type, int min_growth, int output, const std::vector<int> &level_limits) {
        setAnisotropicRefinement(type, min_growth, output, level_limits);
        return getNeededPoints();
    }

    std::vector<int> estimateAnisotropicCoefficients(TypeDepth type, int output) const{ std::vector<int> w; estimateAnisotropicCoefficients(type, output, w); return w; }
    void estimateAnisotropicCoefficients(TypeDepth type, int output, std::vector<int> &weights) const;

    void setSurplusRefinement(double tolerance, int output, std::vector<int> const &level_limits);
    void setSurplusRefinement(double tolerance, int output, const int *level_limits = nullptr);
    std::vector<double> getSurplusRefinement(double tolerance, int output, std::vector<int> const &level_limits) {
        setSurplusRefinement(tolerance, output, level_limits);
        return getNeededPoints();
    }

    void setSurplusRefinement(double tolerance, TypeRefinement criteria, int output,
                              std::vector<int> const &level_limits, std::vector<double> const &scale_correction = std::vector<double>()); // -1 indicates using all outputs
    void setSurplusRefinement(double tolerance, TypeRefinement criteria, int output = -1, const int *level_limits = nullptr, const double *scale_correction = nullptr);
    std::vector<double> getSurplusRefinement(double tolerance, TypeRefinement criteria, int output, std::vector<int> const &level_limits,
                                             std::vector<double> const &scale_correction = std::vector<double>()) {
        setSurplusRefinement(tolerance, criteria, output, level_limits, scale_correction);
        return getNeededPoints();
    }

    void clearRefinement();
    void mergeRefinement(); // merges all points but resets all loaded values to 0.0

    void beginConstruction();
    bool isUsingConstruction() const{ return using_dynamic_construction; }
    std::vector<double> getCandidateConstructionPoints(TypeDepth type, std::vector<int> const &anisotropic_weights = std::vector<int>(),
                                                       std::vector<int> const &level_limits = std::vector<int>());
    std::vector<double> getCandidateConstructionPoints(TypeDepth type, int output, std::vector<int> const &level_limits = std::vector<int>());

    std::vector<double> getCandidateConstructionPoints(double tolerance, TypeRefinement criteria, int output = -1,
                                                       std::vector<int> const &level_limits = std::vector<int>(),
                                                       std::vector<double> const &scale_correction = std::vector<double>());
    void loadConstructedPoints(std::vector<double> const &x, std::vector<double> const &y);
    void loadConstructedPoints(const double x[], int numx, const double y[]);
    void finishConstruction();

    const double* getHierarchicalCoefficients() const;

    void getHierarchicalCoefficientsStatic(double *coeff) const{
        std::copy(getHierarchicalCoefficients(), getHierarchicalCoefficients()
                  + ((isFourier()) ? 2 : 1) * getNumOutputs() * getNumLoaded(), coeff);
    }

    void setHierarchicalCoefficients(const std::vector<double> &c);
    void setHierarchicalCoefficients(const double c[]);

    void evaluateHierarchicalFunctions(std::vector<double> const &x, std::vector<double> &y) const;
    std::vector<double> evaluateHierarchicalFunctions(std::vector<double> const &x) const{
        std::vector<double> y;
        evaluateHierarchicalFunctions(x, y);
        return y;
    }

    void evaluateHierarchicalFunctions(const double x[], int num_x, double y[]) const;

    void evaluateSparseHierarchicalFunctions(const std::vector<double> &x, std::vector<int> &pntr, std::vector<int> &indx, std::vector<double> &vals) const;

    std::vector<double> getHierarchicalSupport() const;

    std::vector<double> integrateHierarchicalFunctions() const{
        std::vector<double> integrals;
        integrateHierarchicalFunctions(integrals);
        return integrals;
    }

    void integrateHierarchicalFunctions(std::vector<double> &integrals) const{
        integrals.resize(getNumPoints());
        integrateHierarchicalFunctions(integrals.data());
    }

    void integrateHierarchicalFunctions(double integrals[]) const;

    std::vector<int> getGlobalPolynomialSpace(bool interpolation) const;

    void printStats(std::ostream &os = std::cout) const;

    void enableAcceleration(TypeAcceleration acc);
    void enableAcceleration(TypeAcceleration acc, int new_gpu_id);
    void favorSparseAcceleration(bool favor);
    void setCuBlasHandle(void *handle);
    void setCuSparseHandle(void *handle);
    void setCuSolverHandle(void *handle);
    void setRocBlasHandle(void *handle);
    void setRocSparseHandle(void *handle);
    void setSycleQueue(void *queue);
    TypeAcceleration getAccelerationType() const{ return acceleration->mode; }
    static bool isAccelerationAvailable(TypeAcceleration acc);

    void setGPUID(int new_gpu_id);
    int getGPUID() const{ return acceleration->device; }
    static int getNumGPUs(){ return AccelerationMeta::getNumGpuDevices(); }
    static int getGPUMemory(int gpu);
    static std::string getGPUName(int gpu);

    template<typename FloatType>
    void evaluateHierarchicalFunctionsGPU(const FloatType gpu_x[], int cpu_num_x, FloatType gpu_y[]) const;
    template<typename FloatType>
    void evaluateSparseHierarchicalFunctionsGPU(const FloatType gpu_x[], int cpu_num_x, int* &gpu_pntr, int* &gpu_indx, FloatType* &gpu_vals, int &num_nz) const;

    using EvaluateCallable = std::function<void(std::vector<double> const&, std::vector<double>&)>;

    operator EvaluateCallable() const{
        return [&](std::vector<double> const &x, std::vector<double> &y)->void{ evaluateBatch(x, y); };
    }

    using DomainInsideSignature = std::function<bool(std::vector<double> const &)>; // trick for Doxygen formatting purposes

     DomainInsideSignature getDomainInside() const{
        auto rule = getRule();
        if ((rule == TasGrid::rule_gausshermite) || (rule == TasGrid::rule_gausshermiteodd)){ // unbounded domain
            return [](std::vector<double> const &)->bool{ return true; };
        }else if ((rule == TasGrid::rule_gausslaguerre) || (rule == TasGrid::rule_gausslaguerreodd)){ // bounded from below
            if (domain_transform_a.empty()){
                return [](std::vector<double> const &x)->bool{
                    for(auto const &v : x) if (v < 0.0) return false;
                    return true;
                };
            }else{
                size_t dims = (size_t) getNumDimensions();
                return [dta=domain_transform_a, dims](std::vector<double> const &x)->bool{
                    for(size_t i=0; i<dims; i++) if (x[i] < dta[i]) return false;
                    return true;
                };
            }
        }else{
            if (domain_transform_a.empty()){
                if (isFourier()){
                    return [](std::vector<double> const &x)->bool{
                        for(auto const &v : x) if ((v < 0.0) || (v > 1.0)) return false;
                        return true;
                    };
                }else{
                    return [](std::vector<double> const &x)->bool{
                        for(auto const &v : x) if ((v < -1.0) || (v > 1.0)) return false;
                        return true;
                    };
                }
            }else{
                size_t dims = (size_t) getNumDimensions();
                return [dta=domain_transform_a, dtb=domain_transform_b, dims](std::vector<double> const &x)->bool{
                    for(size_t i=0; i<dims; i++)
                        if (x[i] < dta[i] or x[i] > dtb[i]) return false;
                    return true;
                };
            }
        }
    }

    void removePointsByHierarchicalCoefficient(double tolerance, int output = -1, const double *scale_correction = nullptr);
    void removePointsByHierarchicalCoefficient(int num_new_points, int output = -1, const double *scale_correction = nullptr);
 
    #ifndef __TASMANIAN_DOXYGEN_SKIP_INTERNAL
    int evaluateSparseHierarchicalFunctionsGetNZ(const double x[], int num_x) const;
    void evaluateSparseHierarchicalFunctionsStatic(const double x[], int num_x, int pntr[], int indx[], double vals[]) const;
    const int* getPointsIndexes() const;
    const int* getNeededIndexes() const;

    template<class T> inline T* get(){
        return std::get_if<T>(std::addressof(base));
    }
    template<class T> inline T const* get() const{
        return std::get_if<T>(std::addressof(base));
    }
    AccelerationContext const* getAccelerationContext() const{ return acceleration.get(); }
    #endif // __TASMANIAN_DOXYGEN_SKIP_INTERNAL
 
protected:
    #ifndef __TASMANIAN_DOXYGEN_SKIP_INTERNAL
    enum class BaseMethod {
        getNumDimensions,
        getNumOutputs,
        getNumLoaded,
        getNumNeeded,
        getNumPoints
    };
    template<BaseMethod method>
    int makeCall() const{
        return std::visit([](auto const &g) -> int {
            if constexpr (not std::is_same_v<std::decay_t<decltype(g)>, std::monostate>) {
                if constexpr (method == BaseMethod::getNumDimensions)
                    return g.getNumDimensions();
                else if constexpr (method == BaseMethod::getNumOutputs)
                    return g.getNumOutputs();
                else if constexpr (method == BaseMethod::getNumLoaded)
                    return g.getNumLoaded();
                else if constexpr (method == BaseMethod::getNumNeeded)
                    return g.getNumNeeded();
                else if constexpr (method == BaseMethod::getNumPoints)
                    return g.getNumPoints();
            } else {
                return int{0};
            }
        }, base);
    }

    void clear();

    void mapCanonicalToTransformed(int num_dimensions, int num_points, TypeOneDRule rule, double x[]) const;
    template<typename FloatType> void mapTransformedToCanonical(int num_dimensions, int num_points, TypeOneDRule rule, FloatType x[]) const;
    double getQuadratureScale(int num_dimensions, TypeOneDRule rule) const;

    void mapConformalCanonicalToTransformed(int num_dimensions, int num_points, double x[]) const;
    template<typename FloatType> void mapConformalTransformedToCanonical(int num_dimensions, int num_points, Data2D<FloatType> &x) const;
    void mapConformalWeights(int num_dimensions, int num_points, double weights[]) const;

    template<typename FloatType> const FloatType* formCanonicalPoints(const FloatType *x, Data2D<FloatType> &x_temp, int num_x) const;
    template<typename T>
    const T* formCanonicalPointsGPU(const T *gpu_x, int num_x, GpuVector<T> &gpu_x_temp) const;
    template<typename FloatType> std::vector<double> diffCanonicalTransform() const;
    void formTransformedPoints(int num_points, double x[]) const; // when calling get***Points()

    void writeAscii(std::ostream &ofs) const;
    void readAscii(std::istream &ifs);
    void writeBinary(std::ostream &ofs) const;
    void readBinary(std::istream &ifs);
    #endif // __TASMANIAN_DOXYGEN_SKIP_INTERNAL
 
private:
    // must not invalidate aliases when the grid is moved, must be destroyed last for sycl
    std::unique_ptr<AccelerationContext> acceleration;
 
    // contains all possible sparse grids with a null-grid std::monostate as the default
    using grid_variant = std::variant<std::monostate,
                                      GridGlobal,           // index 1
                                      GridSequence,         // index 2
                                      GridFourier,          // index 3
                                      GridLocalPolynomial,  // index 4
                                      GridWavelet           // index 5
                                      >;
    grid_variant base;
 
    std::vector<double> domain_transform_a, domain_transform_b;
    std::vector<int> conformal_asin_power;
    std::vector<int> llimits;
 
    bool using_dynamic_construction = false;
 
    mutable AccelerationDomainTransform acc_domain;
};

inline TasmanianSparseGrid makeEmpty(){ return TasmanianSparseGrid(); }

inline TasmanianSparseGrid
makeGlobalGrid(int dimensions, int outputs, int depth, TypeDepth type, TypeOneDRule rule,
               std::vector<int> const &anisotropic_weights = std::vector<int>(), double alpha = 0.0, double beta = 0.0,
               const char* custom_filename = nullptr, std::vector<int> const &level_limits = std::vector<int>()){
    TasmanianSparseGrid grid;
    grid.makeGlobalGrid(dimensions, outputs, depth, type, rule, anisotropic_weights, alpha, beta, custom_filename, level_limits);
    return grid;
}

inline TasmanianSparseGrid
makeSequenceGrid(int dimensions, int outputs, int depth, TypeDepth type, TypeOneDRule rule,
                 std::vector<int> const &anisotropic_weights = std::vector<int>(), std::vector<int> const &level_limits = std::vector<int>()){
    TasmanianSparseGrid grid;
    grid.makeSequenceGrid(dimensions, outputs, depth, type, rule, anisotropic_weights, level_limits);
    return grid;
}

inline TasmanianSparseGrid
makeLocalPolynomialGrid(int dimensions, int outputs, int depth, int order = 1, TypeOneDRule rule = rule_localp, std::vector<int> const &level_limits = std::vector<int>()){
    TasmanianSparseGrid grid;
    grid.makeLocalPolynomialGrid(dimensions, outputs, depth, order, rule, level_limits);
    return grid;
}

inline TasmanianSparseGrid
makeWaveletGrid(int dimensions, int outputs, int depth, int order = 1, std::vector<int> const &level_limits = std::vector<int>()){
    TasmanianSparseGrid grid;
    grid.makeWaveletGrid(dimensions, outputs, depth, order, level_limits);
    return grid;
}

inline TasmanianSparseGrid
makeFourierGrid(int dimensions, int outputs, int depth, TypeDepth type,
                std::vector<int> const &anisotropic_weights = std::vector<int>(), std::vector<int> const &level_limits = std::vector<int>()){
    TasmanianSparseGrid grid;
    grid.makeFourierGrid(dimensions, outputs, depth, type, anisotropic_weights, level_limits);
    return grid;
}

inline TasmanianSparseGrid readGrid(const char *filename){
    TasmanianSparseGrid grid;
    grid.read(filename);
    return grid;
}

inline TasmanianSparseGrid readGrid(std::string const &filename){ return readGrid(filename.c_str()); }

inline TasmanianSparseGrid copyGrid(TasmanianSparseGrid const &source, int outputs_begin = 0, int outputs_end = -1){
    TasmanianSparseGrid grid;
    grid.copyGrid(source, outputs_begin, outputs_end);
    return grid;
}


enum class NumDimensions : int {};

enum class NumOutputs : int {};

enum class DepthLevel : int {};

enum class BasisOrder : int {};

struct IsotropicWeights {};

struct AnsitropicWeights {
    AnsitropicWeights() = default;
    AnsitropicWeights(IsotropicWeights) {};
    AnsitropicWeights(std::initializer_list<int> list) : weights(list) {}
    explicit AnsitropicWeights(std::vector<int> const &w) : weights(w) {}
    std::vector<int> weights;
};

struct AlphaValue {
    AlphaValue() = default;
    explicit AlphaValue(double a) : alpha(a) {}
    double alpha = 0;
};

struct BetaValue {
    BetaValue() = default;
    explicit BetaValue(double b) : beta(b) {}
    double beta = 0;
};

struct CustomRuleFilename {
    CustomRuleFilename() = default;
    explicit CustomRuleFilename(std::string_view str) : name(str) {}
    const char *cname() const { return name.c_str(); }
    std::string name;
};

struct NoLevelLimits {};

struct LevelLimits {
    LevelLimits() = default;
    LevelLimits(NoLevelLimits) {};
    LevelLimits(std::initializer_list<int> list) : limits(list) {}
    explicit LevelLimits(std::vector<int> const &ll) : limits(ll) {}
    std::vector<int> limits;
};

inline TasmanianSparseGrid
makeGlobalGrid(NumDimensions dims, NumOutputs outs, DepthLevel depth, TypeDepth type, TypeOneDRule rule,
               AnsitropicWeights const &anisotropic_weights = IsotropicWeights{},
               AlphaValue alpha = AlphaValue{}, BetaValue beta = BetaValue{},
               CustomRuleFilename const &custom_filename = CustomRuleFilename{},
               LevelLimits const &level_limits = NoLevelLimits{}){
    TasmanianSparseGrid grid;
    grid.makeGlobalGrid(static_cast<int>(dims), static_cast<int>(outs), static_cast<int>(depth), type, rule, anisotropic_weights.weights, alpha.alpha, beta.beta, custom_filename.cname(), level_limits.limits);
    return grid;
}

inline TasmanianSparseGrid
makeSequenceGrid(NumDimensions dims, NumOutputs outs, DepthLevel depth, TypeDepth type, TypeOneDRule rule,
                 AnsitropicWeights const &anisotropic_weights = IsotropicWeights{},
                 LevelLimits const &level_limits = NoLevelLimits{}){
    TasmanianSparseGrid grid;
    grid.makeSequenceGrid(static_cast<int>(dims), static_cast<int>(outs), static_cast<int>(depth), type, rule, anisotropic_weights.weights, level_limits.limits);
    return grid;
}

inline TasmanianSparseGrid
makeLocalPolynomialGrid(NumDimensions dims, NumOutputs outs, DepthLevel depth,
                        BasisOrder order = BasisOrder{1}, TypeOneDRule rule = rule_localp,
                        LevelLimits const &level_limits = NoLevelLimits{}){
    TasmanianSparseGrid grid;
    grid.makeLocalPolynomialGrid(static_cast<int>(dims), static_cast<int>(outs), static_cast<int>(depth), static_cast<int>(order), rule, level_limits.limits);
    return grid;
}

inline TasmanianSparseGrid
makeWaveletGrid(NumDimensions dims, NumOutputs outs, DepthLevel depth,
                BasisOrder order = BasisOrder{1},
                LevelLimits const &level_limits = NoLevelLimits{}){
    TasmanianSparseGrid grid;
    grid.makeWaveletGrid(static_cast<int>(dims), static_cast<int>(outs), static_cast<int>(depth), static_cast<int>(order), level_limits.limits);
    return grid;
}

inline TasmanianSparseGrid
makeFourierGrid(NumDimensions dims, NumOutputs outs, DepthLevel depth, TypeDepth type,
                AnsitropicWeights const &anisotropic_weights = IsotropicWeights{},
                LevelLimits const &level_limits = NoLevelLimits{}){
    TasmanianSparseGrid grid;
    grid.makeFourierGrid(static_cast<int>(dims), static_cast<int>(outs), static_cast<int>(depth), type, anisotropic_weights.weights, level_limits.limits);
    return grid;
}
 
}
 
#endif