Fix random cone test

Fixes #870 (closed).

In a recent MR, the function random_direction_in_cone that generates a random unit vector within a given cone was modified. It no longer takes both the ID of a BH and a gas particle as inputs, but only the former. However, this means that the testRandomCone.h test was less random than we wanted, and was thus failing.

I have modified that test in a very simple way. It still does exactly the same thing, but we now generate a random id id_random and pass id_BH + id_random instead of id_BH to the random_direction_in_cone function. The test again passes with this change.

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);
         }
       }
     }