From 106bdce2ac871e9c5ac91b83e23f89a91bf7e934 Mon Sep 17 00:00:00 2001 From: Miroslav Stoyanov Date: Fri, 22 Mar 2024 17:40:44 -0400 Subject: [PATCH] new templates --- SparseGrids/tsgDConstructGridGlobal.hpp | 15 + SparseGrids/tsgGridLocalPolynomial.cpp | 709 +++++++--- SparseGrids/tsgGridLocalPolynomial.hpp | 143 +- SparseGrids/tsgHierarchyManipulator.cpp | 121 +- SparseGrids/tsgHierarchyManipulator.hpp | 60 +- SparseGrids/tsgRuleLocalPolynomial.hpp | 1602 +++++++++++++---------- 6 files changed, 1716 insertions(+), 934 deletions(-) diff --git a/SparseGrids/tsgDConstructGridGlobal.hpp b/SparseGrids/tsgDConstructGridGlobal.hpp index 5bcacc1b4..c3f4ebaac 100644 --- a/SparseGrids/tsgDConstructGridGlobal.hpp +++ b/SparseGrids/tsgDConstructGridGlobal.hpp @@ -201,6 +201,21 @@ std::vector listToNodes(std::forward_list const &node_list, si ix = MultiIndexManipulations::indexesToNodes(t.point, rule, ix); return result; } +/*! + * \internal + * \ingroup TasmanianRefinement + * \brief Using MultiIndexManipulations::indexesToNodes() convert the \b node_list to actual points according to the rule. + * + * \endinternal + */ +template +std::vector listToLocalNodes(std::forward_list const &node_list, size_t num_dimensions, callable_method rule){ + std::vector result(Utils::size_mult(std::distance(node_list.begin(), node_list.end()), num_dimensions)); + auto ix = result.begin(); + for(auto const &t : node_list) + ix = std::transform(t.point.begin(), t.point.end(), ix, rule); + return result; +} /*! * \internal diff --git a/SparseGrids/tsgGridLocalPolynomial.cpp b/SparseGrids/tsgGridLocalPolynomial.cpp index e01e0d3f0..c3ac21163 100644 --- a/SparseGrids/tsgGridLocalPolynomial.cpp +++ b/SparseGrids/tsgGridLocalPolynomial.cpp @@ -39,7 +39,7 @@ namespace TasGrid{ template void GridLocalPolynomial::write(std::ostream &os) const{ if (iomode == mode_ascii){ os << std::scientific; os.precision(17); } IO::writeNumbers(os, num_dimensions, num_outputs, order, top_level); - IO::writeRule(rule->getType(), os); + IO::writeRule(RuleLocal::getRule(effective_rule), os); IO::writeFlag(!points.empty(), os); if (!points.empty()) points.write(os); if (iomode == mode_ascii){ // backwards compatible: surpluses and needed, or needed and surpluses @@ -72,19 +72,40 @@ template void GridLocalPolynomial::write(std::ostream &) const; GridLocalPolynomial::GridLocalPolynomial(AccelerationContext const *acc, int cnum_dimensions, int cnum_outputs, int depth, int corder, TypeOneDRule crule, const std::vector &level_limits) : BaseCanonicalGrid(acc, cnum_dimensions, cnum_outputs, MultiIndexSet(), MultiIndexSet(), StorageSet()), order(corder), - rule(makeRuleLocalPolynomial(((crule == rule_semilocalp) && (order < 2)) ? rule_localp : crule, corder)) - { + effective_rule(RuleLocal::getEffectiveRule(order, (crule == rule_semilocalp and order < 2) ? rule_localp : crule)) + { MultiIndexSet tensors = MultiIndexManipulations::selectTensors((size_t) num_dimensions, depth, type_level, [&](int i) -> int{ return i; }, std::vector(), level_limits); - needed = MultiIndexManipulations::generateNestedPoints(tensors, [&](int l) -> int{ return rule->getNumPoints(l); }); + switch(effective_rule) { + case RuleLocal::erule::pwc: + needed = MultiIndexManipulations::generateNestedPoints(tensors, + [&](int l) -> int{ return RuleLocal::getNumPoints(l); }); + break; + case RuleLocal::erule::localp: + needed = MultiIndexManipulations::generateNestedPoints(tensors, + [&](int l) -> int{ return RuleLocal::getNumPoints(l); }); + break; + case RuleLocal::erule::semilocalp: + needed = MultiIndexManipulations::generateNestedPoints(tensors, + [&](int l) -> int{ return RuleLocal::getNumPoints(l); }); + break; + case RuleLocal::erule::localp0: + needed = MultiIndexManipulations::generateNestedPoints(tensors, + [&](int l) -> int{ return RuleLocal::getNumPoints(l); }); + break; + default: // case RuleLocal::erule::localpb: + needed = MultiIndexManipulations::generateNestedPoints(tensors, + [&](int l) -> int{ return RuleLocal::getNumPoints(l); }); + break; + }; buildTree(); if (num_outputs == 0){ points = std::move(needed); needed = MultiIndexSet(); - parents = HierarchyManipulations::computeDAGup(points, rule.get()); + parents = HierarchyManipulations::computeDAGup(points, effective_rule); }else{ values.resize(num_outputs, needed.getNumIndexes()); } @@ -99,7 +120,7 @@ GridLocalPolynomial::GridLocalPolynomial(AccelerationContext const *acc, GridLoc roots(pwpoly->roots), pntr(pwpoly->pntr), indx(pwpoly->indx), - rule(makeRuleLocalPolynomial(pwpoly->rule->getType(), pwpoly->order)){ + effective_rule(pwpoly->effective_rule){ if (pwpoly->dynamic_values){ dynamic_values = Utils::make_unique(*pwpoly->dynamic_values); @@ -113,24 +134,65 @@ GridLocalPolynomial::GridLocalPolynomial(AccelerationContext const *acc, int cnu StorageSet(cnum_outputs, static_cast(vals.size() / cnum_outputs), std::move(vals))), order(corder), surpluses(Data2D(cnum_outputs, points.getNumIndexes(), std::move(surps))), - rule(makeRuleLocalPolynomial(crule, corder)){ + effective_rule(RuleLocal::getEffectiveRule(order, crule)){ buildTree(); } -void GridLocalPolynomial::getLoadedPoints(double *x) const{ - int num_points = points.getNumIndexes(); +struct localp_loaded{}; +struct localp_needed{}; + +template +void GridLocalPolynomial::getPoints(double *x) const { + int num_points = (std::is_same::value) ? points.getNumIndexes() : needed.getNumIndexes(); Utils::Wrapper2D split(num_dimensions, x); #pragma omp parallel for schedule(static) - for(int i=0; i::value) ? points.getIndex(i) : needed.getIndex(i); + double *s = split.getStrip(i); + for(int j=0; j(p[j]); + } +} +void GridLocalPolynomial::getLoadedPoints(double *x) const{ + using pmode = localp_loaded; + switch(effective_rule) { + case RuleLocal::erule::pwc: + getPoints(x); + break; + case RuleLocal::erule::localp: + getPoints(x); + break; + case RuleLocal::erule::semilocalp: + getPoints(x); + break; + case RuleLocal::erule::localp0: + getPoints(x); + break; + default: // case RuleLocal::erule::localpb: + getPoints(x); + break; + }; } void GridLocalPolynomial::getNeededPoints(double *x) const{ - int num_points = needed.getNumIndexes(); - Utils::Wrapper2D split(num_dimensions, x); - #pragma omp parallel for schedule(static) - for(int i=0; i(x); + break; + case RuleLocal::erule::localp: + getPoints(x); + break; + case RuleLocal::erule::semilocalp: + getPoints(x); + break; + case RuleLocal::erule::localp0: + getPoints(x); + break; + default: // case RuleLocal::erule::localpb: + getPoints(x); + break; + }; } void GridLocalPolynomial::getPoints(double *x) const{ if (points.empty()){ getNeededPoints(x); }else{ getLoadedPoints(x); } @@ -220,7 +282,7 @@ void GridLocalPolynomial::evaluateBatch(const double x[], int num_x, double y[]) void GridLocalPolynomial::loadNeededValuesGPU(const double *vals){ updateValues(vals); - std::vector levels = HierarchyManipulations::computeLevels(points, rule.get()); + std::vector levels = HierarchyManipulations::computeLevels(points, effective_rule); std::vector> lpnts = HierarchyManipulations::splitByLevels(points, levels); std::vector> lvals = HierarchyManipulations::splitByLevels(values, levels); @@ -308,23 +370,23 @@ void GridLocalPolynomial::evaluateBatchGPU(const float gpu_x[], int cpu_num_x, f } void GridLocalPolynomial::evaluateHierarchicalFunctionsGPU(const double gpu_x[], int cpu_num_x, double *gpu_y) const{ loadGpuBasis(); - TasGpu::devalpwpoly(acceleration, order, rule->getType(), num_dimensions, cpu_num_x, getNumPoints(), gpu_x, gpu_cache->nodes.data(), gpu_cache->support.data(), gpu_y); + TasGpu::devalpwpoly(acceleration, order, RuleLocal::getRule(effective_rule), num_dimensions, cpu_num_x, getNumPoints(), gpu_x, gpu_cache->nodes.data(), gpu_cache->support.data(), gpu_y); } void GridLocalPolynomial::buildSparseBasisMatrixGPU(const double gpu_x[], int cpu_num_x, GpuVector &gpu_spntr, GpuVector &gpu_sindx, GpuVector &gpu_svals) const{ loadGpuBasis(); loadGpuHierarchy(); - TasGpu::devalpwpoly_sparse(acceleration, order, rule->getType(), num_dimensions, cpu_num_x, gpu_x, + TasGpu::devalpwpoly_sparse(acceleration, order, RuleLocal::getRule(effective_rule), num_dimensions, cpu_num_x, gpu_x, gpu_cache->nodes, gpu_cache->support, gpu_cache->hpntr, gpu_cache->hindx, gpu_cache->hroots, gpu_spntr, gpu_sindx, gpu_svals); } void GridLocalPolynomial::evaluateHierarchicalFunctionsGPU(const float gpu_x[], int cpu_num_x, float *gpu_y) const{ loadGpuBasis(); - TasGpu::devalpwpoly(acceleration, order, rule->getType(), num_dimensions, cpu_num_x, getNumPoints(), gpu_x, gpu_cachef->nodes.data(), gpu_cachef->support.data(), gpu_y); + TasGpu::devalpwpoly(acceleration, order, RuleLocal::getRule(effective_rule), num_dimensions, cpu_num_x, getNumPoints(), gpu_x, gpu_cachef->nodes.data(), gpu_cachef->support.data(), gpu_y); } void GridLocalPolynomial::buildSparseBasisMatrixGPU(const float gpu_x[], int cpu_num_x, GpuVector &gpu_spntr, GpuVector &gpu_sindx, GpuVector &gpu_svals) const{ loadGpuBasis(); loadGpuHierarchy(); - TasGpu::devalpwpoly_sparse(acceleration, order, rule->getType(), num_dimensions, cpu_num_x, gpu_x, + TasGpu::devalpwpoly_sparse(acceleration, order, RuleLocal::getRule(effective_rule), num_dimensions, cpu_num_x, gpu_x, gpu_cachef->nodes, gpu_cachef->support, gpu_cachef->hpntr, gpu_cachef->hindx, gpu_cachef->hroots, gpu_spntr, gpu_sindx, gpu_svals); } @@ -339,28 +401,22 @@ template void GridLocalPolynomial::loadGpuBasis() const{ Data2D cpu_support = [&](void)->Data2D{ const MultiIndexSet &work = (points.empty()) ? needed : points; - if (rule->getType() == rule_localp){ - switch(order){ - case 0: return encodeSupportForGPU<0, rule_localp, T>(work); - case 2: return encodeSupportForGPU<2, rule_localp, T>(work); - default: - return encodeSupportForGPU<1, rule_localp, T>(work); - } - }else if (rule->getType() == rule_semilocalp){ - return encodeSupportForGPU<2, rule_semilocalp, T>(work); - }else if (rule->getType() == rule_localpb){ - switch(order){ - case 2: return encodeSupportForGPU<2, rule_localpb, T>(work); - default: - return encodeSupportForGPU<1, rule_localpb, T>(work); - } - }else{ - switch(order){ - case 2: return encodeSupportForGPU<2, rule_localp0, T>(work); - default: - return encodeSupportForGPU<1, rule_localp0, T>(work); - } - } + switch(effective_rule) { + case RuleLocal::erule::pwc: + return encodeSupportForGPU<0, rule_localp, T>(work); + case RuleLocal::erule::localp: + return (order == 1) ? encodeSupportForGPU<1, rule_localp, T>(work) + : encodeSupportForGPU<2, rule_localp, T>(work); + case RuleLocal::erule::semilocalp: + return (order == 1) ? encodeSupportForGPU<1, rule_semilocalp, T>(work) + : encodeSupportForGPU<2, rule_semilocalp, T>(work); + case RuleLocal::erule::localp0: + return (order == 1) ? encodeSupportForGPU<1, rule_localp0, T>(work) + : encodeSupportForGPU<2, rule_localp0, T>(work); + default: // case RuleLocal::erule::localpb: + return (order == 1) ? encodeSupportForGPU<1, rule_localpb, T>(work) + : encodeSupportForGPU<2, rule_localpb, T>(work); + }; }(); ccache->support.load(acceleration, cpu_support.begin(), cpu_support.end()); } @@ -452,10 +508,11 @@ void GridLocalPolynomial::readConstructionData(std::istream &is, bool iomode){ else dynamic_values = Utils::make_unique(is, num_dimensions, num_outputs, IO::mode_binary_type()); } +template std::vector GridLocalPolynomial::getCandidateConstructionPoints(double tolerance, TypeRefinement criteria, int output, std::vector const &level_limits, double const *scale_correction){ // combine the initial points with negative weights and the refinement candidates with surplus weights (no need to normalize, the sort uses relative values) - MultiIndexSet refine_candidates = getRefinementCanidates(tolerance, criteria, output, level_limits, scale_correction); + MultiIndexSet refine_candidates = getRefinementCanidates(tolerance, criteria, output, level_limits, scale_correction); MultiIndexSet new_points = (dynamic_values->initial_points.empty()) ? std::move(refine_candidates) : refine_candidates - dynamic_values->initial_points; // compute the weights for the new_points points @@ -488,14 +545,14 @@ std::vector GridLocalPolynomial::getCandidateConstructionPoints(double t double weight = 0.0; std::vector p = new_points.copyIndex(i); - HierarchyManipulations::touchAllImmediateRelatives(p, points, rule.get(), - [&](int relative)->void{ weight = std::max(weight, getDominantSurplus(relative)); }); + HierarchyManipulations::touchAllImmediateRelatives(p, points, + [&](int relative)->void{ weight = std::max(weight, getDominantSurplus(relative)); }); refine_weights[i] = weight; // those will be inverted } // if using stable refinement, ensure the weight of the parents is never less than the children if (!new_points.empty() && ((criteria == refine_parents_first) || (criteria == refine_fds))){ - auto rlevels = HierarchyManipulations::computeLevels(new_points, rule.get()); + auto rlevels = HierarchyManipulations::computeLevels(new_points); auto split = HierarchyManipulations::splitByLevels(new_points, rlevels); for(auto is = split.rbegin(); is != split.rend(); is++){ for(int i=0; igetNumStrips(); i++){ @@ -503,10 +560,10 @@ std::vector GridLocalPolynomial::getCandidateConstructionPoints(double t double correction = refine_weights[new_points.getSlot(parent)]; // will never be missing for(auto &p : parent){ int r = p; - p = rule->getParent(r); + p = RuleLocal::getParent(r); int ip = (p == -1) ? -1 : new_points.getSlot(parent); // if parent is among the refined if (ip != -1) refine_weights[ip] += correction; - p = rule->getStepParent(r); + p = RuleLocal::getStepParent(r); ip = (p == -1) ? -1 : new_points.getSlot(parent); // if parent is among the refined if (ip != -1) refine_weights[ip] += correction; p = r; @@ -515,7 +572,7 @@ std::vector GridLocalPolynomial::getCandidateConstructionPoints(double t } }else if (!new_points.empty() && (criteria == refine_stable)){ // stable refinement, ensure that if level[i] < level[j] then weight[i] > weight[j] - auto rlevels = HierarchyManipulations::computeLevels(new_points, rule.get()); + auto rlevels = HierarchyManipulations::computeLevels(new_points); auto split = HierarchyManipulations::splitByLevels(new_points, rlevels); double max_weight = 0.0; for(auto is = split.rbegin(); is != split.rend(); is++){ // loop backwards in levels @@ -529,7 +586,7 @@ std::vector GridLocalPolynomial::getCandidateConstructionPoints(double t } // compute the weights for the initial points - std::vector initial_levels = HierarchyManipulations::computeLevels(dynamic_values->initial_points, rule.get()); + std::vector initial_levels = HierarchyManipulations::computeLevels(dynamic_values->initial_points); std::forward_list weighted_points; for(int i=0; iinitial_points.getNumIndexes(); i++) @@ -540,35 +597,68 @@ std::vector GridLocalPolynomial::getCandidateConstructionPoints(double t // sort and return the sorted list weighted_points.sort([&](const NodeData &a, const NodeData &b)->bool{ return (a.value[0] < b.value[0]); }); - return listToNodes(weighted_points, num_dimensions, *rule); + return listToLocalNodes(weighted_points, num_dimensions, [&](int i)->double{ return RuleLocal::getNode(i); }); } +std::vector GridLocalPolynomial::getCandidateConstructionPoints(double tolerance, TypeRefinement criteria, int output, + std::vector const &level_limits, double const *scale_correction) { + switch(effective_rule) { + case RuleLocal::erule::pwc: + return getCandidateConstructionPoints(tolerance, criteria, output, level_limits, scale_correction); + case RuleLocal::erule::localp: + return getCandidateConstructionPoints(tolerance, criteria, output, level_limits, scale_correction); + case RuleLocal::erule::semilocalp: + return getCandidateConstructionPoints(tolerance, criteria, output, level_limits, scale_correction); + case RuleLocal::erule::localp0: + return getCandidateConstructionPoints(tolerance, criteria, output, level_limits, scale_correction); + default: // case RuleLocal::erule::localpb: + return getCandidateConstructionPoints(tolerance, criteria, output, level_limits, scale_correction); + }; +} + +template std::vector GridLocalPolynomial::getMultiIndex(const double x[]){ std::vector p(num_dimensions); // convert x to p, maybe expensive for(int j=0; jgetNode(i) - x[j]) > Maths::num_tol) i++; + while(std::abs(RuleLocal::getNode(i) - x[j]) > Maths::num_tol) i++; p[j] = i; } return p; } + +template void GridLocalPolynomial::loadConstructedPoint(const double x[], const std::vector &y){ - auto p = getMultiIndex(x); + auto p = getMultiIndex(x); dynamic_values->initial_points.removeIndex(p); bool isConnected = false; - HierarchyManipulations::touchAllImmediateRelatives(p, points, rule.get(), - [&](int)->void{ isConnected = true; }); - int lvl = rule->getLevel(p[0]); - for(int j=1; jgetLevel(p[j]); + HierarchyManipulations::touchAllImmediateRelatives(p, points, [&](int)->void{ isConnected = true; }); + int lvl = RuleLocal::getLevel(p[0]); + for(int j=1; j(p[j]); if (isConnected || (lvl == 0)){ - expandGrid(p, y); - loadConstructedPoints(); + expandGrid(p, y); + loadConstructedPoints(); }else{ dynamic_values->data.push_front({p, y}); } } +void GridLocalPolynomial::loadConstructedPoint(const double x[], const std::vector &y){ + switch(effective_rule) { + case RuleLocal::erule::pwc: + return loadConstructedPoint(x, y); + case RuleLocal::erule::localp: + return loadConstructedPoint(x, y); + case RuleLocal::erule::semilocalp: + return loadConstructedPoint(x, y); + case RuleLocal::erule::localp0: + return loadConstructedPoint(x, y); + default: // case RuleLocal::erule::localpb: + return loadConstructedPoint(x, y); + }; +} +template void GridLocalPolynomial::expandGrid(const std::vector &point, const std::vector &value){ if (points.empty()){ // only one point points = MultiIndexSet((size_t) num_dimensions, std::vector(point)); @@ -576,12 +666,15 @@ void GridLocalPolynomial::expandGrid(const std::vector &point, const std::v surpluses = Data2D(num_outputs, 1, std::vector(value)); // one value is its own surplus }else{ // merge with existing points // compute the surplus for the point - std::vector xnode = MultiIndexManipulations::getIndexesToNodes(point, *rule); + std::vector xnode(num_dimensions); + for(int j=0; j(point[j]); + std::vector approximation(num_outputs), surp(num_outputs); evaluate(xnode.data(), approximation.data()); std::transform(approximation.begin(), approximation.end(), value.begin(), surp.begin(), [&](double e, double v)->double{ return v - e; }); - std::vector graph = getSubGraph(point); // get the descendant nodes that must be updated later + std::vector graph = getSubGraph(point); // get the descendant nodes that must be updated later values.addValues(points, MultiIndexSet(num_dimensions, std::vector(point)), value.data()); // added the value @@ -594,24 +687,25 @@ void GridLocalPolynomial::expandGrid(const std::vector &point, const std::v std::vector levels(points.getNumIndexes(), 0); // compute the levels, but only for the new indexes for(auto &g : graph){ int const *pnt = points.getIndex(g); - int l = rule->getLevel(pnt[0]); - for(int j=1; jgetLevel(pnt[j]); + int l = RuleLocal::getLevel(pnt[0]); + for(int j=1; j(pnt[j]); levels[g] = l; std::copy_n(values.getValues(g), num_outputs, surpluses.getStrip(g)); // reset the surpluses to the values (will be updated) } // compute the current DAG and update the surplused for the descendants - updateSurpluses(points, top_level + 1, levels, HierarchyManipulations::computeDAGup(points, rule.get())); + updateSurpluses(points, top_level + 1, levels, HierarchyManipulations::computeDAGup(points)); } buildTree(); // the tree is needed for evaluate(), must be rebuild every time the points set is updated } +template void GridLocalPolynomial::loadConstructedPoint(const double x[], int numx, const double y[]){ Utils::Wrapper2D wrapx(num_dimensions, x); std::vector> pnts(numx); #pragma omp parallel for for(int i=0; i(wrapx.getStrip(i)); if (!dynamic_values->initial_points.empty()){ Data2D combined_pnts(num_dimensions, numx); @@ -624,15 +718,30 @@ void GridLocalPolynomial::loadConstructedPoint(const double x[], int numx, const for(int i=0; idata.push_front({std::move(pnts[i]), std::vector(wrapy.getStrip(i), wrapy.getStrip(i) + num_outputs)}); - loadConstructedPoints(); + loadConstructedPoints(); +} +void GridLocalPolynomial::loadConstructedPoint(const double x[], int numx, const double y[]){ + switch(effective_rule) { + case RuleLocal::erule::pwc: + return loadConstructedPoint(x, numx, y); + case RuleLocal::erule::localp: + return loadConstructedPoint(x, numx, y); + case RuleLocal::erule::semilocalp: + return loadConstructedPoint(x, numx, y); + case RuleLocal::erule::localp0: + return loadConstructedPoint(x, numx, y); + default: // case RuleLocal::erule::localpb: + return loadConstructedPoint(x, numx, y); + }; } +template void GridLocalPolynomial::loadConstructedPoints(){ Data2D candidates(num_dimensions, (int) std::distance(dynamic_values->data.begin(), dynamic_values->data.end())); for(struct{ int i; std::forward_list::iterator d; } p = {0, dynamic_values->data.begin()}; p.d != dynamic_values->data.end(); p.i++, p.d++){ std::copy_n(p.d->point.begin(), num_dimensions, candidates.getIStrip(p.i)); } - auto new_points = HierarchyManipulations::getLargestConnected(points, MultiIndexSet(candidates), rule.get()); + auto new_points = HierarchyManipulations::getLargestConnected(points, MultiIndexSet(candidates)); if (new_points.empty()) return; clearGpuBasisHierarchy(); // the points will change, clear the cache @@ -647,14 +756,15 @@ void GridLocalPolynomial::loadConstructedPoints(){ points += new_points; } buildTree(); - recomputeSurpluses(); // costly, but the only option under the circumstances + recomputeSurpluses(); // costly, but the only option under the circumstances } void GridLocalPolynomial::finishConstruction(){ dynamic_values.reset(); } +template std::vector GridLocalPolynomial::getSubGraph(std::vector const &point) const{ std::vector graph, p = point; std::vector used(points.getNumIndexes(), false); - int max_1d_kids = rule->getMaxNumKids(); + int max_1d_kids = RuleLocal::getMaxNumKids(); int max_kids = max_1d_kids * num_dimensions; std::vector monkey_count(1, 0), monkey_tail; @@ -663,7 +773,7 @@ std::vector GridLocalPolynomial::getSubGraph(std::vector const &point) if (monkey_count.back() < max_kids){ int dim = monkey_count.back() / max_1d_kids; monkey_tail.push_back(p[dim]); - p[dim] = rule->getKid(monkey_tail.back(), monkey_count.back() % max_1d_kids); + p[dim] = RuleLocal::getKid(monkey_tail.back(), monkey_count.back() % max_1d_kids); int slot = points.getSlot(p); if ((slot == -1) || used[slot]){ // this kid is missing p[dim] = monkey_tail.back(); @@ -721,17 +831,62 @@ void GridLocalPolynomial::evaluateHierarchicalFunctions(const double x[], int nu int num_points = work.getNumIndexes(); Utils::Wrapper2D xwrap(num_dimensions, x); Utils::Wrapper2D ywrap(num_points, y); - #pragma omp parallel for - for(int i=0; i(work.getIndex(j), this_x, dummy); + } + break; + case RuleLocal::erule::localp: + #pragma omp parallel for + for(int i=0; i(work.getIndex(j), this_x, dummy); + } + break; + case RuleLocal::erule::semilocalp: + #pragma omp parallel for + for(int i=0; i(work.getIndex(j), this_x, dummy); + } + break; + case RuleLocal::erule::localp0: + #pragma omp parallel for + for(int i=0; i(work.getIndex(j), this_x, dummy); + } + break; + default: // case RuleLocal::erule::localpb: + #pragma omp parallel for + for(int i=0; i(work.getIndex(j), this_x, dummy); + } + break; + }; } -void GridLocalPolynomial::recomputeSurpluses(){ +template +void GridLocalPolynomial::recomputeSurpluses() { surpluses = Data2D(num_outputs, points.getNumIndexes(), std::vector(values.begin(), values.end())); // There are two available algorithms here: @@ -763,19 +918,19 @@ void GridLocalPolynomial::recomputeSurpluses(){ // higher dimensions will favor (kron) even more if (num_dimensions <= 2 or (num_dimensions == 3 and points.getNumIndexes() > 2000000)) { - Data2D dagUp = HierarchyManipulations::computeDAGup(points, rule.get()); - std::vector level = HierarchyManipulations::computeLevels(points, rule.get()); - updateSurpluses(points, top_level, level, dagUp); + Data2D dagUp = HierarchyManipulations::computeDAGup(points); + std::vector level = HierarchyManipulations::computeLevels(points); + updateSurpluses(points, top_level, level, dagUp); return; } bool is_complete = true; - Data2D dagUp = HierarchyManipulations::computeDAGup(points, rule.get(), is_complete); + Data2D dagUp = HierarchyManipulations::computeDAGup(points, is_complete); if (not is_complete) { // incomplete hierarchy, must use the slow algorithm - std::vector level = HierarchyManipulations::computeLevels(points, rule.get()); - updateSurpluses(points, top_level, level, dagUp); + std::vector level = HierarchyManipulations::computeLevels(points); + updateSurpluses(points, top_level, level, dagUp); return; } @@ -783,7 +938,7 @@ void GridLocalPolynomial::recomputeSurpluses(){ std::vector vpntr, vindx; std::vector vvals; - rule->van_matrix(num_nodes, vpntr, vindx, vvals); + RuleLocal::van_matrix(order, num_nodes, vpntr, vindx, vvals); std::vector> map; std::vector> lines1d; @@ -823,9 +978,30 @@ void GridLocalPolynomial::recomputeSurpluses(){ } } +void GridLocalPolynomial::recomputeSurpluses() { + switch(effective_rule) { + case RuleLocal::erule::pwc: + recomputeSurpluses(); + break; + case RuleLocal::erule::localp: + recomputeSurpluses(); + break; + case RuleLocal::erule::semilocalp: + recomputeSurpluses(); + break; + case RuleLocal::erule::localp0: + recomputeSurpluses(); + break; + default: // case RuleLocal::erule::localpb: + recomputeSurpluses(); + break; + }; +} + +template void GridLocalPolynomial::updateSurpluses(MultiIndexSet const &work, int max_level, std::vector const &level, Data2D const &dagUp){ int num_points = work.getNumIndexes(); - int max_parents = num_dimensions * rule->getMaxNumParents(); + int max_parents = num_dimensions * RuleLocal::getMaxNumParents(); std::vector> indexses_for_levels((size_t) max_level+1); for(int i=0; i x = MultiIndexManipulations::getIndexesToNodes(work.getIndex(i), num_dimensions, *rule); - double *surpi = surpluses.getStrip(i); + #pragma omp parallel + { + std::vector x(num_dimensions); std::vector monkey_count(max_level + 1); std::vector monkey_tail(max_level + 1); std::vector used(num_points, false); - int current = 0; + #pragma omp for schedule(dynamic) + for(int s=0; s(pnt[j]); - while(monkey_count[0] < max_parents){ - if (monkey_count[current] < max_parents){ - int branch = dagUp.getStrip(monkey_tail[current])[monkey_count[current]]; - if ((branch == -1) || (used[branch])){ - monkey_count[current]++; - }else{ - const double *branch_surp = surpluses.getStrip(branch); - double basis_value = evalBasisRaw(work.getIndex(branch), x.data()); - for(int k=0; k(work.getIndex(branch), x.data()); + for(int k=0; k +template void GridLocalPolynomial::applyTransformationTransposed(double weights[], const MultiIndexSet &work, const std::vector &active_points) const { Data2D lparents = (parents.getNumStrips() != work.getNumIndexes()) ? // if the current dag loaded in parents does not reflect the indexes in work - HierarchyManipulations::computeDAGup(work, rule.get()) : + HierarchyManipulations::computeDAGup(work) : Data2D(); const Data2D &dagUp = (parents.getNumStrips() != work.getNumIndexes()) ? lparents : parents; @@ -884,9 +1070,9 @@ void GridLocalPolynomial::applyTransformationTransposed(double weights[], const int active_top_level = 0; for(size_t i=0; igetLevel(p[0]); + int current_level = RuleLocal::getLevel(p[0]); for(int j=1; jgetLevel(p[j]); + current_level += RuleLocal::getLevel(p[j]); } if (active_top_level < current_level) active_top_level = current_level; level[i] = current_level; @@ -895,12 +1081,16 @@ void GridLocalPolynomial::applyTransformationTransposed(double weights[], const std::vector monkey_count(top_level+1); std::vector monkey_tail(top_level+1); std::vector used(work.getNumIndexes()); - int max_parents = rule->getMaxNumParents() * num_dimensions; + int max_parents = RuleLocal::getMaxNumParents() * num_dimensions; + + std::vector node(num_dimensions); for(int l=active_top_level; l>0; l--){ for(size_t i=0; i node = MultiIndexManipulations::getIndexesToNodes(work.getIndex(active_points[i]), num_dimensions, *rule); + int const *pnt = work.getIndex(active_points[i]); + for(int j=0; j(pnt[j]); std::fill(used.begin(), used.end(), false); @@ -915,8 +1105,8 @@ void GridLocalPolynomial::applyTransformationTransposed(double weights[], const monkey_count[current]++; }else{ const int *func = work.getIndex(branch); - double basis_value = rule->evalRaw(func[0], node[0]); - for(int j=1; jevalRaw(func[j], node[j]); + double basis_value = RuleLocal::evalRaw(order, func[0], node[0]); + for(int j=1; j(order, func[j], node[j]); if (mode == 0) { weights[branch] -= weights[active_points[i]] * basis_value; @@ -939,35 +1129,25 @@ void GridLocalPolynomial::applyTransformationTransposed(double weights[], const } } -double GridLocalPolynomial::evalBasisRaw(const int point[], const double x[]) const{ - double f = rule->evalRaw(point[0], x[0]); - for(int j=1; jevalRaw(point[j], x[j]); - return f; -} -double GridLocalPolynomial::evalBasisSupported(const int point[], const double x[], bool &isSupported) const{ - double f = rule->evalSupport(point[0], x[0], isSupported); - if (!isSupported) return 0.0; - for(int j=1; jevalSupport(point[j], x[j], isSupported); - if (!isSupported) return 0.0; - } - return f; -} - -void GridLocalPolynomial::diffBasisSupported(const int point[], const double x[], double diff_values[], bool &isSupported) const{ - isSupported = false; - for(int i=0; ievalSupport(point[k], x[k], isDimSupported); - isSupported = isDimSupported or isSupported; - for(int j=0; jdiffSupport(point[k], x[k], isDimSupported); - isSupported = isDimSupported or isSupported; - } +template +void GridLocalPolynomial::applyTransformationTransposed(double weights[], const MultiIndexSet &work, const std::vector &active_points) const { + switch(effective_rule) { + case RuleLocal::erule::pwc: + applyTransformationTransposed(weights, work, active_points); + break; + case RuleLocal::erule::localp: + applyTransformationTransposed(weights, work, active_points); + break; + case RuleLocal::erule::semilocalp: + applyTransformationTransposed(weights, work, active_points); + break; + case RuleLocal::erule::localp0: + applyTransformationTransposed(weights, work, active_points); + break; + default: // case RuleLocal::erule::localpb: + applyTransformationTransposed(weights, work, active_points); + break; + }; } void GridLocalPolynomial::buildSpareBasisMatrix(const double x[], int num_x, int num_chunk, std::vector &spntr, std::vector &sindx, std::vector &svals) const{ @@ -1058,10 +1238,32 @@ void GridLocalPolynomial::buildTree(){ const MultiIndexSet &work = (points.empty()) ? needed : points; int num_points = work.getNumIndexes(); - std::vector level = HierarchyManipulations::computeLevels(work, rule.get()); - top_level = *std::max_element(level.begin(), level.end()); + Data2D kids; + std::vector level; - Data2D kids = HierarchyManipulations::computeDAGDown(work, rule.get()); + switch(effective_rule) { + case RuleLocal::erule::pwc: + kids = HierarchyManipulations::computeDAGDown(work); + level = HierarchyManipulations::computeLevels(work); + break; + case RuleLocal::erule::localp: + kids = HierarchyManipulations::computeDAGDown(work); + level = HierarchyManipulations::computeLevels(work); + break; + case RuleLocal::erule::semilocalp: + kids = HierarchyManipulations::computeDAGDown(work); + level = HierarchyManipulations::computeLevels(work); + break; + case RuleLocal::erule::localp0: + kids = HierarchyManipulations::computeDAGDown(work); + level = HierarchyManipulations::computeLevels(work); + break; + default: // case RuleLocal::erule::localpb: + kids = HierarchyManipulations::computeDAGDown(work); + level = HierarchyManipulations::computeLevels(work); + break; + }; + top_level = *std::max_element(level.begin(), level.end()); int max_kids = (int) kids.getStride(); @@ -1120,16 +1322,58 @@ void GridLocalPolynomial::integrateHierarchicalFunctions(double integrals[]) con const MultiIndexSet &work = (points.empty()) ? needed : points; std::vector w, x; - if ((rule->getMaxOrder() == -1) || (rule->getMaxOrder() > 3) ) - OneDimensionalNodes::getGaussLegendre(((rule->getMaxOrder() == -1) ? top_level : rule->getMaxOrder()) / 2 + 1, w, x); - - for(int i=0; igetArea(p[0], w, x); - for(int j=1; jgetArea(p[j], w, x); - } + if (order == -1 or order > 3) + OneDimensionalNodes::getGaussLegendre(((order == -1) ? top_level : order) / 2 + 1, w, x); + + switch(effective_rule) { + case RuleLocal::erule::pwc: + #pragma omp parallel for + for(int i=0; i(order, p[0], w, x); + for(int j=1; j(order, p[j], w, x); + } + break; + case RuleLocal::erule::localp: + #pragma omp parallel for + for(int i=0; i(order, p[0], w, x); + for(int j=1; j(order, p[j], w, x); + } + break; + case RuleLocal::erule::semilocalp: + #pragma omp parallel for + for(int i=0; i(order, p[0], w, x); + for(int j=1; j(order, p[j], w, x); + } + break; + case RuleLocal::erule::localp0: + #pragma omp parallel for + for(int i=0; i(order, p[0], w, x); + for(int j=1; j(order, p[j], w, x); + } + break; + default: // case RuleLocal::erule::localpb: + #pragma omp parallel for + for(int i=0; i(order, p[0], w, x); + for(int j=1; j(order, p[j], w, x); + } + break; + }; } +template void GridLocalPolynomial::getQuadratureWeights(double *weights) const{ const MultiIndexSet &work = (points.empty()) ? needed : points; integrateHierarchicalFunctions(weights); @@ -1141,21 +1385,26 @@ void GridLocalPolynomial::getQuadratureWeights(double *weights) const{ Data2D lparents; if (parents.getNumStrips() != work.getNumIndexes()) - lparents = HierarchyManipulations::computeDAGup(work, rule.get()); + lparents = HierarchyManipulations::computeDAGup(work); const Data2D &dagUp = (parents.getNumStrips() != work.getNumIndexes()) ? lparents : parents; int num_points = work.getNumIndexes(); - std::vector level = HierarchyManipulations::computeLevels(work, rule.get()); + std::vector level = HierarchyManipulations::computeLevels(work); + + int max_parents = RuleLocal::getMaxNumParents() * num_dimensions; - int max_parents = rule->getMaxNumParents() * num_dimensions; + std::vector node(num_dimensions); + std::vector used(work.getNumIndexes(), false); for(int l=top_level; l>0; l--){ for(int i=0; i node = MultiIndexManipulations::getIndexesToNodes(work.getIndex(i), num_dimensions, *rule); + int const *pnt = work.getIndex(i); + for(int j=0; j(pnt[j]); - std::vector used(work.getNumIndexes(), false); + std::fill(used.begin(), used.end(), false); monkey_count[0] = 0; monkey_tail[0] = i; @@ -1168,8 +1417,8 @@ void GridLocalPolynomial::getQuadratureWeights(double *weights) const{ monkey_count[current]++; }else{ const int *func = work.getIndex(branch); - basis_value = rule->evalRaw(func[0], node[0]); - for(int j=1; jevalRaw(func[j], node[j]); + basis_value = RuleLocal::evalRaw(order, func[0], node[0]); + for(int j=1; j(order, func[j], node[j]); weights[branch] -= weights[i] * basis_value; used[branch] = true; @@ -1184,6 +1433,25 @@ void GridLocalPolynomial::getQuadratureWeights(double *weights) const{ } } } +void GridLocalPolynomial::getQuadratureWeights(double *weights) const{ + switch(effective_rule) { + case RuleLocal::erule::pwc: + getQuadratureWeights(weights); + break; + case RuleLocal::erule::localp: + getQuadratureWeights(weights); + break; + case RuleLocal::erule::semilocalp: + getQuadratureWeights(weights); + break; + case RuleLocal::erule::localp0: + getQuadratureWeights(weights); + break; + default: // case RuleLocal::erule::localpb: + getQuadratureWeights(weights); + break; + }; +} void GridLocalPolynomial::integrate(double q[], double *conformal_correction) const{ int num_points = points.getNumIndexes(); @@ -1227,6 +1495,7 @@ std::vector GridLocalPolynomial::getNormalization() const{ return norms; } +template Data2D GridLocalPolynomial::buildUpdateMap(double tolerance, TypeRefinement criteria, int output, const double *scale_correction) const{ int num_points = points.getNumIndexes(); Data2D pmap(num_dimensions, num_points, @@ -1261,9 +1530,9 @@ Data2D GridLocalPolynomial::buildUpdateMap(double tolerance, TypeRefinement } }else{ // construct a series of 1D interpolants and use a refinement criteria that is a combination of the two hierarchical coefficients - Data2D dagUp = HierarchyManipulations::computeDAGup(points, rule.get()); + Data2D dagUp = HierarchyManipulations::computeDAGup(points); - int max_1D_parents = rule->getMaxNumParents(); + int max_1D_parents = RuleLocal::getMaxNumParents(); HierarchyManipulations::SplitDirections split(points); @@ -1296,14 +1565,14 @@ Data2D GridLocalPolynomial::buildUpdateMap(double tolerance, TypeRefinement for(int i=0; igetLevel(points.getIndex(pnts[i])[d]); + levels[i] = RuleLocal::getLevel(points.getIndex(pnts[i])[d]); if (max_level < levels[i]) max_level = levels[i]; } } else { for(int i=0; igetLevel(points.getIndex(pnts[i])[d]); + levels[i] = RuleLocal::getLevel(points.getIndex(pnts[i])[d]); if (max_level < levels[i]) max_level = levels[i]; } } @@ -1312,13 +1581,13 @@ Data2D GridLocalPolynomial::buildUpdateMap(double tolerance, TypeRefinement for(int l=1; l<=max_level; l++){ for(int i=0; igetNode(points.getIndex(pnts[i])[d]); + double x = RuleLocal::getNode(points.getIndex(pnts[i])[d]); double *valsi = vals.getStrip(i); int branch = dagUp.getStrip(pnts[i])[d]; while(branch != -1) { const int *branch_point = points.getIndex(branch); - double basis_value = rule->evalRaw(branch_point[d], x); + double basis_value = RuleLocal::evalRaw(order, branch_point[d], x); const double *branch_vals = vals.getStrip(global_to_pnts[branch]); for(int k=0; k GridLocalPolynomial::buildUpdateMap(double tolerance, TypeRefinement for(int l=1; l<=max_level; l++){ for(int i=0; igetNode(points.getIndex(pnts[i])[d]); + double x = RuleLocal::getNode(points.getIndex(pnts[i])[d]); double *valsi = vals.getStrip(i); int current = 0; @@ -1345,7 +1614,7 @@ Data2D GridLocalPolynomial::buildUpdateMap(double tolerance, TypeRefinement monkey_count[current]++; }else{ const int *branch_point = points.getIndex(branch); - double basis_value = rule->evalRaw(branch_point[d], x); + double basis_value = RuleLocal::evalRaw(order, branch_point[d], x); const double *branch_vals = vals.getStrip(global_to_pnts[branch]); for(int k=0; k GridLocalPolynomial::buildUpdateMap(double tolerance, TypeRefinement } return pmap; } +template MultiIndexSet GridLocalPolynomial::getRefinementCanidates(double tolerance, TypeRefinement criteria, int output, const std::vector &level_limits, const double *scale_correction) const{ - Data2D pmap = buildUpdateMap(tolerance, criteria, output, scale_correction); + Data2D pmap = buildUpdateMap(tolerance, criteria, output, scale_correction); bool useParents = (criteria == refine_fds) || (criteria == refine_parents_first); @@ -1402,8 +1672,8 @@ MultiIndexSet GridLocalPolynomial::getRefinementCanidates(double tolerance, Type const int *map = pmap.getStrip(i); for(int j=0; j(points.getIndex(i), j, points, lrefined))){ + addChild(points.getIndex(i), j, points, lrefined); } } } @@ -1414,8 +1684,8 @@ MultiIndexSet GridLocalPolynomial::getRefinementCanidates(double tolerance, Type const int *map = pmap.getStrip(i); for(int j=0; j(points.getIndex(i), j, points, lrefined))){ + addChildLimited(points.getIndex(i), j, points, level_limits, lrefined); } } } @@ -1433,8 +1703,8 @@ MultiIndexSet GridLocalPolynomial::getRefinementCanidates(double tolerance, Type const int *map = pmap.getStrip(i); for(int j=0; j(points.getIndex(i), j, points, refined))){ + addChild(points.getIndex(i), j, points, refined); } } } @@ -1444,8 +1714,8 @@ MultiIndexSet GridLocalPolynomial::getRefinementCanidates(double tolerance, Type const int *map = pmap.getStrip(i); for(int j=0; j(points.getIndex(i), j, points, refined))){ + addChildLimited(points.getIndex(i), j, points, level_limits, refined); } } } @@ -1455,42 +1725,46 @@ MultiIndexSet GridLocalPolynomial::getRefinementCanidates(double tolerance, Type MultiIndexSet result(refined); if (criteria == refine_stable) - HierarchyManipulations::completeToLower(points, result, rule.get()); + HierarchyManipulations::completeToLower(points, result); + return result; } +template bool GridLocalPolynomial::addParent(const int point[], int direction, const MultiIndexSet &exclude, Data2D &destination) const{ std::vector dad(point, point + num_dimensions); bool added = false; - dad[direction] = rule->getParent(point[direction]); + dad[direction] = RuleLocal::getParent(point[direction]); if ((dad[direction] != -1) && exclude.missing(dad)){ destination.appendStrip(dad); added = true; } - dad[direction] = rule->getStepParent(point[direction]); + dad[direction] = RuleLocal::getStepParent(point[direction]); if ((dad[direction] != -1) && exclude.missing(dad)){ destination.appendStrip(dad); added = true; } return added; } +template void GridLocalPolynomial::addChild(const int point[], int direction, const MultiIndexSet &exclude, Data2D &destination) const{ std::vector kid(point, point + num_dimensions); - int max_1d_kids = rule->getMaxNumKids(); + int max_1d_kids = RuleLocal::getMaxNumKids(); for(int i=0; igetKid(point[direction], i); + kid[direction] = RuleLocal::getKid(point[direction], i); if ((kid[direction] != -1) && exclude.missing(kid)){ destination.appendStrip(kid); } } } +template void GridLocalPolynomial::addChildLimited(const int point[], int direction, const MultiIndexSet &exclude, const std::vector &level_limits, Data2D &destination) const{ std::vector kid(point, point + num_dimensions); - int max_1d_kids = rule->getMaxNumKids(); + int max_1d_kids = RuleLocal::getMaxNumKids(); for(int i=0; igetKid(point[direction], i); + kid[direction] = RuleLocal::getKid(point[direction], i); if ((kid[direction] != -1) - && ((level_limits[direction] == -1) || (rule->getLevel(kid[direction]) <= level_limits[direction])) + && ((level_limits[direction] == -1) || (RuleLocal::getLevel(kid[direction]) <= level_limits[direction])) && exclude.missing(kid)){ destination.appendStrip(kid); } @@ -1508,14 +1782,53 @@ std::vector GridLocalPolynomial::getSupport() const{ MultiIndexSet const &work = (points.empty()) ? needed : points; std::vector support(Utils::size_mult(work.getNumIndexes(), work.getNumDimensions())); - std::transform(work.begin(), work.end(), support.begin(), [&](int p)->double{ return rule->getSupport(p); }); + + switch(effective_rule) { + case RuleLocal::erule::pwc: + std::transform(work.begin(), work.end(), support.begin(), + [&](int p)->double{ return RuleLocal::getSupport(p); }); + break; + case RuleLocal::erule::localp: + std::transform(work.begin(), work.end(), support.begin(), + [&](int p)->double{ return RuleLocal::getSupport(p); }); + break; + case RuleLocal::erule::semilocalp: + std::transform(work.begin(), work.end(), support.begin(), + [&](int p)->double{ return RuleLocal::getSupport(p); }); + break; + case RuleLocal::erule::localp0: + std::transform(work.begin(), work.end(), support.begin(), + [&](int p)->double{ return RuleLocal::getSupport(p); }); + break; + default: // case RuleLocal::erule::localpb: + std::transform(work.begin(), work.end(), support.begin(), + [&](int p)->double{ return RuleLocal::getSupport(p); }); + break; + }; return support; } void GridLocalPolynomial::setSurplusRefinement(double tolerance, TypeRefinement criteria, int output, const std::vector &level_limits, const double *scale_correction){ clearRefinement(); - needed = getRefinementCanidates(tolerance, criteria, output, level_limits, scale_correction); + switch(effective_rule) { + case RuleLocal::erule::pwc: + needed = getRefinementCanidates(tolerance, criteria, output, level_limits, scale_correction); + break; + case RuleLocal::erule::localp: + needed = getRefinementCanidates(tolerance, criteria, output, level_limits, scale_correction); + break; + case RuleLocal::erule::semilocalp: + needed = getRefinementCanidates(tolerance, criteria, output, level_limits, scale_correction); + break; + case RuleLocal::erule::localp0: + needed = getRefinementCanidates(tolerance, criteria, output, level_limits, scale_correction); + break; + default: // case RuleLocal::erule::localpb: + needed = getRefinementCanidates(tolerance, criteria, output, level_limits, scale_correction); + break; + }; + } std::vector GridLocalPolynomial::getScaledCoefficients(int output, const double *scale_correction){ int num_points = points.getNumIndexes(); diff --git a/SparseGrids/tsgGridLocalPolynomial.hpp b/SparseGrids/tsgGridLocalPolynomial.hpp index e1a4542da..cb32e6265 100644 --- a/SparseGrids/tsgGridLocalPolynomial.hpp +++ b/SparseGrids/tsgGridLocalPolynomial.hpp @@ -50,13 +50,17 @@ class GridLocalPolynomial : public BaseCanonicalGrid{ template void write(std::ostream &os) const; - TypeOneDRule getRule() const override{ return rule->getType(); } + TypeOneDRule getRule() const override{ return RuleLocal::getRule(effective_rule); } int getOrder() const{ return order; } + template + void getPoints(double *x) const; void getLoadedPoints(double *x) const override; void getNeededPoints(double *x) const override; void getPoints(double *x) const override; // returns the loaded points unless no points are loaded, then returns the needed points + template + void getQuadratureWeights(double weights[]) const; void getQuadratureWeights(double weights[]) const override; void getInterpolationWeights(const double x[], double weights[]) const override; void getDifferentiationWeights(const double x[], double weights[]) const override; @@ -102,8 +106,14 @@ class GridLocalPolynomial : public BaseCanonicalGrid{ void beginConstruction() override; void writeConstructionData(std::ostream &os, bool) const override; void readConstructionData(std::istream &is, bool) override; + template std::vector getCandidateConstructionPoints(double tolerance, TypeRefinement criteria, int output, std::vector const &level_limits, double const *scale_correction); + std::vector getCandidateConstructionPoints(double tolerance, TypeRefinement criteria, int output, std::vector const &level_limits, double const *scale_correction); + template + void loadConstructedPoint(const double x[], const std::vector &y); void loadConstructedPoint(const double x[], const std::vector &y) override; + template + void loadConstructedPoint(const double x[], int numx, const double y[]); void loadConstructedPoint(const double x[], int numx, const double y[]) override; void finishConstruction() override; @@ -131,18 +141,25 @@ class GridLocalPolynomial : public BaseCanonicalGrid{ void buildTree(); //! \brief Returns a list of indexes of the nodes in \b points that are descendants of the \b point. + template std::vector getSubGraph(std::vector const &point) const; //! \brief Add the \b point to the grid using the \b values. + template void expandGrid(std::vector const &point, std::vector const &value); //! \brief Return the multi-index of canonical point \b x. + template std::vector getMultiIndex(const double x[]); //! \brief Looks for a batch of constructed points and processes all that will result in a connected graph. + template void loadConstructedPoints(); //! \brief Fast algorithm, uses global Kronecker algorithm to recompute all surpluses + template + void recomputeSurpluses(); + void recomputeSurpluses(); /*! @@ -160,14 +177,45 @@ class GridLocalPolynomial : public BaseCanonicalGrid{ * Note: see the comments inside recomputeSurpluses() for the performance comparison between different algorithms * also note that this method can be used to partially update, i.e., update the surpluses for some of the indexes */ + template void updateSurpluses(MultiIndexSet const &work, int max_level, std::vector const &level, Data2D const &dagUp); // Same idea as in GridSequence::applyTransformationTransposed(). - template void applyTransformationTransposed(double weights[], const MultiIndexSet &work, const std::vector &active_points) const; + template + void applyTransformationTransposed(double weights[], const MultiIndexSet &work, const std::vector &active_points) const; + template + void applyTransformationTransposed(double weights[], const MultiIndexSet &work, const std::vector &active_points) const; void buildSparseMatrixBlockForm(const double x[], int num_x, int num_chunk, std::vector &numnz, std::vector> &tindx, std::vector> &tvals) const; + template + double evalBasisSupported(const int point[], const double x[], bool &isSupported) const{ + double f = RuleLocal::evalSupport(order, point[0], x[0], isSupported); + if (!isSupported) return 0.0; + for(int j=1; j(order, point[j], x[j], isSupported); + if (!isSupported) return 0.0; + } + return f; + } + template + void diffBasisSupported(const int point[], const double x[], double diff_values[], bool &isSupported) const{ + isSupported = false; + for(int i=0; i(order, point[k], x[k], isDimSupported); + isSupported = isDimSupported or isSupported; + for(int j=0; j(order, point[k], x[k], isDimSupported); + isSupported = isDimSupported or isSupported; + } + } + /*! * \brief Walk through all the nodes of the tree and touches only the nodes supported at \b x. * @@ -183,17 +231,21 @@ class GridLocalPolynomial : public BaseCanonicalGrid{ * * In all cases, \b work is the \b points or \b needed set that has been used to construct the tree. */ - template + template void walkTree(const MultiIndexSet &work, const double x[], std::vector &sindx, std::vector &svals, double *y) const{ std::vector monkey_count(top_level+1); // traverse the tree, counts the branches of the current node std::vector monkey_tail(top_level+1); // traverse the tree, keeps track of the previous node (history) + bool isSupported; + double basis_value; + std::vector basis_derivative(num_dimensions); + for(const auto &r : roots){ - bool isSupported; - double basis_value = (mode == 3 or mode == 4) ? 0.0 : evalBasisSupported(work.getIndex(r), x, isSupported); - std::vector basis_derivative(num_dimensions); - if (mode == 3 or mode == 4) - diffBasisSupported(work.getIndex(r), x, basis_derivative.data(), isSupported); + if (mode == 3 or mode == 4) { + diffBasisSupported(work.getIndex(r), x, basis_derivative.data(), isSupported); + } else { + basis_value = evalBasisSupported(work.getIndex(r), x, isSupported); + } if (isSupported){ if (mode == 0){ @@ -222,9 +274,9 @@ class GridLocalPolynomial : public BaseCanonicalGrid{ if (monkey_count[current] < pntr[monkey_tail[current]+1]){ int p = indx[monkey_count[current]]; if (mode == 3 or mode == 4){ - diffBasisSupported(work.getIndex(p), x, basis_derivative.data(), isSupported); + diffBasisSupported(work.getIndex(p), x, basis_derivative.data(), isSupported); }else{ - basis_value = evalBasisSupported(work.getIndex(p), x, isSupported); + basis_value = evalBasisSupported(work.getIndex(p), x, isSupported); } if (isSupported){ if (mode == 0){ @@ -271,18 +323,47 @@ class GridLocalPolynomial : public BaseCanonicalGrid{ std::transform(map.begin(), map.end(), svals.begin(), [&](int i)->double{ return vls[i]; }); } } + // Explicitly instantiates based on the effective_rule + template + void walkTree(const MultiIndexSet &work, const double x[], std::vector &sindx, std::vector &svals, double *y) const{ + switch(effective_rule) { + case RuleLocal::erule::pwc: + walkTree(work,x, sindx, svals, y); + break; + case RuleLocal::erule::localp: + walkTree(work,x, sindx, svals, y); + break; + case RuleLocal::erule::semilocalp: + walkTree(work,x, sindx, svals, y); + break; + case RuleLocal::erule::localp0: + walkTree(work,x, sindx, svals, y); + break; + default: // case RuleLocal::erule::localpb: + walkTree(work,x, sindx, svals, y); + break; + }; + } - double evalBasisRaw(const int point[], const double x[]) const; - double evalBasisSupported(const int point[], const double x[], bool &isSupported) const; - void diffBasisSupported(const int point[], const double x[], double diff_values[], bool &isSupported) const; + template + double evalBasisRaw(const int point[], const double x[]) const { + double f = RuleLocal::evalRaw(order, point[0], x[0]); + for(int j=1; j(order, point[j], x[j]); + return f; + } std::vector getNormalization() const; + template Data2D buildUpdateMap(double tolerance, TypeRefinement criteria, int output, const double *scale_correction) const; + template MultiIndexSet getRefinementCanidates(double tolerance, TypeRefinement criteria, int output, const std::vector &level_limits, const double *scale_correction) const; + template bool addParent(const int point[], int direction, const MultiIndexSet &exclude, Data2D &destination) const; + template void addChild(const int point[], int direction, const MultiIndexSet &exclude, Data2D &destination) const; + template void addChildLimited(const int point[], int direction, const MultiIndexSet &exclude, const std::vector &level_limits, Data2D &destination) const; void clearGpuSurpluses(); @@ -300,7 +381,7 @@ class GridLocalPolynomial : public BaseCanonicalGrid{ std::vector pntr; std::vector indx; - std::unique_ptr rule; + RuleLocal::erule effective_rule; std::unique_ptr dynamic_values; @@ -312,8 +393,23 @@ class GridLocalPolynomial : public BaseCanonicalGrid{ const int* p = work.getIndex(i); T *s = cpu_support.getStrip(i); for(int j=0; j(rule->getSupport(p[j])); - if (ord != 0){ + if (ord == 0){ + s[j] = static_cast(RuleLocal::getSupport(p[j])); + } else { + switch(crule) { + case rule_localp: + s[j] = static_cast(RuleLocal::getSupport(p[j])); + break; + case rule_semilocalp: + s[j] = static_cast(RuleLocal::getSupport(p[j])); + break; + case rule_localp0: + s[j] = static_cast(RuleLocal::getSupport(p[j])); + break; + case rule_localpb: + s[j] = static_cast(RuleLocal::getSupport(p[j])); + break; + }; if (ord == 2) s[j] *= s[j]; if ((crule == rule_localp) || (crule == rule_semilocalp)) if (p[j] == 0) s[j] = static_cast(-1.0); // constant function if ((crule == rule_localp) && (ord == 2)){ @@ -354,7 +450,7 @@ template<> struct GridReaderVersion5{ grid->order = IO::readNumber(is); grid->top_level = IO::readNumber(is); TypeOneDRule rule = IO::readRule(is); - grid->rule = makeRuleLocalPolynomial(rule, grid->order); + grid->effective_rule = RuleLocal::getEffectiveRule(grid->order, rule); if (IO::readFlag(is)) grid->points = MultiIndexSet(is, iomode()); if (std::is_same::value){ // backwards compatible: surpluses and needed, or needed and surpluses @@ -366,8 +462,19 @@ template<> struct GridReaderVersion5{ if (IO::readFlag(is)) grid->surpluses = IO::readData2D(is, grid->num_outputs, grid->points.getNumIndexes()); } + int max_parents = [&]()->int { + switch(grid->effective_rule) { + case RuleLocal::erule::pwc: return RuleLocal::getMaxNumParents(); + case RuleLocal::erule::localp: return RuleLocal::getMaxNumParents(); + case RuleLocal::erule::semilocalp: return RuleLocal::getMaxNumParents(); + case RuleLocal::erule::localp0: return RuleLocal::getMaxNumParents(); + default: // case RuleLocal::erule::localpb: + return RuleLocal::getMaxNumParents(); + }; + }(); + if (IO::readFlag(is)) - grid->parents = IO::readData2D(is, grid->rule->getMaxNumParents() * grid->num_dimensions, grid->points.getNumIndexes()); + grid->parents = IO::readData2D(is, max_parents * grid->num_dimensions, grid->points.getNumIndexes()); size_t num_points = (size_t) ((grid->points.empty()) ? grid->needed.getNumIndexes() : grid->points.getNumIndexes()); grid->roots = std::vector((size_t) IO::readNumber(is)); diff --git a/SparseGrids/tsgHierarchyManipulator.cpp b/SparseGrids/tsgHierarchyManipulator.cpp index b3803ba19..a9aa637ce 100644 --- a/SparseGrids/tsgHierarchyManipulator.cpp +++ b/SparseGrids/tsgHierarchyManipulator.cpp @@ -37,13 +37,14 @@ namespace TasGrid{ namespace HierarchyManipulations{ -Data2D computeDAGup(MultiIndexSet const &mset, const BaseRuleLocalPolynomial *rule){ +template +Data2D computeDAGup(MultiIndexSet const &mset){ size_t num_dimensions = mset.getNumDimensions(); int num_points = mset.getNumIndexes(); - if (rule->getMaxNumParents() > 1){ // allow for multiple parents and level 0 may have more than one node - int max_parents = rule->getMaxNumParents() * (int) num_dimensions; + if (RuleLocal::getMaxNumParents() > 1){ // allow for multiple parents and level 0 may have more than one node + int max_parents = RuleLocal::getMaxNumParents() * (int) num_dimensions; Data2D parents(max_parents, num_points, -1); - int level0_offset = rule->getNumPoints(0); + int level0_offset = RuleLocal::getNumPoints(0); #pragma omp parallel for schedule(static) for(int i=0; i computeDAGup(MultiIndexSet const &mset, const BaseRuleLocalPolynomia for(size_t j=0; j= level0_offset){ int current = p[j]; - dad[j] = rule->getParent(current); + dad[j] = RuleLocal::getParent(current); pp[2*j] = mset.getSlot(dad); while ((dad[j] >= level0_offset) && (pp[2*j] == -1)){ current = dad[j]; - dad[j] = rule->getParent(current); + dad[j] = RuleLocal::getParent(current); pp[2*j] = mset.getSlot(dad); } - dad[j] = rule->getStepParent(current); + dad[j] = RuleLocal::getStepParent(current); if (dad[j] != -1){ pp[2*j + 1] = mset.getSlot(dad); } @@ -81,10 +82,10 @@ Data2D computeDAGup(MultiIndexSet const &mset, const BaseRuleLocalPolynomia if (dad[j] == 0){ pp[j] = -1; }else{ - dad[j] = rule->getParent(dad[j]); + dad[j] = RuleLocal::getParent(dad[j]); pp[j] = mset.getSlot(dad.data()); while((dad[j] != 0) && (pp[j] == -1)){ - dad[j] = rule->getParent(dad[j]); + dad[j] = RuleLocal::getParent(dad[j]); pp[j] = mset.getSlot(dad); } dad[j] = p[j]; @@ -94,13 +95,36 @@ Data2D computeDAGup(MultiIndexSet const &mset, const BaseRuleLocalPolynomia return parents; } } -Data2D computeDAGup(MultiIndexSet const &mset, const BaseRuleLocalPolynomial *rule, bool &is_complete){ + +template Data2D computeDAGup(MultiIndexSet const &mset); +template Data2D computeDAGup(MultiIndexSet const &mset); +template Data2D computeDAGup(MultiIndexSet const &mset); +template Data2D computeDAGup(MultiIndexSet const &mset); +template Data2D computeDAGup(MultiIndexSet const &mset); + +Data2D computeDAGup(MultiIndexSet const &mset, RuleLocal::erule effrule) { + switch(effrule) { + case RuleLocal::erule::pwc: + return computeDAGup(mset); + case RuleLocal::erule::localp: + return computeDAGup(mset); + case RuleLocal::erule::semilocalp: + return computeDAGup(mset); + case RuleLocal::erule::localp0: + return computeDAGup(mset); + default: // case RuleLocal::erule::localpb: + return computeDAGup(mset); + }; +} + +template +Data2D computeDAGup(MultiIndexSet const &mset, bool &is_complete){ size_t num_dimensions = mset.getNumDimensions(); int num_points = mset.getNumIndexes(); - if (rule->getMaxNumParents() > 1){ // allow for multiple parents and level 0 may have more than one node - int max_parents = rule->getMaxNumParents() * (int) num_dimensions; + if (RuleLocal::getMaxNumParents() > 1){ // allow for multiple parents and level 0 may have more than one node + int max_parents = RuleLocal::getMaxNumParents() * (int) num_dimensions; Data2D parents(max_parents, num_points, -1); - int level0_offset = rule->getNumPoints(0); + int level0_offset = RuleLocal::getNumPoints(0); int any_fail = 0; // count if there are failures #pragma omp parallel { @@ -115,16 +139,16 @@ Data2D computeDAGup(MultiIndexSet const &mset, const BaseRuleLocalPolynomia for(size_t j=0; j= level0_offset){ int current = p[j]; - dad[j] = rule->getParent(current); + dad[j] = RuleLocal::getParent(current); pp[2*j] = mset.getSlot(dad); if (pp[2*j] == -1) fail = 1; while ((dad[j] >= level0_offset) && (pp[2*j] == -1)){ current = dad[j]; - dad[j] = rule->getParent(current); + dad[j] = RuleLocal::getParent(current); pp[2*j] = mset.getSlot(dad); } - dad[j] = rule->getStepParent(current); + dad[j] = RuleLocal::getStepParent(current); if (dad[j] != -1){ pp[2*j + 1] = mset.getSlot(dad); if (pp[2*j + 1] == -1) @@ -157,12 +181,12 @@ Data2D computeDAGup(MultiIndexSet const &mset, const BaseRuleLocalPolynomia if (dad[j] == 0){ pp[j] = -1; }else{ - dad[j] = rule->getParent(dad[j]); + dad[j] = RuleLocal::getParent(dad[j]); pp[j] = mset.getSlot(dad.data()); if (pp[j] == -1) fail = 1; while((dad[j] != 0) && (pp[j] == -1)){ - dad[j] = rule->getParent(dad[j]); + dad[j] = RuleLocal::getParent(dad[j]); pp[j] = mset.getSlot(dad); } dad[j] = p[j]; @@ -178,9 +202,16 @@ Data2D computeDAGup(MultiIndexSet const &mset, const BaseRuleLocalPolynomia } } -Data2D computeDAGDown(MultiIndexSet const &mset, const BaseRuleLocalPolynomial *rule){ +template Data2D computeDAGup(MultiIndexSet const &mset, bool &is_complete); +template Data2D computeDAGup(MultiIndexSet const &mset, bool &is_complete); +template Data2D computeDAGup(MultiIndexSet const &mset, bool &is_complete); +template Data2D computeDAGup(MultiIndexSet const &mset, bool &is_complete); +template Data2D computeDAGup(MultiIndexSet const &mset, bool &is_complete); + +template +Data2D computeDAGDown(MultiIndexSet const &mset){ size_t num_dimensions = mset.getNumDimensions(); - int max_1d_kids = rule->getMaxNumKids(); + int max_1d_kids = RuleLocal::getMaxNumKids(); int num_points = mset.getNumIndexes(); Data2D kids(Utils::size_mult(max_1d_kids, num_dimensions), num_points); @@ -193,7 +224,7 @@ Data2D computeDAGDown(MultiIndexSet const &mset, const BaseRuleLocalPolynom for(size_t j=0; jgetKid(current, k); + kid[j] = RuleLocal::getKid(current, k); *family++ = (kid[j] == -1) ? -1 : mset.getSlot(kid); } kid[j] = current; @@ -203,23 +234,53 @@ Data2D computeDAGDown(MultiIndexSet const &mset, const BaseRuleLocalPolynom return kids; } -std::vector computeLevels(MultiIndexSet const &mset, BaseRuleLocalPolynomial const *rule){ +template Data2D computeDAGDown(MultiIndexSet const &mset); +template Data2D computeDAGDown(MultiIndexSet const &mset); +template Data2D computeDAGDown(MultiIndexSet const &mset); +template Data2D computeDAGDown(MultiIndexSet const &mset); +template Data2D computeDAGDown(MultiIndexSet const &mset); + + +template +std::vector computeLevels(MultiIndexSet const &mset){ size_t num_dimensions = mset.getNumDimensions(); int num_points = mset.getNumIndexes(); std::vector level((size_t) num_points); #pragma omp parallel for schedule(static) for(int i=0; igetLevel(p[0]); + int current_level = RuleLocal::getLevel(p[0]); for(size_t j=1; jgetLevel(p[j]); + current_level += RuleLocal::getLevel(p[j]); } level[i] = current_level; } return level; } -void completeToLower(MultiIndexSet const &mset, MultiIndexSet &refined, BaseRuleLocalPolynomial const *rule){ +template std::vector computeLevels(MultiIndexSet const &mset); +template std::vector computeLevels(MultiIndexSet const &mset); +template std::vector computeLevels(MultiIndexSet const &mset); +template std::vector computeLevels(MultiIndexSet const &mset); +template std::vector computeLevels(MultiIndexSet const &mset); + +std::vector computeLevels(MultiIndexSet const &mset, RuleLocal::erule effrule) { + switch(effrule) { + case RuleLocal::erule::pwc: + return computeLevels(mset); + case RuleLocal::erule::localp: + return computeLevels(mset); + case RuleLocal::erule::semilocalp: + return computeLevels(mset); + case RuleLocal::erule::localp0: + return computeLevels(mset); + default: // case RuleLocal::erule::localpb: + return computeLevels(mset); + }; +} + +template +void completeToLower(MultiIndexSet const &mset, MultiIndexSet &refined){ size_t num_dimensions = mset.getNumDimensions(); size_t num_added = 1; // set to 1 to start the loop while(num_added > 0){ @@ -230,10 +291,10 @@ void completeToLower(MultiIndexSet const &mset, MultiIndexSet &refined, BaseRule std::vector parent(refined.getIndex(i), refined.getIndex(i) + num_dimensions); for(auto &p : parent){ int r = p; - p = rule->getParent(r); + p = RuleLocal::getParent(r); if ((p != -1) && refined.missing(parent) && mset.missing(parent)) addons.appendStrip(parent); - p = rule->getStepParent(r); + p = RuleLocal::getStepParent(r); if ((p != -1) && refined.missing(parent) && mset.missing(parent)) addons.appendStrip(parent); p = r; @@ -245,6 +306,12 @@ void completeToLower(MultiIndexSet const &mset, MultiIndexSet &refined, BaseRule } } +template void completeToLower(MultiIndexSet const &mset, MultiIndexSet &refined); +template void completeToLower(MultiIndexSet const &mset, MultiIndexSet &refined); +template void completeToLower(MultiIndexSet const &mset, MultiIndexSet &refined); +template void completeToLower(MultiIndexSet const &mset, MultiIndexSet &refined); +template void completeToLower(MultiIndexSet const &mset, MultiIndexSet &refined); + SplitDirections::SplitDirections(const MultiIndexSet &points){ // split the points into "jobs", where each job represents a batch of // points that lay on a line in some direction diff --git a/SparseGrids/tsgHierarchyManipulator.hpp b/SparseGrids/tsgHierarchyManipulator.hpp index 3c6f4eb72..fb41a27a3 100644 --- a/SparseGrids/tsgHierarchyManipulator.hpp +++ b/SparseGrids/tsgHierarchyManipulator.hpp @@ -84,7 +84,17 @@ namespace HierarchyManipulations{ * (or -1 if the parent is missing from \b mset). * \endinternal */ -Data2D computeDAGup(MultiIndexSet const &mset, const BaseRuleLocalPolynomial *rule); +template +Data2D computeDAGup(MultiIndexSet const &mset); +/*! + * \internal + * \ingroup TasmanianHierarchyManipulations + * \brief Cache the indexes slot numbers of the parents of the multi-indexes in \b mset. + * + * The effective rule is passed in at runtime and a switch statement finds the correct template. + * \endinternal + */ +Data2D computeDAGup(MultiIndexSet const &mset, RuleLocal::erule effrule); /*! * \internal * \ingroup TasmanianHierarchyManipulations @@ -95,7 +105,8 @@ Data2D computeDAGup(MultiIndexSet const &mset, const BaseRuleLocalPolynomia * On exit, \b is_complete will indicate whether there are points with missing parents. * \endinternal */ -Data2D computeDAGup(MultiIndexSet const &mset, const BaseRuleLocalPolynomial *rule, bool &is_complete); +template +Data2D computeDAGup(MultiIndexSet const &mset, bool &is_complete); /*! * \internal @@ -108,7 +119,8 @@ Data2D computeDAGup(MultiIndexSet const &mset, const BaseRuleLocalPolynomia * (or -1 if the kid is missing from \b mset). * \endinternal */ -Data2D computeDAGDown(MultiIndexSet const &mset, const BaseRuleLocalPolynomial *rule); +template +Data2D computeDAGDown(MultiIndexSet const &mset); /*! * \internal @@ -116,7 +128,15 @@ Data2D computeDAGDown(MultiIndexSet const &mset, const BaseRuleLocalPolynom * \ingroup TasmanianHierarchyManipulations * \endinternal */ -std::vector computeLevels(MultiIndexSet const &mset, BaseRuleLocalPolynomial const *rule); +template +std::vector computeLevels(MultiIndexSet const &mset); +/*! + * \internal + * \brief Overload that turns switch statement into template instantiations + * \ingroup TasmanianHierarchyManipulations + * \endinternal + */ +std::vector computeLevels(MultiIndexSet const &mset, RuleLocal::erule effrule); /*! * \internal @@ -124,7 +144,8 @@ std::vector computeLevels(MultiIndexSet const &mset, BaseRuleLocalPolynomia * \brief Complete \b refined so that the union of \b refined and \b mset is lower w.r.t. the \b rule. * \endinternal */ -void completeToLower(MultiIndexSet const &mset, MultiIndexSet &refined, BaseRuleLocalPolynomial const *rule); +template +void completeToLower(MultiIndexSet const &mset, MultiIndexSet &refined); /*! * \internal @@ -132,20 +153,21 @@ void completeToLower(MultiIndexSet const &mset, MultiIndexSet &refined, BaseRule * \ingroup TasmanianHierarchyManipulations * \endinternal */ -inline void touchAllImmediateRelatives(std::vector &point, MultiIndexSet const &mset, BaseRuleLocalPolynomial const *rule, std::function apply){ - int max_kids = rule->getMaxNumKids(); +template +void touchAllImmediateRelatives(std::vector &point, MultiIndexSet const &mset, callable_method apply){ + int max_kids = RuleLocal::getMaxNumKids(); for(auto &v : point){ int save = v; // replace one by one each index of p with either parent or kid // check the parents - v = rule->getParent(save); + v = RuleLocal::getParent(save); if (v > -1){ int parent_index = mset.getSlot(point); if (parent_index > -1) apply(parent_index); } - v = rule->getStepParent(save); + v = RuleLocal::getStepParent(save); if (v > -1){ int parent_index = mset.getSlot(point); if (parent_index > -1) @@ -153,7 +175,7 @@ inline void touchAllImmediateRelatives(std::vector &point, MultiIndexSet co } for(int k=0; kgetKid(save, k); + v = RuleLocal::getKid(save, k); if (v > -1){ int kid_index = mset.getSlot(point); if (kid_index > -1) @@ -172,9 +194,10 @@ inline void touchAllImmediateRelatives(std::vector &point, MultiIndexSet co * * \endinternal */ -inline MultiIndexSet getLevelZeroPoints(size_t num_dimensions, BaseRuleLocalPolynomial const *rule){ +template +MultiIndexSet getLevelZeroPoints(size_t num_dimensions){ int num_parents = 0; - while(rule->getParent(num_parents) == -1) num_parents++; + while(RuleLocal::getParent(num_parents) == -1) num_parents++; return MultiIndexManipulations::generateFullTensorSet(std::vector(num_dimensions, num_parents)); } @@ -185,12 +208,13 @@ inline MultiIndexSet getLevelZeroPoints(size_t num_dimensions, BaseRuleLocalPoly * * \endinternal */ -inline MultiIndexSet getLargestConnected(MultiIndexSet const ¤t, MultiIndexSet const &candidates, BaseRuleLocalPolynomial const *rule){ +template +MultiIndexSet getLargestConnected(MultiIndexSet const ¤t, MultiIndexSet const &candidates){ if (candidates.empty()) return MultiIndexSet(); auto num_dimensions = candidates.getNumDimensions(); // always consider the points without parents - MultiIndexSet level_zero = getLevelZeroPoints(num_dimensions, rule); + MultiIndexSet level_zero = getLevelZeroPoints(num_dimensions); MultiIndexSet result; // keep track of the cumulative result MultiIndexSet total = current; // forms a working copy of the entire merged graph @@ -212,8 +236,8 @@ inline MultiIndexSet getLargestConnected(MultiIndexSet const ¤t, MultiInde if (total.empty()) return MultiIndexSet(); // current was empty and no roots could be added - int max_kids = rule->getMaxNumKids(); - int max_relatives = rule->getMaxNumParents() + max_kids; + int max_kids = RuleLocal::getMaxNumKids(); + int max_relatives = RuleLocal::getMaxNumParents() + max_kids; Data2D update; do{ update = Data2D(num_dimensions, 0); @@ -223,7 +247,8 @@ inline MultiIndexSet getLargestConnected(MultiIndexSet const ¤t, MultiInde for(auto &r : relative){ int k = r; // save the value for(int j=0; jgetKid(k, j) : ((j - max_kids == 0) ? rule->getParent(k) : rule->getStepParent(k)); + r = (j < max_kids) ? RuleLocal::getKid(k, j) + : ((j - max_kids == 0) ? RuleLocal::getParent(k) : RuleLocal::getStepParent(k)); if ((r != -1) && !candidates.missing(relative) && total.missing(relative)) update.appendStrip(relative); } @@ -241,6 +266,7 @@ inline MultiIndexSet getLargestConnected(MultiIndexSet const ¤t, MultiInde return result; } + /*! * \internal * \ingroup TasmanianHierarchyManipulations diff --git a/SparseGrids/tsgRuleLocalPolynomial.hpp b/SparseGrids/tsgRuleLocalPolynomial.hpp index 62644f75a..fe416776b 100644 --- a/SparseGrids/tsgRuleLocalPolynomial.hpp +++ b/SparseGrids/tsgRuleLocalPolynomial.hpp @@ -36,679 +36,322 @@ namespace TasGrid{ #ifndef __TASMANIAN_DOXYGEN_SKIP -class BaseRuleLocalPolynomial{ -public: - BaseRuleLocalPolynomial() : max_order(0){} - virtual ~BaseRuleLocalPolynomial() = default; - virtual int getMaxOrder() const = 0; - virtual void setMaxOrder(int order) = 0; +namespace RuleLocal { + // effective rule type + enum class erule { + pwc, localp, semilocalp, localp0, localpb + }; - virtual TypeOneDRule getType() const = 0; - - virtual int getNumPoints(int level) const = 0; // for building the initial grid - - virtual double getNode(int point) const = 0; - - virtual int getLevel(int point) const = 0; - virtual double getSupport(int point) const = 0; - - virtual int getMaxNumKids() const = 0; - virtual int getMaxNumParents() const = 0; - - virtual int getParent(int point) const = 0; - virtual int getStepParent(int point) const = 0; - virtual int getKid(int point, int kid_number) const = 0; - - virtual double evalRaw(int point, double x) const = 0; // normalizes x (i.e., (x-node) / support), but it does not check the support - virtual double evalSupport(int point, double x, bool &isSupported) const = 0; // // normalizes x (i.e., (x-node) / support) and checks if x is within the support - - // Analogous to the eval functions, but differentiates at x instead. - virtual double diffRaw(int point, double x) const = 0; - virtual double diffSupport(int point, double x, bool &isSupported) const = 0; - - virtual double getArea(int point, std::vector const &w, std::vector const &x) const = 0; - // integrate the function associated with the point, constant to cubic are known analytically, higher order need a 1-D quadrature rule - - // constructs a Vandemond matrix for the first num_rows points - virtual void van_matrix(int num_rows, std::vector &pntr, std::vector &indx, std::vector &vals) = 0; - -protected: - int max_order; -}; - - -template -class templRuleLocalPolynomial : public BaseRuleLocalPolynomial{ -public: - templRuleLocalPolynomial() = default; - ~templRuleLocalPolynomial() = default; - - int getMaxOrder() const override{ return max_order; } - void setMaxOrder(int order) override{ max_order = order; } - - TypeOneDRule getType() const override{ return rule; } - - int getNumPoints(int level) const override{ - if (isZeroOrder){ - int n = 1; - while (level-- > 0) n *= 3; - return n; - }else{ - if ((rule == rule_localp) || (rule == rule_semilocalp)){ - return (level == 0) ? 1 : ((1 << level) + 1); - }else if (rule == rule_localp0){ - return (1 << (level+1)) -1; - }else{ // rule == rule_localpb - return ((1 << level) + 1); - } + inline erule getEffectiveRule(int order, TypeOneDRule rule) { + if (order == 0) return erule::pwc; + switch(rule) { + case rule_semilocalp: return erule::semilocalp; + case rule_localp0: return erule::localp0; + case rule_localpb: return erule::localpb; + default: + return erule::localp; } } - - double getNode(int point) const override{ - if (isZeroOrder){ - return -2.0 + (1.0 / ((double) Maths::int3log3(point))) * (3*point + 2 - point % 2); - }else{ - if ((rule == rule_localp) || (rule == rule_semilocalp)){ - if (point == 0) return 0.0; - if (point == 1) return -1.0; - if (point == 2) return 1.0; - return ((double)(2*point - 1)) / ((double) Maths::int2log2(point - 1)) - 3.0; - }else if (rule == rule_localp0){ - return ((double)(2*point +3) ) / ((double) Maths::int2log2(point + 1) ) - 3.0; - }else{ // rule == rule_localpb - if (point == 0) return -1.0; - if (point == 1) return 1.0; - if (point == 2) return 0.0; - return ((double)(2*point - 1)) / ((double) Maths::int2log2(point - 1)) - 3.0; - } - } + inline TypeOneDRule getRule(erule effective_rule) { + switch(effective_rule) { + case erule::pwc: + case erule::localp: return rule_localp; + case erule::semilocalp: return rule_semilocalp; + case erule::localp0: return rule_localp0; + default: //case erule::localpb: + return rule_localpb; + }; } - int getLevel(int point) const override{ - if (isZeroOrder){ - int level = 0; - while(point >= 1){ point /= 3; level += 1; } - return level; - }else{ - if ((rule == rule_localp) || (rule == rule_semilocalp)){ - return (point == 0) ? 0 : (point == 1) ? 1 : (Maths::intlog2(point - 1) + 1); - }else if (rule == rule_localp0){ - return Maths::intlog2(point + 1); - }else{ // rule == rule_localpb - return ((point == 0) || (point == 1)) ? 0 : (Maths::intlog2(point - 1) + 1); - } - } - } - double getSupport(int point) const override{ - if (isZeroOrder){ - return 1.0 / (double) Maths::int3log3(point); - }else{ - if ((rule == rule_localp) || (rule == rule_semilocalp)){ - return (point == 0) ? 1.0 : 1.0 / ((double) Maths::int2log2(point - 1)); - }else if (rule == rule_localp0){ - return 1.0 / ((double) Maths::int2log2(point + 1)); - }else{ // rule == rule_localpb - return ((point == 0) || (point == 1)) ? 2.0 : 1.0 / ((double) Maths::int2log2(point - 1)); + template + int getNumPoints(int level) { + switch(effective_rule) { + case erule::pwc: { + int n = 1; + while (level-- > 0) n *= 3; + return n; } - } + case erule::localp: + case erule::semilocalp: + return (level == 0) ? 1 : ((1 << level) + 1); + case erule::localp0: + return (1 << (level+1)) -1; + default: // case erule::localpb: + return ((1 << level) + 1); + }; } - int getMaxNumKids() const override{ return (isZeroOrder) ? 4 : 2; } - int getMaxNumParents() const override{ return ((isZeroOrder || (rule == rule_semilocalp) || (rule == rule_localpb)) ? 2 : 1); } + template + int getMaxNumKids() { return (effective_rule == erule::pwc) ? 4 : 2; } + template + int getMaxNumParents() { + return ((effrule == erule::pwc or effrule == erule::semilocalp or effrule == erule::localpb) ? 2 : 1); + } - int getParent(int point) const override{ - if (isZeroOrder){ - return (point == 0) ? -1 : point / 3; - }else{ - if ((rule == rule_localp) || (rule == rule_semilocalp)){ + template + int getParent(int point) { + switch(effective_rule) { + case erule::pwc: + return (point == 0) ? -1 : point / 3; + case erule::localp: + case erule::semilocalp: { int dad = (point + 1) / 2; if (point < 4) dad--; return dad; - }else if (rule == rule_localp0){ - if (point == 0) return -1; - return (point - 1) / 2; - }else{ // rule == rule_localpb - return (point < 2) ? -1 : ((point + 1) / 2); } - } + case erule::localp0: + return (point == 0) ? -1 : (point - 1) / 2; + default: // case erule::localpb: + return (point < 2) ? -1 : ((point + 1) / 2); + }; } - int getStepParent(int point) const override{ - if (isZeroOrder){ + template + int getStepParent(int point) { + if (effective_rule == erule::pwc){ int i3l3 = Maths::int3log3(point); if (point == i3l3/3) return -1; if (point == i3l3-1) return -1; - if ((point % 3 == 2) && (point % 2 == 0)) return point / 3 + 1; - if ((point % 3 == 0) && (point % 2 == 1)) return point / 3 - 1; + int mod3 = point % 3; + int mod2 = point % 2; + if (mod3 == 2 and mod2 == 0) return point / 3 + 1; + if (mod3 == 0 and mod2 == 1) return point / 3 - 1; return -1; }else{ - if (rule == rule_semilocalp){ - if (point == 3) return 2; - if (point == 4) return 1; - }else if (rule == rule_localpb){ - if (point == 2) return 0; + if (effective_rule == erule::semilocalp){ + switch(point) { + case 3: return 2; + case 4: return 1; + default: + return -1; + }; + }else if (effective_rule == erule::localpb){ + return (point == 2) ? 0 : -1; } return -1; } } - int getKid(int point, int kid_number) const override{ - if (isZeroOrder){ - if (point == 0) return (kid_number == 0) ? 1 : (kid_number==1) ? 2 : -1; - if (kid_number == 3){ - int i3l3 = Maths::int3log3(point); - if (point == i3l3/3) return -1; - if (point == i3l3-1) return -1; - return (point % 2 == 0) ? 3*point + 3 : 3*point - 1; - } - return 3*point + kid_number; - }else if ((rule == rule_localp) || (rule == rule_semilocalp)){ - if (kid_number == 0){ - if (point == 0) return 1; - if (point == 1) return 3; - if (point == 2) return 4; - return 2*point-1; - }else{ - if (point == 0) return 2; - if ((point == 1) || (point == 2)) return -1; - return 2*point; - } - }else if (rule == rule_localp0){ - return 2*point + ((kid_number == 0) ? 1 : 2); - }else{ // localpb - if ((point == 0) || (point == 1)) return (kid_number == 0) ? 2 : -1; - return 2*point - ((kid_number == 0) ? 1 : 0); - } - } - - double evalRaw(int point, double x) const override{ - if (isZeroOrder){ - if (std::abs(x - getNode(point)) > getSupport(point)) return 0.0; // maybe can speed this up - return 1.0; - }else{ - if ((rule == rule_localp) || (rule == rule_semilocalp)){ - if (point == 0) return 1.0; - if (rule == rule_semilocalp){ - if (point == 1) return 0.5 * x * (x - 1.0); - if (point == 2) return 0.5 * x * (x + 1.0); + template + int getKid(int point, int kid_number) { + switch(effective_rule) { + case erule::pwc: + if (point == 0) return (kid_number == 0) ? 1 : (kid_number==1) ? 2 : -1; + if (kid_number == 3){ + int i3l3 = Maths::int3log3(point); + if (point == i3l3/3) return -1; + if (point == i3l3-1) return -1; + return (point % 2 == 0) ? 3*point + 3 : 3*point - 1; } - } - double xn = scaleX(point, x); - if (rule != rule_semilocalp) if (max_order == 1) return 1.0 - std::abs(xn); - if (max_order == 2) return evalPWQuadratic(point, xn); - if (max_order == 3) return evalPWCubic(point, xn); - return evalPWPower(point, xn); - } - } - - double evalSupport(int point, double x, bool &isSupported) const override{ - if (isZeroOrder){ - double distance = std::abs(x - getNode(point)), support = getSupport(point); // can speed this up, see above - isSupported = (distance <= (2.0) * support); - return (distance > support) ? 0.0 : 1.0; - }else{ - isSupported = true; - if ((rule == rule_localp) || (rule == rule_semilocalp)){ - if (point == 0) return 1.0; - if (rule == rule_semilocalp){ - if (point == 1) return 0.5 * x * (x - 1.0); - if (point == 2) return 0.5 * x * (x + 1.0); + return 3*point + kid_number; + case erule::localp: + case erule::semilocalp: + if (kid_number == 0){ + switch(point) { + case 0: return 1; + case 1: return 3; + case 2: return 4; + default: + return 2 * point - 1; + }; + } else { + switch(point) { + case 0: return 2; + case 1: + case 2: return -1; + default: + return 2 * point; + }; } - } - double xn = scaleX(point, x); - if (std::abs(xn) <= 1.0){ - if (rule != rule_semilocalp) if (max_order == 1) return 1.0 - std::abs(xn); - if (max_order == 2) return evalPWQuadratic(point, xn); - if (max_order == 3) return evalPWCubic(point, xn); - return evalPWPower(point, xn); - }else{ - isSupported = false; - return 0.0; - } - } + case erule::localp0: + return 2 * point + ((kid_number == 0) ? 1 : 2); + default: // case erule::localpb: + switch(point) { + case 0: + case 1: + return (kid_number == 0) ? 2 : -1; + default: + return 2*point - ((kid_number == 0) ? 1 : 0); + }; + }; } - double diffRaw(int point, double x) const override{ - if (isZeroOrder) { - return 0.0; - } else { - if (rule == rule_localp or rule == rule_semilocalp) { - if (point == 0) return 0.0; - if (rule == rule_semilocalp) { - if (point == 1) return x - 0.5; - if (point == 2) return x + 0.5; - } - } - double xn = scaleX(point, x); - double an = scaleDiffX(point); - if (rule != rule_semilocalp and max_order == 1) { - if (x == 1.0) { - if (rule == rule_localp0 and point == 0) - return -1.0; - else if (rule == rule_localpb and point == 1) - return 0.5; - else if (rule == rule_localpb and point == 2) - return -1.0; - else if (point == 2) - return an; - } - return (xn >= 0 ? -1.0 : 1.0) * an; - } - if (max_order == 2) return an * diffPWQuadratic(point, xn); - if (max_order == 3) return an * diffPWCubic(point, xn); - return an * diffPWPower(point, xn); - } + template + double getNode(int point) { + switch(effective_rule) { + case erule::pwc: + return -2.0 + (1.0 / ((double) Maths::int3log3(point))) * (3*point + 2 - point % 2); + case erule::localp: + case erule::semilocalp: + switch(point) { + case 0: return 0.0; + case 1: return -1.0; + case 2: return 1.0; + default: + return ((double)(2*point - 1)) / ((double) Maths::int2log2(point - 1)) - 3.0; + }; + case erule::localp0: + return ((double)(2*point + 3) ) / ((double) Maths::int2log2(point + 1) ) - 3.0; + default: // case erule::localpb: + switch(point) { + case 0: return -1.0; + case 1: return 1.0; + case 2: return 0.0; + default: + return ((double)(2*point - 1)) / ((double) Maths::int2log2(point - 1)) - 3.0; + }; + }; } - double diffSupport(int point, double x, bool &isSupported) const override{ - if (isZeroOrder) { - isSupported = false; - return 0.0; - } else { - // We still need (isSupported == true) at point = 0 so that TasGrid::GridLocalPolynomial::walkTree<3>() - // walks through the children of this point. - isSupported = true; - if (rule == rule_localp or rule == rule_semilocalp) { - if (point == 0) { - return 0.0; - } - if (rule == rule_semilocalp) { - if (point == 1) return x - 0.5; - if (point == 2) return x + 0.5; - } - } - // Note that the support will be of the form [a, b), except for x = +1.0, i.e., the rightmost node in the domain. This - // is to avoid double-counting derivatives at points of discontinuity. - double xn = scaleX(point, x); - isSupported = (-1.0 <= xn and xn < 1.0) or (x == 1.0 and xn == 1.0); - if (isSupported) { - double an = scaleDiffX(point); - if (rule != rule_semilocalp and max_order == 1) { - // Edge case due to the logic of avoiding double-counting. - if (x == 1.0 and rule == rule_localp and point == 2) return an; - return (xn >= 0 ? -1.0 : 1.0) * an; - } - if (max_order == 2) return an * diffPWQuadratic(point, xn); - if (max_order == 3) return an * diffPWCubic(point, xn); - return an * diffPWPower(point, xn); - } else { - return 0.0; + template + int getLevel(int point) { + switch(effective_rule) { + case erule::pwc: { + int level = 0; + while(point >= 1){ point /= 3; level += 1; } + return level; } - } + case erule::localp: + case erule::semilocalp: + return (point == 0) ? 0 : (point == 1) ? 1 : (Maths::intlog2(point - 1) + 1); + case erule::localp0: + return Maths::intlog2(point + 1); + default: // case erule::localpb: + return (point <= 1) ? 0 : (Maths::intlog2(point - 1) + 1); + }; } - double getArea(int point, std::vector const &w, std::vector const &x) const override{ - if (isZeroOrder){ - return 2.0 * getSupport(point); - }else{ - if (rule == rule_localp){ - if (point == 0) return 2.0; - if ((point == 1) || (point == 2)) return 0.5; - if (max_order == 1) return getSupport(point); - if ((max_order == 2) || (max_order == 3) || (point <= 8)) return (4.0/3.0) * getSupport(point); - }else if (rule == rule_semilocalp){ - if (point == 0) return 2.0; - if ((point == 1) || (point == 2)) return 1.0/3.0; - if ((max_order == 2) || (max_order == 3) || (point <= 4)) return (4.0/3.0) * getSupport(point); - }else if (rule == rule_localpb){ - if ((point == 0) || (point == 1)) return 1.0; - if (max_order == 1) return getSupport(point); - if ((max_order == 2) || (max_order == 3) || (point <= 4)) return (4.0/3.0) * getSupport(point); - }else if (rule == rule_localp0){ - if (max_order == 1) return getSupport(point); - if ((max_order == 2) || (max_order == 3) || (point <= 2)) return (4.0/3.0) * getSupport(point); - } - double sum = 0.0; - for(size_t i=0; i + double getSupport(int point) { + switch(effective_rule) { + case erule::pwc: + return 1.0 / (double) Maths::int3log3(point); + case erule::localp: + case erule::semilocalp: + return (point == 0) ? 1.0 : 1.0 / ((double) Maths::int2log2(point - 1)); + case erule::localp0: + return 1.0 / ((double) Maths::int2log2(point + 1)); + default: // case erule::localpb: + return (point <= 1) ? 2.0 : 1.0 / ((double) Maths::int2log2(point - 1)); + }; } - void van_matrix(int num_rows, std::vector &pntr, std::vector &indx, std::vector &vals) override{ - int max_level = getLevel(num_rows); - - if (isZeroOrder) { - switch(num_rows) { - case 0: - indx = {0,}; - vals = {1.0,}; - pntr = {0, 1}; - break; - case 1: - indx = {0, 0, 1}; - vals = {1.0, 1.0, 1.0}; - pntr = {0, 1, 3}; - break; - case 2: - indx = {0, 0, 1, 0, 2}; - vals = {1.0, 1.0, 1.0, 1.0, 1.0}; - pntr = {0, 1, 3, 5}; - break; - default: - indx.clear(); - indx.reserve(num_rows * max_level); - vals.clear(); - vals.reserve(num_rows * max_level); - pntr = std::vector(num_rows + 1, 0); - for(auto i : std::array{0, 0, 1, 0, 2}) - indx.push_back(i); - for(size_t i=0; i<5; i++) - vals.push_back(1.0); - pntr[1] = 1; - pntr[2] = 3; - { - std::vector ancestors; - ancestors.reserve(max_level); - for(int r=3; r(indx.size()); - ancestors.clear(); - int kid = r; - int dad = kid / 3; - while(dad != 0) { - if (dad % 2 == 0) { - if (kid != 3 * dad) - ancestors.push_back(dad); - } else { - if (kid != 3 * dad + 2) - ancestors.push_back(dad); - } - kid = dad; - dad = kid / 3; - } - indx.push_back(0); - indx.insert(indx.end(), ancestors.rbegin(), ancestors.rend()); - indx.push_back(r); - for(size_t i=0; i(indx.size()); - } - break; - }; - } else { - if (rule == rule_localp) { - switch(num_rows) { - case 0: - indx = {0,}; - vals = {1.0,}; - pntr = {0, 1}; - break; - default: - indx.clear(); - indx.reserve(num_rows * max_level); - vals.clear(); - vals.reserve(num_rows * max_level); - pntr = std::vector(num_rows + 1, 0); - indx.push_back(0); - vals.push_back(1.0); - { - std::vector ancestors; - std::vector ancestors_vals; - ancestors.reserve(max_level); - ancestors_vals.reserve(max_level); - for(int r=1; r(indx.size()); - ancestors.clear(); - ancestors_vals.clear(); - int kid = r; - int dad = (kid + 1) / 2; - if (kid < 4) dad--; - while(dad != 0) { - ancestors.push_back(dad); - ancestors_vals.push_back( evalRaw(dad, x) ); - kid = dad; - dad = (kid + 1) / 2; - if (kid < 4) dad--; - } - indx.push_back(0); - indx.insert(indx.end(), ancestors.rbegin(), ancestors.rend()); - indx.push_back(r); - vals.push_back(1.0); - vals.insert(vals.end(), ancestors_vals.rbegin(), ancestors_vals.rend()); - vals.push_back(1.0); - } - pntr.back() = static_cast(indx.size()); - } - break; - }; - } else if (rule == rule_semilocalp) { - switch(num_rows) { + template + double scaleDiffX(int point) { + switch(effective_rule) { + case erule::pwc: + return 0.0; + case erule::localp: + return (point <= 2) ? 1.0 : static_cast(Maths::int2log2(point - 1)); + case erule::semilocalp: + return static_cast(Maths::int2log2(point - 1)); + case erule::localp0: + return (point == 0) ? 1.0 : static_cast(Maths::int2log2(point + 1)); + default: // case erule::localpb: + switch(point) { case 0: - indx = {0,}; - vals = {1.0,}; - pntr = {0, 1}; - break; - case 1: - indx = {0, 0, 1}; - vals = {1.0, 1.0, 1.0}; - pntr = {0, 1, 3}; - break; - case 2: - indx = {0, 0, 1, 0, 2}; - vals = {1.0, 1.0, 1.0, 1.0, 1.0}; - pntr = {0, 1, 3, 5}; - break; + case 1: return 0.5; + case 2: return 1.0; default: - indx.clear(); - indx.reserve(num_rows * max_level); - vals.clear(); - vals.reserve(num_rows * max_level); - pntr = std::vector(num_rows + 1, 0); - for(auto i : std::array{0, 0, 1, 0, 2}) - indx.push_back(i); - for(int i=0; i<5; i++) - vals.push_back(1.0); - pntr[1] = 1; - pntr[2] = 3; - { - std::vector ancestors; - std::vector ancestors_vals; - ancestors.reserve(max_level); - ancestors_vals.reserve(max_level); - for(int r=3; r(indx.size()); - double x = getNode(r); - ancestors.clear(); - ancestors_vals.clear(); - int kid = r; - int dad = (kid + 1) / 2; - while(dad > 2) { - ancestors.push_back(dad); - ancestors_vals.push_back( evalRaw(dad, x) ); - kid = dad; - dad = (kid + 1) / 2; - } - indx.push_back(0); - indx.push_back(1); - indx.push_back(2); - indx.insert(indx.end(), ancestors.rbegin(), ancestors.rend()); - indx.push_back(r); - vals.push_back(1.0); - vals.push_back( evalRaw(1, x) ); - vals.push_back( evalRaw(2, x) ); - vals.insert(vals.end(), ancestors_vals.rbegin(), ancestors_vals.rend()); - vals.push_back(1.0); - } - pntr.back() = static_cast(indx.size()); - } - break; + return static_cast(Maths::int2log2(point - 1)); }; - } else if (rule == rule_localp0) { - switch(num_rows) { - case 0: - indx = {0,}; - vals = {1.0,}; - pntr = {0, 1}; - break; + }; + } + + template + double scaleX(int point, double x) { + switch(effective_rule) { + case erule::pwc: + return 0.0; // should never be called + case erule::localp: + switch(point) { + case 0: return x; + case 1: return (x + 1.0); + case 2: return (x - 1.0); default: - indx.clear(); - indx.reserve(num_rows * max_level); - vals.clear(); - vals.reserve(num_rows * max_level); - pntr = std::vector(num_rows + 1, 0); - indx.push_back(0); - vals.push_back(1.0); - { - std::vector ancestors; - std::vector ancestors_vals; - ancestors.reserve(max_level); - ancestors_vals.reserve(max_level); - for(int r=1; r(indx.size()); - ancestors.clear(); - ancestors_vals.clear(); - int kid = r; - int dad = (kid - 1) / 2; - while(dad != 0) { - ancestors.push_back(dad); - ancestors_vals.push_back( evalRaw(dad, x) ); - kid = dad; - dad = (kid - 1) / 2; - } - indx.push_back(0); - indx.insert(indx.end(), ancestors.rbegin(), ancestors.rend()); - indx.push_back(r); - vals.push_back( evalRaw(0, x) ); - vals.insert(vals.end(), ancestors_vals.rbegin(), ancestors_vals.rend()); - vals.push_back(1.0); - } - pntr.back() = static_cast(indx.size()); - } - break; + return ((double) Maths::int2log2(point - 1) * (x + 3.0) + 1.0 - (double) (2*point)); }; - } else { // if (rule == rule_localpb) { - switch(num_rows) { - case 0: - indx = {0,}; - vals = {1.0,}; - pntr = {0, 1}; - break; - case 1: - indx = {0, 1}; - vals = {1.0, 1.0}; - pntr = {0, 1, 2}; - break; + case erule::semilocalp: + return ((double) Maths::int2log2(point - 1) * (x + 3.0) + 1.0 - (double) (2*point)); + case erule::localp0: + return ((double) Maths::int2log2(point + 1) * (x + 3.0) - 3.0 - (double) (2*point)); + default: // case erule::localpb: + switch(point) { + case 0: return (x + 1.0) / 2.0; + case 1: return (x - 1.0) / 2.0; + case 2: return x; default: - indx.clear(); - indx.reserve(num_rows * max_level); - vals.clear(); - vals.reserve(num_rows * max_level); - pntr = std::vector(num_rows + 1, 0); - indx.push_back(0); - indx.push_back(1); - vals.push_back(1.0); - vals.push_back(1.0); - pntr[1] = 1; - { - std::vector ancestors; - std::vector ancestors_vals; - ancestors.reserve(max_level); - ancestors_vals.reserve(max_level); - for(int r=2; r(indx.size()); - double x = getNode(r); - ancestors.clear(); - ancestors_vals.clear(); - int kid = r; - int dad = (kid + 1) / 2; - while(dad > 1) { - ancestors.push_back(dad); - ancestors_vals.push_back( evalRaw(dad, x) ); - kid = dad; - dad = (kid + 1) / 2; - } - indx.push_back(0); - indx.push_back(1); - indx.insert(indx.end(), ancestors.rbegin(), ancestors.rend()); - indx.push_back(r); - vals.push_back( evalRaw(0, x) ); - vals.push_back( evalRaw(1, x) ); - vals.insert(vals.end(), ancestors_vals.rbegin(), ancestors_vals.rend()); - vals.push_back(1.0); - } - pntr.back() = static_cast(indx.size()); - } - break; + return ((double) Maths::int2log2(point - 1) * (x + 3.0) + 1.0 - (double) (2*point)); }; - } - } - } - -protected: - double scaleX(int point, double x) const{ - if (rule == rule_localp0){ - if (point == 0) return x; - return ((double) Maths::int2log2(point + 1) * (x + 3.0) - 3.0 - (double) (2*point)); - }else{ - if (rule == rule_localp){ - if (point == 0) return x; - if (point == 1) return (x + 1.0); - if (point == 2) return (x - 1.0); - }else if (rule == rule_localpb){ - if (point == 0) return (x + 1.0) / 2.0; - if (point == 1) return (x - 1.0) / 2.0; - if (point == 2) return x; - } - return ((double) Maths::int2log2(point - 1) * (x + 3.0) + 1.0 - (double) (2*point)); - } - } - - double scaleDiffX(int point) const{ - if (rule == rule_localp0) { - if (point == 0) return 1.0; - return (double) Maths::int2log2(point + 1); - } else { - if (rule == rule_localp) { - if (point == 0 or point == 1 or point == 2) return 1.0; - } else if (rule == rule_localpb) { - if (point == 0 or point == 1) return 0.5; - if (point == 2) return 1.0; - } - return (double) Maths::int2log2(point - 1); - } + }; } - double evalPWQuadratic(int point, double x) const{ - if (rule == rule_localp){ - if (point == 1) return 1.0 - x; - if (point == 2) return 1.0 + x; - }else if (rule == rule_localpb){ - if (point == 0) return 1.0 - x; - if (point == 1) return 1.0 + x; + template + double evalPWQuadratic(int point, double x) { + if (effective_rule == erule::localp){ + switch(point) { + case 1: return 1.0 - x; + case 2: return 1.0 + x; + default: + return (1.0 - x) * (1.0 + x); + }; + }else if (effective_rule == erule::localpb){ + switch(point) { + case 0: return 1.0 - x; + case 1: return 1.0 + x; + default: + return (1.0 - x) * (1.0 + x); + }; } return (1.0 - x) * (1.0 + x); } - double evalPWCubic(int point, double x) const{ - if (rule == rule_localp){ - if (point == 0) return 1.0; - if (point == 1) return 1.0 - x; - if (point == 2) return 1.0 + x; - if (point <= 4) return (1.0 - x) * (1.0 + x); - }else if (rule == rule_localpb){ - if (point == 0) return 1.0 - x; - if (point == 1) return 1.0 + x; - if (point == 2) return (1.0 - x) * (1.0 + x); - }else if (rule == rule_localp0){ + template + double evalPWCubic(int point, double x) { + if (effective_rule == erule::localp){ + switch(point) { + case 0: return 1.0; + case 1: return 1.0 - x; + case 2: return 1.0 + x; + case 3: + case 4: return (1.0 - x) * (1.0 + x); + default: + return (point % 2 == 0) ? (1.0 - x) * (1.0 + x) * (3.0 + x) / 3.0 : (1.0 - x) * (1.0 + x) * (3.0 - x) / 3.0; + }; + }else if (effective_rule == erule::localp0){ if (point == 0) return (1.0 - x) * (1.0 + x); + }else if (effective_rule == erule::localpb){ + switch(point) { + case 0: return 1.0 - x; + case 1: return 1.0 + x; + case 2: return (1.0 - x) * (1.0 + x); + default: + return (point % 2 == 0) ? (1.0 - x) * (1.0 + x) * (3.0 + x) / 3.0 : (1.0 - x) * (1.0 + x) * (3.0 - x) / 3.0; + }; } return (point % 2 == 0) ? (1.0 - x) * (1.0 + x) * (3.0 + x) / 3.0 : (1.0 - x) * (1.0 + x) * (3.0 - x) / 3.0; } - double evalPWPower(int point, double x) const{ - if (rule == rule_localp) if (point <= 8) return evalPWCubic(point, x); // if order is cubic or less, use the hard-coded functions - if (rule == rule_semilocalp) if (point <= 4) return evalPWCubic(point, x); - if (rule == rule_localpb) if (point <= 4) return evalPWCubic(point, x); - if (rule == rule_localp0) if (point <= 2) return evalPWCubic(point, x); - int level = getLevel(point); + + template + double evalPWPower(int max_order, int point, double x) { + // use the cubic implementation until we have enough points + if (effective_rule == erule::localp) if (point <= 8) return evalPWCubic(point, x); // if order is cubic or less, use the hard-coded functions + if (effective_rule == erule::semilocalp) if (point <= 4) return evalPWCubic(point, x); + if (effective_rule == erule::localpb) if (point <= 4) return evalPWCubic(point, x); + if (effective_rule == erule::localp0) if (point <= 2) return evalPWCubic(point, x); + int level = getLevel(point); int most_turns = 1; - double value = (1.0 - x)*(1.0 + x), phantom_distance = 1.0; - int max_ancestors; // maximum number of ancestors to consider, first set the possible max then constrain by max_order - if (rule == rule_localp) max_ancestors = level-2; // number of ancestors to consider (two counted in the initialization of value) - if (rule == rule_semilocalp) max_ancestors = level-1; // semi-localp and localpb add one more ancestor due to the global support at levels 1 and 0 (respectively) - if (rule == rule_localpb) max_ancestors = level-1; - if (rule == rule_localp0) max_ancestors = level; // localp0 adds two ancestors, the assumed zero nodes at the boundary + double value = (1.0 - x)*(1.0 + x); + double phantom_distance = 1.0; + int max_ancestors = [&]()->int{ // maximum number of ancestors to consider, first set the possible max then constrain by max_order + switch(effective_rule) { + case erule::pwc: return 0; // should never happen + case erule::localp: return max_ancestors = level-2; + case erule::semilocalp: return max_ancestors = level-1; + case erule::localpb: return max_ancestors = level-1; + default: return max_ancestors = level; + }; + }(); if (max_order > 0) max_ancestors = std::min(max_ancestors, max_order - 2); // use the minimum of the available ancestors or the order restriction for(int j=0; j < max_ancestors; j++){ @@ -722,50 +365,206 @@ class templRuleLocalPolynomial : public BaseRuleLocalPolynomial{ // Every time we turn, we backtrack and we lose 2 units, the phantom node is at maximum distance minus 2 time the number of turns (to the left or right) most_turns *= 2; phantom_distance = 2.0 * phantom_distance + 1.0; - int turns = (rule == rule_localp0) ? ((point+1) % most_turns) : ((point-1) % most_turns); + int turns = (effective_rule == erule::localp0) ? ((point+1) % most_turns) : ((point-1) % most_turns); double node = (turns < most_turns / 2) ? (phantom_distance - 2.0 * ((double) turns)) : (-phantom_distance + 2.0 * ((double) (most_turns - 1 - turns))); - value *= ( x - node ) / ( - node); + value *= - ( x - node ) / node; } return value; } - // Same logic as the corresponding eval functions. - double diffPWQuadratic(int point, double x) const{ - if (rule == rule_localp) { - if (point == 1) return -1.0; - if (point == 2) return 1.0; - } else if (rule == rule_localpb) { - if (point == 0) return -1.0; - if (point == 1) return 1.0; + template + double evalSupport(int max_order, int point, double x, bool &isSupported) { + switch(effective_rule) { + case erule::pwc: { + double distance = std::abs(x - getNode(point)); + double support = getSupport(point); + isSupported = (distance <= 2.0 * support); + return (distance <= support) ? 1.0 : 0.0; + } + case erule::localp: + isSupported = true; + if (point == 0) { + return 1.0; + } else { + double xn = scaleX(point, x); + if (std::abs(xn) <= 1.0) { + switch(max_order) { + case 1: return 1.0 - std::abs(xn); + case 2: return evalPWQuadratic(point, xn); + case 3: return evalPWCubic(point, xn); + default: + return evalPWPower(max_order, point, xn); + }; + } else { + isSupported = false; + return 0.0; + } + } + case erule::semilocalp: + isSupported = true; + switch(point) { + case 0: return 1.0; + case 1: return 0.5 * x * (x - 1.0); + case 2: return 0.5 * x * (x + 1.0); + default: { + double xn = scaleX(point, x); + if (std::abs(xn) <= 1.0) { + switch(max_order) { + case 1: return 1.0 - std::abs(xn); + case 2: return evalPWQuadratic(point, xn); + case 3: return evalPWCubic(point, xn); + default: + return evalPWPower(max_order, point, xn); + }; + } else { + isSupported = false; + return 0.0; + } + } + }; + case erule::localp0: + default: { // case erule::localpb: + double xn = scaleX(point, x); + if (std::abs(xn) <= 1.0) { + isSupported = true; + switch(max_order) { + case 1: return 1.0 - std::abs(xn); + case 2: return evalPWQuadratic(point, xn); + case 3: return evalPWCubic(point, xn); + default: + return evalPWPower(max_order, point, xn); + }; + } else { + isSupported = false; + return 0.0; + } + } + }; + } + + template + double evalRaw(int max_order, int point, double x) { + switch(effective_rule) { + case erule::pwc: + return (std::abs(x - getNode(point)) <= getSupport(point)) ? 1.0 : 0.0; + case erule::localp: + if (point == 0) { + return 1.0; + } else { + double xn = scaleX(point, x); + if (std::abs(xn) <= 1.0) { + switch(max_order) { + case 1: return 1.0 - std::abs(xn); + case 2: return evalPWQuadratic(point, xn); + case 3: return evalPWCubic(point, xn); + default: + return evalPWPower(max_order, point, xn); + }; + } else { + return 0.0; + } + } + case erule::semilocalp: + switch(point) { + case 0: return 1.0; + case 1: return 0.5 * x * (x - 1.0); + case 2: return 0.5 * x * (x + 1.0); + default: { + double xn = scaleX(point, x); + if (std::abs(xn) <= 1.0) { + switch(max_order) { + case 1: return 1.0 - std::abs(xn); + case 2: return evalPWQuadratic(point, xn); + case 3: return evalPWCubic(point, xn); + default: + return evalPWPower(max_order, point, xn); + }; + } else { + return 0.0; + } + } + }; + case erule::localp0: + default: { // case erule::localpb: + double xn = scaleX(point, x); + if (std::abs(xn) <= 1.0) { + switch(max_order) { + case 1: return 1.0 - std::abs(xn); + case 2: return evalPWQuadratic(point, xn); + case 3: return evalPWCubic(point, xn); + default: + return evalPWPower(max_order, point, xn); + }; + } else { + return 0.0; + } + } + }; + } + + template + double diffPWQuadratic(int point, double x) { + if (effective_rule == erule::localp) { + switch(point) { + case 1: return -1.0; + case 2: return 1.0; + default: + return -2.0 * x; + }; + } else if (effective_rule == erule::localpb) { + switch(point) { + case 0: return -1.0; + case 1: return 1.0; + default: + return -2.0 * x; + }; } return -2.0 * x; } - double diffPWCubic(int point, double x) const{ - if (rule == rule_localp) { - if (point == 0) return 0.0; - if (point == 1) return -1.0; - if (point == 2) return 1.0; - if (point <= 4) return -2.0 * x; - } else if (rule == rule_localpb) { - if (point == 0) return -1.0; - if (point == 1) return 1.0; - if (point == 2) return -2.0 * x; - } else if (rule == rule_localp0) { - if (point == 0) return -2.0 * x; + template + double diffPWCubic(int point, double x) { + if (effective_rule == erule::localp) { + switch(point) { + case 0: return 0.0; + case 1: return -1.0; + case 2: return 1.0; + case 3: + case 4: return -2.0 * x; + default: + return (point % 2 == 0) ? 1.0 / 3.0 - x * (x + 2.0) : -1.0 / 3.0 + x * (x - 2.0); + }; + } else if (effective_rule == erule::localpb) { + switch(point) { + case 0: return -1.0; + case 1: return 1.0; + case 2: return -2.0 * x; + default: + return (point % 2 == 0) ? 1.0 / 3.0 - x * (x + 2.0) : -1.0 / 3.0 + x * (x - 2.0); + }; + } else if (effective_rule == erule::localp0) { + return (point == 0) ? -2.0 * x : ((point % 2 == 0) ? 1.0 / 3.0 - x * (x + 2.0) : -1.0 / 3.0 + x * (x - 2.0)); } return (point % 2 == 0) ? 1.0 / 3.0 - x * (x + 2.0) : -1.0 / 3.0 + x * (x - 2.0); } - double diffPWPower(int point, double x) const{ - if (rule == rule_localp and point <= 8) return diffPWCubic(point, x); - if (rule == rule_semilocalp and point <= 4) return diffPWCubic(point, x); - if (rule == rule_localpb and point <= 4) return diffPWCubic(point, x); - if (rule == rule_localp0 and point <= 2) return diffPWCubic(point, x); - int level = getLevel(point); - int max_ancestors; - if (rule == rule_localp) max_ancestors = level-2; - if (rule == rule_semilocalp) max_ancestors = level-1; - if (rule == rule_localpb) max_ancestors = level-1; - if (rule == rule_localp0) max_ancestors = level; + + template + double diffPWPower(int max_order, int point, double x) { + // use the cubic implementation until we have enough points + if (effrule == erule::localp and point <= 8) return diffPWCubic(point, x); + if (effrule == erule::semilocalp and point <= 4) return diffPWCubic(point, x); + if (effrule == erule::localpb and point <= 4) return diffPWCubic(point, x); + if (effrule == erule::localp0 and point <= 2) return diffPWCubic(point, x); + int level = getLevel(point); + int max_ancestors = [&]()->int { + switch(effrule) { + case erule::pwc: return 0; + case erule::localp: return level - 2; + case erule::semilocalp: + case erule::localpb: return level - 1; + default: return level; + }; + }(); + if (max_order > 0) max_ancestors = std::min(max_ancestors, max_order - 2); // This lambda captures most_turns and phantom_distance by reference, uses those as internal state variables, and on each @@ -775,7 +574,7 @@ class templRuleLocalPolynomial : public BaseRuleLocalPolynomial{ auto update_and_get_next_node = [&]() { most_turns *= 2; phantom_distance = 2.0 * phantom_distance + 1.0; - int turns = (rule == rule_localp0) ? ((point+1) % most_turns) : ((point-1) % most_turns); + int turns = (effrule == erule::localp0) ? ((point+1) % most_turns) : ((point-1) % most_turns); return (turns < most_turns / 2) ? (phantom_distance - 2.0 * ((double) turns)) : (-phantom_distance + 2.0 * ((double) (most_turns - 1 - turns))); @@ -785,7 +584,7 @@ class templRuleLocalPolynomial : public BaseRuleLocalPolynomial{ auto rollback_and_get_prev_node = [&]() { most_turns /= 2; phantom_distance = 0.5 * (phantom_distance - 1.0); - int turns = (rule == rule_localp0) ? ((point+1) % most_turns) : ((point-1) % most_turns); + int turns = (effrule == erule::localp0) ? ((point+1) % most_turns) : ((point-1) % most_turns); return (turns < most_turns / 2) ? (phantom_distance - 2.0 * ((double) turns)) : (-phantom_distance + 2.0 * ((double) (most_turns - 1 - turns))); @@ -815,30 +614,485 @@ class templRuleLocalPolynomial : public BaseRuleLocalPolynomial{ return derivative; } -}; - -inline std::unique_ptr makeRuleLocalPolynomial(TypeOneDRule rule, int order){ - if (order == 0){ - return Utils::make_unique>(); - }else{ - std::unique_ptr rule_pntr = [=]()->std::unique_ptr{ - if (rule == rule_localp){ - return Utils::make_unique>(); - }else if (rule == rule_semilocalp){ - return Utils::make_unique>(); - }else if (rule == rule_localp0){ - return Utils::make_unique>(); - }else{ // must be (rule == rule_localpb) - return Utils::make_unique>(); + + template + double diffRaw(int max_order, int point, double x) { + switch(effrule) { + case erule::pwc: return 0.0; + case erule::localp: { + if (point == 0) return 0.0; + double xn = scaleX(point, x); + double an = scaleDiffX(point); + switch(max_order) { + case 1: return ((xn >= 0 ? -1.0 : 1.0) * an); + case 2: return an * diffPWQuadratic(point, xn); + case 3: return an * diffPWCubic(point, xn); + default: + return an * diffPWPower(point, xn); + }; } - }(); - rule_pntr->setMaxOrder(order); - return rule_pntr; + case erule::semilocalp: + switch(point) { + case 0: return 0.0; + case 1: return x - 0.5; + case 2: return x + 0.5; + default: { + double xn = scaleX(point, x); + double an = scaleDiffX(point); + switch(max_order) { + case 2: return an * diffPWQuadratic(point, xn); + case 3: return an * diffPWCubic(point, xn); + default: + return an * diffPWPower(point, xn); + }; + } + }; + case erule::localp0: { + double xn = scaleX(point, x); + double an = scaleDiffX(point); + switch(max_order) { + case 1: return (x == 1.0 and point == 0) ? -1.0 : ((xn >= 0 ? -1.0 : 1.0) * an); + case 2: return an * diffPWQuadratic(point, xn); + case 3: return an * diffPWCubic(point, xn); + default: + return an * diffPWPower(point, xn); + }; + } + default: { // case erule::localpb: + double xn = scaleX(point, x); + double an = scaleDiffX(point); + switch(max_order) { + case 1: return ((xn >= 0 ? -1.0 : 1.0) * an); + case 2: return an * diffPWQuadratic(point, xn); + case 3: return an * diffPWCubic(point, xn); + default: + return an * diffPWPower(point, xn); + }; + } + }; } -} -#endif // __TASMANIAN_DOXYGEN_SKIP + template + double diffSupport(int max_order, int point, double x, bool &isSupported) { + switch(effrule) { + case erule::pwc: + isSupported = false; + return 0.0; + case erule::localp: { + if (point == 0) { isSupported = true; return 0.0; } + double xn = scaleX(point, x); + isSupported = (-1.0 <= xn and xn < 1.0) or (x == 1.0 and xn == 1.0); + if (isSupported) { + double an = scaleDiffX(point); + switch(max_order) { + case 1: return (x == 1.0 and point == 2) ? an : (((xn >= 0 ? -1.0 : 1.0) * an)); + case 2: return an * diffPWQuadratic(point, xn); + case 3: return an * diffPWCubic(point, xn); + default: + return an * diffPWPower(max_order, point, xn); + }; + } else { + return 0.0; + } + } + case erule::semilocalp: + switch(point) { + case 0: { isSupported = true; return 0.0; } + case 1: { isSupported = true; return x - 0.5; } + case 2: { isSupported = true; return x + 0.5; } + default: { + double xn = scaleX(point, x); + isSupported = (-1.0 <= xn and xn < 1.0) or (x == 1.0 and xn == 1.0); + if (isSupported) { + double an = scaleDiffX(point); + switch(max_order) { + case 2: return an * diffPWQuadratic(point, xn); + case 3: return an * diffPWCubic(point, xn); + default: + return an * diffPWPower(max_order, point, xn); + }; + } else { + return 0.0; + } + } + }; + case erule::localp0: { + double xn = scaleX(point, x); + isSupported = (-1.0 <= xn and xn < 1.0) or (x == 1.0 and xn == 1.0); + if (isSupported) { + double an = scaleDiffX(point); + switch(max_order) { + case 1: return (x == 1.0 and point == 0) ? -1.0 : ((xn >= 0 ? -1.0 : 1.0) * an); + case 2: return an * diffPWQuadratic(point, xn); + case 3: return an * diffPWCubic(point, xn); + default: + return an * diffPWPower(max_order, point, xn); + }; + } else { + return 0.0; + } + } + default: { // case erule::localpb: + double xn = scaleX(point, x); + isSupported = (-1.0 <= xn and xn < 1.0) or (x == 1.0 and xn == 1.0); + if (isSupported) { + double an = scaleDiffX(point); + switch(max_order) { + case 1: return ((xn >= 0 ? -1.0 : 1.0) * an); + case 2: return an * diffPWQuadratic(point, xn); + case 3: return an * diffPWCubic(point, xn); + default: + return an * diffPWPower(max_order, point, xn); + }; + } else { + return 0.0; + } + } + }; + } + + template + double getArea(int max_order, int point, std::vector const &w, std::vector const &x) { + switch(effrule) { + case erule::pwc: + return 2.0 * getSupport(point); + case erule::localp: + switch(point) { + case 0: return 2.0; + case 1: + case 2: return 0.5; + default: + switch(max_order) { + case 1: return getSupport(point); + case 2: + case 3: return (4.0/3.0) * getSupport(point); + default: + if (point <= 8) return (4.0/3.0) * getSupport(point); + break; + }; + break; + }; + break; + case erule::semilocalp: + switch(point) { + case 0: return 2.0; + case 1: + case 2: return 1.0/3.0; + default: + switch(max_order) { + case 2: + case 3: return (4.0/3.0) * getSupport(point); + default: + if (point <= 4) return (4.0/3.0) * getSupport(point); + break; + }; + break; + }; + break; + case erule::localp0: + switch(max_order) { + case 1: return getSupport(point); + case 2: + case 3: return (4.0/3.0) * getSupport(point); + default: + if (point <= 2) return (4.0/3.0) * getSupport(point); + break; + }; + break; + default: // case erule::localpb: + if (point <= 1) { + return 1.0; + } else { + switch(max_order) { + case 1: return getSupport(point); + case 2: + case 3: return (4.0/3.0) * getSupport(point); + default: + if (point <= 4) return (4.0/3.0) * getSupport(point); + break; + }; + } + break; + }; + double sum = 0.0; + for(size_t i=0; i(max_order, point, x[i]); + return sum * getSupport(point); + } + template + void van_matrix(int max_order, int num_rows, std::vector &pntr, std::vector &indx, std::vector &vals) { + int max_level = getLevel(num_rows); + + if (effrule == erule::pwc) { + switch(num_rows) { + case 0: + indx = {0,}; + vals = {1.0,}; + pntr = {0, 1}; + break; + case 1: + indx = {0, 0, 1}; + vals = {1.0, 1.0, 1.0}; + pntr = {0, 1, 3}; + break; + case 2: + indx = {0, 0, 1, 0, 2}; + vals = {1.0, 1.0, 1.0, 1.0, 1.0}; + pntr = {0, 1, 3, 5}; + break; + default: + indx.clear(); + indx.reserve(num_rows * max_level); + vals.clear(); + vals.reserve(num_rows * max_level); + pntr = std::vector(num_rows + 1, 0); + for(auto i : std::array{0, 0, 1, 0, 2}) + indx.push_back(i); + for(size_t i=0; i<5; i++) + vals.push_back(1.0); + pntr[1] = 1; + pntr[2] = 3; + { + std::vector ancestors; + ancestors.reserve(max_level); + for(int r=3; r(indx.size()); + ancestors.clear(); + int kid = r; + int dad = kid / 3; + while(dad != 0) { + if (dad % 2 == 0) { + if (kid != 3 * dad) + ancestors.push_back(dad); + } else { + if (kid != 3 * dad + 2) + ancestors.push_back(dad); + } + kid = dad; + dad = kid / 3; + } + indx.push_back(0); + indx.insert(indx.end(), ancestors.rbegin(), ancestors.rend()); + indx.push_back(r); + for(size_t i=0; i(indx.size()); + } + break; + }; + } else if (effrule == erule::localp) { + switch(num_rows) { + case 0: + indx = {0,}; + vals = {1.0,}; + pntr = {0, 1}; + break; + default: + indx.clear(); + indx.reserve(num_rows * max_level); + vals.clear(); + vals.reserve(num_rows * max_level); + pntr = std::vector(num_rows + 1, 0); + indx.push_back(0); + vals.push_back(1.0); + { + std::vector ancestors; + std::vector ancestors_vals; + ancestors.reserve(max_level); + ancestors_vals.reserve(max_level); + for(int r=1; r(r); + pntr[r] = static_cast(indx.size()); + ancestors.clear(); + ancestors_vals.clear(); + int kid = r; + int dad = (kid + 1) / 2; + if (kid < 4) dad--; + while(dad != 0) { + ancestors.push_back(dad); + ancestors_vals.push_back( evalRaw(max_order, dad, x) ); + kid = dad; + dad = (kid + 1) / 2; + if (kid < 4) dad--; + } + indx.push_back(0); + indx.insert(indx.end(), ancestors.rbegin(), ancestors.rend()); + indx.push_back(r); + vals.push_back(1.0); + vals.insert(vals.end(), ancestors_vals.rbegin(), ancestors_vals.rend()); + vals.push_back(1.0); + } + pntr.back() = static_cast(indx.size()); + } + break; + }; + } else if (effrule == erule::semilocalp) { + switch(num_rows) { + case 0: + indx = {0,}; + vals = {1.0,}; + pntr = {0, 1}; + break; + case 1: + indx = {0, 0, 1}; + vals = {1.0, 1.0, 1.0}; + pntr = {0, 1, 3}; + break; + case 2: + indx = {0, 0, 1, 0, 2}; + vals = {1.0, 1.0, 1.0, 1.0, 1.0}; + pntr = {0, 1, 3, 5}; + break; + default: + indx.clear(); + indx.reserve(num_rows * max_level); + vals.clear(); + vals.reserve(num_rows * max_level); + pntr = std::vector(num_rows + 1, 0); + for(auto i : std::array{0, 0, 1, 0, 2}) + indx.push_back(i); + for(int i=0; i<5; i++) + vals.push_back(1.0); + pntr[1] = 1; + pntr[2] = 3; + { + std::vector ancestors; + std::vector ancestors_vals; + ancestors.reserve(max_level); + ancestors_vals.reserve(max_level); + for(int r=3; r(indx.size()); + double x = getNode(r); + ancestors.clear(); + ancestors_vals.clear(); + int kid = r; + int dad = (kid + 1) / 2; + while(dad > 2) { + ancestors.push_back(dad); + ancestors_vals.push_back( evalRaw(max_order, dad, x) ); + kid = dad; + dad = (kid + 1) / 2; + } + indx.push_back(0); + indx.push_back(1); + indx.push_back(2); + indx.insert(indx.end(), ancestors.rbegin(), ancestors.rend()); + indx.push_back(r); + vals.push_back(1.0); + vals.push_back( evalRaw(max_order, 1, x) ); + vals.push_back( evalRaw(max_order, 2, x) ); + vals.insert(vals.end(), ancestors_vals.rbegin(), ancestors_vals.rend()); + vals.push_back(1.0); + } + pntr.back() = static_cast(indx.size()); + } + break; + }; + } else if (effrule == erule::localp0) { + switch(num_rows) { + case 0: + indx = {0,}; + vals = {1.0,}; + pntr = {0, 1}; + break; + default: + indx.clear(); + indx.reserve(num_rows * max_level); + vals.clear(); + vals.reserve(num_rows * max_level); + pntr = std::vector(num_rows + 1, 0); + indx.push_back(0); + vals.push_back(1.0); + { + std::vector ancestors; + std::vector ancestors_vals; + ancestors.reserve(max_level); + ancestors_vals.reserve(max_level); + for(int r=1; r(r); + pntr[r] = static_cast(indx.size()); + ancestors.clear(); + ancestors_vals.clear(); + int kid = r; + int dad = (kid - 1) / 2; + while(dad != 0) { + ancestors.push_back(dad); + ancestors_vals.push_back( evalRaw(max_order, dad, x) ); + kid = dad; + dad = (kid - 1) / 2; + } + indx.push_back(0); + indx.insert(indx.end(), ancestors.rbegin(), ancestors.rend()); + indx.push_back(r); + vals.push_back( evalRaw(max_order, 0, x) ); + vals.insert(vals.end(), ancestors_vals.rbegin(), ancestors_vals.rend()); + vals.push_back(1.0); + } + pntr.back() = static_cast(indx.size()); + } + break; + }; + } else { // if (effrule == erule::localpb) { + switch(num_rows) { + case 0: + indx = {0,}; + vals = {1.0,}; + pntr = {0, 1}; + break; + case 1: + indx = {0, 1}; + vals = {1.0, 1.0}; + pntr = {0, 1, 2}; + break; + default: + indx.clear(); + indx.reserve(num_rows * max_level); + vals.clear(); + vals.reserve(num_rows * max_level); + pntr = std::vector(num_rows + 1, 0); + indx.push_back(0); + indx.push_back(1); + vals.push_back(1.0); + vals.push_back(1.0); + pntr[1] = 1; + { + std::vector ancestors; + std::vector ancestors_vals; + ancestors.reserve(max_level); + ancestors_vals.reserve(max_level); + for(int r=2; r(indx.size()); + double x = getNode(r); + ancestors.clear(); + ancestors_vals.clear(); + int kid = r; + int dad = (kid + 1) / 2; + while(dad > 1) { + ancestors.push_back(dad); + ancestors_vals.push_back( evalRaw(max_order, dad, x) ); + kid = dad; + dad = (kid + 1) / 2; + } + indx.push_back(0); + indx.push_back(1); + indx.insert(indx.end(), ancestors.rbegin(), ancestors.rend()); + indx.push_back(r); + vals.push_back( evalRaw(max_order, 0, x) ); + vals.push_back( evalRaw(max_order, 1, x) ); + vals.insert(vals.end(), ancestors_vals.rbegin(), ancestors_vals.rend()); + vals.push_back(1.0); + } + pntr.back() = static_cast(indx.size()); + } + break; + }; + } + } + +} // namespace RuleLocal + +#endif // __TASMANIAN_DOXYGEN_SKIP }