Skip to content
Snippets Groups Projects
Commit d5480423 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Merge branch 'fix_random_cone_test' into 'master' 

Fix random cone test

Closes #870

See merge request !1793
parents f5da7045 e82d2ee9
Branches
Tags
1 merge request!1793Fix random cone test
Hide whitespace changes
Inline Side-by-side
Showing
with 95 additions and 61 deletions
tests/testRandomCone.c
+ 95
− 61
View file @ d5480423
......@@ -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);
}
}
}
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment