Solved bug in pdf estimation. Added pdf russian roulette.

main/compilation-tests/multiple-parameters-2d.cc

@@ -100,7 +100,7 @@ int main() {
  std::vector<std::vector<float>> sol(bins,std::vector<float>(bins,0.0f));
  integrate(integrator_adaptive_fubini_variance_reduction<2>(
  nested(simpson,trapezoidal),error_heuristic_default(error_metric_absolute()),128,10,
- rr_integral_region(),cv_fixed_weight(0),region_sampling_importance(),samples),
+ rr_pdf_region(),cv_fixed_weight(0),region_sampling_importance(),samples),
  sol,f,range_primary<3>(),logger);
  for (const auto& vv : sol) {
  for (float v : vv)

src/control-variates/region-russian-roulette.h

@@ -107,6 +107,47 @@ class rr_error_region {
  }
 };
 
+template<typename Norm = NormDefault>
+class rr_pdf_region {
+ Norm norm;
+public:
+ rr_pdf_region(const Norm& n = Norm()) : norm(n) {}
+ class RR {
+ private:
+ //This is the MonteCarlo residual, RR among regions
+ std::discrete_distribution<std::size_t> rr;
+ std::vector<double> weights;
+ Norm norm;
+
+ template<typename Rs>
+ RR(const Rs& regions, const Norm& n) : weights(regions.size()), norm(n) {
+ std::size_t i = 0; 
+ for (const auto& [r,region_bin_range] : regions) {
+ weights[i++] = r->pdf_integral_subrange(region_bin_range,norm);
+ }
+ //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());
+ }
+
+ public:
+ template<typename RNG>
+ std::tuple<std::size_t,double> choose(RNG& rng) { 
+ std::size_t choice = rr(rng);
+ return std::tuple(choice,1.0/rr.probabilities()[choice]); 
+ }
+ friend class rr_pdf_region;
+ };
+ template<typename Regions>
+ RR russian_roulette(const Regions& regions) const {
+ return RR(regions,norm);
+ }
+};
+
+
 
 
 } 

src/newton-cotes/region.h

@@ -250,21 +250,25 @@ class Region {
  return range().volume()*pdf_sub_last(data,range().pos_in_range(a),range().pos_in_range(b),norm);
  }
 
- 
+public:
  template<std::size_t DIMSUB, typename Norm = NormDefault>
  Float pdf_integral_subrange(const std::array<Float,DIMSUB>& a, 
  const std::array<Float,DIMSUB>& b, const Norm& norm = Norm()) const {
  return pdf_integral_subrange_last(a,b,norm);
  }
 
+ template<std::size_t DIMSUB,typename Norm = NormDefault>
+ Float pdf_integral_subrange(const Range<Float,DIMSUB>& range, const Norm& norm = Norm()) const {
+ return this->pdf_integral_subrange(range.min(),range.max(),norm);
+ } 
 public:
  template<std::size_t DIMSUB,typename Norm = NormDefault>
  Float pdf_subrange(const std::array<Float,DIMSUB>& pos, const std::array<Float,DIMSUB>& a, 
  const std::array<Float,DIMSUB>& b,const Norm& norm = Norm()) const {
  static_assert(DIM>=DIMSUB,"Cannot calculate the subrange pdf for that many dimensions, as the region has less dimensions");
  Float den = this->pdf_integral_subrange(a,b,norm);
- if (den < 1.e-10) return Float(0);
- else return Range(a,b).volume()*norm(this->approximation_at(pos))/den;
+ if (den < 1.e-10) den=1;
+ return Range(a,b).volume()*norm(this->approximation_at(pos))/den;
  }
 
  template<std::size_t DIMSUB,typename Norm = NormDefault>