Added several options to define integrands: either an array parameter…

… or a multiple parameter function.

CMakeLists.txt

@@ -13,8 +13,11 @@ include(Compiler)
 ###########################################################################################
 find_package(TBB REQUIRED)
 
-add_executable(test main/test.cc)
-target_link_libraries(test PRIVATE TBB::tbb)
+add_executable(test-array-parameter main/compilation-tests/array-parameter.cc)
+target_link_libraries(test-array-parameter PRIVATE TBB::tbb)
+add_executable(test-multiple-parameters main/compilation-tests/multiple-parameters.cc)
+target_link_libraries(test-multiple-parameters PRIVATE TBB::tbb)
+add_executable(test-random-sequence main/compilation-tests/random-sequence.cc)
 
 ##########
 # FOR DOCUMENTATION

main/test.cc → main/compilation-tests/array-parameter.cc

@@ -1,4 +1,4 @@
-#include "../viltrum.h"
+#include "../../viltrum.h"
 #include <iostream>
 #include <iomanip>
 
@@ -7,10 +7,12 @@ using namespace viltrum;
 int main() {
  const std::size_t bins = 24;
  const unsigned long samples = 1000000;
+
  auto f =[] (const std::array<float,2>& x) -> double {
-  if ((x[0]+x[1])<1) return 1.0;
-  else return 0.0;
+ if ((x[0]+x[1])<1) return 1.0; 
+ else return 0.0;
  };
+
  //With floats, due to numerica stability, for a very large number of samples, monte_carlo fails to converge?
  {
  LoggerProgress logger("Simple");

main/compilation-tests/multiple-parameters.cc

@@ -0,0 +1,44 @@
+#include "../../viltrum.h"
+#include <iostream>
+#include <iomanip>
+
+using namespace viltrum;
+
+int main() {
+ const std::size_t bins = 24;
+ const unsigned long samples = 1000000;
+
+ auto f =[] (float x, float y) -> double {
+ if ((x+y)<1) return 1.0;
+ else return 0.0;
+ };
+
+ //With floats, due to numerica stability, for a very large number of samples, monte_carlo fails to converge?
+ {
+ LoggerProgress logger("Simple");
+ std::cout<<integrate(monte_carlo(samples*bins),f,range_primary<2>(),logger)<<std::endl; 
+ }
+
+ {
+ LoggerProgress logger("Monte-Carlo");
+ std::vector<float> sol(bins,0.0f);
+ integrate(monte_carlo(bins*samples),sol,f,range_primary<2>(),logger);
+ for (float v : sol) std::cout<<std::fixed<<std::setprecision(2)<<std::setw(4)<<v<<" ";
+ std::cout<<std::endl;
+ }
+ {
+ LoggerProgress logger("Per bin");
+ std::vector<float> sol(bins,0.0f); 
+ integrate(integrator_per_bin(monte_carlo(samples)),sol,f,range_primary<2>(),logger);
+ for (float v : sol) std::cout<<std::fixed<<std::setprecision(2)<<std::setw(4)<<v<<" ";
+ std::cout<<std::endl;
+ }
+ {
+ 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_primary<2>(),logger);
+ for (float v : sol) std::cout<<std::fixed<<std::setprecision(2)<<std::setw(4)<<v<<" ";
+ std::cout<<std::endl;
+ }
+
+} 

main/compilation-tests/random-sequence.cc

@@ -0,0 +1,16 @@
+#include "../../src/monte-carlo/random-sequence.h"
+#include <iostream>
+
+
+
+using namespace viltrum;
+
+int main() {
+ auto l = random_sequence(0.0,1.0);
+ int i = 0;
+ for (double d : l) {
+ std::cout<<d<<" ";
+ if (++i>=10) break;
+ } 
+ std::cout<<std::endl;
+} 

src/integrate.h

@@ -7,11 +7,57 @@
 
 namespace viltrum {
 
+namespace detail {
+ template<typename P, typename F, typename Enable = void>
+ struct Integrand {
+ static const F& adapt(const F& f) { return f; } 
+ };
+
+ template<typename P, typename F>
+ struct Integrand<P,F,typename std::enable_if_t<std::is_invocable_v<F,P>> > {
+ static auto adapt(const F& f) { 
+ return [&f] (const std::array<P,1>& x) { return f(x[0]);}; 
+ } 
+ };
+
+ template<typename P, typename F>
+ struct Integrand<P,F,typename std::enable_if_t<std::is_invocable_v<F,P,P>> > {
+ static auto adapt(const F& f) { 
+ return [&f] (const std::array<P,2>& x) { return f(x[0],x[1]);}; 
+ } 
+ };
+
+ template<typename P, typename F>
+ struct Integrand<P,F,typename std::enable_if_t<std::is_invocable_v<F,P,P,P>> > {
+ static auto adapt(const F& f) { 
+ return [&f] (const std::array<P,3>& x) { return f(x[0],x[1],x[2]);}; 
+ } 
+ };
+
+ template<typename P, typename F>
+ struct Integrand<P,F,typename std::enable_if_t<std::is_invocable_v<F,P,P,P,P>> > {
+ static auto adapt(const F& f) { 
+ return [&f] (const std::array<P,4>& x) { return f(x[0],x[1],x[2],x[3]);}; 
+ } 
+ };
+
+ template<typename P, typename F>
+ auto adapt(const F& f) {
+ return Integrand<P,std::decay_t<F>>::adapt(f);
+ } 
+} 
+
 template<typename Integrator, typename Bins, std::size_t DIMBINS, typename F, typename Float, std::size_t DIM, typename Logger>
-void integrate(const Integrator& integrator, Bins& bins, const std::array<std::size_t,DIMBINS>& resolution, const F& function, const Range<Float,DIM>& range, Logger& logger) {
- integrator.integrate(bins,resolution, function, range, logger);
+void integrate(const Integrator& integrator, Bins& bins, const std::array<std::size_t,DIMBINS>& resolution, 
+ const F& function, const Range<Float,DIM>& range, Logger& logger,
+ std::enable_if_t<std::is_convertible_v<
+ std::decay_t<decltype(detail::adapt<Float>(function)(std::array<Float,DIM>()))>,
+ std::decay_t<decltype(bins(std::array<std::size_t,DIMBINS>()))>>,int> dummy=0) {
+ integrator.integrate(bins,resolution, detail::adapt<Float>(function), range, logger);
 }
 
+
+
 template<typename Integrator, typename Bins, std::size_t DIMBINS, typename F, typename Float, std::size_t DIM>
 void integrate(const Integrator& integrator, Bins& bins, const std::array<std::size_t,DIMBINS>& resolution, const F& function, const Range<Float,DIM>& range) {
  LoggerNull log;
@@ -20,8 +66,8 @@ void integrate(const Integrator& integrator, Bins& bins, const std::array<std::s
 
 template<typename Integrator, typename F, typename R, typename Logger>
 auto integrate(const Integrator& integrator, const F& function, const R& range, Logger& logger) {
- using T = decltype(function(range.min()));
- T sol(0);
+ using T = decltype(detail::adapt<typename R::value_type>(function)(range.min()));
+ T sol(0.0);
  auto bins = [&sol] (const std::array<std::size_t,1>& i) -> T& { return sol; };
  std::array<std::size_t,1> res{1};
  integrate(integrator,bins,res,function,range,logger);

src/monte-carlo/random-sequence.h

@@ -0,0 +1,39 @@
+#pragma once
+#include <random>
+
+namespace viltrum {
+
+//I am implementing it so it is predictable once we initialize the RNG.
+template<typename RNG, typename Number>
+class RandomSequence {
+ RNG rng;
+ std::uniform_real_distribution<Number> dis;
+public:
+ RandomSequence(RNG&& r, const Number& lower = Number(0), const Number& upper = Number(1)) : 
+ rng(std::forward<RNG>(r)), dis(lower,upper) {}
+
+ class const_iterator {
+ private:
+ RNG rng;
+ std::uniform_real_distribution<Number> dis;
+ Number n;
+ const_iterator(const RNG& r, const std::uniform_real_distribution<Number>& d) :
+ rng(r), dis(d), n(dis(rng)) {} 
+ friend class RandomSequence<RNG,Number>;
+ public:
+ const Number& operator*() const { return n; }
+ const_iterator& operator++() { n = dis(rng); return (*this); }
+ bool operator==(const const_iterator& that) const { return false; } //Infinite list
+ bool operator!=(const const_iterator& that) const { return true; } //Infinite list
+ }; 
+
+ const_iterator begin() const { return const_iterator(rng,dis); }
+ const_iterator end() const { return const_iterator(rng,dis); }
+};
+
+template<typename Number,typename RNG = std::mt19937>
+auto random_sequence(const Number& lower, const Number& upper, std::size_t seed = std::random_device()()) {
+ return RandomSequence<RNG,Number>(RNG(seed),lower,upper);
+} 
+
+}

src/range.h

@@ -22,6 +22,9 @@ 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]);
  }
+
+ using value_type = T;
+ static constexpr std::size_t size = DIM;
 
  const std::array<T,DIM>& min() const { return (*this)[0]; }
  T min(std::size_t i) const { return min()[i]; }