Skip to content

Commit

Permalink
WIP
Browse files Browse the repository at this point in the history
  • Loading branch information
shutsch committed Mar 28, 2024
1 parent c4ae22a commit 5a9d1e3
Show file tree
Hide file tree
Showing 5 changed files with 135 additions and 41 deletions.
2 changes: 1 addition & 1 deletion c_library/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -66,7 +66,7 @@ if (IMAGINEMODELS_COMPILE_TESTS)
set(TESTSOURCES uniform helix regular)

if (IMAGINEMODELS_USE_FFTW)
set(TESTSOURCES_RANDOM) #gaussianscalar)
set(TESTSOURCES_RANDOM gaussianscalar)
endif()

add_executable(test_main "${PROJECT_SOURCE_DIR}/test/test_main.cc")
Expand Down
4 changes: 2 additions & 2 deletions c_library/headers/ImagineModelsRandom/GaussianScalar.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,8 +12,8 @@ class GaussianScalarField : public RandomScalarField {
public:
using RandomScalarField :: RandomScalarField;

double mean = 0;
double rms = 1;
double mean = 0.;
double rms = 1.;

//bool apply_spectrum = true;
double spectral_offset = .001;
Expand Down
42 changes: 24 additions & 18 deletions c_library/headers/ImagineModelsRandom/random.tpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -58,25 +58,22 @@ void RandomField<POSTYPE, GRIDTYPE>::seed_complex_random_numbers(fftw_complex* v
float nyquist_x = shp[0]/2.;
float nyquist_y = shp[1]/2.;
float nyquist_z = shp[2]/2.;

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];
double kz = (double)l / lz;

for (int i = 0; i < shp[0]; ++i) {
const int idx_lv1 = i * shp[1] * size_z;
double kx = (double)i / lx;
if (i > nyquist_x)
kx -= 1./ inc[0];
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 * 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 * shp[1] * size_z;
for (int l = 0; l < size_z; ++l) {
double kz = (double)l / lz;
const int idx = idx_lv2 + l;

//if (debug_random) {
// std::cout << "At Index (i, j, k): (" << i << j << l << ")" << std::endl;
Expand Down Expand Up @@ -148,7 +145,7 @@ void RandomField<POSTYPE, GRIDTYPE>::seed_complex_random_numbers(fftw_complex* v
}
}
else { // complex conjugate on the line, mirroring the x coordinate
cg_idx = (shp[0] - i) + j * shp[0] + l * shp[1] * shp[0];
cg_idx = (shp[0] - i) * shp[1] * size_z + j * size_z + l;
vec[idx][0] = vec[cg_idx][0];
vec[idx][1] = - vec[cg_idx][1];
}
Expand All @@ -163,11 +160,20 @@ void RandomField<POSTYPE, GRIDTYPE>::seed_complex_random_numbers(fftw_complex* v
no_of_free += 1;
}
}
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 {
if (i_is_zero_or_nyquist) {
cg_idx = i * shp[1] * size_z + (shp[1] - j) * size_z + l ;
vec[idx][0] = vec[cg_idx][0];
vec[idx][1] = - vec[cg_idx][1];
}
else {
// mirrored complex conjugate of (shp[1] - j) line
//cg_idx = (shp[0] - i - 1) + (shp[1] - j) * shp[0] + l * shp[1] * shp[0];
cg_idx = (shp[0] - i - 1) * shp[1] * size_z + (shp[1] - j) * size_z + l ;
//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];
}
}
}
}
Expand Down
83 changes: 69 additions & 14 deletions c_library/source/randomscalarfield.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -141,6 +141,56 @@ double* RandomScalarField::random_numbers_on_grid(const std::array<int, 3> &shp,
return val;
}

fftw_complex* RandomScalarField::complex_random_numbers_on_grid(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];

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;

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

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

mean = mean / (gs - 1);

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


var = var /(gs - 2);

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


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);
Expand All @@ -154,30 +204,35 @@ fftw_complex* RandomScalarField::test_random_numbers(const std::array<int, 3> &s
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) {
int size_z = static_cast<int>(shp[2]/2.) + 1;

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


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 << "\n" << std::endl;

for (int i = 0; i < padded_shp[0]; ++i) {
const int idx_lv1 = i * padded_shp[1] * size_z;
for (int j = 0; j < padded_shp[1]; ++j) {
const int idx_lv2 = idx_lv1 + j * size_z;
for (int k = 0; k < size_z; ++k) {
const int idx = idx_lv2 + k;
std::cout << i << ", " << j << ", " << k << ": " << val_comp[idx][0] << ", " << val_comp[idx][1] << std::endl;

}
}
}

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

Expand Down
45 changes: 39 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,12 +14,12 @@ 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, 600}};
const std::array<int, 3> shape {{4, 4, 3}};
const std::array<int, 3> shape {{150, 200, 60}};
//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 All @@ -29,7 +29,6 @@ TEST_CASE("GaussianScalarField") {
const std::vector<double> grid_y {{4. , 6. , -0.1, 0., 0.2}};
const std::vector<double> grid_z {{-0.2, 0.8, 0.2, 0., 1.}};


GaussianScalarField gauss_plain = GaussianScalarField();

gauss_plain.apply_spectrum = false;
Expand All @@ -45,6 +44,39 @@ TEST_CASE("GaussianScalarField") {
std::cout << "variance_of_the_sample_mean: " << variance_of_the_sample_mean << std::endl;
std::cout << "variance_of_the_sample_variance: " << variance_of_the_sample_variance << std::endl;

for (int i = 0; i < n_seeds; i++) {
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++) {
sample_mean = sample_mean + ev_random[j];
//std::cout << j << ": " << ev_random[j] << ", " << sample_mean << std::endl;
}

sample_mean = sample_mean / (size - 1);

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] - gauss_plain.mean;
sample_variance = sample_variance + diff*diff;
}
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/2, Catch::Matchers::WithinAbs(model_var, 5.*std::sqrt(variance_of_the_sample_variance)) );
}

/*
gauss_plain.mean = -3.21;
gauss_plain.rms = 0.632;
for (int i = 0; i < n_seeds; i++) {
double sample_mean = 0.;
double sample_variance = 0.;
Expand All @@ -70,8 +102,9 @@ TEST_CASE("GaussianScalarField") {
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.*std::sqrt(variance_of_the_sample_variance)) );
CHECK_THAT(sample_variance/2, Catch::Matchers::WithinAbs(gauss_plain.rms * gauss_plain.rms, 5.*std::sqrt(variance_of_the_sample_variance)) );
}

*/

}

0 comments on commit 5a9d1e3

Please sign in to comment.