Skip to content

Commit

Permalink
ongoing work on random number generation
Browse files Browse the repository at this point in the history
  • Loading branch information
shutsch committed Feb 27, 2024
1 parent 7962ccc commit ff9865c
Show file tree
Hide file tree
Showing 4 changed files with 120 additions and 36 deletions.
1 change: 1 addition & 0 deletions c_library/headers/ImagineModelsRandom/RandomScalarField.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -45,6 +45,7 @@ class RandomScalarField : public RandomField<number, double*> {
virtual double spatial_profile(const double &x, const double &y, const double &z) const = 0;

double* random_numbers_on_grid(const std::array<int, 3> &shp, const std::array<double, 3> &inc, const int seed);
fftw_complex* test_random_numbers(const std::array<int, 3> &shp, const std::array<double, 3> &inc, const int seed);

double* on_grid(const std::array<int, 3> &shp, const std::array<double, 3> &rpt, const std::array<double, 3> &inc, const int seed);

Expand Down
62 changes: 38 additions & 24 deletions c_library/headers/ImagineModelsRandom/random.tpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -32,7 +32,11 @@ void RandomField<POSTYPE, GRIDTYPE>::remove_padding(double* val, const std::arra
int sz = shp[2];
int n = shp[0]*shp[1];
for (int i = 1; i<n; i++) {
std::copy(val + i*(sz + pad), val + i*(sz + pad) + sz, val + i*sz);
auto start_new = val + i*sz;
auto start_old = start_new + i*pad;
auto end_old = start_old + sz;

std::copy(start_old, end_old, start_new);
}
}

Expand All @@ -57,25 +61,27 @@ void RandomField<POSTYPE, GRIDTYPE>::seed_complex_random_numbers(fftw_complex* v

int size_z = static_cast<int>(nyquist_z) + 1;

// In the following, we walk FORTRAN style through the array. The reason for this is that implementing the complex conjugates is easier this way

for (int l = 0; l < size_z; ++l) {
const int idx_lv1 = l * shp[1] * shp[0];
const int idx_lv1 = l; //* shp[1] * shp[0];
double kz = (double)l / lz;

for (int j = 0; j < shp[1]; ++j) {
double ky = (double)j / ly;
if (j > nyquist_y)
ky -= 1./ inc[1];
const int idx_lv2 = idx_lv1 + j * shp[0];
const int idx_lv2 = idx_lv1 + j * size_z; //* shp[0];
for (int i = 0; i < shp[0]; ++i) {
double kx = (double)i / lx;
if (i > nyquist_x)
kx -= 1./ inc[0];
const int idx = idx_lv2 + i;
const int idx = idx_lv2 + i * shp[1] * size_z;

if (debug_random) {
std::cout << "At Index (i, j, k): (" << i << j << l << ")" << std::endl;
//if (debug_random) {
// std::cout << "At Index (i, j, k): (" << i << j << l << ")" << std::endl;
//std::cout << "flattened array index " << idx << std::endl;
}
//}
if (l == 0 and j == 0 and i == 0) {
// Full Monopole is set to zero, dealt with seperately in the outer scope
vec[0][0] = 0.;
Expand All @@ -100,29 +106,34 @@ void RandomField<POSTYPE, GRIDTYPE>::seed_complex_random_numbers(fftw_complex* v
bool i_is_zero_or_nyquist = (i == 0 or i == nyquist_x);
int cg_idx;
/*
Below, the random numbers are dostributed on the grid. Since we require real output, certain parts of the grid must obey hermitian symmetry. To illustrate this, below is an example with (x, y, z ) )= (4,4, 3) dimensions.
Below, the random numbers are distributed on the grid. Since we require real output, certain parts of the grid must obey hermitian symmetry. To illustrate this, below is an example of such a complex grid with (x, y, z) = (4, 4, 3) dimensions.
RX means real type, CX complex type, and C*X complex conjugate w.r.t to CX.
We number the values by type.
Just as in the script, the respective indices are i, j, l.
The l = 2 case is just the complex conjugate of l = 1 and dealt with autmomatically by FFTW.
l = 0 l = 1 (just complex numbers)
l = 0 l = 1 (just complex numbers,
no symmetry, labels omitted)
i -->
i -->
R1 C1 R2 C*1 -- real y line after x-directed trafo C C C C
C2 C3 C4 C5 -- complex y line after x-directed trafo C C C C
R3 C6 R4 C*6 -- real y line after x-directed trafo C C C C
C*5 C*4 C*3 C*2 -- complex conjugate to j=1 line after x-directed trafo C C C C
j R1 C1 R2 C*1 -- real y line after x-directed trafo C C C C
C2 C3 C4 C5 -- complex y line after x-directed trafo C C C C
| R3 C6 R4 C*6 -- real y line after x-directed trafo C C C C
v C*5 C*4 C*3 C*2 -- complex conjugate to j=1 line after x-directed trafo C C C C
Most of the logic below deals with the l=0 case and the symmetries in there.
*/

if (l_is_zero_or_nyquist) { // real z planes (l = 0 )
if (j_is_zero_or_nyquist) { // real y_lines
if (l_is_zero_or_nyquist) { // real z planes (l = 0, l=nyq_z)
if (j_is_zero_or_nyquist) { // real y_lines (j = 0, j=nyq_y)
if (i_is_zero_or_nyquist) { // line monopole or nyquist, ->draw real numbers (global monopole is dealt with earlier)
vec[idx][0] = nd(gen);
vec[idx][1] = 0.;
if (debug_random) {
var += vec[idx][0]*vec[idx][0];
no_of_rand += 1;
no_of_real += 1;
}
Expand All @@ -131,6 +142,7 @@ void RandomField<POSTYPE, GRIDTYPE>::seed_complex_random_numbers(fftw_complex* v
vec[idx][0] = nd(gen);
vec[idx][1] = nd(gen);
if (debug_random) {
var += vec[idx][0]*vec[idx][0] + vec[idx][1] * vec[idx][1];
no_of_rand += 2;
no_of_free += 1;
}
Expand All @@ -141,34 +153,36 @@ void RandomField<POSTYPE, GRIDTYPE>::seed_complex_random_numbers(fftw_complex* v
vec[idx][1] = - vec[cg_idx][1];
}
}
else { // complex lines
else { // complex y lines
if (j < nyquist_y) { // just draw numbers
vec[idx][0] = nd(gen);
vec[idx][1] = nd(gen);
if (debug_random) {
var += vec[idx][0]*vec[idx][0] + vec[idx][1] * vec[idx][1];
no_of_rand += 2;
no_of_free += 1;
}
}
else { // mirrored complex conjugate of above
cg_idx = (shp[0] - i + 1) + (shp[1] - j + 1) * shp[0] + l * shp[1] * shp[0];
else { // mirrored complex conjugate of above (note that the -1 in the mirrored i index reflects the different symmetry in x here, as we also need to consider the x monopole of the line.)
cg_idx = (shp[0] - i - 1) + (shp[1] - j) * shp[0] + l * shp[1] * shp[0];
//std::cout << " update index: " << (shp[0] - i - 1) << ", " << (shp[1] - j - 1) << std::endl;
vec[idx][0] = vec[cg_idx][0];
vec[idx][1] = - vec[cg_idx][1];
}
}
}
else { // complex z planes
else { // complex z planes, just draw random numbers. Complex conjugate dealt with by FFTW.
vec[idx][0] = nd(gen);
vec[idx][1] = nd(gen);
if (debug_random) {
var += vec[idx][0]*vec[idx][0] + vec[idx][1] * vec[idx][1];
no_of_rand += 2;
no_of_free += 1;
}
}
if (debug_random) {
var += vec[idx][0]*vec[idx][0] + vec[idx][1] * vec[idx][1];
std::cout << "Array val (real, imag): " << vec[idx][0] << ", " << vec[idx][1] << std::endl;
}
//if (debug_random) {
// std::cout << "Array val (real, imag): " << vec[idx][0] << ", " << vec[idx][1] << std::endl;
// }

}
}
Expand Down
80 changes: 74 additions & 6 deletions c_library/source/randomscalarfield.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -102,8 +102,18 @@ double* RandomScalarField::random_numbers_on_grid(const std::array<int, 3> &shp,

seed_complex_random_numbers(val_comp, shp, inc, seed);
fftw_execute(c2r);


//std::cout << "val[s] before padding: \n " << std::endl;
//for (int s = 0; s < padded_size; s++) {
// std::cout << s << ": " << val[s]/sqrt_gs << std::endl;
//}
remove_padding(val, shp, pad);

//std::cout << "val[s] after padding: \n " << std::endl;
//for (int s = 0; s < padded_size; s++) {
// std::cout << s << ": " << val[s]/sqrt_gs << std::endl;
//}


std::cout << "gs: " << gs << std::endl;

Expand All @@ -112,17 +122,17 @@ double* RandomScalarField::random_numbers_on_grid(const std::array<int, 3> &shp,

for (int s = 0; s < gs; s++) {
val[s] /= sqrt_gs;
mean = mean + val[s];
if (s!=gs-1) mean = mean + val[s];
}

mean = mean / gs;
mean = mean / (gs - 1);

for (int s = 0; s < gs; s++) {
var = var + val[s]*val[s];
for (int s = 0; s < gs - 1; s++) {
var = var + (val[s] -mean)*(val[s] - mean);
}


var = var /(gs - 1);
var = var /(gs - 2);

std::cout << "MEAN: " << mean << std::endl;
std::cout << "VAR: " << var << std::endl;
Expand All @@ -131,6 +141,64 @@ double* RandomScalarField::random_numbers_on_grid(const std::array<int, 3> &shp,
return val;
}


fftw_complex* RandomScalarField::test_random_numbers(const std::array<int, 3> &shp, const std::array<double, 3> &inc, const int seed) {
double* val = allocate_memory(shp);
fftw_complex* val_comp = construct_plans(val, shp);
int gs = grid_size(shp);
double sqrt_gs = std::sqrt(gs);
std::array<int, 3> padded_shp = {shp[0], shp[1], 2*(shp[2]/2 + 1)};
int padded_size = grid_size(padded_shp);
int pad = padded_shp[2] - shp[2];

auto gen = std::mt19937(seed);
std::normal_distribution<double> nd{0., 1.};

for (int i = 0; i < shp[0]; ++i) {
const int idx_lv1 = i * shp[1] * padded_shp[2];
for (int j = 0; j < shp[1]; ++j) {
const int idx_lv2 = idx_lv1 + j * padded_shp[2];
for (int k = 0; k < shp[2]; ++k) {
const int idx = idx_lv2 + k;
val[idx] = nd(gen);
}
}
}


fftw_execute(r2c);

std::cout << "test: val[s]: \n " << std::endl;
for (int s = 0; s < padded_size; s++) {
std::cout << s << ": " << val[s] << std::endl;
}

std::cout << "test val_comp[s]: \n " << std::endl;
for (int s = 0; s < padded_size; s++) {
std::cout << s << ": " << val_comp[s][0]/sqrt_gs << ", " << val_comp[s][1]/sqrt_gs << std::endl;
}


std::cout << "gs: " << gs << std::endl;

double var = 0.;
double mean = 0.;

for (int s = 0; s < gs; s++) {
mean = mean + val[s];
var = var + val[s]*val[s];
}

mean = mean /gs;
var = var / gs;

std::cout << "MEAN: " << mean << std::endl;
std::cout << "VAR: " << var << std::endl;


return val_comp;
}

void RandomScalarField::_on_grid(double* val, const std::array<int, 3> &shp, const std::array<double, 3> &rpt, const std::array<double, 3> &inc, const int seed) {

fftw_complex* val_comp = construct_plans(val, shp);
Expand Down
13 changes: 7 additions & 6 deletions c_library/test/random/test_gaussianscalar.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,11 +14,11 @@ TEST_CASE("GaussianScalarField") {

const int n_seeds = 10;

std::array<int, n_seeds> seeds{23, 35, 645, 1, 75, 968876 , 323424, 89987,86786, 2342};
std::array<int, n_seeds> seeds{23, 35, 645, 1, 75, 968876 , 323424, 89987, 86786, 2342};

UNSCOPED_INFO("Start testing gaussian scalar test case");

//const std::array<int, 3> shape {{150, 200, 60}};
//const std::array<int, 3> shape {{150, 200, 600}};
const std::array<int, 3> shape {{4, 4, 3}};
const std::array<double, 3> refpoint {{-4., -0.1, -3.2}};
const std::array<double, 3> increment {{0.5, 0.01, 1.6}};
Expand Down Expand Up @@ -49,6 +49,7 @@ TEST_CASE("GaussianScalarField") {
double sample_mean = 0.;
double sample_variance = 0.;

fftw_complex* ev_random_c = gauss_plain.test_random_numbers(shape, increment, seeds[i]);
double* ev_random = gauss_plain.random_numbers_on_grid(shape, increment, seeds[i]);

for (int j = 0; j < size - 1; j++) {
Expand All @@ -58,18 +59,18 @@ TEST_CASE("GaussianScalarField") {

sample_mean = sample_mean / (size - 1);

CHECK_THAT(sample_mean, Catch::Matchers::WithinAbs(gauss_plain.mean, 2.* variance_of_the_sample_mean) );
CHECK_THAT(sample_mean, Catch::Matchers::WithinAbs(gauss_plain.mean, 2.* std::sqrt(variance_of_the_sample_mean)) );

for (int j = 0; j < size - 1; j++) {
double diff = ev_random[j] - sample_mean;
double diff = ev_random[j] - gauss_plain.mean;
sample_variance = sample_variance + diff*diff;
}
sample_variance = sample_variance / (size - 2) ; /// M_PI;
sample_variance = sample_variance / (size - 1) ; /// M_PI;

std::cout << "sample_mean: " << sample_mean << std::endl;
std::cout << "sample_variance: " << sample_variance << std::endl;

CHECK_THAT(sample_variance, Catch::Matchers::WithinAbs(model_var, 2.*variance_of_the_sample_variance) );
CHECK_THAT(sample_variance, Catch::Matchers::WithinAbs(model_var, 2.*std::sqrt(variance_of_the_sample_variance)) );
}


Expand Down

0 comments on commit ff9865c

Please sign in to comment.