GSI_Toolbox 1.0.0 A toolbox for Gas-Surface Interaction simulations

Loading...
Searching...
No Matches
Random_Sampler.h
1#include <iostream>
2#include <cmath>
3#include <random>
4#include <vector>
5#include <tuple>
6#include <functional>
7
8template <typename T>
9T random_sampler_1D(T (*pdf)(T, std::vector<T>), std::vector<T> parameters, T x0, T size, T domain_x[]) {
10
11 unsigned long steps = 100;
12
13 T x = x0;
14
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);
20
21 for(auto i = 0; i < steps; i ++) {
22
23 T step = normDist(rng) * size;
24 T u_sample = uniformDist(rng);
25 T x_p = x + step;
26 T acceptance = 0.0;
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));
29 }
30 if(acceptance >= u_sample) {
31 x = x_p;
32 }
33 }
34
35 return x;
36}
37
38template <typename T>
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[]) {
40
41 unsigned long steps = 100;
42
43 T x = x0;
44 T y = y0;
45
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);
51
52 for(auto i = 0; i < steps; i ++) {
53
54 T step_x = normDist(rng) * size_x;
55 T step_y = normDist(rng) * size_y;
56 T u_sample = uniformDist(rng);
57 T x_p = x + step_x;
58 T y_p = y + step_y;
59 T acceptance = 0.0;
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));
62 }
63 if(acceptance >= u_sample) {
64 x = x_p;
65 y = y_p;
66 }
67 }
68
69 return std::make_tuple(x, y);
70}
71
72template <typename T>
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[]) {
74
75 unsigned long steps = 100;
76
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);
82
83 for(auto i = 0; i < steps; i ++) {
84
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;
90 T acceptance = 0.0;
91 if(theta_r1_p < 0.0) {
92 theta_r1_p = std::abs(theta_r1_p);
93 theta_r2_p += M_PI;
94 }
95 if(theta_r2_p < 0.0) {
96 theta_r2_p += 2.0 * M_PI;
97 }
98 if(theta_r2_p > 2.0 * M_PI) {
99 theta_r2_p -= 2.0 * M_PI;
100 }
101
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));
104 }
105
106 if(acceptance >= u_sample) {
107 theta_r1 = theta_r1_p;
108 theta_r2 = theta_r2_p;
109 }
110 }
111
112 return std::make_tuple(theta_r1, theta_r2);
113}