Solved stupid bug that happened when the pdf integral was zero.

adolfomunoz · May 7, 2024 · 5c12064 · 5c12064
1 parent d27017d
commit 5c12064
Show file tree

Hide file tree

Showing 3 changed files with 27 additions and 18 deletions.
diff --git a/src/control-variates/regions-integrator-variance-reduction.h b/src/control-variates/regions-integrator-variance-reduction.h
@@ -76,8 +76,10 @@ class RegionsIntegratorVarianceReduction {
  f_regdim(sample)*double(factor)*rrfactor*sfactor,
  r->approximation_at(sample)*double(factor)*rrfactor*sfactor
  );
+// std::cerr<<"Usage : "<<f_regdim(sample)<<" "<<r->approximation_at(sample)<<std::endl;
  }
- bins(pos) = accumulator.integral(approximation); 
+ bins(pos) = accumulator.integral(approximation);
+// std::cerr<<"Integral : "<<accumulator.integral(approximation); 
  } ,logger_control_variates);
  }
 };

diff --git a/src/newton-cotes/region.h b/src/newton-cotes/region.h
@@ -227,37 +227,44 @@ class Region {
  } 
 
 private:
- template<typename MA, std::size_t DIMSUB>
- value_type pdf_sub_last(const MA& ma, 
- const std::array<Float,DIMSUB>& a, const std::array<Float,DIMSUB>& b) const {
+ template<typename MA, std::size_t DIMSUB, typename Norm = NormDefault>
+ Float pdf_sub_last(const MA& ma, 
+ const std::array<Float,DIMSUB>& a, const std::array<Float,DIMSUB>& b,
+ const Norm& norm = Norm()) const {
  if constexpr (MA::dimensions > DIMSUB)
  return pdf_sub_last(ma.fold(quadrature,DIMSUB),a,b);
  else if constexpr (MA::dimensions > 1)
- return pdf_sub_last(ma.fold([&] (const auto& v) { return quadrature.pdf_integral_subrange(a[MA::dimensions-1],b[MA::dimensions-1],v); },MA::dimensions-1),a,b);
+ return pdf_sub_last(ma.fold([&] (const auto& v) { 
+ return quadrature.pdf_integral_subrange(a[MA::dimensions-1],b[MA::dimensions-1],v,norm); 
+ },MA::dimensions-1),a,b);
  else
- return ma.fold([&] (const auto& v) { return quadrature.pdf_integral_subrange(a[0],b[0],v); }).value();
+ return ma.fold([&] (const auto& v) { 
+ return quadrature.pdf_integral_subrange(a[0],b[0],v,norm); 
+ }).value();
  }
 
- template<std::size_t DIMSUB>
- value_type pdf_integral_subrange_last(const std::array<Float,DIMSUB>& a, 
- const std::array<Float,DIMSUB>& b) const {
+ template<std::size_t DIMSUB,typename Norm = NormDefault>
+ Float pdf_integral_subrange_last(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 integral for that many dimensions, as the region has less dimensions");
- return range().volume()*pdf_sub_last(data,range().pos_in_range(a),range().pos_in_range(b));
+ return range().volume()*pdf_sub_last(data,range().pos_in_range(a),range().pos_in_range(b),norm);
  }
 
 
- template<std::size_t DIMSUB>
+ template<std::size_t DIMSUB, typename Norm = NormDefault>
  value_type pdf_integral_subrange(const std::array<Float,DIMSUB>& a, 
- const std::array<Float,DIMSUB>& b) const {
- return pdf_integral_subrange_last(a,b);
+ const std::array<Float,DIMSUB>& b, const Norm& norm = Norm()) const {
+ return pdf_integral_subrange_last(a,b,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");
- return Range(a,b).volume()*norm(this->approximation_at(pos))/(norm(this->pdf_integral_subrange(a,b)));
+ 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;
  }
 
  template<std::size_t DIMSUB,typename Norm = NormDefault>

diff --git a/src/newton-cotes/rules.h b/src/newton-cotes/rules.h
@@ -189,7 +189,7 @@ struct Simpson {
  //In case you are curious, this is Cardano's method for one solution:
  auto p = c/a - pow(b,2.)/(3.*pow(a,2.));
  auto q = 2*pow(b,3.)/(27.*pow(a,3.)) - b*c/(3.*pow(a,2.)) + d/a;
- std::cerr<<a<<" "<<b<<" "<<c<<" "<<d<<" "<<p<<" "<<q<<std::endl; 
+// std::cerr<<a<<" "<<b<<" "<<c<<" "<<d<<" "<<p<<" "<<q<<std::endl; 
  auto sqr = pow(q/2.,2.) + pow(p/3.,3.);
  if(sqr >= 0.){ 
  //One solution
@@ -214,9 +214,9 @@ struct Simpson {
  if ((r>=t0) && (r<=t1) && (!std::isnan(r))) return r; 
  } 
  //Uniform sampling if not found a solution before
- std::cout<<"Warning, no solutions for range "<<t0<<" - "<<t1;
- for (Float r : res) std::cout<<" "<<r;
- std::cout<<std::endl; 
+// std::cout<<"Warning, no solutions for range "<<t0<<" - "<<t1;
+// for (Float r : res) std::cout<<" "<<r;
+// std::cout<<std::endl; 
  return s*(t1-t0) + t0;
  }