diff --git a/tests/testRandomCone.c b/tests/testRandomCone.c index 08c953c00acf2adb6997be010b44ed045229b5ab..a7e4568de213770cbd6ce4a7afd4c41cb6244100 100644 --- a/tests/testRandomCone.c +++ b/tests/testRandomCone.c @@ -18,6 +18,9 @@ ******************************************************************************/ #include <config.h> +/* System includes. */ +#include <fenv.h> + /* Local headers. */ #include "swift.h" @@ -47,10 +50,9 @@ const int N_cube = 5; * the uniformity of the distribution. * @param tolerance The tolerance of each bin relative to the expected value. */ -void test_cone(int64_t id_bh, const integertime_t ti_current, - const enum random_number_type type, double opening_angle, - float unit_vector[3], const int N_test, const int N_bins, - const double tolerance) { +float test_cone(int64_t id_bh, const integertime_t ti_current, + const enum random_number_type type, double opening_angle, + float unit_vector[3]) { /* Compute cosine that corresponds to the maximum opening angle */ const double cos_theta_max = cos(opening_angle); @@ -58,56 +60,41 @@ void test_cone(int64_t id_bh, const integertime_t ti_current, /* Initialize an array that will hold a random vector every step */ float rand_vector[3]; - /* Initialize an array that will hold the binned number of drawn cosines, - i.e. this is the probability density function that we wish to test. */ - double binned_cosines[N_bins]; - for (int j = 0; j < N_bins; ++j) { - binned_cosines[j] = 0.; + /* Generate a random unit vector within a cone around unit_vector */ + random_direction_in_cone(id_bh, ti_current, type, opening_angle, unit_vector, + rand_vector); + + /* Check that this vector is actually within the cone we want */ + const double cos_rand_unit = rand_vector[0] * unit_vector[0] + + rand_vector[1] * unit_vector[1] + + rand_vector[2] * unit_vector[2]; + if (cos_rand_unit < 0.99999 * cos_theta_max) { + printf("Cos_opening_angle is: %f, Random cos is: %f\n", cos_theta_max, + cos_rand_unit); + error("Generated random unit vector is outside cone."); } - for (int k = 0; k < N_test; ++k) { + return cos_rand_unit; +} - /* Generate a random unit vector within a cone around unit_vector */ - random_direction_in_cone(id_bh, ti_current, type, opening_angle, - unit_vector, rand_vector); - - /* Check that this vector is actually within the cone we want */ - const double cos_rand_unit = rand_vector[0] * unit_vector[0] + - rand_vector[1] * unit_vector[1] + - rand_vector[2] * unit_vector[2]; - if (cos_rand_unit < 0.99999 * cos_theta_max) { - printf("Cos_opening_angle is: %f, Random cos is: %f\n", cos_theta_max, - cos_rand_unit); - error("Generated random unit vector is outside cone."); - } +int main(int argc, char *argv[]) { - /* Add the unit vector to the probability density function array. The solid - * angle subtended by some angle theta grows as (1-cos(theta)). Furthermore, - * we are limited to the spherical cap defined by the angles [0, theta_max]. - * Therefore the variable which we expect to be uniformly distributed is (1 - * - cos(theta)) / (1 - cos(theta_max)). */ - double uniform_variable = (1. - cos_rand_unit) / (1 - cos_theta_max); - for (int j = 0; j < N_bins; ++j) { - if ((uniform_variable > (double)j / (double)N_bins) && - (uniform_variable < (double)(j + 1) / (double)N_bins)) { - binned_cosines[j] = binned_cosines[j] + 1. / (double)N_test; - } - } - } + /* Initialize CPU frequency, this also starts time. */ + unsigned long long cpufreq = 0; + clocks_set_cpufreq(cpufreq); - /* Check whether the binned quantity is really uniformly distributed. If it - * is, the density (value) of each bin should be 1/N_bin. */ - for (int j = 0; j < N_bins; ++j) { - if ((binned_cosines[j] < (1. - tolerance) / (double)N_bins) || - (binned_cosines[j] > (1. + tolerance) / (double)N_bins)) { - error( - "Generated distribution of random unit vectors within a cone exceeds " - "the limit imposed by the tolerance."); - } - } -} +/* Choke on FPEs */ +#ifdef HAVE_FE_ENABLE_EXCEPT + feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW); +#endif -int main(int argc, char *argv[]) { + /* Get some randomness going */ + const int seed = time(NULL); + message("Seed = %d", seed); + srand(seed); + + /* Log the swift random seed */ + message("SWIFT random seed = %d", SWIFT_RANDOM_SEED_XOR); /* Test the random-vector-in-cone function, for different values of opening * angle from 0 to pi/2 (in radians). For each of these opening angles we draw @@ -132,15 +119,16 @@ int main(int argc, char *argv[]) { const double cos_unit = random_unit_interval(id_bh, ti_current, random_number_BH_kick); const double sin_unit = sqrtf(max(0., (1. - cos_unit) * (1. + cos_unit))); - const double phi_unit = random_unit_interval(id_bh * id_bh, ti_current, - random_number_BH_kick); + const double phi_unit = + (2. * M_PI) * random_unit_interval(id_bh * id_bh, ti_current, + random_number_BH_kick); unit_vector[0] = sin_unit * cos(phi_unit); unit_vector[1] = sin_unit * sin(phi_unit); unit_vector[2] = cos_unit; /* Do the test. */ test_cone(id_bh, ti_current, random_number_BH_kick, opening_angle, - unit_vector, 100000, 10, 0.1); + unit_vector); } } @@ -148,26 +136,72 @@ int main(int argc, char *argv[]) { * bins, but for just one opening angle and one randomly generated cone */ const double opening_angle_0 = 0.2; - /* Generate an id for the bh and a time. We do this for every opening - * angle and every cone. */ - const long long id_bh_0 = rand() * (1LL << 31) + rand(); - const integertime_t ti_current_0 = rand() * (1LL << 31) + rand(); + /* Compute cosine that corresponds to the maximum opening angle */ + const double cos_theta_max = cos(opening_angle_0); /* Generate a random unit vector that defines a cone, along with the * opening angle. */ + const long long id_bh_0 = rand() * (1LL << 31) + rand(); + const integertime_t ti_current_0 = rand() * (1LL << 31) + rand(); + float unit_vector_0[3]; const double cos_unit = random_unit_interval(id_bh_0, ti_current_0, random_number_BH_kick); const double sin_unit = sqrtf(max(0., (1. - cos_unit) * (1. + cos_unit))); - const double phi_unit = random_unit_interval(id_bh_0 * id_bh_0, ti_current_0, - random_number_BH_kick); + const double phi_unit = + (2. * M_PI) * random_unit_interval(id_bh_0 * id_bh_0, ti_current_0, + random_number_BH_kick); unit_vector_0[0] = sin_unit * cos(phi_unit); unit_vector_0[1] = sin_unit * sin(phi_unit); unit_vector_0[2] = cos_unit; - /* Do the test. */ - test_cone(id_bh_0, ti_current_0, random_number_BH_kick, opening_angle_0, - unit_vector_0, 100000000, 100, 0.01); + /* Some parameters to test the uniformity of drawn vectors */ + int N_test = 10000000; + int N_bins = 100; + float tolerance = 0.05; + + /* Initialize an array that will hold the binned number of drawn cosines, + i.e. this is the probability density function that we wish to test. */ + double binned_cosines[N_bins]; + for (int j = 0; j < N_bins; ++j) { + binned_cosines[j] = 0.; + } + + /* Draw N_test vectors and bin them to test uniformity */ + for (int k = 0; k < N_test; ++k) { + + const long long id_bh = rand() * (1LL << 31) + rand(); + const integertime_t ti_current = rand() * (1LL << 31) + rand(); + + /* Do the test, with a newly generated BH id and time */ + const float cos_rand_unit = + test_cone(id_bh, ti_current, random_number_BH_kick, opening_angle_0, + unit_vector_0); + + /* Add the unit vector to the probability density function array. The solid + * angle subtended by some angle theta grows as (1-cos(theta)). Furthermore, + * we are limited to the spherical cap defined by the angles [0, theta_max]. + * Therefore the variable which we expect to be uniformly distributed is (1 + * - cos(theta)) / (1 - cos(theta_max)). */ + double uniform_variable = (1. - cos_rand_unit) / (1 - cos_theta_max); + for (int j = 0; j < N_bins; ++j) { + if ((uniform_variable > (double)j / (double)N_bins) && + (uniform_variable < (double)(j + 1) / (double)N_bins)) { + binned_cosines[j] = binned_cosines[j] + 1. / (double)N_test; + } + } + } + + /* Check whether the binned quantity is really uniformly distributed. If it + * is, the density (value) of each bin should be 1/N_bin. */ + for (int j = 0; j < N_bins; ++j) { + if ((binned_cosines[j] < (1. - tolerance) / (double)N_bins) || + (binned_cosines[j] > (1. + tolerance) / (double)N_bins)) { + error( + "Generated distribution of random unit vectors within a cone exceeds " + "the limit imposed by the tolerance."); + } + } /* We now repeat the same process, but we do not generate random unit vectors * to define the cones. Instead, we sample unit vectors along the grid @@ -204,7 +238,7 @@ int main(int argc, char *argv[]) { /* Do the test. */ test_cone(id_bh, ti_current, random_number_BH_kick, opening_angle, - unit_vector, 100000, 10, 0.1); + unit_vector); } } }