Skip to content

Commit

Permalink
Merge pull request #760 from mkstoyanov/improve_updatemap
Browse files Browse the repository at this point in the history 
* adding a few improvements to the refinement algorithm
  • Loading branch information
@mkstoyanov
mkstoyanov authored Mar 21, 2024
2 parents 98b7d12 + 2206070 commit 1de7210
Show file tree
Hide file tree
Showing 6 changed files with 150 additions and 88 deletions.
21 changes: 12 additions & 9 deletions Addons/testLoadUnstructured.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -88,9 +88,10 @@ inline double coefficientDifference(TasmanianSparseGrid const &gA, TasmanianSpar
size_t num_coeffs = Utils::size_mult(gA.getNumOutputs(), gA.getNumPoints());
double const *c1 = gA.getHierarchicalCoefficients();
double const *c2 = gB.getHierarchicalCoefficients();
return std::inner_product(c1, c1 + num_coeffs, c2, 0.0,
[](double a, double b)->double{ return std::max(a, b); },
[](double a, double b)->double{ return std::abs(a - b); });
double err = 0.0;
for(size_t i=0; i<num_coeffs; i++)
err = std::max(std::abs(c1[i] - c2[i]), err);
return err;
}

/*!
Expand All @@ -102,16 +103,18 @@ inline double evalDifference(std::vector<double> const&x, TasmanianSparseGrid co
gA.evaluateBatch(x, yA);
gB.evaluateBatch(x, yB);

return std::inner_product(yA.begin(), yA.end(), yB.begin(), 0.0,
[](double a, double b)->double{ return std::max(a, b); },
[](double a, double b)->double{ return std::abs(a - b); });
double err = 0.0;
for(size_t i=0; i<yA.size(); i++)
err = std::max(std::abs(yA[i] - yB[i]), err);
return err;
}

inline double vecDifference(std::vector<double> const &x, std::vector<double> const &y){
if (x.size() != y.size()) return 1.E+20;
return std::inner_product(x.begin(), x.end(), y.begin(), 0.0,
[](double a, double b)->double{ return std::max(a, b); },
[](double a, double b)->double{ return std::abs(a - b); });
double err = 0.0;
for(size_t i=0; i<x.size(); i++)
err = std::max(std::abs(x[i] - y[i]), err);
return err;
}

/*!
Expand Down
9 changes: 5 additions & 4 deletions SparseGrids/gridtestExternalTests.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1041,10 +1041,11 @@ bool ExternalTester::testDynamicRefinement(const BaseFunction *f, TasmanianSpars
std::vector<double> xpnts = genRandom(10, grid.getNumDimensions());
grid.evaluateBatch(xpnts, res1);
grid2.evaluateBatch(xpnts, res2);
double err = std::inner_product(res1.begin(), res1.end(), res2.begin(), 0.0,
[](double a, double b)->double{ return std::max(a, b); },
[](double a, double b)->double{ return std::abs(a - b); });
if (err > Maths::num_tol){ cout << "ERROR: failed evaluate after loading batch points." << endl; return false; }
double err = 0.0;
for(size_t i=0; i<res1.size(); i++)
err = std::max(std::abs(res1[i] - res2[i]), err);

if (err > Maths::num_tol){ cout << "ERROR: failed evaluate after loading batch points. " << endl; return false; }
return true;
}

Expand Down
192 changes: 119 additions & 73 deletions SparseGrids/tsgGridLocalPolynomial.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1266,83 +1266,117 @@ Data2D<int> GridLocalPolynomial::buildUpdateMap(double tolerance, TypeRefinement
int max_1D_parents = rule->getMaxNumParents();

HierarchyManipulations::SplitDirections split(points);
std::vector<int> global_to_pnts(num_points);
#pragma omp parallel for firstprivate(global_to_pnts)
for(int j=0; j<split.getNumJobs(); j++){ // split.getNumJobs() gives the number of 1D interpolants to construct
int d = split.getJobDirection(j);
int nump = split.getJobNumPoints(j);
const int *pnts = split.getJobPoints(j);

std::vector<int> levels(nump);
#pragma omp parallel
{
int max_nump = split.getMaxNumPoints();
std::vector<int> global_to_pnts(num_points);
std::vector<int> levels(max_nump);

std::vector<int> monkey_count;
std::vector<int> monkey_tail;
std::vector<bool> used;
if (max_1D_parents > 1) {
monkey_count.resize(top_level + 1);
monkey_tail.resize(top_level + 1);
used.resize(max_nump, false);
}

int max_level = 0;
#pragma omp for
for(int j=0; j<split.getNumJobs(); j++){ // split.getNumJobs() gives the number of 1D interpolants to construct
int d = split.getJobDirection(j);
int nump = split.getJobNumPoints(j);
const int *pnts = split.getJobPoints(j);

Data2D<double> vals(active_outputs, nump);
int max_level = 0;

for(int i=0; i<nump; i++){
const double* v = values.getValues(pnts[i]);
const int *p = points.getIndex(pnts[i]);
if (output == -1){
std::copy(v, v + num_outputs, vals.getStrip(i));
}else{
vals.getStrip(i)[0] = v[output];
Data2D<double> vals(active_outputs, nump);

if (output == -1) {
for(int i=0; i<nump; i++){
std::copy_n(values.getValues(pnts[i]), num_outputs, vals.getStrip(i));
global_to_pnts[pnts[i]] = i;
levels[i] = rule->getLevel(points.getIndex(pnts[i])[d]);
if (max_level < levels[i]) max_level = levels[i];
}
} else {
for(int i=0; i<nump; i++){
vals.getStrip(i)[0] = values.getValues(pnts[i])[output];
global_to_pnts[pnts[i]] = i;
levels[i] = rule->getLevel(points.getIndex(pnts[i])[d]);
if (max_level < levels[i]) max_level = levels[i];
}
}
global_to_pnts[pnts[i]] = i;
levels[i] = rule->getLevel(p[d]);
if (max_level < levels[i]) max_level = levels[i];
}

std::vector<int> monkey_count(max_level + 1);
std::vector<int> monkey_tail(max_level + 1);
if (max_1D_parents == 1) {
for(int l=1; l<=max_level; l++){
for(int i=0; i<nump; i++){
if (levels[i] == l){
double x = rule->getNode(points.getIndex(pnts[i])[d]);
double *valsi = vals.getStrip(i);

for(int l=1; l<=max_level; l++){
for(int i=0; i<nump; i++){
if (levels[i] == l){
double x = rule->getNode(points.getIndex(pnts[i])[d]);
double *valsi = vals.getStrip(i);

int current = 0;
monkey_count[0] = d * max_1D_parents;
monkey_tail[0] = pnts[i]; // uses the global indexes
std::vector<bool> used(nump, false);

while(monkey_count[0] < (d+1) * max_1D_parents){
if (monkey_count[current] < (d+1) * max_1D_parents){
int branch = dagUp.getStrip(monkey_tail[current])[monkey_count[current]];
if ((branch == -1) || (used[global_to_pnts[branch]])){
monkey_count[current]++;
}else{
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);
const double *branch_vals = vals.getStrip(global_to_pnts[branch]);
for(int k=0; k<active_outputs; k++) valsi[k] -= basis_value * branch_vals[k];

used[global_to_pnts[branch]] = true;
monkey_count[++current] = d * max_1D_parents;
monkey_tail[current] = branch;
branch = dagUp.getStrip(branch)[d];
}
}
}
}
} else {
for(int l=1; l<=max_level; l++){
for(int i=0; i<nump; i++){
if (levels[i] == l){
double x = rule->getNode(points.getIndex(pnts[i])[d]);
double *valsi = vals.getStrip(i);

int current = 0;
monkey_count[0] = d * max_1D_parents;
monkey_tail[0] = pnts[i]; // uses the global indexes
std::fill_n(used.begin(), nump, false);

while(monkey_count[0] < (d+1) * max_1D_parents){
if (monkey_count[current] < (d+1) * max_1D_parents){
int branch = dagUp.getStrip(monkey_tail[current])[monkey_count[current]];
if ((branch == -1) || (used[global_to_pnts[branch]])){
monkey_count[current]++;
}else{
const int *branch_point = points.getIndex(branch);
double basis_value = rule->evalRaw(branch_point[d], x);
const double *branch_vals = vals.getStrip(global_to_pnts[branch]);
for(int k=0; k<active_outputs; k++) valsi[k] -= basis_value * branch_vals[k];

used[global_to_pnts[branch]] = true;
monkey_count[++current] = d * max_1D_parents;
monkey_tail[current] = branch;
}
}else{
monkey_count[--current]++;
}
}
}else{
monkey_count[--current]++;
}
}
}
}
}

// at this point, vals contains the one directional surpluses
for(int i=0; i<nump; i++){
const double *s = surpluses.getStrip(pnts[i]);
const double *c = scale.getStrip(pnts[i]);
const double *v = vals.getStrip(i);
bool small = true;
if (output == -1){
for(int k=0; k<num_outputs; k++){
small = small && (((c[k] * std::abs(s[k]) / norm[k]) <= tolerance) || ((c[k] * std::abs(v[k]) / norm[k]) <= tolerance));
// at this point, vals contains the one directional surpluses
for(int i=0; i<nump; i++){
const double *s = surpluses.getStrip(pnts[i]);
const double *c = scale.getStrip(pnts[i]);
const double *v = vals.getStrip(i);
bool small = true;
if (output == -1){
for(int k=0; k<num_outputs; k++){
small = small && (((c[k] * std::abs(s[k]) / norm[k]) <= tolerance) || ((c[k] * std::abs(v[k]) / norm[k]) <= tolerance));
}
}else{
small = ((c[0] * std::abs(s[output]) / norm[output]) <= tolerance) || ((c[0] * std::abs(v[0]) / norm[output]) <= tolerance);
}
}else{
small = ((c[0] * std::abs(s[output]) / norm[output]) <= tolerance) || ((c[0] * std::abs(v[0]) / norm[output]) <= tolerance);
pmap.getStrip(pnts[i])[d] = (small) ? 0 : 1;;
}
pmap.getStrip(pnts[i])[d] = (small) ? 0 : 1;;
}
}
}
Expand All @@ -1357,28 +1391,40 @@ MultiIndexSet GridLocalPolynomial::getRefinementCanidates(double tolerance, Type

int num_points = points.getNumIndexes();

if (level_limits.empty()){
for(int i=0; i<num_points; i++){
const int *map = pmap.getStrip(i);
for(int j=0; j<num_dimensions; j++){
if (map[j] == 1){ // if this dimension needs to be refined
if (!(useParents && addParent(points.getIndex(i), j, points, refined))){
addChild(points.getIndex(i), j, points, refined);
#pragma omp parallel
{
Data2D<int> lrefined(num_dimensions, 0);

if (level_limits.empty()){
#pragma omp for
for(int i=0; i<num_points; i++){
const int *map = pmap.getStrip(i);
for(int j=0; j<num_dimensions; j++){
if (map[j] == 1){ // if this dimension needs to be refined
if (!(useParents && addParent(points.getIndex(i), j, points, lrefined))){
addChild(points.getIndex(i), j, points, lrefined);
}
}
}
}
}
}else{
for(int i=0; i<num_points; i++){
const int *map = pmap.getStrip(i);
for(int j=0; j<num_dimensions; j++){
if (map[j] == 1){ // if this dimension needs to be refined
if (!(useParents && addParent(points.getIndex(i), j, points, refined))){
addChildLimited(points.getIndex(i), j, points, level_limits, refined);
}else{
#pragma omp for
for(int i=0; i<num_points; i++){
const int *map = pmap.getStrip(i);
for(int j=0; j<num_dimensions; j++){
if (map[j] == 1){ // if this dimension needs to be refined
if (!(useParents && addParent(points.getIndex(i), j, points, lrefined))){
addChildLimited(points.getIndex(i), j, points, level_limits, lrefined);
}
}
}
}
}

#pragma omp critical
{
refined.append(lrefined);
}
}

MultiIndexSet result(refined);
Expand Down
3 changes: 3 additions & 0 deletions SparseGrids/tsgHierarchyManipulator.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -309,6 +309,9 @@ class SplitDirections{
int getJobNumPoints(int job) const{ return (int) job_pnts[job].size(); }
//! \brief Return the indexes of the points associated with the job.
const int* getJobPoints(int job) const{ return job_pnts[job].data(); }
//! \brief Return the max number of points for any job.
int getMaxNumPoints() const { return (job_pnts.size() > 0) ? (int) std::max_element(job_pnts.begin(), job_pnts.end(),
[&](std::vector<int> const &a, std::vector<int> const& b)->bool{ return (a.size() < b.size()); })->size() : 0; }

private:
std::vector<int> job_directions;
Expand Down
6 changes: 6 additions & 0 deletions SparseGrids/tsgIndexSets.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -192,6 +192,12 @@ class Data2D{
num_strips++;
}

//! \brief Uses std::vector::insert to append all the data from the other to this
void append(Data2D<T> const &other) {
vec.insert(vec.end(), other.vec.begin(), other.vec.end());
num_strips += other.num_strips;
}

private:
size_t stride, num_strips;
std::vector<T> vec;
Expand Down
7 changes: 5 additions & 2 deletions SparseGrids/tsgSequenceOptimizer.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -346,8 +346,11 @@ template<> double getValue<rule_mindeltaodd>(CurrentNodes<rule_mindeltaodd> cons
auto lag = evalLagrange(current.nodes, current.coeff, x);
auto lag_less1 = evalLagrange(current.nodes_less1, current.coeff_less1, x);

return std::abs(lag.back()) + std::inner_product(lag_less1.begin(), lag_less1.end(), lag.begin(), 0.0,
std::plus<double>(), [](double a, double b)->double{ return std::abs(a - b); });
double result = std::abs(lag.back());
for(size_t i=0; i<lag_less1.size(); i++) {
result += std::abs(lag_less1[i] - lag[i]);
}
return result;
}

/*!
Expand Down

0 comments on commit 1de7210

Please sign in to comment.