diff --git a/c_library/headers/ImagineModelsRandom/RandomScalarField.h b/c_library/headers/ImagineModelsRandom/RandomScalarField.h index 19039d1..48a32d5 100644 --- a/c_library/headers/ImagineModelsRandom/RandomScalarField.h +++ b/c_library/headers/ImagineModelsRandom/RandomScalarField.h @@ -45,6 +45,7 @@ class RandomScalarField : public RandomField { virtual double spatial_profile(const double &x, const double &y, const double &z) const = 0; double* random_numbers_on_grid(const std::array &shp, const std::array &inc, const int seed); + fftw_complex* test_random_numbers(const std::array &shp, const std::array &inc, const int seed); double* on_grid(const std::array &shp, const std::array &rpt, const std::array &inc, const int seed); diff --git a/c_library/headers/ImagineModelsRandom/random.tpp b/c_library/headers/ImagineModelsRandom/random.tpp index 519c876..25cbedb 100644 --- a/c_library/headers/ImagineModelsRandom/random.tpp +++ b/c_library/headers/ImagineModelsRandom/random.tpp @@ -32,7 +32,11 @@ void RandomField::remove_padding(double* val, const std::arra int sz = shp[2]; int n = shp[0]*shp[1]; for (int i = 1; i::seed_complex_random_numbers(fftw_complex* v int size_z = static_cast(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.; @@ -100,29 +106,34 @@ void RandomField::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; } @@ -131,6 +142,7 @@ void RandomField::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; } @@ -141,34 +153,36 @@ void RandomField::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; + // } } } diff --git a/c_library/source/randomscalarfield.cc b/c_library/source/randomscalarfield.cc index f0f78f0..fafbbd4 100644 --- a/c_library/source/randomscalarfield.cc +++ b/c_library/source/randomscalarfield.cc @@ -102,8 +102,18 @@ double* RandomScalarField::random_numbers_on_grid(const std::array &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; @@ -112,17 +122,17 @@ double* RandomScalarField::random_numbers_on_grid(const std::array &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; @@ -131,6 +141,64 @@ double* RandomScalarField::random_numbers_on_grid(const std::array &shp, return val; } + +fftw_complex* RandomScalarField::test_random_numbers(const std::array &shp, const std::array &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 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 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 &shp, const std::array &rpt, const std::array &inc, const int seed) { fftw_complex* val_comp = construct_plans(val, shp); diff --git a/c_library/test/random/test_gaussianscalar.cc b/c_library/test/random/test_gaussianscalar.cc index 6b8d8f8..f62ff6b 100644 --- a/c_library/test/random/test_gaussianscalar.cc +++ b/c_library/test/random/test_gaussianscalar.cc @@ -14,11 +14,11 @@ TEST_CASE("GaussianScalarField") { const int n_seeds = 10; - std::array seeds{23, 35, 645, 1, 75, 968876 , 323424, 89987,86786, 2342}; + std::array seeds{23, 35, 645, 1, 75, 968876 , 323424, 89987, 86786, 2342}; UNSCOPED_INFO("Start testing gaussian scalar test case"); - //const std::array shape {{150, 200, 60}}; + //const std::array shape {{150, 200, 600}}; const std::array shape {{4, 4, 3}}; const std::array refpoint {{-4., -0.1, -3.2}}; const std::array increment {{0.5, 0.01, 1.6}}; @@ -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++) { @@ -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)) ); }