Added new integrator for variance reduction that does not share MC sa…

…mples among control variate and residual.

main/compilation-tests/covariance.cc

@@ -1,4 +1,4 @@
-#include "../../src/control-variates/regions-integrator-parallel-variance-reduction.h"
+#include "../../src/control-variates/weight-strategy.h"
 #include <random>
 #include <list>
 #include <string>

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

@@ -81,7 +81,10 @@ int main() {
     {
         LoggerProgress logger("New Fubini");
         std::vector<std::vector<float>> sol(bins,std::vector<float>(bins,0.0f));
-        integrate(integrator_fubini({3,2},integrator_adaptive_variance_reduction_parallel(nested(simpson,trapezoidal),128,rr_integral_region(),cv_optimize_weight(),samples),monte_carlo(10)),sol,f,range_primary<3>(),logger);
+        integrate(integrator_adaptive_fubini_variance_reduction_parallel<2>(
+                    nested(simpson,trapezoidal),error_heuristic_default(error_metric_absolute()),128,10,
+                    rr_integral_region(),cv_optimize_weight(),region_sampling_uniform(),samples),
+                sol,f,range_primary<3>(),logger);
         for (const auto& vv : sol) {
             for (float v : vv)
                 std::cout<<std::fixed<<std::setprecision(2)<<std::setw(4)<<v<<" ";

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

@@ -84,5 +84,16 @@ int main() {
         for (float v : sol) std::cout<<std::fixed<<std::setprecision(2)<<std::setw(4)<<v<<" ";
         std::cout<<std::endl;
     }
+
+    {
+        LoggerProgress logger("New Fubini");
+        std::vector<float> sol(bins,0.0f); 
+        integrate(integrator_adaptive_fubini_variance_reduction_parallel<4>(
+                    nested(simpson,trapezoidal),error_heuristic_default(error_metric_absolute()),128,10,
+                    rr_integral_region(),cv_optimize_weight(),region_sampling_uniform(),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;
+    }
 
 } 

src/combination/fubini.h

@@ -3,13 +3,75 @@
 #include "concat.h"
 #include "reorder.h"
 #include <tuple>
+#include "../range.h"
+#include "../range-infinite.h"
+#include "../array.h"
+
+
 
 /* Fubini combines different integration techniques, one of them for the first N dimensions and the other
  * for the rest of the dimensions, even if infinite.
  */
 
 namespace viltrum {
 
+template<std::size_t N, typename Float, std::size_t DIM>
+std::tuple<Range<Float,N>, Range<Float,DIM-N>> range_split_at(const Range<Float,DIM>& range) {
+    std::array<Float,N> range_first_min, range_first_max;
+    for (std::size_t i = 0; i<N; ++i) {
+        range_first_min[i] = range.min(i);
+        range_first_max[i] = range.max(i);
+    }
+    std::array<Float,DIM-N> range_rest_min, range_rest_max;
+    for (std::size_t i = N; i<DIM; ++i) {
+        range_rest_min[i-N] = range.min(i);
+        range_rest_max[i-N] = range.max(i);
+    }      
+    return std::tuple(Range<Float,N>(range_first_min, range_first_max),
+                      Range<Float,DIM-N>(range_rest_min, range_rest_max));
+}
+
+template<std::size_t N, typename Float>
+std::tuple<Range<Float,N>, RangeInfinite<Float>> range_split_at(const RangeInfinite<Float>& range) {
+    std::array<Float,N> range_first_min, range_first_max;
+    for (std::size_t i = 0; i<N; ++i) {
+        range_first_min[i] = range.min(i);
+        range_first_max[i] = range.max(i);
+    }
+        std::vector<Float> range_rest_min = range.min();
+        std::vector<Float> range_rest_max = range.max();
+        for (std::size_t i = 0; i<N; ++i) {
+            if (!range_rest_min.empty()) range_rest_min.erase(range_rest_min.begin());
+            if (!range_rest_max.empty()) range_rest_max.erase(range_rest_max.begin());
+        }    
+        return std::tuple(Range<Float,N>(range_first_min, range_first_max),
+                     RangeInfinite<Float>(range_rest_min, range_rest_max));
+}
+
+template<std::size_t N, typename F, typename Integrator, typename Float>
+auto function_split_and_integrate_at(const F& f, const Integrator& integrator, const Range<Float,0>& range) {
+    return f; 
+} 
+
+template<std::size_t N, typename F, typename Integrator, typename Float, std::size_t DIM>
+auto function_split_and_integrate_at(const F& f, const Integrator& integrator, const Range<Float,DIM>& range) {
+    return [=] (const std::array<Float,N>& x) {
+        return viltrum::integrate(integrator,[&f,&x] (const std::array<float,DIM>& xr) {
+            return f( (x | xr) );
+        }, range);
+    };  
+}
+
+template<std::size_t N, typename F, typename Integrator, typename Float>
+auto function_split_and_integrate_at(const F& f, const Integrator& integrator, const RangeInfinite<Float>& range) {
+    return [=] (const std::array<Float,N>& x) {
+        return viltrum::integrate(integrator,[&f,&x] (const auto& xr) {
+            return f(concat(x,xr));
+        }, range);
+    }; 
+} 
+
+
 template<typename IntegratorFirst, typename IntegratorRest, std::size_t N>
 class IntegratorFubini {
     IntegratorFirst integrator_first;
@@ -18,61 +80,18 @@ class IntegratorFubini {
     IntegratorFubini(const IntegratorFirst& i_f, const IntegratorRest& i_r) :
         integrator_first(i_f), integrator_rest(i_r) {}
 
-	template<typename Bins, std::size_t DIMBINS, typename F, typename Float, std::size_t DIM, typename Logger>
+	template<typename Bins, std::size_t DIMBINS, typename F, typename Range, typename Logger>
 	void integrate(Bins& bins, const std::array<std::size_t,DIMBINS>& bin_resolution,
-		const F& f, const Range<Float,DIM>& range, Logger& logger) const {
+		const F& f, const Range& range, Logger& logger) const {
 
         static_assert(N>=DIMBINS,"Fubini does not work with that many dimensions on bin resolution");
 
-        std::array<Float,N> range_first_min, range_first_max;
-        for (std::size_t i = 0; i<N; ++i) {
-            range_first_min[i] = range.min(i);
-            range_first_max[i] = range.max(i);
-        }
-        Range<Float,N> range_first(range_first_min, range_first_max);
-
-
-        std::array<Float,DIM-N> range_rest_min, range_rest_max;
-        for (std::size_t i = N; i<DIM; ++i) {
-            range_rest_min[i-N] = range.min(i);
-            range_rest_max[i-N] = range.max(i);
-        }      
-        Range<Float,DIM-N> range_rest(range_rest_min, range_rest_max);
+        auto [range_first, range_rest] = range_split_at<N>(range);
 
         viltrum::integrate(integrator_first,bins,bin_resolution,
-            [&] (const std::array<float,N>& x) {
-                return viltrum::integrate(integrator_rest,[&f,&x] (const std::array<float,DIM-N>& xr) {
-                    return f(x | xr);
-                }, range_rest);
-            }, range_first, logger);
+            function_split_and_integrate_at<N>(f,integrator_rest,range_rest),
+            range_first, logger);
 	}
-
-	template<typename Bins, std::size_t DIMBINS, typename F, typename Float, typename Logger>
-	void integrate(Bins& bins, const std::array<std::size_t,DIMBINS>& bin_resolution,
-		const F& f, const RangeInfinite<Float>& range, Logger& logger) const {
-
-        static_assert(N>=DIMBINS,"Fubini does not work with that many dimensions on bin resolution");
-
-        std::array<Float,N> first_min, first_max;
-        for (std::size_t i = 0; i<N; ++i) {
-            first_min[i] = range.min(i);
-            first_max[i] = range.max(i);
-        }
-
-        std::vector<Float> rest_min = range.min();
-        std::vector<Float> rest_max = range.max();
-        for (std::size_t i = 0; i<N; ++i) {
-            if (!rest_min.empty()) rest_min.erase(rest_min.begin());
-            if (!rest_max.empty()) rest_max.erase(rest_max.begin());
-        }
-
-        viltrum::integrate(integrator_first,bins,bin_resolution,
-            [&] (const std::array<Float,N>& x) {
-                return viltrum::integrate(integrator_rest,[&f,&x] (const auto& xr) {
-                    return f(concat(x,xr));
-                }, range_infinite(rest_min,rest_max));
-            }, Range<Float,N>(first_min,first_max), logger);
-    }
 };
 
 template<std::size_t N, typename IntegratorFirst, typename IntegratorRest>

src/combination/regions-generator-fubini.h

@@ -0,0 +1,31 @@
+#pragma once
+
+#include "fubini.h"
+
+namespace viltrum {
+
+template<typename BaseGenerator, typename RestIntegrator, std::size_t N>
+class RegionsGeneratorFubini {
+    BaseGenerator base_generator;
+    RestIntegrator rest_integrator;
+
+public:
+    RegionsGeneratorFubini(const BaseGenerator& bg, const RestIntegrator& rs) 
+        : base_generator(bg), rest_integrator(rs) {}
+
+    template<std::size_t DIMBINS, typename F, typename IntegrationRange, typename Logger>
+    auto generate(const std::array<std::size_t,DIMBINS>& bin_resolution,
+		const F& f, const IntegrationRange& range, Logger& logger) const {
+            auto [range_first, range_rest] = range_split_at<N>(range);
+            auto f_gen = function_split_and_integrate_at<N>(f,rest_integrator,range_rest);
+            return base_generator.generate(bin_resolution, f_gen, range_first, logger); 
+        }
+};
+
+template<std::size_t N, typename BaseGenerator, typename RestIntegrator>
+auto regions_generator_fubini(const BaseGenerator& bg, const RestIntegrator& ri) {
+    return RegionsGeneratorFubini<BaseGenerator,RestIntegrator,N>(bg,ri);
+}
+
+
+}

src/control-variates/integrator-adaptive-fubini-variance-reduction.h

@@ -0,0 +1,32 @@
+#pragma once
+#include "integrator-adaptive-variance-reduction.h"
+#include "../combination/regions-generator-fubini.h"
+#include "../monte-carlo/monte-carlo.h"
+
+
+
+
+
+
+
+
+namespace viltrum {
+
+template<std::size_t N, typename RR, typename CV, typename RS, typename R, typename EH, typename RNG>
+auto integrator_adaptive_fubini_variance_reduction_parallel(const R& rule, const EH& error_heuristic, std::size_t iterations, unsigned long mc_samples, RR&& rr, CV&& cv, RS&& region_sampler, RNG&& rng, unsigned long spp, std::size_t nmutexes = 16) {
+    RNG rn = rng;
+    return integrator_region_based(
+        regions_generator_fubini<N>(regions_generator_adaptive_heap(rule,error_heuristic,iterations),monte_carlo(rn,mc_samples)),
+        regions_integrator_parallel_variance_reduction(std::forward<RR>(rr),std::forward<CV>(cv), std::forward<RS>(region_sampler), rn,spp,nmutexes));
+}
+
+template<std::size_t N, typename RR, typename CV, typename RS, typename R, typename EH>
+auto integrator_adaptive_fubini_variance_reduction_parallel(const R& rule, const EH& error_heuristic, std::size_t iterations, unsigned long mc_samples, 
+        RR&& rr, CV&& cv, RS&& region_sampler, unsigned long spp, std::size_t seed = std::random_device()(), std::size_t nmutexes = 16) {
+    return integrator_region_based(
+        regions_generator_fubini<N>(regions_generator_adaptive_heap(rule,error_heuristic,iterations),monte_carlo(std::mt19937(seed),mc_samples)),
+        regions_integrator_parallel_variance_reduction(std::forward<RR>(rr),std::forward<CV>(cv), std::forward<RS>(region_sampler), std::mt19937(seed+1),spp,nmutexes));
+}
+
+
+}

src/control-variates/integrator-adaptive-variance-reduction.h

@@ -8,6 +8,11 @@
 
 namespace viltrum {
 
+template<typename RR, typename CV, typename RS, typename R, typename EH, typename RNG>
+auto integrator_adaptive_variance_reduction_parallel(const R& rule, const EH& error_heuristic, std::size_t iterations, RR& rr, CV&& cv, RS&& region_sampler, RNG&& rng, unsigned long spp, std::size_t nmutexes = 16) {
+    return integrator_region_based(regions_generator_adaptive_heap(rule,error_heuristic,iterations),regions_integrator_parallel_variance_reduction(std::forward<RR>(rr),std::forward<CV>(cv), std::forward<RS>(region_sampler), std::forward<RNG>(rng),spp,nmutexes));
+}
+
 template<typename RR, typename CV, typename R, typename EH, typename RNG>
 auto integrator_adaptive_variance_reduction_parallel(const R& rule, const EH& error_heuristic, std::size_t iterations, RR& rr, CV&& cv, RNG&& rng, unsigned long spp, std::size_t nmutexes = 16) {
     return integrator_region_based(regions_generator_adaptive_heap(rule,error_heuristic,iterations),regions_integrator_parallel_variance_reduction(std::forward<RR>(rr),std::forward<CV>(cv), std::forward<RNG>(rng),spp,nmutexes));
@@ -18,6 +23,18 @@ auto integrator_adaptive_variance_reduction_parallel(const R& rule, const EH& er
     return integrator_region_based(regions_generator_adaptive_heap(rule,error_heuristic,iterations),regions_integrator_parallel_variance_reduction(std::forward<RR>(rr),std::forward<CV>(cv), spp,seed,nmutexes));
 }
 
+template<typename RR, typename CV, typename RS, typename R, typename EH>
+auto integrator_adaptive_variance_reduction_parallel(const R& rule, const EH& error_heuristic, std::size_t iterations, RR&& rr, CV&& cv, RS&& region_sampler, unsigned long spp, std::size_t seed = std::random_device()(), std::size_t nmutexes = 16) {
+    return integrator_region_based(regions_generator_adaptive_heap(rule,error_heuristic,iterations),regions_integrator_parallel_variance_reduction(std::forward<RR>(rr),std::forward<CV>(cv), std::forward<RS>(region_sampler), spp,seed,nmutexes));
+}
+
+template<typename RR, typename CV, typename RS, typename R>
+auto integrator_adaptive_variance_reduction_parallel(const R& rule, std::size_t iterations, RR&& rr, CV&& cv, RS&& region_sampler, unsigned long spp, std::size_t seed = std::random_device()(), std::size_t nmutexes = 16) {
+    return integrator_adaptive_variance_reduction_parallel(rule,error_heuristic_default(error_metric_absolute()),iterations,std::forward<RR>(rr),std::forward<CV>(cv), std::forward<RS>(region_sampler),spp,seed,nmutexes);
+}
+
+
+
 template<typename RR, typename CV, typename R>
 auto integrator_adaptive_variance_reduction_parallel(const R& rule, std::size_t iterations, RR&& rr, CV&& cv, unsigned long spp, std::size_t seed = std::random_device()(), std::size_t nmutexes = 16) {
     return integrator_adaptive_variance_reduction_parallel(rule,error_heuristic_default(error_metric_absolute()),iterations,std::forward<RR>(rr),std::forward<CV>(cv), spp,seed,nmutexes);

src/control-variates/region-sampling.h

@@ -0,0 +1,21 @@
+#pragma once
+#include "../range.h"
+#include <array>
+#include <random>
+
+namespace viltrum {
+
+class region_sampling_uniform {
+public:
+    template<typename R, typename Float, std::size_t DIM, typename RNG>
+    std::tuple<std::array<Float,DIM>,Float> sample(const R& reg, const Range<Float,DIM>& range, RNG& rng) const {
+        std::array<Float,DIM> sample;
+        for (std::size_t i=0;i<DIM;++i) {
+            std::uniform_real_distribution<Float> dis(range.min(i),range.max(i));
+            sample[i] = dis(rng);
+        }
+        return std::tuple<std::array<Float,DIM>,Float>(sample,range.volume());
+    } 
+};
+
+} 

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

@@ -6,27 +6,34 @@
 #include "../foreach.h"
 #include "region-russian-roulette.h"
 #include "weight-strategy.h"
-
-
+#include "region-sampling.h"
+#include "../combination/fubini.h"
+#include "../monte-carlo/monte-carlo.h"
 
 namespace viltrum {
 
-template<typename RR, typename CV, typename RNG>
+template<typename RR, typename CV, typename RS, typename RNG>
 class RegionsIntegratorParallelVarianceReduction {
     RR rr;
     CV cv;
+    RS region_sampler;
     mutable RNG rng;
     unsigned long samples;
     std::size_t nmutexes;
 
 public:
-    RegionsIntegratorParallelVarianceReduction(RR&& rr, CV&& cv, RNG&& r, unsigned long s,std::size_t n = 16) 
-        : rr(rr), cv(cv), rng(r), samples(s), nmutexes(n) {}
-	template<typename Bins, std::size_t DIMBINS, typename SeqRegions, typename F, typename Float, std::size_t DIM, typename Logger>
+    RegionsIntegratorParallelVarianceReduction(RR&& rr, CV&& cv, RS&& rs, RNG&& r, unsigned long s,std::size_t n = 16) 
+        : rr(rr), cv(cv), region_sampler(rs), rng(r), samples(s), nmutexes(n) {}
+	template<typename Bins, std::size_t DIMBINS, typename SeqRegions, typename F, typename IntegrationRange, typename Logger>
 	void integrate_regions(Bins& bins, const std::array<std::size_t,DIMBINS>& bin_resolution,
-		const SeqRegions& seq_regions, const F& f, const Range<Float,DIM>& range, Logger& logger) const {
+		const SeqRegions& seq_regions, const F& f, const IntegrationRange& range, Logger& logger) const {
 
         using Reg = typename SeqRegions::value_type;
+        using RegRange = std::decay_t<decltype(std::declval<Reg>().range())>;
+        constexpr std::size_t DIM = RegRange::size;
+        using Float = typename RegRange::value_type;
+        auto [range_first,range_rest] = range_split_at<DIM>(range); 
+        auto f_regdim = function_split_and_integrate_at<DIM>(f,monte_carlo(rng,1),range_rest);
 
         std::array<Float,DIMBINS> drange;
         for (std::size_t i=0;i<DIMBINS;++i) drange[i] = (range.max(i) - range.min(i))/Float(bin_resolution[i]);
@@ -37,16 +44,16 @@ class RegionsIntegratorParallelVarianceReduction {
         auto logger_bins = logger_step(logger, "region bin stratification");
         logger_bins.log_progress(progress,final_progress);
         std::for_each(std::execution::par_unseq, seq_regions.begin(),seq_regions.end(),[&] (const auto& r) {
-            for (auto pos : pixels_in_region(r,bin_resolution,range)) {
+            for (auto pos : pixels_in_region(r,bin_resolution,range_first)) {
                 perbin.push_at(pos,&r);
             } 
         });  
         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),
             [&] (const std::array<std::size_t, DIMBINS>& pos) {
-                using Sample = decltype(f(std::declval<std::array<Float,DIM>>()));
-                Range<Float,DIM> binrange = range;
+                using Sample = std::decay_t<decltype(f_regdim(std::declval<std::array<Float,DIM>>()))>;
+                Range<Float,DIM> binrange = range_first;
                 for (std::size_t i=0;i<DIMBINS;++i)
                     binrange = binrange.subrange_dimension(i,range.min(i)+pos[i]*drange[i],range.min(i)+(pos[i]+1)*drange[i]);
                 std::vector<std::tuple<const Reg*,Range<Float,DIM>>> regions_ranges(perbin[pos].size(),std::tuple((const Reg*)nullptr,range_primary<DIM,Float>()));
@@ -64,24 +71,25 @@ class RegionsIntegratorParallelVarianceReduction {
                 for (unsigned long i=0;i<samples;++i) {
                     auto [chosen,rrfactor] = roulette.choose(rng);
                     const auto& [r, region_bin_range]  = regions_ranges[chosen];
-                    std::array<Float,DIM> sample;
-	                for (std::size_t i=0;i<DIM;++i) {
-		                std::uniform_real_distribution<Float> dis(region_bin_range.min(i),region_bin_range.max(i));
-		                sample[i] = dis(rng);
-                    }
+                    const auto& [sample, sfactor] = region_sampler.sample(r,region_bin_range,rng); 
                     accumulator.push(
-                        f(sample)*double(factor)*rrfactor*region_bin_range.volume(),
-                        r->approximation_at(sample)*double(factor)*rrfactor*region_bin_range.volume()
+                        f_regdim(sample)*double(factor)*rrfactor*sfactor,
+                        r->approximation_at(sample)*double(factor)*rrfactor*sfactor
                     );
                 }
                 bins(pos) = accumulator.integral(approximation);    
         }   ,logger_control_variates);
     }
 };
 
+template<typename RR, typename CV, typename RS, typename RNG>
+auto regions_integrator_parallel_variance_reduction(RR&& rr, CV&& cv, RS&& rs, RNG&& rng, unsigned long samples, std::size_t nmutexes = 16, std::enable_if_t<!std::is_integral_v<RNG>,int> dummy = 0) {
+    return RegionsIntegratorParallelVarianceReduction<RR,CV,RS,RNG>(std::forward<RR>(rr), std::forward<CV>(cv), std::forward<RS>(rs), std::forward<RNG>(rng),samples,nmutexes);
+}
+
 template<typename RR, typename CV, typename RNG>
 auto regions_integrator_parallel_variance_reduction(RR&& rr, CV&& cv, RNG&& rng, unsigned long samples, std::size_t nmutexes = 16, std::enable_if_t<!std::is_integral_v<RNG>,int> dummy = 0) {
-    return RegionsIntegratorParallelVarianceReduction<RR,CV,RNG>(std::forward<RR>(rr), std::forward<CV>(cv), std::forward<RNG>(rng),samples,nmutexes);
+    return regions_integrator_parallel_variance_reduction(std::forward<RR>(rr), std::forward<CV>(cv), region_sampling_uniform(), std::forward<RNG>(rng),samples,nmutexes);
 }
 
 template<typename RR, typename CV>