Found a couple of corner cases with issues with alpha optimization an…

…d russian roulette.

main/compilation-tests/multiple-parameters-sequence.cc

@@ -46,7 +46,7 @@ int main() {
  {
  LoggerProgress logger("Parallel"); //I'm suprised this works with a RNG shared among threads
  std::vector<float> sol(bins,0.0f); 
- integrate(integrator_per_bin_parallel(monte_carlo(samples)),sol,f,range_infinite(-1.0,-1.0,1.0,1.0),logger);
+ integrate(integrator_per_bin_parallel(monte_carlo(100*samples)),sol,f,range_infinite(-1.0,-1.0,1.0,1.0),logger);
  for (float v : sol) std::cout<<std::fixed<<std::setprecision(2)<<std::setw(4)<<v<<" ";
  std::cout<<std::endl;
  }
@@ -59,19 +59,28 @@ int main() {
  for (float v : sol) std::cout<<std::fixed<<std::setprecision(2)<<std::setw(4)<<v<<" ";
  std::cout<<std::endl;
  }
-*/
+
  {
  LoggerProgress logger("Fubini adaptive"); 
  std::vector<float> sol(bins,0.0f); 
  integrate(integrator_fubini<4>(integrator_adaptive_tolerance(nested(simpson,trapezoidal),1.e-4f),monte_carlo(16)),sol,f,range_infinite(-1.0,-1.0,1.0,1.0),logger);
  for (float v : sol) std::cout<<std::fixed<<std::setprecision(2)<<std::setw(4)<<v<<" ";
  std::cout<<std::endl;
  }
+*/
+
+ {
+ LoggerProgress logger("PPA IS");
+ std::vector<float> sol(bins,0.0f); 
+ integrate(integrator_fubini({0,1,4,5},integrator_adaptive_variance_reduction_parallel(nested(simpson,trapezoidal),128,rr_integral_region(),cv_fixed_weight(0.0),samples),monte_carlo(10)),sol,f,range_infinite(-1.0,-1.0,1.0,1.0),logger);
+ for (float v : sol) std::cout<<std::fixed<<std::setprecision(2)<<std::setw(4)<<v<<" ";
+ std::cout<<std::endl;
+ }
 
  {
- LoggerProgress logger("New Fubini");
+ LoggerProgress logger("PPA CV");
  std::vector<float> sol(bins,0.0f); 
- integrate(integrator_fubini({0,1,4,5},integrator_adaptive_variance_reduction_parallel(nested(simpson,trapezoidal),128,rr_integral_region(),cv_optimize_weight(),samples),monte_carlo(10)),sol,f,range_infinite(-1.0,-1.0,1.0,1.0),logger);
+ integrate(integrator_fubini({0,1,4,5},integrator_adaptive_variance_reduction_parallel(nested(simpson,trapezoidal),128,rr_uniform_region(),cv_fixed_weight(1.0),samples),monte_carlo(10)),sol,f,range_infinite(-1.0,-1.0,1.0,1.0),logger);
  for (float v : sol) std::cout<<std::fixed<<std::setprecision(2)<<std::setw(4)<<v<<" ";
  std::cout<<std::endl;
  }

src/control-variates/regions-integrator-parallel-variance-reduction.h

@@ -58,10 +58,15 @@ class rr_integral_region {
 
  template<typename Rs>
  RR(const Rs& regions, const Norm& n) : weights(regions.size()), norm(n) {
- std::size_t i = 0;
+ std::size_t i = 0; 
  for (const auto& [r,region_bin_range] : regions) {
  weights[i++] = norm(r->integral_subrange(region_bin_range));
  }
+ //This is to avoid zeroes and negative numbers: put a minimum relative weight
+ double sum = 0.0;
+ for (double w : weights) sum += w;
+ if (sum<=0.0) for(double& w : weights) w = 1.0;
+ else for (double& w : weights) w = std::max(w,0.01*sum/double(weights.size()));
  rr = std::discrete_distribution<std::size_t>(weights.begin(),weights.end());
  }
 
@@ -148,13 +153,16 @@ class cv_optimize_weight {
  } 
 
  double alpha() const {
+ if (size < 2) return 1;
  auto c = std::max(0.0,covariance());
+ if (c<=0.0) return 0.0;
  auto v = std::max(c,variance());
  return c/v;
  } 
 
  Sample integral(const Sample& approximation) const {
  auto a = alpha();
+// std::cerr<<"Alpha = "<<a<<std::endl;
  return (sum_f - a*sum_app)/double(size) + a*approximation;
  }
  friend class cv_optimize_weight;
@@ -200,7 +208,7 @@ class RegionsIntegratorParallelVarianceReduction {
  }); 
  logger_bins.log_progress(final_progress,final_progress);
  auto logger_control_variates = logger_step(logger,"residual and variance reduction");
- for_each(parallel,multidimensional_range(bin_resolution),
+ for_each(sequential,multidimensional_range(bin_resolution),
  [&] (const std::array<std::size_t, DIMBINS>& pos) {
  using Sample = decltype(f(std::declval<std::array<Float,DIM>>()));
  Range<Float,DIM> binrange = range;