From 63f0e5fedbec612614b6c23332d6bb685e705f7a Mon Sep 17 00:00:00 2001 From: Folkert Nobels Date: Tue, 19 Feb 2019 11:33:34 +0100 Subject: [PATCH 01/35] Update the random number generator to have MWC -> Xorshift -> Quadratic congruential generator --- src/random.h | 39 ++++++++++++++++++++++++++++----------- 1 file changed, 28 insertions(+), 11 deletions(-) diff --git a/src/random.h b/src/random.h index 4d665a269..1fd8bd25e 100644 --- a/src/random.h +++ b/src/random.h @@ -33,10 +33,10 @@ * doing! */ enum random_number_type { - random_number_star_formation = 7, - random_number_stellar_feedback = 53, - random_number_stellar_enrichment = 197, - random_number_BH_feedback = 491 + random_number_star_formation = 4294957665, + random_number_stellar_feedback = 3947008974, + random_number_stellar_enrichment = 2936881968, + random_number_BH_feedback = 1640531364 }; /** @@ -59,16 +59,33 @@ INLINE static double random_unit_interval(const long long int id, /* Range used for the seeds. Best if prime */ static const long long seed_range = RAND_MAX; static const double RAND_MAX_inv = 1. / ((double)RAND_MAX); + static const long long mwc_number =(int)pow(2,32) - 1; /* Calculate the seed */ /* WARNING: Only change the math if you really know what you are doing! - The numbers are carefully chosen prime numbers that prevent correlation - with either the current integer time or the particle IDs. - The calculation overflows on purpose. */ - unsigned int seed = ((937LL * id + 1109LL) % 2147987LL + - (ti_current - 1LL) % 1514917LL + (long long)type) % - seed_range; - + * The numbers are carefully chosen prime numbers that prevent correlation + * with either the current integer time or the particle IDs. The current + * method also prevents any correlation between different random number + * types. + * The calculation overflows on purpose. + * 1. The first step is calculating the seed by using a multiply with carry + * (MWC) method, this method depends on the type of random number and + * this therefore also prevents that there is any correlation between + * the different types of random numbers. + * 2. After this we use the 64 bit Xorshift method to randomize the seeds + * even more. + * 3. We calculate a prime multiplication for the id with a quadratic + * term. + * 4. We calculate the seed by using a Quadratic congruential generator, + * in which we use the id part and the current time step bin. + */ + unsigned long long number = ti_current; + number = type * ( number & (mwc_number)) + (number >> 32); + number ^= number << 21; + number ^= number >> 35; + number ^= number << 4; + const unsigned long long idpart = 3457LL * id + 593LL * id; + unsigned int seed = (937LL * number + 5171LL * number * number + idpart + 1109LL)%9996361LL % seed_range; /* Generate a random number between 0 and 1. */ return rand_r(&seed) * RAND_MAX_inv; } -- GitLab From d207cba7c41021a2e2666cbbe894380564604911 Mon Sep 17 00:00:00 2001 From: Folkert Nobels Date: Tue, 19 Feb 2019 14:30:52 +0100 Subject: [PATCH 02/35] Update the random number generator to have a nonlinear term of the ID and the ti_current, preventing correlation between IDs --- src/random.h | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/random.h b/src/random.h index 1fd8bd25e..3acdaf3fc 100644 --- a/src/random.h +++ b/src/random.h @@ -1,6 +1,7 @@ /******************************************************************************* * This file is part of SWIFT. * Copyright (c) 2018 Matthieu Schaller (schaller@strw.leidenuniv.nl) + * 2019 Folkert Nobels (nobels@strw.leidenuniv.nl) * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published @@ -84,7 +85,7 @@ INLINE static double random_unit_interval(const long long int id, number ^= number << 21; number ^= number >> 35; number ^= number << 4; - const unsigned long long idpart = 3457LL * id + 593LL * id; + const unsigned long long idpart = 3457LL * id + 593LL * id * ti_current; unsigned int seed = (937LL * number + 5171LL * number * number + idpart + 1109LL)%9996361LL % seed_range; /* Generate a random number between 0 and 1. */ return rand_r(&seed) * RAND_MAX_inv; -- GitLab From 58a4eec1cdc2a289eccd3b0d82c559826ddc86d5 Mon Sep 17 00:00:00 2001 From: Folkert Nobels Date: Tue, 19 Feb 2019 14:46:36 +0100 Subject: [PATCH 03/35] First draft of testRandom modified but also with writing to file --- tests/testRandom.c | 33 +++++++++++++++++++++++++++++---- 1 file changed, 29 insertions(+), 4 deletions(-) diff --git a/tests/testRandom.c b/tests/testRandom.c index 4ac323070..dc5a9812a 100644 --- a/tests/testRandom.c +++ b/tests/testRandom.c @@ -41,6 +41,8 @@ int main(int argc, char* argv[]) { message("Seed = %d", seed); srand(seed); + + /* Time-step size */ const int time_bin = 29; @@ -51,10 +53,21 @@ int main(int argc, char* argv[]) { const integertime_t increment = (1LL << time_bin); message("Testing id=%lld time_bin=%d", id, time_bin); + char buffer[32]; + snprintf(buffer, sizeof(char)*32, "fileII%i.txt", i); + + FILE *fp; + fp = fopen(buffer,"w"); double total = 0., total2 = 0.; int count = 0; + /* Pearson correlation variables */ + double sum_previous_current = 0.; + double previous = 0.; + + message("Max nr timesteps = %lld",max_nr_timesteps); + /* Check that the numbers are uniform over the full-range of useful * time-steps */ for (integertime_t ti_current = 0LL; ti_current < max_nr_timesteps; @@ -68,18 +81,30 @@ int main(int argc, char* argv[]) { total += r; total2 += r * r; count++; + const unsigned int test = 127LL*(ti_current - 1LL) + 124429LL; + fprintf(fp, "%f %lld %lld\n", r, (test) % 1514917LL, ti_current ); + + /* For the pearson correlation */ + sum_previous_current += r * previous; + previous = r; } + fclose(fp); const double mean = total / (double)count; const double var = total2 / (double)count - mean * mean; + double mean_xy = sum_previous_current / ( (double)count -1.f); + double correlation = (mean_xy-mean*mean)/var; + message("Correlation = %f", correlation); + /* Verify that the mean and variance match the expected values for a uniform * distribution */ - if ((fabs(mean - 0.5) / 0.5 > 1e-4) || - (fabs(var - 1. / 12.) / (1. / 12.) > 1e-4)) { + if ((fabs(mean - 0.5) / 0.5 > 2e-4) || + (fabs(var - 1. / 12.) / (1. / 12.) > 1e-3) || + (fabs(correlation) > 3e-4)) { message("Test failed!"); - message("Result: count=%d mean=%f var=%f", count, mean, var); - message("Expected: count=%d mean=%f var=%f", count, 0.5f, 1. / 12.); + message("Result: count=%d mean=%f var=%f, correlation=%f", count, mean, var, correlation); + message("Expected: count=%d mean=%f var=%f, correlation=%f", count, 0.5f, 1. / 12., 0.); return 1; } } -- GitLab From aa97c92e73c403991723bc0310a120ae2f372346 Mon Sep 17 00:00:00 2001 From: Folkert Nobels Date: Wed, 20 Feb 2019 13:30:47 +0100 Subject: [PATCH 04/35] Add the ID and adjacent ID check --- tests/testRandom.c | 57 ++++++++++++++++++++++++++++++++++------------ 1 file changed, 42 insertions(+), 15 deletions(-) diff --git a/tests/testRandom.c b/tests/testRandom.c index dc5a9812a..4e1b71458 100644 --- a/tests/testRandom.c +++ b/tests/testRandom.c @@ -51,22 +51,28 @@ int main(int argc, char* argv[]) { const long long id = rand() * (1LL << 31) + rand(); const integertime_t increment = (1LL << time_bin); + const long long idoffset = id + 2; message("Testing id=%lld time_bin=%d", id, time_bin); - char buffer[32]; - snprintf(buffer, sizeof(char)*32, "fileII%i.txt", i); + //char buffer[32]; + //snprintf(buffer, sizeof(char)*32, "fileII%i.txt", i); - FILE *fp; - fp = fopen(buffer,"w"); + //FILE *fp; + //fp = fopen(buffer,"w"); double total = 0., total2 = 0.; int count = 0; - /* Pearson correlation variables */ + /* Pearson correlation variables for different times */ double sum_previous_current = 0.; double previous = 0.; + + /* Pearson correlation for two different IDs */ + double pearsonIDs = 0.; + double totalID = 0.; + double total2ID = 0.; - message("Max nr timesteps = %lld",max_nr_timesteps); + //message("Max nr timesteps = %lld",max_nr_timesteps); /* Check that the numbers are uniform over the full-range of useful * time-steps */ @@ -78,33 +84,54 @@ int main(int argc, char* argv[]) { const double r = random_unit_interval(id, ti_current, random_number_star_formation); + const double r_2ndid = + random_unit_interval(idoffset, ti_current, random_number_star_formation); + total += r; total2 += r * r; count++; - const unsigned int test = 127LL*(ti_current - 1LL) + 124429LL; - fprintf(fp, "%f %lld %lld\n", r, (test) % 1514917LL, ti_current ); + //const unsigned int test = 127LL*(ti_current - 1LL) + 124429LL; + //fprintf(fp, "%f %lld %lld\n", r, (test) % 1514917LL, ti_current ); - /* For the pearson correlation */ + /* For the pearson correlation of time i and i-1 */ sum_previous_current += r * previous; previous = r; + + /* Pearson correlation for small different IDs */ + pearsonIDs += r * r_2ndid; + totalID += r_2ndid; + total2ID += r_2ndid * r_2ndid; + } - fclose(fp); + //fclose(fp); const double mean = total / (double)count; const double var = total2 / (double)count - mean * mean; - double mean_xy = sum_previous_current / ( (double)count -1.f); - double correlation = (mean_xy-mean*mean)/var; + /* Pearson correlation calculation for different times */ + const double mean_xy = sum_previous_current / ( (double)count -1.f); + const double correlation = (mean_xy-mean*mean)/var; message("Correlation = %f", correlation); + + /* Pearson correlation for different IDs */ + const double meanID = totalID / (double)count; + const double varID = total2ID / (double)count - meanID * meanID; + + const double meanID_xy = pearsonIDs / (double)count; + const double correlationID = (meanID_xy - mean*meanID) / pow(var * varID, .5f); + + message("Correlation ID = %f", correlationID); + /* Verify that the mean and variance match the expected values for a uniform * distribution */ if ((fabs(mean - 0.5) / 0.5 > 2e-4) || (fabs(var - 1. / 12.) / (1. / 12.) > 1e-3) || - (fabs(correlation) > 3e-4)) { + (fabs(correlation) > 3e-4) || + (fabs(correlationID) > 3e-4) ) { message("Test failed!"); - message("Result: count=%d mean=%f var=%f, correlation=%f", count, mean, var, correlation); - message("Expected: count=%d mean=%f var=%f, correlation=%f", count, 0.5f, 1. / 12., 0.); + message("Result: count=%d mean=%f var=%f, correlation=%f ID correlation=%f", count, mean, var, correlation, correlationID); + message("Expected: count=%d mean=%f var=%f, correlation=%f ID correlation=%f", count, 0.5f, 1. / 12., 0., 0.); return 1; } } -- GitLab From e54523a1227c192a7be9d89e8cde7fd19eddd488 Mon Sep 17 00:00:00 2001 From: Folkert Nobels Date: Wed, 20 Feb 2019 13:33:02 +0100 Subject: [PATCH 05/35] Clean the testRandom.c file --- tests/testRandom.c | 17 ----------------- 1 file changed, 17 deletions(-) diff --git a/tests/testRandom.c b/tests/testRandom.c index 4e1b71458..417d8a3de 100644 --- a/tests/testRandom.c +++ b/tests/testRandom.c @@ -41,8 +41,6 @@ int main(int argc, char* argv[]) { message("Seed = %d", seed); srand(seed); - - /* Time-step size */ const int time_bin = 29; @@ -54,11 +52,6 @@ int main(int argc, char* argv[]) { const long long idoffset = id + 2; message("Testing id=%lld time_bin=%d", id, time_bin); - //char buffer[32]; - //snprintf(buffer, sizeof(char)*32, "fileII%i.txt", i); - - //FILE *fp; - //fp = fopen(buffer,"w"); double total = 0., total2 = 0.; int count = 0; @@ -71,8 +64,6 @@ int main(int argc, char* argv[]) { double pearsonIDs = 0.; double totalID = 0.; double total2ID = 0.; - - //message("Max nr timesteps = %lld",max_nr_timesteps); /* Check that the numbers are uniform over the full-range of useful * time-steps */ @@ -90,8 +81,6 @@ int main(int argc, char* argv[]) { total += r; total2 += r * r; count++; - //const unsigned int test = 127LL*(ti_current - 1LL) + 124429LL; - //fprintf(fp, "%f %lld %lld\n", r, (test) % 1514917LL, ti_current ); /* For the pearson correlation of time i and i-1 */ sum_previous_current += r * previous; @@ -101,17 +90,14 @@ int main(int argc, char* argv[]) { pearsonIDs += r * r_2ndid; totalID += r_2ndid; total2ID += r_2ndid * r_2ndid; - } - //fclose(fp); const double mean = total / (double)count; const double var = total2 / (double)count - mean * mean; /* Pearson correlation calculation for different times */ const double mean_xy = sum_previous_current / ( (double)count -1.f); const double correlation = (mean_xy-mean*mean)/var; - message("Correlation = %f", correlation); /* Pearson correlation for different IDs */ const double meanID = totalID / (double)count; @@ -120,9 +106,6 @@ int main(int argc, char* argv[]) { const double meanID_xy = pearsonIDs / (double)count; const double correlationID = (meanID_xy - mean*meanID) / pow(var * varID, .5f); - message("Correlation ID = %f", correlationID); - - /* Verify that the mean and variance match the expected values for a uniform * distribution */ if ((fabs(mean - 0.5) / 0.5 > 2e-4) || -- GitLab From 7862afa27f7f3c36af36bd46b4d0c57ba1d3343f Mon Sep 17 00:00:00 2001 From: Folkert Nobels Date: Wed, 20 Feb 2019 13:39:42 +0100 Subject: [PATCH 06/35] Add 2 file types to gitignore --- .gitignore | 2 ++ 1 file changed, 2 insertions(+) diff --git a/.gitignore b/.gitignore index e61f87609..74d17aa77 100644 --- a/.gitignore +++ b/.gitignore @@ -150,6 +150,8 @@ tests/testOutputList tests/testCbrt tests/testFormat.sh tests/testCooling +tests/*.png +tests/*.txt theory/latex/swift.pdf theory/SPH/Kernels/kernels.pdf -- GitLab From eee112e242318d93f1ff7a1bc27c02e6d73f23f6 Mon Sep 17 00:00:00 2001 From: Folkert Nobels Date: Wed, 20 Feb 2019 13:42:51 +0100 Subject: [PATCH 07/35] Format the testRandom code --- tests/testRandom.c | 27 +++++++++++++++++---------- 1 file changed, 17 insertions(+), 10 deletions(-) diff --git a/tests/testRandom.c b/tests/testRandom.c index 417d8a3de..a0f6490a2 100644 --- a/tests/testRandom.c +++ b/tests/testRandom.c @@ -1,6 +1,7 @@ /******************************************************************************* * This file is part of SWIFT. * Copyright (C) 2019 Matthieu Schaller (schaller@strw.leidenuniv.nl) + * 2019 Folkert Nobels (nobels@strw.leidenuniv.nl) * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published @@ -75,8 +76,8 @@ int main(int argc, char* argv[]) { const double r = random_unit_interval(id, ti_current, random_number_star_formation); - const double r_2ndid = - random_unit_interval(idoffset, ti_current, random_number_star_formation); + const double r_2ndid = random_unit_interval(idoffset, ti_current, + random_number_star_formation); total += r; total2 += r * r; @@ -96,25 +97,31 @@ int main(int argc, char* argv[]) { const double var = total2 / (double)count - mean * mean; /* Pearson correlation calculation for different times */ - const double mean_xy = sum_previous_current / ( (double)count -1.f); - const double correlation = (mean_xy-mean*mean)/var; - + const double mean_xy = sum_previous_current / ((double)count - 1.f); + const double correlation = (mean_xy - mean * mean) / var; + /* Pearson correlation for different IDs */ const double meanID = totalID / (double)count; const double varID = total2ID / (double)count - meanID * meanID; const double meanID_xy = pearsonIDs / (double)count; - const double correlationID = (meanID_xy - mean*meanID) / pow(var * varID, .5f); + const double correlationID = + (meanID_xy - mean * meanID) / pow(var * varID, .5f); /* Verify that the mean and variance match the expected values for a uniform * distribution */ if ((fabs(mean - 0.5) / 0.5 > 2e-4) || (fabs(var - 1. / 12.) / (1. / 12.) > 1e-3) || - (fabs(correlation) > 3e-4) || - (fabs(correlationID) > 3e-4) ) { + (fabs(correlation) > 3e-4) || (fabs(correlationID) > 3e-4)) { message("Test failed!"); - message("Result: count=%d mean=%f var=%f, correlation=%f ID correlation=%f", count, mean, var, correlation, correlationID); - message("Expected: count=%d mean=%f var=%f, correlation=%f ID correlation=%f", count, 0.5f, 1. / 12., 0., 0.); + message( + "Result: count=%d mean=%f var=%f, correlation=%f ID " + "correlation=%f", + count, mean, var, correlation, correlationID); + message( + "Expected: count=%d mean=%f var=%f, correlation=%f ID " + "correlation=%f", + count, 0.5f, 1. / 12., 0., 0.); return 1; } } -- GitLab From eca4bea786fd9ce989af04aaafd42d9ab8c05144 Mon Sep 17 00:00:00 2001 From: Folkert Nobels Date: Wed, 20 Feb 2019 13:43:59 +0100 Subject: [PATCH 08/35] code formatting of random.h --- src/random.h | 20 +++++++++++--------- 1 file changed, 11 insertions(+), 9 deletions(-) diff --git a/src/random.h b/src/random.h index 3acdaf3fc..17ab4e2de 100644 --- a/src/random.h +++ b/src/random.h @@ -60,33 +60,35 @@ INLINE static double random_unit_interval(const long long int id, /* Range used for the seeds. Best if prime */ static const long long seed_range = RAND_MAX; static const double RAND_MAX_inv = 1. / ((double)RAND_MAX); - static const long long mwc_number =(int)pow(2,32) - 1; + static const long long mwc_number = (int)pow(2, 32) - 1; /* Calculate the seed */ /* WARNING: Only change the math if you really know what you are doing! * The numbers are carefully chosen prime numbers that prevent correlation - * with either the current integer time or the particle IDs. The current - * method also prevents any correlation between different random number + * with either the current integer time or the particle IDs. The current + * method also prevents any correlation between different random number * types. - * The calculation overflows on purpose. + * The calculation overflows on purpose. * 1. The first step is calculating the seed by using a multiply with carry - * (MWC) method, this method depends on the type of random number and - * this therefore also prevents that there is any correlation between + * (MWC) method, this method depends on the type of random number and + * this therefore also prevents that there is any correlation between * the different types of random numbers. * 2. After this we use the 64 bit Xorshift method to randomize the seeds - * even more. + * even more. * 3. We calculate a prime multiplication for the id with a quadratic * term. * 4. We calculate the seed by using a Quadratic congruential generator, * in which we use the id part and the current time step bin. */ unsigned long long number = ti_current; - number = type * ( number & (mwc_number)) + (number >> 32); + number = type * (number & (mwc_number)) + (number >> 32); number ^= number << 21; number ^= number >> 35; number ^= number << 4; const unsigned long long idpart = 3457LL * id + 593LL * id * ti_current; - unsigned int seed = (937LL * number + 5171LL * number * number + idpart + 1109LL)%9996361LL % seed_range; + unsigned int seed = + (937LL * number + 5171LL * number * number + idpart + 1109LL) % + 9996361LL % seed_range; /* Generate a random number between 0 and 1. */ return rand_r(&seed) * RAND_MAX_inv; } -- GitLab From 7755f7edd8be98b79b03e7850352af6a61cec9aa Mon Sep 17 00:00:00 2001 From: Folkert Nobels Date: Tue, 19 Feb 2019 11:33:34 +0100 Subject: [PATCH 09/35] Update the random number generator to have MWC -> Xorshift -> Quadratic congruential generator --- src/random.h | 39 ++++++++++++++++++++++++++++----------- 1 file changed, 28 insertions(+), 11 deletions(-) diff --git a/src/random.h b/src/random.h index 4d665a269..1fd8bd25e 100644 --- a/src/random.h +++ b/src/random.h @@ -33,10 +33,10 @@ * doing! */ enum random_number_type { - random_number_star_formation = 7, - random_number_stellar_feedback = 53, - random_number_stellar_enrichment = 197, - random_number_BH_feedback = 491 + random_number_star_formation = 4294957665, + random_number_stellar_feedback = 3947008974, + random_number_stellar_enrichment = 2936881968, + random_number_BH_feedback = 1640531364 }; /** @@ -59,16 +59,33 @@ INLINE static double random_unit_interval(const long long int id, /* Range used for the seeds. Best if prime */ static const long long seed_range = RAND_MAX; static const double RAND_MAX_inv = 1. / ((double)RAND_MAX); + static const long long mwc_number =(int)pow(2,32) - 1; /* Calculate the seed */ /* WARNING: Only change the math if you really know what you are doing! - The numbers are carefully chosen prime numbers that prevent correlation - with either the current integer time or the particle IDs. - The calculation overflows on purpose. */ - unsigned int seed = ((937LL * id + 1109LL) % 2147987LL + - (ti_current - 1LL) % 1514917LL + (long long)type) % - seed_range; - + * The numbers are carefully chosen prime numbers that prevent correlation + * with either the current integer time or the particle IDs. The current + * method also prevents any correlation between different random number + * types. + * The calculation overflows on purpose. + * 1. The first step is calculating the seed by using a multiply with carry + * (MWC) method, this method depends on the type of random number and + * this therefore also prevents that there is any correlation between + * the different types of random numbers. + * 2. After this we use the 64 bit Xorshift method to randomize the seeds + * even more. + * 3. We calculate a prime multiplication for the id with a quadratic + * term. + * 4. We calculate the seed by using a Quadratic congruential generator, + * in which we use the id part and the current time step bin. + */ + unsigned long long number = ti_current; + number = type * ( number & (mwc_number)) + (number >> 32); + number ^= number << 21; + number ^= number >> 35; + number ^= number << 4; + const unsigned long long idpart = 3457LL * id + 593LL * id; + unsigned int seed = (937LL * number + 5171LL * number * number + idpart + 1109LL)%9996361LL % seed_range; /* Generate a random number between 0 and 1. */ return rand_r(&seed) * RAND_MAX_inv; } -- GitLab From c31113a8e3a81865d2366b4faa3ad1d0bf769d73 Mon Sep 17 00:00:00 2001 From: Folkert Nobels Date: Tue, 19 Feb 2019 14:30:52 +0100 Subject: [PATCH 10/35] Update the random number generator to have a nonlinear term of the ID and the ti_current, preventing correlation between IDs --- src/random.h | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/random.h b/src/random.h index 1fd8bd25e..3acdaf3fc 100644 --- a/src/random.h +++ b/src/random.h @@ -1,6 +1,7 @@ /******************************************************************************* * This file is part of SWIFT. * Copyright (c) 2018 Matthieu Schaller (schaller@strw.leidenuniv.nl) + * 2019 Folkert Nobels (nobels@strw.leidenuniv.nl) * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published @@ -84,7 +85,7 @@ INLINE static double random_unit_interval(const long long int id, number ^= number << 21; number ^= number >> 35; number ^= number << 4; - const unsigned long long idpart = 3457LL * id + 593LL * id; + const unsigned long long idpart = 3457LL * id + 593LL * id * ti_current; unsigned int seed = (937LL * number + 5171LL * number * number + idpart + 1109LL)%9996361LL % seed_range; /* Generate a random number between 0 and 1. */ return rand_r(&seed) * RAND_MAX_inv; -- GitLab From c157975adc9bf7e51669da132ffa6f01f44545b0 Mon Sep 17 00:00:00 2001 From: Folkert Nobels Date: Tue, 19 Feb 2019 14:46:36 +0100 Subject: [PATCH 11/35] First draft of testRandom modified but also with writing to file --- tests/testRandom.c | 33 +++++++++++++++++++++++++++++---- 1 file changed, 29 insertions(+), 4 deletions(-) diff --git a/tests/testRandom.c b/tests/testRandom.c index 4ac323070..dc5a9812a 100644 --- a/tests/testRandom.c +++ b/tests/testRandom.c @@ -41,6 +41,8 @@ int main(int argc, char* argv[]) { message("Seed = %d", seed); srand(seed); + + /* Time-step size */ const int time_bin = 29; @@ -51,10 +53,21 @@ int main(int argc, char* argv[]) { const integertime_t increment = (1LL << time_bin); message("Testing id=%lld time_bin=%d", id, time_bin); + char buffer[32]; + snprintf(buffer, sizeof(char)*32, "fileII%i.txt", i); + + FILE *fp; + fp = fopen(buffer,"w"); double total = 0., total2 = 0.; int count = 0; + /* Pearson correlation variables */ + double sum_previous_current = 0.; + double previous = 0.; + + message("Max nr timesteps = %lld",max_nr_timesteps); + /* Check that the numbers are uniform over the full-range of useful * time-steps */ for (integertime_t ti_current = 0LL; ti_current < max_nr_timesteps; @@ -68,18 +81,30 @@ int main(int argc, char* argv[]) { total += r; total2 += r * r; count++; + const unsigned int test = 127LL*(ti_current - 1LL) + 124429LL; + fprintf(fp, "%f %lld %lld\n", r, (test) % 1514917LL, ti_current ); + + /* For the pearson correlation */ + sum_previous_current += r * previous; + previous = r; } + fclose(fp); const double mean = total / (double)count; const double var = total2 / (double)count - mean * mean; + double mean_xy = sum_previous_current / ( (double)count -1.f); + double correlation = (mean_xy-mean*mean)/var; + message("Correlation = %f", correlation); + /* Verify that the mean and variance match the expected values for a uniform * distribution */ - if ((fabs(mean - 0.5) / 0.5 > 1e-4) || - (fabs(var - 1. / 12.) / (1. / 12.) > 1e-4)) { + if ((fabs(mean - 0.5) / 0.5 > 2e-4) || + (fabs(var - 1. / 12.) / (1. / 12.) > 1e-3) || + (fabs(correlation) > 3e-4)) { message("Test failed!"); - message("Result: count=%d mean=%f var=%f", count, mean, var); - message("Expected: count=%d mean=%f var=%f", count, 0.5f, 1. / 12.); + message("Result: count=%d mean=%f var=%f, correlation=%f", count, mean, var, correlation); + message("Expected: count=%d mean=%f var=%f, correlation=%f", count, 0.5f, 1. / 12., 0.); return 1; } } -- GitLab From cc7087853b14d3b2bee291046835601029df1d3d Mon Sep 17 00:00:00 2001 From: Folkert Nobels Date: Wed, 20 Feb 2019 13:30:47 +0100 Subject: [PATCH 12/35] Add the ID and adjacent ID check --- tests/testRandom.c | 57 ++++++++++++++++++++++++++++++++++------------ 1 file changed, 42 insertions(+), 15 deletions(-) diff --git a/tests/testRandom.c b/tests/testRandom.c index dc5a9812a..4e1b71458 100644 --- a/tests/testRandom.c +++ b/tests/testRandom.c @@ -51,22 +51,28 @@ int main(int argc, char* argv[]) { const long long id = rand() * (1LL << 31) + rand(); const integertime_t increment = (1LL << time_bin); + const long long idoffset = id + 2; message("Testing id=%lld time_bin=%d", id, time_bin); - char buffer[32]; - snprintf(buffer, sizeof(char)*32, "fileII%i.txt", i); + //char buffer[32]; + //snprintf(buffer, sizeof(char)*32, "fileII%i.txt", i); - FILE *fp; - fp = fopen(buffer,"w"); + //FILE *fp; + //fp = fopen(buffer,"w"); double total = 0., total2 = 0.; int count = 0; - /* Pearson correlation variables */ + /* Pearson correlation variables for different times */ double sum_previous_current = 0.; double previous = 0.; + + /* Pearson correlation for two different IDs */ + double pearsonIDs = 0.; + double totalID = 0.; + double total2ID = 0.; - message("Max nr timesteps = %lld",max_nr_timesteps); + //message("Max nr timesteps = %lld",max_nr_timesteps); /* Check that the numbers are uniform over the full-range of useful * time-steps */ @@ -78,33 +84,54 @@ int main(int argc, char* argv[]) { const double r = random_unit_interval(id, ti_current, random_number_star_formation); + const double r_2ndid = + random_unit_interval(idoffset, ti_current, random_number_star_formation); + total += r; total2 += r * r; count++; - const unsigned int test = 127LL*(ti_current - 1LL) + 124429LL; - fprintf(fp, "%f %lld %lld\n", r, (test) % 1514917LL, ti_current ); + //const unsigned int test = 127LL*(ti_current - 1LL) + 124429LL; + //fprintf(fp, "%f %lld %lld\n", r, (test) % 1514917LL, ti_current ); - /* For the pearson correlation */ + /* For the pearson correlation of time i and i-1 */ sum_previous_current += r * previous; previous = r; + + /* Pearson correlation for small different IDs */ + pearsonIDs += r * r_2ndid; + totalID += r_2ndid; + total2ID += r_2ndid * r_2ndid; + } - fclose(fp); + //fclose(fp); const double mean = total / (double)count; const double var = total2 / (double)count - mean * mean; - double mean_xy = sum_previous_current / ( (double)count -1.f); - double correlation = (mean_xy-mean*mean)/var; + /* Pearson correlation calculation for different times */ + const double mean_xy = sum_previous_current / ( (double)count -1.f); + const double correlation = (mean_xy-mean*mean)/var; message("Correlation = %f", correlation); + + /* Pearson correlation for different IDs */ + const double meanID = totalID / (double)count; + const double varID = total2ID / (double)count - meanID * meanID; + + const double meanID_xy = pearsonIDs / (double)count; + const double correlationID = (meanID_xy - mean*meanID) / pow(var * varID, .5f); + + message("Correlation ID = %f", correlationID); + /* Verify that the mean and variance match the expected values for a uniform * distribution */ if ((fabs(mean - 0.5) / 0.5 > 2e-4) || (fabs(var - 1. / 12.) / (1. / 12.) > 1e-3) || - (fabs(correlation) > 3e-4)) { + (fabs(correlation) > 3e-4) || + (fabs(correlationID) > 3e-4) ) { message("Test failed!"); - message("Result: count=%d mean=%f var=%f, correlation=%f", count, mean, var, correlation); - message("Expected: count=%d mean=%f var=%f, correlation=%f", count, 0.5f, 1. / 12., 0.); + message("Result: count=%d mean=%f var=%f, correlation=%f ID correlation=%f", count, mean, var, correlation, correlationID); + message("Expected: count=%d mean=%f var=%f, correlation=%f ID correlation=%f", count, 0.5f, 1. / 12., 0., 0.); return 1; } } -- GitLab From f4fafbfe3c8ca5b7cb1d25226401044822241402 Mon Sep 17 00:00:00 2001 From: Folkert Nobels Date: Wed, 20 Feb 2019 13:33:02 +0100 Subject: [PATCH 13/35] Clean the testRandom.c file --- tests/testRandom.c | 17 ----------------- 1 file changed, 17 deletions(-) diff --git a/tests/testRandom.c b/tests/testRandom.c index 4e1b71458..417d8a3de 100644 --- a/tests/testRandom.c +++ b/tests/testRandom.c @@ -41,8 +41,6 @@ int main(int argc, char* argv[]) { message("Seed = %d", seed); srand(seed); - - /* Time-step size */ const int time_bin = 29; @@ -54,11 +52,6 @@ int main(int argc, char* argv[]) { const long long idoffset = id + 2; message("Testing id=%lld time_bin=%d", id, time_bin); - //char buffer[32]; - //snprintf(buffer, sizeof(char)*32, "fileII%i.txt", i); - - //FILE *fp; - //fp = fopen(buffer,"w"); double total = 0., total2 = 0.; int count = 0; @@ -71,8 +64,6 @@ int main(int argc, char* argv[]) { double pearsonIDs = 0.; double totalID = 0.; double total2ID = 0.; - - //message("Max nr timesteps = %lld",max_nr_timesteps); /* Check that the numbers are uniform over the full-range of useful * time-steps */ @@ -90,8 +81,6 @@ int main(int argc, char* argv[]) { total += r; total2 += r * r; count++; - //const unsigned int test = 127LL*(ti_current - 1LL) + 124429LL; - //fprintf(fp, "%f %lld %lld\n", r, (test) % 1514917LL, ti_current ); /* For the pearson correlation of time i and i-1 */ sum_previous_current += r * previous; @@ -101,17 +90,14 @@ int main(int argc, char* argv[]) { pearsonIDs += r * r_2ndid; totalID += r_2ndid; total2ID += r_2ndid * r_2ndid; - } - //fclose(fp); const double mean = total / (double)count; const double var = total2 / (double)count - mean * mean; /* Pearson correlation calculation for different times */ const double mean_xy = sum_previous_current / ( (double)count -1.f); const double correlation = (mean_xy-mean*mean)/var; - message("Correlation = %f", correlation); /* Pearson correlation for different IDs */ const double meanID = totalID / (double)count; @@ -120,9 +106,6 @@ int main(int argc, char* argv[]) { const double meanID_xy = pearsonIDs / (double)count; const double correlationID = (meanID_xy - mean*meanID) / pow(var * varID, .5f); - message("Correlation ID = %f", correlationID); - - /* Verify that the mean and variance match the expected values for a uniform * distribution */ if ((fabs(mean - 0.5) / 0.5 > 2e-4) || -- GitLab From ac87133a3791833d3689d51ed2f8bb7700044ee5 Mon Sep 17 00:00:00 2001 From: Folkert Nobels Date: Wed, 20 Feb 2019 13:39:42 +0100 Subject: [PATCH 14/35] Add 2 file types to gitignore --- .gitignore | 2 ++ 1 file changed, 2 insertions(+) diff --git a/.gitignore b/.gitignore index 2066ba7d5..f85ebc57e 100644 --- a/.gitignore +++ b/.gitignore @@ -136,6 +136,8 @@ tests/testOutputList tests/testCbrt tests/testFormat.sh tests/testCooling +tests/*.png +tests/*.txt theory/latex/swift.pdf theory/SPH/Kernels/kernels.pdf -- GitLab From fd045b0eb5b4011f5692d29e9b17e78c8a2aefcd Mon Sep 17 00:00:00 2001 From: Folkert Nobels Date: Wed, 20 Feb 2019 13:42:51 +0100 Subject: [PATCH 15/35] Format the testRandom code --- tests/testRandom.c | 27 +++++++++++++++++---------- 1 file changed, 17 insertions(+), 10 deletions(-) diff --git a/tests/testRandom.c b/tests/testRandom.c index 417d8a3de..a0f6490a2 100644 --- a/tests/testRandom.c +++ b/tests/testRandom.c @@ -1,6 +1,7 @@ /******************************************************************************* * This file is part of SWIFT. * Copyright (C) 2019 Matthieu Schaller (schaller@strw.leidenuniv.nl) + * 2019 Folkert Nobels (nobels@strw.leidenuniv.nl) * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published @@ -75,8 +76,8 @@ int main(int argc, char* argv[]) { const double r = random_unit_interval(id, ti_current, random_number_star_formation); - const double r_2ndid = - random_unit_interval(idoffset, ti_current, random_number_star_formation); + const double r_2ndid = random_unit_interval(idoffset, ti_current, + random_number_star_formation); total += r; total2 += r * r; @@ -96,25 +97,31 @@ int main(int argc, char* argv[]) { const double var = total2 / (double)count - mean * mean; /* Pearson correlation calculation for different times */ - const double mean_xy = sum_previous_current / ( (double)count -1.f); - const double correlation = (mean_xy-mean*mean)/var; - + const double mean_xy = sum_previous_current / ((double)count - 1.f); + const double correlation = (mean_xy - mean * mean) / var; + /* Pearson correlation for different IDs */ const double meanID = totalID / (double)count; const double varID = total2ID / (double)count - meanID * meanID; const double meanID_xy = pearsonIDs / (double)count; - const double correlationID = (meanID_xy - mean*meanID) / pow(var * varID, .5f); + const double correlationID = + (meanID_xy - mean * meanID) / pow(var * varID, .5f); /* Verify that the mean and variance match the expected values for a uniform * distribution */ if ((fabs(mean - 0.5) / 0.5 > 2e-4) || (fabs(var - 1. / 12.) / (1. / 12.) > 1e-3) || - (fabs(correlation) > 3e-4) || - (fabs(correlationID) > 3e-4) ) { + (fabs(correlation) > 3e-4) || (fabs(correlationID) > 3e-4)) { message("Test failed!"); - message("Result: count=%d mean=%f var=%f, correlation=%f ID correlation=%f", count, mean, var, correlation, correlationID); - message("Expected: count=%d mean=%f var=%f, correlation=%f ID correlation=%f", count, 0.5f, 1. / 12., 0., 0.); + message( + "Result: count=%d mean=%f var=%f, correlation=%f ID " + "correlation=%f", + count, mean, var, correlation, correlationID); + message( + "Expected: count=%d mean=%f var=%f, correlation=%f ID " + "correlation=%f", + count, 0.5f, 1. / 12., 0., 0.); return 1; } } -- GitLab From 5ac11dbb58db32ff338536afa77fba00204f3d14 Mon Sep 17 00:00:00 2001 From: Folkert Nobels Date: Wed, 20 Feb 2019 13:43:59 +0100 Subject: [PATCH 16/35] code formatting of random.h --- src/random.h | 20 +++++++++++--------- 1 file changed, 11 insertions(+), 9 deletions(-) diff --git a/src/random.h b/src/random.h index 3acdaf3fc..17ab4e2de 100644 --- a/src/random.h +++ b/src/random.h @@ -60,33 +60,35 @@ INLINE static double random_unit_interval(const long long int id, /* Range used for the seeds. Best if prime */ static const long long seed_range = RAND_MAX; static const double RAND_MAX_inv = 1. / ((double)RAND_MAX); - static const long long mwc_number =(int)pow(2,32) - 1; + static const long long mwc_number = (int)pow(2, 32) - 1; /* Calculate the seed */ /* WARNING: Only change the math if you really know what you are doing! * The numbers are carefully chosen prime numbers that prevent correlation - * with either the current integer time or the particle IDs. The current - * method also prevents any correlation between different random number + * with either the current integer time or the particle IDs. The current + * method also prevents any correlation between different random number * types. - * The calculation overflows on purpose. + * The calculation overflows on purpose. * 1. The first step is calculating the seed by using a multiply with carry - * (MWC) method, this method depends on the type of random number and - * this therefore also prevents that there is any correlation between + * (MWC) method, this method depends on the type of random number and + * this therefore also prevents that there is any correlation between * the different types of random numbers. * 2. After this we use the 64 bit Xorshift method to randomize the seeds - * even more. + * even more. * 3. We calculate a prime multiplication for the id with a quadratic * term. * 4. We calculate the seed by using a Quadratic congruential generator, * in which we use the id part and the current time step bin. */ unsigned long long number = ti_current; - number = type * ( number & (mwc_number)) + (number >> 32); + number = type * (number & (mwc_number)) + (number >> 32); number ^= number << 21; number ^= number >> 35; number ^= number << 4; const unsigned long long idpart = 3457LL * id + 593LL * id * ti_current; - unsigned int seed = (937LL * number + 5171LL * number * number + idpart + 1109LL)%9996361LL % seed_range; + unsigned int seed = + (937LL * number + 5171LL * number * number + idpart + 1109LL) % + 9996361LL % seed_range; /* Generate a random number between 0 and 1. */ return rand_r(&seed) * RAND_MAX_inv; } -- GitLab From d3444f20ad3d23f6afc0bf6182d6fc81a10d4d62 Mon Sep 17 00:00:00 2001 From: Folkert Nobels Date: Thu, 21 Feb 2019 09:43:41 +0100 Subject: [PATCH 17/35] Remove the pow() and replace by bit shift --- src/random.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/random.h b/src/random.h index 17ab4e2de..5d27c0373 100644 --- a/src/random.h +++ b/src/random.h @@ -60,7 +60,7 @@ INLINE static double random_unit_interval(const long long int id, /* Range used for the seeds. Best if prime */ static const long long seed_range = RAND_MAX; static const double RAND_MAX_inv = 1. / ((double)RAND_MAX); - static const long long mwc_number = (int)pow(2, 32) - 1; + static const long long mwc_number = (1LL << 32) - 1LL; /* Calculate the seed */ /* WARNING: Only change the math if you really know what you are doing! -- GitLab From a65150da9802ab00ddc15a19c9524c4fe845525b Mon Sep 17 00:00:00 2001 From: Folkert Nobels Date: Mon, 25 Feb 2019 17:22:37 +0100 Subject: [PATCH 18/35] Add improved random number generator --- src/random.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/random.h b/src/random.h index 5d27c0373..b74e16845 100644 --- a/src/random.h +++ b/src/random.h @@ -85,7 +85,7 @@ INLINE static double random_unit_interval(const long long int id, number ^= number << 21; number ^= number >> 35; number ^= number << 4; - const unsigned long long idpart = 3457LL * id + 593LL * id * ti_current; + const unsigned long long idpart = 3457LL * id + 593LL * id * ti_current + 5417LL * id * id; unsigned int seed = (937LL * number + 5171LL * number * number + idpart + 1109LL) % 9996361LL % seed_range; -- GitLab From a35dbcdc70fb1439d053b53a161cf0ffc93540ed Mon Sep 17 00:00:00 2001 From: Folkert Nobels Date: Tue, 26 Feb 2019 10:24:21 +0100 Subject: [PATCH 19/35] Add long long prefix to enum numbers --- src/random.h | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/random.h b/src/random.h index b74e16845..50e20bbb0 100644 --- a/src/random.h +++ b/src/random.h @@ -34,10 +34,10 @@ * doing! */ enum random_number_type { - random_number_star_formation = 4294957665, - random_number_stellar_feedback = 3947008974, - random_number_stellar_enrichment = 2936881968, - random_number_BH_feedback = 1640531364 + random_number_star_formation = 4294957665LL, + random_number_stellar_feedback = 3947008974LL, + random_number_stellar_enrichment = 2936881968LL, + random_number_BH_feedback = 1640531364LL }; /** -- GitLab From 2d2601c5c497ee44144405821cf210d964bc63d5 Mon Sep 17 00:00:00 2001 From: Folkert Nobels Date: Wed, 27 Feb 2019 16:22:11 +0100 Subject: [PATCH 20/35] Add test for different types of processes --- tests/testRandom.c | 83 +++++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 78 insertions(+), 5 deletions(-) diff --git a/tests/testRandom.c b/tests/testRandom.c index a0f6490a2..524f99b04 100644 --- a/tests/testRandom.c +++ b/tests/testRandom.c @@ -26,6 +26,13 @@ /* Local headers. */ #include "swift.h" +double pearsonfunc(double mean1, double mean2, double total12, double var1, double var2, int counter) { + + const double mean12 = total12 / (double)counter; + const double correlation = (mean12 - mean1 * mean2)/ pow(var1 * var2, .5f); + return correlation; +} + int main(int argc, char* argv[]) { /* Initialize CPU frequency, this also starts time. */ @@ -66,6 +73,23 @@ int main(int argc, char* argv[]) { double totalID = 0.; double total2ID = 0.; + /* Pearson correlation for different processes */ + double pearson_star_sf = 0.; + double pearson_star_se = 0.; + double pearson_star_bh = 0.; + double pearson_sf_se = 0.; + double pearson_sf_bh = 0.; + double pearson_se_bh = 0.; + + /* Calculate the mean and for these processes */ + double total_sf = 0.; + double total_se = 0.; + double total_bh = 0.; + + double total2_sf = 0.; + double total2_se = 0.; + double total2_bh = 0.; + /* Check that the numbers are uniform over the full-range of useful * time-steps */ for (integertime_t ti_current = 0LL; ti_current < max_nr_timesteps; @@ -76,21 +100,51 @@ int main(int argc, char* argv[]) { const double r = random_unit_interval(id, ti_current, random_number_star_formation); - const double r_2ndid = random_unit_interval(idoffset, ti_current, - random_number_star_formation); - total += r; total2 += r * r; count++; - /* For the pearson correlation of time i and i-1 */ + /* Calculate for correlation between time. + * For this we use the pearson correlation of time i and i-1 */ sum_previous_current += r * previous; previous = r; + /* Calculate if there is a correlation between different ids */ + const double r_2ndid = random_unit_interval(idoffset, ti_current, + random_number_star_formation); + /* Pearson correlation for small different IDs */ pearsonIDs += r * r_2ndid; totalID += r_2ndid; total2ID += r_2ndid * r_2ndid; + + /* Calculate random numbers for the different processes and check + * that they are uncorrelated */ + + const double r_sf = + random_unit_interval(id, ti_current, random_number_stellar_feedback); + + const double r_se = + random_unit_interval(id, ti_current, random_number_stellar_enrichment); + + const double r_bh = + random_unit_interval(id, ti_current, random_number_BH_feedback); + + /* Calculate the correlation between the different processes */ + total_sf += r_sf; + total_se += r_se; + total_bh += r_bh; + + total2_sf += r_sf * r_sf; + total2_se += r_se * r_se; + total2_bh += r_bh * r_bh; + + pearson_star_sf += r * r_sf; + pearson_star_se += r * r_se; + pearson_star_bh += r * r_bh; + pearson_sf_se += r_sf * r_se; + pearson_sf_bh += r_sf * r_bh; + pearson_se_bh += r_se * r_bh; } const double mean = total / (double)count; @@ -100,14 +154,33 @@ int main(int argc, char* argv[]) { const double mean_xy = sum_previous_current / ((double)count - 1.f); const double correlation = (mean_xy - mean * mean) / var; - /* Pearson correlation for different IDs */ + /* Mean for different IDs */ const double meanID = totalID / (double)count; const double varID = total2ID / (double)count - meanID * meanID; + /* Pearson correlation between different IDs*/ const double meanID_xy = pearsonIDs / (double)count; const double correlationID = (meanID_xy - mean * meanID) / pow(var * varID, .5f); + /* Mean and for different processes */ + const double mean_sf = total_sf / (double)count; + const double mean_se = total_se / (double)count; + const double mean_bh = total_bh / (double)count; + + const double var_sf = total2_sf / (double)count - mean_sf * mean_sf; + const double var_se = total2_se / (double)count - mean_se * mean_se; + const double var_bh = total2_bh / (double)count - mean_bh * mean_bh; + + /* Correlation between different processes */ + const double corr_star_sf = pearsonfunc(mean,mean_sf,pearson_star_sf, var, var_sf, count); + const double corr_star_se = pearsonfunc(mean,mean_se,pearson_star_se, var, var_se, count); + const double corr_star_bh = pearsonfunc(mean,mean_bh,pearson_star_bh, var, var_bh, count); + const double corr_sf_se = pearsonfunc(mean_sf,mean_se,pearson_sf_se, var_sf, var_se, count); + const double corr_sf_bh = pearsonfunc(mean_sf,mean_bh,pearson_sf_bh, var_sf, var_bh, count); + const double corr_se_bh = pearsonfunc(mean_se,mean_bh,pearson_se_bh, var_se, var_bh, count); + message("%e %e %e %e %e %e",corr_star_sf, corr_star_se, corr_star_bh, corr_sf_se, corr_sf_bh, corr_se_bh); + /* Verify that the mean and variance match the expected values for a uniform * distribution */ if ((fabs(mean - 0.5) / 0.5 > 2e-4) || -- GitLab From 0d38331b8d38f0206e84d42dd961310549b85e03 Mon Sep 17 00:00:00 2001 From: Folkert Nobels Date: Wed, 27 Feb 2019 16:55:31 +0100 Subject: [PATCH 21/35] Add more critiria to check if the numbers are random and also complete the correlation between different processes --- tests/testRandom.c | 65 +++++++++++++++++++++++++++++++++++----------- 1 file changed, 50 insertions(+), 15 deletions(-) diff --git a/tests/testRandom.c b/tests/testRandom.c index 524f99b04..309fc8ae9 100644 --- a/tests/testRandom.c +++ b/tests/testRandom.c @@ -30,7 +30,7 @@ double pearsonfunc(double mean1, double mean2, double total12, double var1, doub const double mean12 = total12 / (double)counter; const double correlation = (mean12 - mean1 * mean2)/ pow(var1 * var2, .5f); - return correlation; + return fabs(correlation); } int main(int argc, char* argv[]) { @@ -151,17 +151,16 @@ int main(int argc, char* argv[]) { const double var = total2 / (double)count - mean * mean; /* Pearson correlation calculation for different times */ - const double mean_xy = sum_previous_current / ((double)count - 1.f); - const double correlation = (mean_xy - mean * mean) / var; + //const double mean_xy = sum_previous_current / ((double)count - 1.f); + //const double correlation = (mean_xy - mean * mean) / var; + const double correlation = pearsonfunc(mean,mean, sum_previous_current, var, var, count-1); /* Mean for different IDs */ const double meanID = totalID / (double)count; const double varID = total2ID / (double)count - meanID * meanID; /* Pearson correlation between different IDs*/ - const double meanID_xy = pearsonIDs / (double)count; - const double correlationID = - (meanID_xy - mean * meanID) / pow(var * varID, .5f); + const double correlationID = pearsonfunc(mean, meanID, pearsonIDs, var, varID, count); /* Mean and for different processes */ const double mean_sf = total_sf / (double)count; @@ -183,18 +182,54 @@ int main(int argc, char* argv[]) { /* Verify that the mean and variance match the expected values for a uniform * distribution */ - if ((fabs(mean - 0.5) / 0.5 > 2e-4) || - (fabs(var - 1. / 12.) / (1. / 12.) > 1e-3) || - (fabs(correlation) > 3e-4) || (fabs(correlationID) > 3e-4)) { + const double tolmean = 2e-4; + const double tolvar = 1e-3; + const double tolcorr = 3e-4; + + if ((fabs(mean - 0.5) / 0.5 > tolmean) || + (fabs(var - 1. / 12.) / (1. / 12.) > tolvar) || + (correlation > tolcorr) || (correlationID > tolcorr) || + (fabs(meanID - 0.5) / 0.5 > tolmean) || + (fabs(varID - 1. / 12.) / (1. / 12.) > tolvar) || + (corr_star_sf > tolcorr) || (corr_star_se > tolcorr) || + (corr_star_bh > tolcorr) || (corr_sf_se > tolcorr) || + (corr_sf_bh > tolcorr) || (corr_se_bh > tolcorr) || + (fabs(mean_sf - 0.5) / 0.5 > tolmean) || + (fabs(mean_se - 0.5) / 0.5 > tolmean) || + (fabs(mean_bh - 0.5) / 0.5 > tolmean) || + (fabs(var_sf - 1. / 12.) / (1. / 12.) > tolvar) || + (fabs(var_se - 1. / 12.) / (1. / 12.) > tolvar) || + (fabs(var_bh - 1. / 12.) / (1. / 12.) > tolvar)) { message("Test failed!"); + message("Global result:"); message( - "Result: count=%d mean=%f var=%f, correlation=%f ID " - "correlation=%f", - count, mean, var, correlation, correlationID); + "Result: count=%d mean=%f var=%f, correlation=%f", + count, mean, var, correlation); message( - "Expected: count=%d mean=%f var=%f, correlation=%f ID " - "correlation=%f", - count, 0.5f, 1. / 12., 0., 0.); + "Expected: count=%d mean=%f var=%f, correlation=%f", + count, 0.5f, 1. / 12., 0.); + message("ID part"); + message("Result: count%d mean=%f var=%f" + " correlation=%f", count, meanID, varID, correlationID); + message("Expected: count%d mean=%f var=%f" + " correlation=%f", count, .5f, 1. / 12., 0.); + message("Different physical processes:"); + message("Means: stars=%f stellar feedback=%f stellar " + " enrichement=%f black holes=%f", mean, mean_sf, mean_se, + mean_bh); + message("Expected: stars=%f stellar feedback=%f stellar " + " enrichement=%f black holes=%f", .5f, .5f, .5f, .5f); + message("Var: stars=%f stellar feedback=%f stellar " + " enrichement=%f black holes=%f", var, var_sf, var_se, + var_bh); + message("Expected: stars=%f stellar feedback=%f stellar " + " enrichement=%f black holes=%f", 1./12., 1./12., 1/12., + 1./12.); + message("Correlation: stars-sf=%f stars-se=%f stars-bh=%f" + "sf-se=%f sf-bh=%f se-bh=%f", corr_star_sf, corr_star_se, + corr_star_bh, corr_sf_se, corr_sf_bh, corr_se_bh); + message("Expected: stars-sf=%f stars-se=%f stars-bh=%f" + "sf-se=%f sf-bh=%f se-bh=%f", 0., 0., 0., 0., 0., 0.); return 1; } } -- GitLab From 5bced612a28a0265fce2af5880d7ef87b9beb471 Mon Sep 17 00:00:00 2001 From: Folkert Nobels Date: Thu, 28 Feb 2019 09:41:39 +0100 Subject: [PATCH 22/35] Remove my debugging message, which should not be in the master branch --- tests/testRandom.c | 1 - 1 file changed, 1 deletion(-) diff --git a/tests/testRandom.c b/tests/testRandom.c index 309fc8ae9..d6b623afe 100644 --- a/tests/testRandom.c +++ b/tests/testRandom.c @@ -178,7 +178,6 @@ int main(int argc, char* argv[]) { const double corr_sf_se = pearsonfunc(mean_sf,mean_se,pearson_sf_se, var_sf, var_se, count); const double corr_sf_bh = pearsonfunc(mean_sf,mean_bh,pearson_sf_bh, var_sf, var_bh, count); const double corr_se_bh = pearsonfunc(mean_se,mean_bh,pearson_se_bh, var_se, var_bh, count); - message("%e %e %e %e %e %e",corr_star_sf, corr_star_se, corr_star_bh, corr_sf_se, corr_sf_bh, corr_se_bh); /* Verify that the mean and variance match the expected values for a uniform * distribution */ -- GitLab From ee044fea036602295a44d88e9a054b60116709dc Mon Sep 17 00:00:00 2001 From: Folkert Nobels Date: Thu, 28 Feb 2019 09:44:28 +0100 Subject: [PATCH 23/35] IN random numbers Change pow(arg,.5f) to sqrt(arg) --- tests/testRandom.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/testRandom.c b/tests/testRandom.c index d6b623afe..44d479971 100644 --- a/tests/testRandom.c +++ b/tests/testRandom.c @@ -29,7 +29,7 @@ double pearsonfunc(double mean1, double mean2, double total12, double var1, double var2, int counter) { const double mean12 = total12 / (double)counter; - const double correlation = (mean12 - mean1 * mean2)/ pow(var1 * var2, .5f); + const double correlation = (mean12 - mean1 * mean2)/ sqrt(var1 * var2); return fabs(correlation); } -- GitLab From 0f4d937ead53054a084a1631239c691f17a2b9c5 Mon Sep 17 00:00:00 2001 From: Folkert Nobels Date: Thu, 28 Feb 2019 09:53:10 +0100 Subject: [PATCH 24/35] Update explanation for MWC random number generator a (variable) values --- src/random.h | 14 +++++++++++--- 1 file changed, 11 insertions(+), 3 deletions(-) diff --git a/src/random.h b/src/random.h index 50e20bbb0..c6d1a3362 100644 --- a/src/random.h +++ b/src/random.h @@ -29,9 +29,17 @@ /** * @brief The categories of random number generated. * - * The values of the fields are carefully chose prime - * numbers. Only change them if you know what you are - * doing! + * The values of the fields are carefully chose numbers + * These are values adviced by NR to use on page + * 348 in the first table. We used selected 4 numbers + * and know that they produce no correlation at all + * for the 4 different processes. + * Only change when you know what you are doing, changing + * the numbers to bad values will break the random number + * generator. + * In case new numbers need to be added other possible + * numbers could be: + * 4162943475, 3874257210, 2654432763 */ enum random_number_type { random_number_star_formation = 4294957665LL, -- GitLab From 2e854b7736aa818acf0e79058f4b9b19ad69ac71 Mon Sep 17 00:00:00 2001 From: Folkert Nobels Date: Thu, 28 Feb 2019 09:55:37 +0100 Subject: [PATCH 25/35] Update explanation in random function, to include reference and name of methods used --- src/random.h | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/random.h b/src/random.h index c6d1a3362..afea59f57 100644 --- a/src/random.h +++ b/src/random.h @@ -89,10 +89,13 @@ INLINE static double random_unit_interval(const long long int id, * in which we use the id part and the current time step bin. */ unsigned long long number = ti_current; + /* Multiply with carry (MWC) */ number = type * (number & (mwc_number)) + (number >> 32); + /* 64-bit Xorshift (adviced variables by NR) */ number ^= number << 21; number ^= number >> 35; number ^= number << 4; + /* Nonlinear congruential generator */ const unsigned long long idpart = 3457LL * id + 593LL * id * ti_current + 5417LL * id * id; unsigned int seed = (937LL * number + 5171LL * number * number + idpart + 1109LL) % -- GitLab From 2ebdff41e7d185487dfcb40460a2116556e7789d Mon Sep 17 00:00:00 2001 From: Folkert Nobels Date: Thu, 28 Feb 2019 10:04:50 +0100 Subject: [PATCH 26/35] Add documentation to the pearson correlation function --- tests/testRandom.c | 22 ++++++++++++++++++++++ 1 file changed, 22 insertions(+) diff --git a/tests/testRandom.c b/tests/testRandom.c index 44d479971..2be7cf5fb 100644 --- a/tests/testRandom.c +++ b/tests/testRandom.c @@ -26,6 +26,28 @@ /* Local headers. */ #include "swift.h" +/** + * @brief Compute the Pearson correlation coefficient for two sets of numbers + * + * The pearson correlation coefficient between two sets of numbers can be + * calculated as: + * + * - * + * r_xy = ---------------------- + * (var(x) * var(y))^.5 + * + * In the case that both sets are purely uncorrelated the value of the + * Pearson correlation function is expected to be close to 0. In the case that + * there is positive correlation r_xy > 0 and in the case of negative + * correlation, the function has r_xy < 0. + * + * @param mean1 average of first series of numbers + * @param mean2 average of second series of numbers + * @param total12 sum of x_i * y_i of both series of numbers + * @param var1 variance of the first series of numbers + * @param var2 variance of the second series of numbers + * @param number of elements in both series + * */ double pearsonfunc(double mean1, double mean2, double total12, double var1, double var2, int counter) { const double mean12 = total12 / (double)counter; -- GitLab From 2c4eb2d807ff06b5c3bb0ab8e169160e91007ffd Mon Sep 17 00:00:00 2001 From: Folkert Nobels Date: Thu, 28 Feb 2019 10:18:12 +0100 Subject: [PATCH 27/35] Add documentation to the main function in testRandom.c to explain shortly what the function does --- tests/testRandom.c | 20 +++++++++++++++++++- 1 file changed, 19 insertions(+), 1 deletion(-) diff --git a/tests/testRandom.c b/tests/testRandom.c index 2be7cf5fb..55b26430f 100644 --- a/tests/testRandom.c +++ b/tests/testRandom.c @@ -55,6 +55,24 @@ double pearsonfunc(double mean1, double mean2, double total12, double var1, doub return fabs(correlation); } +/** + * @brief Test to check that the pseodo-random numbers in SWIFT are random + * enough for our purpose. + * + * The test initializes with the current time and than creates 20 ID numbers + * it runs the test using these 20 ID numbers. Using these 20 ID numbers it + * Checks 4 different things: + * 1. The mean and variance are correct for random numbers generated by this + * ID number. + * 2. The random numbers from this ID number do not cause correlation in time. + * 3. A small offset in ID number of 2, doesn't cause correlation between + * the two sets of random numbers and the mean and variance of this set is + * also correct. + * 4. Different physical processes in random.h are also uncorrelated and + * produce the correct mean and variance as expected. + * + * @param none + */ int main(int argc, char* argv[]) { /* Initialize CPU frequency, this also starts time. */ @@ -205,7 +223,7 @@ int main(int argc, char* argv[]) { * distribution */ const double tolmean = 2e-4; const double tolvar = 1e-3; - const double tolcorr = 3e-4; + const double tolcorr = 4e-4; if ((fabs(mean - 0.5) / 0.5 > tolmean) || (fabs(var - 1. / 12.) / (1. / 12.) > tolvar) || -- GitLab From b406d3be5a11b69a3ee0d6f47baa5a2850116891 Mon Sep 17 00:00:00 2001 From: Folkert Nobels Date: Thu, 28 Feb 2019 10:20:48 +0100 Subject: [PATCH 28/35] Update documentation of the main function --- tests/testRandom.c | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/tests/testRandom.c b/tests/testRandom.c index 55b26430f..fd274ff35 100644 --- a/tests/testRandom.c +++ b/tests/testRandom.c @@ -65,11 +65,18 @@ double pearsonfunc(double mean1, double mean2, double total12, double var1, doub * 1. The mean and variance are correct for random numbers generated by this * ID number. * 2. The random numbers from this ID number do not cause correlation in time. + * Correlation is checked using the Pearson correlation coefficient which + * should be sufficiently close to zero. * 3. A small offset in ID number of 2, doesn't cause correlation between - * the two sets of random numbers and the mean and variance of this set is + * the two sets of random numbers (again with the Pearson correlation + * coefficient) and the mean and variance of this set is * also correct. * 4. Different physical processes in random.h are also uncorrelated and - * produce the correct mean and variance as expected. + * produce the correct mean and variance as expected. Again the correlation + * is calculated using the Pearson correlation coefficient. + * + * More information about the Pearson correlation coefficient can be found in + * the function pearsonfunc above this function. * * @param none */ -- GitLab From 13f084febd9568bab694afd66b4ddd0bbfdb2151 Mon Sep 17 00:00:00 2001 From: Folkert Nobels Date: Thu, 28 Feb 2019 10:26:37 +0100 Subject: [PATCH 29/35] Add @return to pearson function --- tests/testRandom.c | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/tests/testRandom.c b/tests/testRandom.c index fd274ff35..cf6d1f883 100644 --- a/tests/testRandom.c +++ b/tests/testRandom.c @@ -26,7 +26,7 @@ /* Local headers. */ #include "swift.h" -/** +/* * @brief Compute the Pearson correlation coefficient for two sets of numbers * * The pearson correlation coefficient between two sets of numbers can be @@ -47,6 +47,7 @@ * @param var1 variance of the first series of numbers * @param var2 variance of the second series of numbers * @param number of elements in both series + * @return the Pearson correlation coefficient * */ double pearsonfunc(double mean1, double mean2, double total12, double var1, double var2, int counter) { -- GitLab From e47790208895a41ac2d03736182ee00dd5940973 Mon Sep 17 00:00:00 2001 From: Folkert Nobels Date: Thu, 28 Feb 2019 10:29:38 +0100 Subject: [PATCH 30/35] @return in the case of main testRandom function --- tests/testRandom.c | 1 + 1 file changed, 1 insertion(+) diff --git a/tests/testRandom.c b/tests/testRandom.c index cf6d1f883..ea65046e6 100644 --- a/tests/testRandom.c +++ b/tests/testRandom.c @@ -80,6 +80,7 @@ double pearsonfunc(double mean1, double mean2, double total12, double var1, doub * the function pearsonfunc above this function. * * @param none + * @return 0 if everything is fine, 1 if random numbers are not random enough. */ int main(int argc, char* argv[]) { -- GitLab From 8874d70f63248de18cd22e424fdcbaccac97e448 Mon Sep 17 00:00:00 2001 From: Folkert Nobels Date: Thu, 28 Feb 2019 10:38:41 +0100 Subject: [PATCH 31/35] Format random.h and testRandom.c --- src/random.h | 13 ++-- tests/testRandom.c | 144 ++++++++++++++++++++++++++------------------- 2 files changed, 89 insertions(+), 68 deletions(-) diff --git a/src/random.h b/src/random.h index afea59f57..35bf47c2c 100644 --- a/src/random.h +++ b/src/random.h @@ -30,14 +30,14 @@ * @brief The categories of random number generated. * * The values of the fields are carefully chose numbers - * These are values adviced by NR to use on page - * 348 in the first table. We used selected 4 numbers - * and know that they produce no correlation at all + * These are values adviced by NR to use on page + * 348 in the first table. We used selected 4 numbers + * and know that they produce no correlation at all * for the 4 different processes. - * Only change when you know what you are doing, changing + * Only change when you know what you are doing, changing * the numbers to bad values will break the random number * generator. - * In case new numbers need to be added other possible + * In case new numbers need to be added other possible * numbers could be: * 4162943475, 3874257210, 2654432763 */ @@ -96,7 +96,8 @@ INLINE static double random_unit_interval(const long long int id, number ^= number >> 35; number ^= number << 4; /* Nonlinear congruential generator */ - const unsigned long long idpart = 3457LL * id + 593LL * id * ti_current + 5417LL * id * id; + const unsigned long long idpart = + 3457LL * id + 593LL * id * ti_current + 5417LL * id * id; unsigned int seed = (937LL * number + 5171LL * number * number + idpart + 1109LL) % 9996361LL % seed_range; diff --git a/tests/testRandom.c b/tests/testRandom.c index ea65046e6..e0dc06a88 100644 --- a/tests/testRandom.c +++ b/tests/testRandom.c @@ -26,19 +26,19 @@ /* Local headers. */ #include "swift.h" -/* +/* * @brief Compute the Pearson correlation coefficient for two sets of numbers * * The pearson correlation coefficient between two sets of numbers can be * calculated as: - * + * * - * * r_xy = ---------------------- * (var(x) * var(y))^.5 * - * In the case that both sets are purely uncorrelated the value of the - * Pearson correlation function is expected to be close to 0. In the case that - * there is positive correlation r_xy > 0 and in the case of negative + * In the case that both sets are purely uncorrelated the value of the + * Pearson correlation function is expected to be close to 0. In the case that + * there is positive correlation r_xy > 0 and in the case of negative * correlation, the function has r_xy < 0. * * @param mean1 average of first series of numbers @@ -49,11 +49,12 @@ * @param number of elements in both series * @return the Pearson correlation coefficient * */ -double pearsonfunc(double mean1, double mean2, double total12, double var1, double var2, int counter) { - +double pearsonfunc(double mean1, double mean2, double total12, double var1, + double var2, int counter) { + const double mean12 = total12 / (double)counter; - const double correlation = (mean12 - mean1 * mean2)/ sqrt(var1 * var2); - return fabs(correlation); + const double correlation = (mean12 - mean1 * mean2) / sqrt(var1 * var2); + return fabs(correlation); } /** @@ -62,21 +63,21 @@ double pearsonfunc(double mean1, double mean2, double total12, double var1, doub * * The test initializes with the current time and than creates 20 ID numbers * it runs the test using these 20 ID numbers. Using these 20 ID numbers it - * Checks 4 different things: + * Checks 4 different things: * 1. The mean and variance are correct for random numbers generated by this * ID number. * 2. The random numbers from this ID number do not cause correlation in time. * Correlation is checked using the Pearson correlation coefficient which * should be sufficiently close to zero. - * 3. A small offset in ID number of 2, doesn't cause correlation between - * the two sets of random numbers (again with the Pearson correlation - * coefficient) and the mean and variance of this set is + * 3. A small offset in ID number of 2, doesn't cause correlation between + * the two sets of random numbers (again with the Pearson correlation + * coefficient) and the mean and variance of this set is * also correct. - * 4. Different physical processes in random.h are also uncorrelated and + * 4. Different physical processes in random.h are also uncorrelated and * produce the correct mean and variance as expected. Again the correlation - * is calculated using the Pearson correlation coefficient. + * is calculated using the Pearson correlation coefficient. * - * More information about the Pearson correlation coefficient can be found in + * More information about the Pearson correlation coefficient can be found in * the function pearsonfunc above this function. * * @param none @@ -169,14 +170,14 @@ int main(int argc, char* argv[]) { /* Calculate random numbers for the different processes and check * that they are uncorrelated */ - - const double r_sf = + + const double r_sf = random_unit_interval(id, ti_current, random_number_stellar_feedback); - const double r_se = - random_unit_interval(id, ti_current, random_number_stellar_enrichment); + const double r_se = random_unit_interval( + id, ti_current, random_number_stellar_enrichment); - const double r_bh = + const double r_bh = random_unit_interval(id, ti_current, random_number_BH_feedback); /* Calculate the correlation between the different processes */ @@ -200,34 +201,42 @@ int main(int argc, char* argv[]) { const double var = total2 / (double)count - mean * mean; /* Pearson correlation calculation for different times */ - //const double mean_xy = sum_previous_current / ((double)count - 1.f); - //const double correlation = (mean_xy - mean * mean) / var; - const double correlation = pearsonfunc(mean,mean, sum_previous_current, var, var, count-1); + // const double mean_xy = sum_previous_current / ((double)count - 1.f); + // const double correlation = (mean_xy - mean * mean) / var; + const double correlation = + pearsonfunc(mean, mean, sum_previous_current, var, var, count - 1); /* Mean for different IDs */ const double meanID = totalID / (double)count; const double varID = total2ID / (double)count - meanID * meanID; /* Pearson correlation between different IDs*/ - const double correlationID = pearsonfunc(mean, meanID, pearsonIDs, var, varID, count); + const double correlationID = + pearsonfunc(mean, meanID, pearsonIDs, var, varID, count); /* Mean and for different processes */ const double mean_sf = total_sf / (double)count; const double mean_se = total_se / (double)count; const double mean_bh = total_bh / (double)count; - + const double var_sf = total2_sf / (double)count - mean_sf * mean_sf; const double var_se = total2_se / (double)count - mean_se * mean_se; const double var_bh = total2_bh / (double)count - mean_bh * mean_bh; /* Correlation between different processes */ - const double corr_star_sf = pearsonfunc(mean,mean_sf,pearson_star_sf, var, var_sf, count); - const double corr_star_se = pearsonfunc(mean,mean_se,pearson_star_se, var, var_se, count); - const double corr_star_bh = pearsonfunc(mean,mean_bh,pearson_star_bh, var, var_bh, count); - const double corr_sf_se = pearsonfunc(mean_sf,mean_se,pearson_sf_se, var_sf, var_se, count); - const double corr_sf_bh = pearsonfunc(mean_sf,mean_bh,pearson_sf_bh, var_sf, var_bh, count); - const double corr_se_bh = pearsonfunc(mean_se,mean_bh,pearson_se_bh, var_se, var_bh, count); - + const double corr_star_sf = + pearsonfunc(mean, mean_sf, pearson_star_sf, var, var_sf, count); + const double corr_star_se = + pearsonfunc(mean, mean_se, pearson_star_se, var, var_se, count); + const double corr_star_bh = + pearsonfunc(mean, mean_bh, pearson_star_bh, var, var_bh, count); + const double corr_sf_se = + pearsonfunc(mean_sf, mean_se, pearson_sf_se, var_sf, var_se, count); + const double corr_sf_bh = + pearsonfunc(mean_sf, mean_bh, pearson_sf_bh, var_sf, var_bh, count); + const double corr_se_bh = + pearsonfunc(mean_se, mean_bh, pearson_se_bh, var_se, var_bh, count); + /* Verify that the mean and variance match the expected values for a uniform * distribution */ const double tolmean = 2e-4; @@ -238,46 +247,57 @@ int main(int argc, char* argv[]) { (fabs(var - 1. / 12.) / (1. / 12.) > tolvar) || (correlation > tolcorr) || (correlationID > tolcorr) || (fabs(meanID - 0.5) / 0.5 > tolmean) || - (fabs(varID - 1. / 12.) / (1. / 12.) > tolvar) || + (fabs(varID - 1. / 12.) / (1. / 12.) > tolvar) || (corr_star_sf > tolcorr) || (corr_star_se > tolcorr) || (corr_star_bh > tolcorr) || (corr_sf_se > tolcorr) || - (corr_sf_bh > tolcorr) || (corr_se_bh > tolcorr) || + (corr_sf_bh > tolcorr) || (corr_se_bh > tolcorr) || (fabs(mean_sf - 0.5) / 0.5 > tolmean) || (fabs(mean_se - 0.5) / 0.5 > tolmean) || (fabs(mean_bh - 0.5) / 0.5 > tolmean) || - (fabs(var_sf - 1. / 12.) / (1. / 12.) > tolvar) || - (fabs(var_se - 1. / 12.) / (1. / 12.) > tolvar) || + (fabs(var_sf - 1. / 12.) / (1. / 12.) > tolvar) || + (fabs(var_se - 1. / 12.) / (1. / 12.) > tolvar) || (fabs(var_bh - 1. / 12.) / (1. / 12.) > tolvar)) { message("Test failed!"); message("Global result:"); + message("Result: count=%d mean=%f var=%f, correlation=%f", count, mean, + var, correlation); + message("Expected: count=%d mean=%f var=%f, correlation=%f", count, 0.5f, + 1. / 12., 0.); + message("ID part"); message( - "Result: count=%d mean=%f var=%f, correlation=%f", - count, mean, var, correlation); + "Result: count%d mean=%f var=%f" + " correlation=%f", + count, meanID, varID, correlationID); message( - "Expected: count=%d mean=%f var=%f, correlation=%f", - count, 0.5f, 1. / 12., 0.); - message("ID part"); - message("Result: count%d mean=%f var=%f" - " correlation=%f", count, meanID, varID, correlationID); - message("Expected: count%d mean=%f var=%f" - " correlation=%f", count, .5f, 1. / 12., 0.); + "Expected: count%d mean=%f var=%f" + " correlation=%f", + count, .5f, 1. / 12., 0.); message("Different physical processes:"); - message("Means: stars=%f stellar feedback=%f stellar " - " enrichement=%f black holes=%f", mean, mean_sf, mean_se, - mean_bh); - message("Expected: stars=%f stellar feedback=%f stellar " - " enrichement=%f black holes=%f", .5f, .5f, .5f, .5f); - message("Var: stars=%f stellar feedback=%f stellar " - " enrichement=%f black holes=%f", var, var_sf, var_se, - var_bh); - message("Expected: stars=%f stellar feedback=%f stellar " - " enrichement=%f black holes=%f", 1./12., 1./12., 1/12., - 1./12.); - message("Correlation: stars-sf=%f stars-se=%f stars-bh=%f" - "sf-se=%f sf-bh=%f se-bh=%f", corr_star_sf, corr_star_se, - corr_star_bh, corr_sf_se, corr_sf_bh, corr_se_bh); - message("Expected: stars-sf=%f stars-se=%f stars-bh=%f" - "sf-se=%f sf-bh=%f se-bh=%f", 0., 0., 0., 0., 0., 0.); + message( + "Means: stars=%f stellar feedback=%f stellar " + " enrichement=%f black holes=%f", + mean, mean_sf, mean_se, mean_bh); + message( + "Expected: stars=%f stellar feedback=%f stellar " + " enrichement=%f black holes=%f", + .5f, .5f, .5f, .5f); + message( + "Var: stars=%f stellar feedback=%f stellar " + " enrichement=%f black holes=%f", + var, var_sf, var_se, var_bh); + message( + "Expected: stars=%f stellar feedback=%f stellar " + " enrichement=%f black holes=%f", + 1. / 12., 1. / 12., 1 / 12., 1. / 12.); + message( + "Correlation: stars-sf=%f stars-se=%f stars-bh=%f" + "sf-se=%f sf-bh=%f se-bh=%f", + corr_star_sf, corr_star_se, corr_star_bh, corr_sf_se, corr_sf_bh, + corr_se_bh); + message( + "Expected: stars-sf=%f stars-se=%f stars-bh=%f" + "sf-se=%f sf-bh=%f se-bh=%f", + 0., 0., 0., 0., 0., 0.); return 1; } } -- GitLab From 957effa5fa14c790bfdf5d55f69da7b48d1ba225 Mon Sep 17 00:00:00 2001 From: Folkert Nobels Date: Thu, 28 Feb 2019 10:41:00 +0100 Subject: [PATCH 32/35] Fix doxygen documentation for pearsonfunc --- tests/testRandom.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/testRandom.c b/tests/testRandom.c index e0dc06a88..1b2a6a5d4 100644 --- a/tests/testRandom.c +++ b/tests/testRandom.c @@ -26,7 +26,7 @@ /* Local headers. */ #include "swift.h" -/* +/** * @brief Compute the Pearson correlation coefficient for two sets of numbers * * The pearson correlation coefficient between two sets of numbers can be -- GitLab From 8b553991753277260fc67131aef1deac23cf2d79 Mon Sep 17 00:00:00 2001 From: Folkert Nobels Date: Fri, 1 Mar 2019 14:26:13 +0100 Subject: [PATCH 33/35] Add random numbers update --- src/random.h | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/src/random.h b/src/random.h index 35bf47c2c..6efc957e0 100644 --- a/src/random.h +++ b/src/random.h @@ -42,7 +42,7 @@ * 4162943475, 3874257210, 2654432763 */ enum random_number_type { - random_number_star_formation = 4294957665LL, + random_number_star_formation = 0LL, //4294957665LL, random_number_stellar_feedback = 3947008974LL, random_number_stellar_enrichment = 2936881968LL, random_number_BH_feedback = 1640531364LL @@ -90,14 +90,16 @@ INLINE static double random_unit_interval(const long long int id, */ unsigned long long number = ti_current; /* Multiply with carry (MWC) */ - number = type * (number & (mwc_number)) + (number >> 32); + number = 4294957665LL * (number & (mwc_number)) + (number >> 32); /* 64-bit Xorshift (adviced variables by NR) */ number ^= number << 21; number ^= number >> 35; number ^= number << 4; + /* Add constant to ID */ + const unsigned long long idt = id + type; /* Nonlinear congruential generator */ const unsigned long long idpart = - 3457LL * id + 593LL * id * ti_current + 5417LL * id * id; + 3457LL * idt + 593LL * idt * ti_current + 5417LL * idt * idt; unsigned int seed = (937LL * number + 5171LL * number * number + idpart + 1109LL) % 9996361LL % seed_range; -- GitLab From 77479d45cad7276a52391f02548c9916e78723e7 Mon Sep 17 00:00:00 2001 From: Folkert Nobels Date: Mon, 4 Mar 2019 10:06:25 +0100 Subject: [PATCH 34/35] Change difference between different kind of random numbers to make it robust --- src/random.h | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/src/random.h b/src/random.h index 6efc957e0..9b894568c 100644 --- a/src/random.h +++ b/src/random.h @@ -30,22 +30,22 @@ * @brief The categories of random number generated. * * The values of the fields are carefully chose numbers - * These are values adviced by NR to use on page - * 348 in the first table. We used selected 4 numbers - * and know that they produce no correlation at all - * for the 4 different processes. + * the numbers are very large primes such that the IDs + * will not have a prime factorization with this coefficient + * this results in a very high period for the random number + * generator. * Only change when you know what you are doing, changing * the numbers to bad values will break the random number * generator. * In case new numbers need to be added other possible * numbers could be: - * 4162943475, 3874257210, 2654432763 + * 4947009007, 5947309451, 6977309513 */ enum random_number_type { - random_number_star_formation = 0LL, //4294957665LL, - random_number_stellar_feedback = 3947008974LL, - random_number_stellar_enrichment = 2936881968LL, - random_number_BH_feedback = 1640531364LL + random_number_star_formation = 0LL, + random_number_stellar_feedback = 3947008991LL, + random_number_stellar_enrichment = 2936881973LL, + random_number_BH_feedback = 1640531371LL }; /** @@ -89,7 +89,7 @@ INLINE static double random_unit_interval(const long long int id, * in which we use the id part and the current time step bin. */ unsigned long long number = ti_current; - /* Multiply with carry (MWC) */ + /* Multiply with carry (MWC), (adviced variables by NR) */ number = 4294957665LL * (number & (mwc_number)) + (number >> 32); /* 64-bit Xorshift (adviced variables by NR) */ number ^= number << 21; -- GitLab From ea4b41f386631c65e53ac5247e571ebf6bfb89aa Mon Sep 17 00:00:00 2001 From: Matthieu Schaller Date: Mon, 4 Mar 2019 16:39:08 +0000 Subject: [PATCH 35/35] Applied code formatting tool. --- src/random.h | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/src/random.h b/src/random.h index 9b894568c..660ae21db 100644 --- a/src/random.h +++ b/src/random.h @@ -30,7 +30,7 @@ * @brief The categories of random number generated. * * The values of the fields are carefully chose numbers - * the numbers are very large primes such that the IDs + * the numbers are very large primes such that the IDs * will not have a prime factorization with this coefficient * this results in a very high period for the random number * generator. @@ -39,10 +39,10 @@ * generator. * In case new numbers need to be added other possible * numbers could be: - * 4947009007, 5947309451, 6977309513 + * 4947009007, 5947309451, 6977309513 */ enum random_number_type { - random_number_star_formation = 0LL, + random_number_star_formation = 0LL, random_number_stellar_feedback = 3947008991LL, random_number_stellar_enrichment = 2936881973LL, random_number_BH_feedback = 1640531371LL @@ -89,20 +89,25 @@ INLINE static double random_unit_interval(const long long int id, * in which we use the id part and the current time step bin. */ unsigned long long number = ti_current; + /* Multiply with carry (MWC), (adviced variables by NR) */ number = 4294957665LL * (number & (mwc_number)) + (number >> 32); + /* 64-bit Xorshift (adviced variables by NR) */ number ^= number << 21; number ^= number >> 35; number ^= number << 4; + /* Add constant to ID */ const unsigned long long idt = id + type; + /* Nonlinear congruential generator */ const unsigned long long idpart = 3457LL * idt + 593LL * idt * ti_current + 5417LL * idt * idt; unsigned int seed = (937LL * number + 5171LL * number * number + idpart + 1109LL) % 9996361LL % seed_range; + /* Generate a random number between 0 and 1. */ return rand_r(&seed) * RAND_MAX_inv; } -- GitLab