Skip to content

Commit

Permalink
A control-variate version that takes less memory.
Browse files Browse the repository at this point in the history
  • Loading branch information
@adolfomunoz
adolfomunoz committed Feb 19, 2024
1 parent ccddd95 commit 7acc1c8
Show file tree
Hide file tree
Showing 4 changed files with 43 additions and 103 deletions.
2 changes: 1 addition & 1 deletion main/compilation-tests/multiple-parameters-2d.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -45,7 +45,7 @@ int main() {
{
LoggerProgress logger("Control variates");
std::vector<std::vector<float>> sol(bins,std::vector<float>(bins,0.0f));
integrate(integrator_adaptive_control_variates_parallel(nested(simpson,trapezoidal),128,samples),sol,f,range_primary<3>(),logger);
integrate(integrator_adaptive_control_variates_parallel_low_memory(nested(simpson,trapezoidal),128,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<<" ";
Expand Down
16 changes: 16 additions & 0 deletions src/control-variates/integrator-adaptive-control-variates.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,7 @@
#include "../nested/nested.h"
#include "../nested/error-heuristic.h"
#include "regions-integrator-parallel-control-variates.h"
#include "regions-integrator-parallel-control-variates-low-memory.h"

namespace viltrum {

Expand All @@ -23,5 +24,20 @@ auto integrator_adaptive_control_variates_parallel(const R& rule, std::size_t it
return integrator_adaptive_control_variates_parallel(rule,error_heuristic_default(error_metric_absolute()),iterations,spp,seed,nmutexes);
}

template<typename R, typename EH, typename RNG>
auto integrator_adaptive_control_variates_parallel_low_memory(const R& rule, const EH& error_heuristic, std::size_t iterations, 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_control_variates_low_memory(std::forward<RNG>(rng),spp,nmutexes));
}

template<typename R, typename EH>
auto integrator_adaptive_control_variates_parallel_low_memory(const R& rule, const EH& error_heuristic, std::size_t iterations, 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_control_variates_low_memory(spp,seed,nmutexes));
}

template<typename R>
auto integrator_adaptive_control_variates_parallel_low_memory(const R& rule, std::size_t iterations, unsigned long spp, std::size_t seed = std::random_device()(), std::size_t nmutexes = 16) {
return integrator_adaptive_control_variates_parallel_low_memory(rule,error_heuristic_default(error_metric_absolute()),iterations,spp,seed,nmutexes);
}


}
74 changes: 0 additions & 74 deletions src/control-variates/regions-integrator-parallel-control-variates-high-memory.h
View file

This file was deleted.

54 changes: 26 additions & 28 deletions src/control-variates/regions-integrator-parallel-control-variates-low-memory.h
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,6 @@
#pragma once
#include <random>
#include "../range.h"
#include "mutexed-tensor-vector.h"

namespace viltrum {

Expand All @@ -12,7 +11,7 @@ class RegionsIntegratorParallelControlVariatesLowMemory {
std::size_t nmutexes;

public:
RegionsIntegratorParallelControlVariates(RNG&& r, unsigned long s,std::size_t n = 16)
RegionsIntegratorParallelControlVariatesLowMemory(RNG&& r, unsigned long s,std::size_t n = 16)
: 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>
void integrate_regions(Bins& bins, const std::array<std::size_t,DIMBINS>& bin_resolution,
Expand All @@ -25,50 +24,49 @@ class RegionsIntegratorParallelControlVariatesLowMemory {
std::size_t factor = 1;
for (std::size_t i=0;i<DIMBINS;++i) factor*=bin_resolution[i];
std::size_t progress = 0; std::size_t final_progress = seq_regions.size();
MutexedTensorVector<std::tuple<const Reg*,Range<Float,DIM>>,DIMBINS> perbin(bin_resolution,nmutexes);
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)) {
Range<Float,DIM> binrange = range;
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]);

Range<Float,DIM> region_bin_range = binrange.intersection(r.range());
if (!region_bin_range.empty()) perbin.push_at(pos,std::tuple(&r,region_bin_range));
}
});
logger_bins.log_progress(final_progress,final_progress);
auto logger_control_variates = logger_step(logger,"control variates");
for_each(parallel,multidimensional_range(bin_resolution),
[&] (const std::array<std::size_t, DIMBINS>& pos) {
//This is the control variate
for (const auto& [r,region_bin_range] : perbin[pos]) {
Range<Float,DIM> binrange = range;
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;
std::for_each(seq_regions.begin(),seq_regions.end(),[&] (const auto& r) {
auto region_bin_range = binrange.intersection(r.range());
if (!region_bin_range.empty()) {
regions_ranges.push_back(std::tuple(&r,region_bin_range));
}
});

//This is for the control variate
for (const auto& [r,region_bin_range] : regions_ranges) {
bins(pos) += double(factor)*r->integral_subrange(region_bin_range);
}
}
//This is the MonteCarlo residual, RR among regions
std::uniform_int_distribution<std::size_t> rr(0,perbin[pos].size()-1);
std::uniform_int_distribution<std::size_t> rr(0,regions_ranges.size()-1);
for (unsigned long i=0;i<samples;++i) {
std::size_t chosen = rr(rng);
const auto& [r,region_bin_range] = perbin[pos][chosen];
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);
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);
}
bins(pos) += (f(sample)-r->approximation_at(sample))*double(factor)*region_bin_range.volume()*double(perbin[pos].size())/double(samples);
bins(pos) += (f(sample)-r->approximation_at(sample))*double(factor)*region_bin_range.volume()*double(regions_ranges.size())/double(samples);
}
} ,logger_control_variates);
}
};

template<typename RNG>
auto regions_integrator_parallel_control_variates(RNG&& rng, unsigned long samples, std::size_t nmutexes = 16, std::enable_if_t<!std::is_integral_v<RNG>,int> dummy = 0) {
return RegionsIntegratorParallelControlVariates<RNG>(std::forward<RNG>(rng),samples,nmutexes);
auto regions_integrator_parallel_control_variates_low_memory(RNG&& rng, unsigned long samples, std::size_t nmutexes = 16, std::enable_if_t<!std::is_integral_v<RNG>,int> dummy = 0) {
return RegionsIntegratorParallelControlVariatesLowMemory<RNG>(std::forward<RNG>(rng),samples,nmutexes);
}

auto regions_integrator_parallel_control_variates(unsigned long samples, std::size_t seed = std::random_device()(), std::size_t nmutexes = 16) {
return regions_integrator_parallel_control_variates(std::mt19937(seed),samples,nmutexes);
auto regions_integrator_parallel_control_variates_low_memory(unsigned long samples, std::size_t seed = std::random_device()(), std::size_t nmutexes = 16) {
return regions_integrator_parallel_control_variates_low_memory(std::mt19937(seed),samples,nmutexes);
}

}

0 comments on commit 7acc1c8

Please sign in to comment.