Fix sum_hills overflows when doing projections

regtest/basic/rt-sum-hills-overflow/HILLS

@@ -0,0 +1,7 @@
+#! FIELDS time d1 vol sigma_d1 sigma_vol height biasf
+#! SET multivariate false
+#! SET kerneltype stretched-gaussian
+      0.050000      1.130546    127.932640      0.100000      0.200000      100.0     10.000000
+      0.100000      1.097928    127.932640      0.100000      0.200000      100.0     10.000000
+      0.150000      1.080244    127.932640      0.100000      0.200000      100.0     10.000000
+      0.200000      1.086855    127.932640      0.100000      0.200000      100.0     10.000000

regtest/basic/rt-sum-hills-overflow/Makefile

@@ -0,0 +1 @@
+include ../../scripts/test.make

regtest/basic/rt-sum-hills-overflow/config

@@ -0,0 +1,3 @@
+type=sum_hills
+# this is to test a different name
+arg=" --hills HILLS --fmt %12.7f --idw d1 --kt 0.1"

regtest/basic/rt-sum-hills-overflow/fes.dat.reference

@@ -0,0 +1,28 @@
+#! FIELDS d1 projection
+#! SET min_d1 0.726691
+#! SET max_d1 1.4841
+#! SET nbins_d1  23
+#! SET periodic_d1 false
+    0.7266910   -0.3044524
+    0.7611187   -0.9426032
+    0.7955464   -3.8206172
+    0.8299740  -11.1560780
+    0.8644017  -26.8774481
+    0.8988294  -56.3034760
+    0.9332571 -104.1750364
+    0.9676848 -171.1615119
+    1.0021125 -250.2626861
+    1.0365401 -325.9693625
+    1.0709678 -378.4433089
+    1.1053955 -391.7742547
+    1.1398232 -361.7430008
+    1.1742509 -297.9659459
+    1.2086785 -218.9493543
+    1.2431062 -143.4879780
+    1.2775339  -83.7886778
+    1.3119616  -43.4845150
+    1.3463893  -19.9096343
+    1.3808170   -7.8582589
+    1.4152446   -2.4672329
+    1.4496723   -0.5443036
+    1.4841000   -0.3044522

src/tools/Grid.cpp

@@ -936,13 +936,7 @@ Grid Grid::project(const std::vector<std::string> & proj, WeightBase *ptr2obj )
     // the vector vhigh now contains at the beginning the index of the low dimension and -1 for the values that need to be integrated
     double initval=0.;
     projectOnLowDimension(initval,vHigh, ptr2obj);
-    smallgrid.setValue(i,initval);
-  }
-  // reset to zero just for biasing (this option can be evtl enabled in a future...)
-  //double vmin;vmin=-smallgrid.getMinValue()+1;
-  for(unsigned i=0; i<smallgrid.getSize(); i++) {
-    double vv=smallgrid.getValue(i);
-    smallgrid.setValue(i,ptr2obj->projectOuterLoop(vv));
+    smallgrid.setValue(i,ptr2obj->projectOuterLoop(initval));
   }
 
   return smallgrid;

src/tools/Grid.h

@@ -44,9 +44,25 @@ class WeightBase {
 class BiasWeight:public WeightBase {
 public:
   double beta,invbeta;
+  double shift=0.0;
   explicit BiasWeight(double v) {beta=v; invbeta=1./beta;}
-  double projectInnerLoop(double &input, double &v) override {return  input+exp(beta*v);}
-  double projectOuterLoop(double &v) override {return -invbeta*std::log(v);}
+  //double projectInnerLoop(double &input, double &v) override {return  input+exp(beta*v);}
+  //double projectOuterLoop(double &v) override {return -invbeta*std::log(v);}
+  double projectInnerLoop(double &input, double &v) override {
+    auto betav=beta*v;
+    auto x=betav-shift;
+    if(x>0) {
+      shift=betav;
+      return input*std::exp(-x)+1;
+    } else {
+      return input+std::exp(x);
+    }
+  }
+  double projectOuterLoop(double &v) override {
+    auto res=-invbeta*(std::log(v)+shift);
+    shift=0.0;
+    return res;
+  }
 };
 
 class ProbWeight:public WeightBase {