Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Improve updatemap #760

Merged
merged 5 commits into from
Mar 21, 2024
Merged

Improve updatemap #760

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
8a1eb19
remove all allocation inside the for-loop
mkstoyanov Mar 20, 2024
9697489
remove if from the loop
mkstoyanov Mar 20, 2024
35eeafc
faster case when having single parent
mkstoyanov Mar 20, 2024
9b7f6da
parallize processing of the update data
mkstoyanov Mar 20, 2024
2206070
remove the use of inner-product with abs
mkstoyanov Mar 20, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
21 changes: 12 additions & 9 deletions Addons/testLoadUnstructured.hpp
View file
Open in desktop
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
Open in desktop
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
Open in desktop
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
Open in desktop
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
Open in desktop
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
Open in desktop
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
Loading