Doxygen 1.17.0
Toolkit for Adaptive Stochastic Modeling and Non-Intrusive ApproximatioN: Tasmanian v8.3
Loading...
Searching...
No Matches
tsgGridWavelet.hpp
1/*
2 * Copyright (c) 2017, Miroslav Stoyanov
3 *
4 * This file is part of
5 * Toolkit for Adaptive Stochastic Modeling And Non-Intrusive ApproximatioN: TASMANIAN
6 *
7 * Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
8 *
9 * 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
10 *
11 * 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions
12 * and the following disclaimer in the documentation and/or other materials provided with the distribution.
13 *
14 * 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse
15 * or promote products derived from this software without specific prior written permission.
16 *
17 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,
18 * INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
19 * IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY,
20 * OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA,
21 * OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
22 * 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.
23 *
24 * UT-BATTELLE, LLC AND THE UNITED STATES GOVERNMENT MAKE NO REPRESENTATIONS AND DISCLAIM ALL WARRANTIES, BOTH EXPRESSED AND IMPLIED.
25 * 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,
26 * 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.
27 * 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,
28 * IN WHOLE OR IN PART THE USE, STORAGE OR DISPOSAL OF THE SOFTWARE.
29 */
30
31#ifndef __TASMANIAN_SPARSE_GRID_WAVELET_HPP
32#define __TASMANIAN_SPARSE_GRID_WAVELET_HPP
33
34#include "tsgRuleWavelet.hpp"
35
36namespace TasGrid{
37
38#ifndef __TASMANIAN_DOXYGEN_SKIP
39class GridWavelet final : public BaseCanonicalGrid{
40public:
41 GridWavelet(AccelerationContext const *acc) : BaseCanonicalGrid(acc), rule1D(1, 10), order(1){}
42 friend struct GridReaderVersion5<GridWavelet>;
43 GridWavelet(AccelerationContext const *acc, const GridWavelet *wav, int ibegin, int iend);
44 GridWavelet(AccelerationContext const *acc, int cnum_dimensions, int cnum_outputs, int depth, int corder, const std::vector<int> &level_limits);
45 GridWavelet(AccelerationContext const *acc, MultiIndexSet &&pset, int cnum_outputs, int corder, Data2D<double> &&vals);
46 ~GridWavelet() = default;
47
48 GridWavelet(GridWavelet &&) = default;
49 GridWavelet &operator = (GridWavelet &&) = default;
50
51 bool isWavelet() const override{ return true; }
52
53 void write(std::ostream &os, bool iomode) const override{ if (iomode == mode_ascii) write<mode_ascii>(os); else write<mode_binary>(os); }
54
55 template<bool iomode> void write(std::ostream &os) const;
56
57 TypeOneDRule getRule() const override{ return rule_wavelet; }
58 int getOrder() const{ return order; }
59
60 void getLoadedPoints(double *x) const override;
61 void getNeededPoints(double *x) const override;
62 void getPoints(double *x) const override; // returns the loaded points unless no points are loaded, then returns the needed points
63
64 void getQuadratureWeights(double weights[]) const override;
65 void getInterpolationWeights(const double x[], double weights[]) const override;
66 void getDifferentiationWeights(const double x[], double weights[]) const override;
67
68 void loadNeededValues(const double *vals) override;
69
70 void evaluate(const double x[], double y[]) const override;
71 void integrate(double q[], double *conformal_correction) const override;
72 void differentiate(const double x[], double jacobian[]) const override;
73
74 void evaluateBatch(const double x[], int num_x, double y[]) const override;
75
76 void evaluateGpuMixed(const double*, int, double[]) const;
77 void evaluateBatchGPU(const double*, int, double[]) const override;
78 void evaluateBatchGPU(const float*, int, float[]) const override;
79 template<typename T> void evaluateBatchGPUtempl(const T*, int, T[]) const;
80 void evaluateHierarchicalFunctionsGPU(const double gpu_x[], int cpu_num_x, double *gpu_y) const override;
81 void evaluateHierarchicalFunctionsGPU(const float gpu_x[], int cpu_num_x, float *gpu_y) const override;
82
83 void setSurplusRefinement(double tolerance, TypeRefinement criteria, int output, const std::vector<int> &level_limits);
84 void clearRefinement() override;
85 void mergeRefinement() override;
86
87 void beginConstruction() override;
88 void writeConstructionData(std::ostream&, bool) const override;
89 void readConstructionData(std::istream&, bool) override;
90 std::vector<double> getCandidateConstructionPoints(double tolerance, TypeRefinement criteria, int output, std::vector<int> const &level_limits);
91 void loadConstructedPoint(const double[], const std::vector<double> &) override;
92 void loadConstructedPoint(const double[], int, const double[]) override;
93 void finishConstruction() override;
94
95 void evaluateHierarchicalFunctions(const double x[], int num_x, double y[]) const override;
96 std::vector<double> getSupport() const override;
97
98 void setHierarchicalCoefficients(const double c[]) override;
99 void integrateHierarchicalFunctions(double integrals[]) const override;
100
101 const double* getSurpluses() const;
102
103 void updateAccelerationData(AccelerationContext::ChangeType change) const override;
104
105protected:
106 double evalBasis(const int p[], const double x[]) const;
107 void buildInterpolationMatrix() const;
108 void recomputeCoefficients();
109 void solveTransposed(double w[]) const;
110 double evalIntegral(const int p[]) const;
111 void evalDiffBasis(const int p[], const double x[], double jacobian[]) const;
112
113 std::vector<double> getNormalization() const;
114 std::vector<int> getMultiIndex(const double x[]);
115
116 Data2D<int> buildUpdateMap(double tolerance, TypeRefinement criteria, int output) const;
117 MultiIndexSet getRefinementCanidates(double tolerance, TypeRefinement criteria, int output, const std::vector<int> &level_limits) const;
118
119 bool addParent(const int point[], int direction, Data2D<int> &destination) const;
120 void addChild(const int point[], int direction, Data2D<int> &destination) const;
121 void addChildLimited(const int point[], int direction, const std::vector<int> &level_limits, Data2D<int> &destination) const;
122
123 void clearGpuCoefficients() const;
124 void clearGpuBasis() const;
125
126private:
127 RuleWavelet rule1D;
128
129 int order;
130
131 Data2D<double> coefficients; // a.k.a., surpluses
132
133 mutable TasSparse::WaveletBasisMatrix inter_matrix;
134
135 std::optional<CudaWaveletData<double>>& getGpuCacheOverload(double) const{ return gpu_cache; }
136 std::optional<CudaWaveletData<float>>& getGpuCacheOverload(float) const{ return gpu_cachef; }
137 template<typename T> std::optional<CudaWaveletData<T>>& getGpuCache() const{
138 return getGpuCacheOverload(static_cast<T>(0.0));
139 }
140 template<typename T> void loadGpuCoefficients() const;
141 template<typename T> void loadGpuBasis() const;
142 mutable std::optional<CudaWaveletData<double>> gpu_cache;
143 mutable std::optional<CudaWaveletData<float>> gpu_cachef;
144
145 std::unique_ptr<SimpleConstructData> dynamic_values;
146};
147
148// Old version reader
149template<> struct GridReaderVersion5<GridWavelet>{
150 template<typename iomode> static void read(std::istream &is, GridWavelet *grid) {
151 grid->num_dimensions = IO::readNumber<iomode, int>(is);
152 grid->num_outputs = IO::readNumber<iomode, int>(is);
153 grid->order = IO::readNumber<iomode, int>(is);
154 grid->rule1D.updateOrder(grid->order);
155
156 if (IO::readFlag<iomode>(is)) grid->points = MultiIndexSet(is, iomode());
157 if (std::is_same<iomode, IO::mode_ascii_type>::value){ // backwards compatible: surpluses and needed, or needed and surpluses
158 if (IO::readFlag<iomode>(is))
159 grid->coefficients = IO::readData2D<iomode, double>(is, grid->num_outputs, grid->points.getNumIndexes());
160 if (IO::readFlag<iomode>(is)) grid->needed = MultiIndexSet(is, iomode());
161 }else{
162 if (IO::readFlag<iomode>(is)) grid->needed = MultiIndexSet(is, iomode());
163 if (IO::readFlag<iomode>(is))
164 grid->coefficients = IO::readData2D<iomode, double>(is, grid->num_outputs, grid->points.getNumIndexes());
165 }
166
167 if (grid->num_outputs > 0) grid->values = StorageSet(is, iomode());
168 grid->buildInterpolationMatrix();
169 }
170};
171#endif // __TASMANIAN_DOXYGEN_SKIP
172
173}
174
175#endif
TasGrid::TypeOneDRule
TypeOneDRule
Used to specify the one dimensional family of rules that induces the sparse grid.
Definition tsgEnumerates.hpp:288
TasGrid::rule_wavelet
@ rule_wavelet
Wavelet basis with uniformly distributed nodes (primarily for internal use).
Definition tsgEnumerates.hpp:373
TasGrid::loadNeededValues
void loadNeededValues(std::function< void(double const x[], double y[], size_t thread_id)> model, TasmanianSparseGrid &grid, size_t num_threads)
Loads the current grid with model values, does not perform any refinement.
Definition tsgLoadNeededValues.hpp:104
TasGrid
Encapsulates the Tasmanian Sparse Grid module.
Definition TasmanianSparseGrid.hpp:68