boost random

main.cpp

#include <iostream>
#include <vector>
#include <algorithm>
#include <boost/random/random_device.hpp>
#include <boost/random/lognormal_distribution.hpp>

#include <tasks/task.h>
#include <utils/distributions.h>

using namespace che::task;
using namespace boost::random;

using value_type = ccd_task::value_type;
using duration_type = ccd_task::duration_type;

ccd_task::lognormal_type make_img_area(ccd_task::value_type mean, ccd_task::value_type sd) {
	const ccd_task::value_type pi = boost::math::constants::pi<ccd_task::value_type>();
	ccd_task::lognormal_type l(che::utils::make_distribution_using_moments<ccd_task::lognormal_type>(mean, sd));

	return ccd_task::lognormal_type(2*l.location() + std::log(pi), 2*l.scale());
}

ccd_task::scaled_lognormal_type make_sky_flux(ccd_task::value_type mean, ccd_task::value_type sd) {
	return ccd_task::scaled_lognormal_type(ccd_task::lognormal_type(-(mean-0.050)/1.086 + std::log(1283112.0), sd / 1.086)
, 0.0);
}

lognormal_distribution<> to_random(const ccd_task::lognormal_type& x) {
	return lognormal_distribution<>(x.location(), x.scale());
}

value_type rel_error2_value(
		const duration_type& exposure,
		const value_type& img_area,
		const value_type& sky_flux,
		const value_type& flux,
		const value_type& area,
		const value_type& f_over_mu,
		const value_type& dark,
		const value_type& readout) {
	const value_type float_exposure = std::chrono::duration_cast<std::chrono::duration<value_type, std::chrono::seconds::period> >(exposure).count();
	const value_type f_over_mu2 = std::pow(f_over_mu, 2);
	const value_type readout2 = std::pow(readout, 2);
	const value_type A = flux * area * float_exposure;
	const value_type B = std::pow(A, 2);
	const value_type const1 = area * float_exposure / B;
	const value_type const2 = f_over_mu2 * (dark * float_exposure + readout2) / B;

	return 1.0/A + const1 * img_area * sky_flux + const2 * img_area;
}

inline value_type solve_exposure_time_equation(const value_type& c0,const value_type& c1) {
	const value_type det = value_type(1) + 4 * c1 / std::pow(c0,2);
	return c0 * (1 + std::sqrt(det)) / 2;
}

value_type exposure_time_value(
		const value_type& rel_error,
		const value_type& img_area,
		const value_type& sky_flux,
		const value_type& flux,
		const value_type& area,
		const value_type& f_over_mu,
		const value_type& dark,
		const value_type& readout) {
        const value_type rel_error2 = std::pow(rel_error,2);
        const value_type f_over_mu2 = std::pow(f_over_mu,2);
        const value_type readout2 = std::pow(readout,2);
        const value_type A0 = 1.0/(flux * area * rel_error2);
        const value_type A1 = A0/flux;
        const value_type A23 = f_over_mu2 * A1 / area;
        const value_type A2 = A23 * dark;
        const value_type A3 = A23 * readout2;
	const value_type C0 = A0 + A1 * sky_flux * img_area + A2 * img_area;
	const value_type C1 = A3 * img_area;
	return solve_exposure_time_equation(C0, C1);
}

template<class G>
std::pair<value_type, value_type> exposure_time_mc(
		G& gen,
		const value_type& rel_error,
		lognormal_distribution<> img_area,
		lognormal_distribution<> sky_flux,
		const value_type& flux,
		const value_type& area,
		const value_type& f_over_mu,
		const value_type& dark,
		const value_type& readout) {
	
	std::vector<double> x;
	x.reserve(10000000);
	
	for(std::size_t i=0;i<x.capacity();++i) {
		x.push_back(exposure_time_value(rel_error,img_area(gen),sky_flux(gen),flux,area,f_over_mu,dark,readout));
	}

	std::sort(x.begin(), x.end());
	auto n = x.size();
	double delta = 1484.372/n;
	std::size_t left = n * 0.95 - n * delta;
	std::size_t right = n * 0.95 + n * delta;
	return std::make_pair(x[left], x[right]);
}
template<class G>
void run_exposure_time_test(G& gen,
	const value_type& rel_error,
	const ccd_task::lognormal_type& img_area,
	const ccd_task::scaled_lognormal_type& sky_flux,
	const value_type& flux,
	const value_type& area,
	const value_type& f_over_mu,
	const value_type& dark,
	const value_type& readout) {

	lognormal_distribution<> img_area1 = to_random(img_area);
	lognormal_distribution<> sky_flux1 = to_random(sky_flux.dist());

	auto pair = exposure_time_mc(gen, rel_error, img_area1, sky_flux1, flux, area, f_over_mu, dark, readout);
	auto q = boost::math::quantile(ccd_task::exposure_time_dist(rel_error, img_area, sky_flux, flux, area, f_over_mu, dark, readout), 0.95);

	std::cout << q << " " << pair.first << " " << pair.second << " " << 2*std::abs(q - (pair.first+pair.second)/2)/(pair.first+pair.second) << std::endl;
}

template<class G>
std::pair<value_type, value_type> rel_error2_mc(
		G& gen,
		const duration_type& exposure,
		lognormal_distribution<> img_area,
		lognormal_distribution<> sky_flux,
		const value_type& flux,
		const value_type& area,
		const value_type& f_over_mu,
		const value_type& dark,
		const value_type& readout) {
	
	std::vector<double> x;
	x.reserve(10000000);
	
	for(std::size_t i=0;i<x.capacity();++i) {
		x.push_back(rel_error2_value(exposure,img_area(gen),sky_flux(gen),flux,area,f_over_mu,dark,readout));
	}

	std::sort(x.begin(), x.end());
	auto n = x.size();
	auto it = std::equal_range(x.begin(), x.end(), 0.01);
	return std::make_pair(value_type(std::distance(it.first, x.end()))/n, value_type(std::distance(it.second, x.end()))/n);
/*	double delta = 1484.372/n;
	std::size_t left = n * 0.95 - n * delta;
	std::size_t right = n * 0.95 + n * delta;
	return std::make_pair(x[left], x[right]);*/
}
template<class G>
void run_rel_error2_test(G& gen,
	const duration_type& exposure,
	const ccd_task::lognormal_type& img_area,
	const ccd_task::scaled_lognormal_type& sky_flux,
	const value_type& flux,
	const value_type& area,
	const value_type& f_over_mu,
	const value_type& dark,
	const value_type& readout) {
	
	lognormal_distribution<> img_area1 = to_random(img_area);
	lognormal_distribution<> sky_flux1 = to_random(sky_flux.dist());

	auto pair = rel_error2_mc(gen, exposure, img_area1, sky_flux1, flux, area, f_over_mu, dark, readout);
	auto q = 1-boost::math::cdf(ccd_task::rel_error2_dist(exposure, img_area, sky_flux, flux, area, f_over_mu, dark, readout), 0.01);

	std::cout << q << " " << pair.first << " " << pair.second << " " << 2*std::abs(q - (pair.first+pair.second)/2)/(pair.first+pair.second) << std::endl;
}


int main(int argc, char** argv) {
	random_device gen;
	auto img_area = make_img_area(1.0,0.03);
	auto sky_flux = make_sky_flux(21.5,0.3);
	ccd_task::duration_type exposure(std::chrono::duration_cast<ccd_task::duration_type>(std::chrono::seconds(10)));

	run_rel_error2_test(gen, exposure, img_area, sky_flux, 0.00325, 40000, 6, 1, 12);
	run_rel_error2_test(gen, exposure, img_area, sky_flux, 0.003375, 40000, 6, 1, 12);
	run_rel_error2_test(gen, exposure, img_area, sky_flux, 0.0035, 40000, 6, 1, 12);
	run_rel_error2_test(gen, exposure, img_area, sky_flux, 0.003625, 40000, 6, 1, 12);
	run_rel_error2_test(gen, exposure, img_area, sky_flux, 0.00375, 40000, 6, 1, 12);
	run_rel_error2_test(gen, exposure, img_area, sky_flux, 0.003875, 40000, 6, 1, 12);
	run_rel_error2_test(gen, exposure, img_area, sky_flux, 0.004, 40000, 6, 1, 12);

	std::cout << "===" << std::endl;

	run_exposure_time_test(gen, 0.1, img_area, sky_flux, 0.00001, 40000, 6, 1, 12);
	run_exposure_time_test(gen, 0.1, img_area, sky_flux, 0.0001, 40000, 6, 1, 12);
	run_exposure_time_test(gen, 0.1, img_area, sky_flux, 0.001, 40000, 6, 1, 12);
	run_exposure_time_test(gen, 0.1, img_area, sky_flux, 0.01, 40000, 6, 1, 12);
	run_exposure_time_test(gen, 0.1, img_area, sky_flux, 0.1, 40000, 6, 1, 12);
	run_exposure_time_test(gen, 0.1, img_area, sky_flux, 1, 40000, 6, 1, 12);
	run_exposure_time_test(gen, 0.1, img_area, sky_flux, 10, 40000, 6, 1, 12);
	
	return 0;
}

Makefile

a.out: main.cpp
	g++ -O3 -DNDEBUG -DEIGEN_NO_DEBUG --std=c++11 main.cpp -o a.out -I/usr/include/eigen3 -I/home/matwey/lab/chelyabinsk/include -L/home/matwey/lab/chelyabinsk/build -lboost_random -lchelyabinsk