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 (base) ? base->getNumDimensions() : 0; }
    int getNumOutputs() const{ return (base) ? base->getNumOutputs() : 0; }
    TypeOneDRule getRule() const;
    const char* getCustomRuleDescription() const; // used only for Global Grids with rule_customtabulated
 
    int getNumLoaded() const{ return (base) ? base->getNumLoaded() : 0; }
    int getNumNeeded() const{ return (base) ? base->getNumNeeded() : 0; }
    int getNumPoints() const{ return (base) ? base->getNumPoints() : 0; }
 
    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{ return (empty()) ? nullptr : base->getLoadedValues(); }
 
    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 isGlobal() const{ return base && base->isGlobal(); }
    bool isSequence() const{ return base && base->isSequence(); }
    bool isLocalPolynomial() const{ return base && base->isLocalPolynomial(); }
    bool isWavelet() const{ return base && base->isWavelet(); }
    bool isFourier() const{ return base && base->isFourier(); }
    bool isEmpty() const{ return !base; }
    bool empty() const{ return !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[]){ base->setHierarchicalCoefficients(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();
                #if __cplusplus > 201103L
                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
                return [=](std::vector<double> const &x)->bool{
                    for(size_t i=0; i<dims; i++) if (x[i] < domain_transform_a[i]) return false;
                    return true;
                };
                #endif
            }
        }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();
                #if __cplusplus > 201103L
                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;
                };
                #else
                return [=](std::vector<double> const &x)->bool{
                    for(size_t i=0; i<dims; i++)
                        if (x[i] < domain_transform_a[i] or x[i] > domain_transform_b[i]) return false;
                    return true;
                };
                #endif
            }
        }
    }
 
    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 dynamic_cast<T*>(base.get()); }
    template<class T> inline T const* get() const{ return dynamic_cast<T const*>(base.get()); }
    AccelerationContext const* getAccelerationContext() const{ return acceleration.get(); }
    #endif // __TASMANIAN_DOXYGEN_SKIP_INTERNAL
 
protected:
    #ifndef __TASMANIAN_DOXYGEN_SKIP_INTERNAL
    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:
    std::unique_ptr<AccelerationContext> acceleration; // must be destroyed last for sycl
 
    std::unique_ptr<BaseCanonicalGrid> base;
 
    std::vector<double> domain_transform_a, domain_transform_b;
    std::vector<int> conformal_asin_power;
    std::vector<int> llimits;
 
    bool using_dynamic_construction;
 
    mutable std::unique_ptr<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;
}
 
}
 
#endif