Added new example that explores different ways of defining integratio…

…n ranges.
adolfomunoz · Jun 15, 2021 · 711a0c9 · 711a0c9
1 parent 612d343
commit 711a0c9
Show file tree

Hide file tree

Showing 4 changed files with 65 additions and 11 deletions.
diff --git a/CMakeLists.txt b/CMakeLists.txt
@@ -25,7 +25,11 @@ add_dependencies(test-quadrature-bins ${function_1d_deps})
 #add_dependencies(test-quadrature-plot svg-cpp-plot ${function_1d_deps})
 add_executable(test-multidimensional-range main/test-multidimensional-range.cc)
 
+##########
+# FOR DOCUMENTATION
+##########
 add_executable(integrands main/doc/integrands.cc)
+add_executable(ranges main/doc/ranges.cc)
 
 #add_executable(convergence-1d main/convergence-1d.cc)
 #add_dependencies(convergence-1d svg-cpp-plot ${function_1d_deps})

diff --git a/main/doc/ranges.cc b/main/doc/ranges.cc
@@ -0,0 +1,46 @@
+#include "../../viltrum.h"
+
+
+class Gaussian {
+public:
+ template<typename real, std::size_t NDIM>
+ real operator()(const std::array<real,NDIM>& x) const {
+ real sum(0);
+ for (real xi : x) sum += xi*xi;
+ return std::exp(-0.5*sum)/std::sqrt(2.0*M_PI);
+ }
+};
+
+int main(int argc, char** argv) {
+ unsigned long samples = 16;
+ for (int i = 0; i<argc-1; ++i) {
+ if (std::string(argv[i])=="-samples") samples = atol(argv[++i]);
+ }
+
+ auto method = viltrum::integrator_monte_carlo_uniform(samples);
+
+ std::array<double,1> min1d{0.0}; std::array<double,1> max1d{1.0};
+ std::cout<<"1D = "<<method.integrate(Gaussian(),viltrum::range(min1d,max1d))<<std::endl;
+ std::array<double,2> min2d{0.0,0.0}; std::array<double,2> max2d{1.0,1.0};
+ std::cout<<"2D = "<<method.integrate(Gaussian(),viltrum::range(min2d,max2d))<<std::endl;
+ std::array<double,3> min3d{0.0,0.0,0.0}; std::array<double,3> max3d{1.0,1.0,1.0};
+ std::cout<<"3D = "<<method.integrate(Gaussian(),viltrum::range(min3d,max3d))<<std::endl;
+
+ std::cout<<"1D = "<<method.integrate(Gaussian(),viltrum::range(0.0,1.0))<<std::endl;
+ std::cout<<"2D = "<<method.integrate(Gaussian(),viltrum::range(0.0,0.0,1.0,1.0))<<std::endl;
+ std::cout<<"3D = "<<method.integrate(Gaussian(),viltrum::range(0.0,0.0,0.0,1.0,1.0,1.0))<<std::endl;
+
+ std::cout<<"1D = "<<method.integrate(Gaussian(),viltrum::range(0.0,1.0))<<std::endl;
+ std::cout<<"2D = "<<method.integrate(Gaussian(),
+ viltrum::range(0.0,1.0) | viltrum::range(0.0,1.0))<<std::endl;
+ std::cout<<"3D = "<<method.integrate(Gaussian(),
+ viltrum::range(0.0,1.0) | viltrum::range(0.0,1.0) | viltrum::range(0.0,1.0))<<std::endl;
+
+ std::cout<<"1D = "<<method.integrate(Gaussian(),viltrum::range_all<1>(0.0,1.0))<<std::endl;
+ std::cout<<"2D = "<<method.integrate(Gaussian(),viltrum::range_all<2>(0.0,1.0))<<std::endl;
+ std::cout<<"3D = "<<method.integrate(Gaussian(),viltrum::range_all<3>(0.0,1.0))<<std::endl;
+
+ std::cout<<"1D = "<<method.integrate(Gaussian(),viltrum::range_primary<1>())<<std::endl;
+ std::cout<<"2D = "<<method.integrate(Gaussian(),viltrum::range_primary<2>())<<std::endl;
+ std::cout<<"3D = "<<method.integrate(Gaussian(),viltrum::range_primary<3>())<<std::endl;
+}
diff --git a/multiarray/array.h b/multiarray/array.h
@@ -21,7 +21,7 @@ auto operator|(const T& t, const std::array<T,N>& a) noexcept -> std::array<T,N+
 }
 
 template<typename T, std::size_t N1, std::size_t N2>
-auto operator|(const std::array<T,N1>& a1,const std::array<T,N1>& a2) noexcept -> std::array<T,N1+N2> {
+auto operator|(const std::array<T,N1>& a1,const std::array<T,N2>& a2) noexcept -> std::array<T,N1+N2> {
  std::array<T,N1+N2> s;
  std::copy(a1.begin(),a1.end(),s.begin());
  std::copy(a2.begin(),a2.end(),s.begin()+N1);

diff --git a/quadrature/range.h b/quadrature/range.h
@@ -1,6 +1,8 @@
 #pragma once
 
-#include <array>
+#include "../multiarray/array.h"
+#include <tuple>
+#include <type_traits>
 
 namespace viltrum {
 
@@ -13,13 +15,13 @@ class Range : public std::array<std::array<T,DIM>,2> {
  _volume = T(1);
  for (std::size_t i = 0; i<DIM; ++i) _volume*=(b[i]-a[i]);
  }
- 
+ 
  const std::array<T,DIM>& min() const { return (*this)[0]; }
  T min(std::size_t i) const { return min()[i]; }
  const std::array<T,DIM>& max() const { return (*this)[1]; }
  T max(std::size_t i) const { return max()[i]; }
  T volume() const { return _volume; }
-
+ 
  bool is_inside(const std::array<T,DIM>& x) const {
  bool is = true;
  for (std::size_t i = 0; (i<DIM) && is; ++i)
@@ -106,13 +108,19 @@ class Range : public std::array<std::array<T,DIM>,2> {
  }
 };
 
+
 template<typename T, std::size_t DIM>
 Range<T,DIM> range(const std::array<T,DIM>& a, const std::array<T,DIM>& b) {
  return Range<T,DIM>(a,b);
 }
 
+template<typename T, std::size_t N1, std::size_t N2>
+auto operator|(const Range<T,N1>& r1,const Range<T,N2>& r2) noexcept -> Range<T,N1+N2> {
+ return range(r1.min() | r2.min(), r1.max() | r2.max());
+}
+
 template<typename T>
-Range<T,1> range(const T& a, const T& b) {
+Range<T,1> range(const T& a, const T& b, std::enable_if_t<std::is_floating_point_v<T>,int> dummy = 0) {
  return Range<T,1>(std::array<T,1>{a},std::array<T,1>{b});
 }
 
@@ -121,6 +129,7 @@ Range<T,2> range(const T& a0, const T& a1, const T& b0, const T& b1) {
  return Range<T,2>(std::array<T,2>{a0,a1},std::array<T,2>{b0,b1});
 }
 
+
 template<typename T>
 Range<T,3> range(const T& a0, const T& a1, const T& a2, const T& b0, const T& b1, const T& b2) {
  return Range<T,3>(std::array<T,3>{a0,a1,a2},std::array<T,3>{b0,b1,b2});
@@ -144,12 +153,7 @@ Range<T,N> range_all(const T& va, const T& vb) {
 
 template<std::size_t N, typename T = float>
 Range<T,N> range_primary() {
- std::array<T,N> a, b;
- for (std::size_t i = 0;i<N;++i) {
- a[i]=T(0);
- b[i]=T(1);
- }
- return range(a,b);
+ return range_all<N>(T(0),T(1));
 }
 
 }