9T random_sampler_1D(T (*pdf)(T, std::vector<T>), std::vector<T> parameters, T x0, T size, T domain_x[]) {
11 unsigned long steps = 100;
15 std::random_device rd;
16 std::seed_seq fullSeed{rd(), rd(), rd(), rd(), rd(), rd(), rd(), rd(), rd(), rd()};
17 std::mt19937 rng(fullSeed);
18 std::uniform_real_distribution<T> uniformDist(0.0f, 1.0f);
19 std::normal_distribution<T> normDist(0.0f, 1.0f);
21 for(
auto i = 0; i < steps; i ++) {
23 T step = normDist(rng) * size;
24 T u_sample = uniformDist(rng);
27 if(x_p >= domain_x[0] && x_p <= domain_x[1]) {
28 T acceptance = std::min(1.0, *pdf(x_p, parameters) / *pdf(x, parameters));
30 if(acceptance >= u_sample) {
39std::tuple<T, T> random_sampler_2D(T (*pdf)(T, T, std::vector<T>), std::vector<T> parameters, T x0, T y0, T size_x, T size_y, T domain_x[], T domain_y[]) {
41 unsigned long steps = 100;
46 std::random_device rd;
47 std::seed_seq fullSeed{rd(), rd(), rd(), rd(), rd(), rd(), rd(), rd(), rd(), rd()};
48 std::mt19937 rng(fullSeed);
49 std::uniform_real_distribution<T> uniformDist(0.0f, 1.0f);
50 std::normal_distribution<T> normDist(0.0f, 1.0f);
52 for(
auto i = 0; i < steps; i ++) {
54 T step_x = normDist(rng) * size_x;
55 T step_y = normDist(rng) * size_y;
56 T u_sample = uniformDist(rng);
60 if(x_p >= domain_x[0] && x_p <= domain_x[1] && y_p >= domain_y[0] && y_p <= domain_y[1]) {
61 T acceptance = std::min(1.0, *pdf(x_p, y_p, parameters) / *pdf(x, y, parameters));
63 if(acceptance >= u_sample) {
69 return std::make_tuple(x, y);
73std::tuple<T, T> random_sampler_angles(T (*pdf)(T, T, std::vector<T>), std::vector<T> parameters, T theta_r1, T theta_r2, T size_x, T size_y, T domain_x[], T domain_y[]) {
75 unsigned long steps = 100;
77 std::random_device rd;
78 std::seed_seq fullSeed{rd(), rd(), rd(), rd(), rd(), rd(), rd(), rd(), rd(), rd()};
79 std::mt19937 rng(fullSeed);
80 std::uniform_real_distribution<T> uniformDist(0.0f, 1.0f);
81 std::normal_distribution<T> normDist(0.0f, 1.0f);
83 for(
auto i = 0; i < steps; i ++) {
85 T step_x = normDist(rng) * size_x;
86 T step_y = normDist(rng) * size_y;
87 T u_sample = uniformDist(rng);
88 T theta_r1_p = theta_r1 + step_x;
89 T theta_r2_p = theta_r2 + step_y;
91 if(theta_r1_p < 0.0) {
92 theta_r1_p = std::abs(theta_r1_p);
95 if(theta_r2_p < 0.0) {
96 theta_r2_p += 2.0 * M_PI;
98 if(theta_r2_p > 2.0 * M_PI) {
99 theta_r2_p -= 2.0 * M_PI;
102 if(theta_r1_p < domain_x[1]) {
103 T acceptance = std::min(1.0, *pdf(theta_r1_p, theta_r2_p, parameters) / *pdf(theta_r1, theta_r2, parameters));
106 if(acceptance >= u_sample) {
107 theta_r1 = theta_r1_p;
108 theta_r2 = theta_r2_p;
112 return std::make_tuple(theta_r1, theta_r2);