RNG overhaul

Previously, we were using random(), which is really really bad.

5 files changed, 44 insertions, 34 deletions

diff --git a/libcrystfel/src/detector.c b/libcrystfel/src/detector.c
index 5cd0ce0f..174abefe 100644
--- a/libcrystfel/src/detector.c
+++ b/libcrystfel/src/detector.c
@@ -311,7 +311,7 @@ double get_tt(struct image *image, double fs, double ss, int *err)
 }
 
 
-void record_image(struct image *image, int do_poisson)
+void record_image(struct image *image, int do_poisson, gsl_rng *rng)
 {
 	int x, y;
 	double total_energy, energy_density;
@@ -371,7 +371,7 @@ void record_image(struct image *image, int do_poisson)
 		sa = proj_area / (dsq + Lsq);
 
 		if ( do_poisson ) {
-			counts = poisson_noise(intensity * ph_per_e * sa);
+			counts = poisson_noise(rng, intensity * ph_per_e * sa);
 		} else {
 			cf = intensity * ph_per_e * sa;
 			counts = cf;
diff --git a/libcrystfel/src/detector.h b/libcrystfel/src/detector.h
index cfe0808a..eaa5618e 100644
--- a/libcrystfel/src/detector.h
+++ b/libcrystfel/src/detector.h
@@ -153,7 +153,7 @@ extern double get_tt(struct image *image, double xs, double ys, int *err);
 
 extern int in_bad_region(struct detector *det, double fs, double ss);
 
-extern void record_image(struct image *image, int do_poisson);
+extern void record_image(struct image *image, int do_poisson, gsl_rng *rng);
 
 extern struct panel *find_panel(struct detector *det, double fs, double ss);
 
diff --git a/libcrystfel/src/reax.c b/libcrystfel/src/reax.c
index 2ab4397a..110d8701 100644
--- a/libcrystfel/src/reax.c
+++ b/libcrystfel/src/reax.c
@@ -793,6 +793,13 @@ static int check_twinning(UnitCell *c1, UnitCell *c2, int verbose)
 	double bx, by, bz;
 	double cx, cy, cz;
 
+	gsl_rng *rng;
+
+	/* This is a rubbish RNG, but it serves for this purpose: nothing more
+	 * than "I couldn't be bothered to think of n_trials sets of random
+	 * indices". */
+	rng = gsl_rng_alloc(gsl_rng_taus2);
+
 	cell_get_reciprocal(c1, &asx, &asy, &asz,
 	                        &bsx, &bsy, &bsz,
 	                        &csx, &csy, &csz);
@@ -807,9 +814,9 @@ static int check_twinning(UnitCell *c1, UnitCell *c2, int verbose)
 		double dev;
 		signed int h2i, k2i, l2i;
 
-		h = flat_noise(0, 10);
-		k = flat_noise(0, 10);
-		l = flat_noise(0, 10);
+		h = gsl_rng_uniform_int(rng, 10);
+		k = gsl_rng_uniform_int(rng, 10);
+		l = gsl_rng_uniform_int(rng, 10);
 
 		/* Position of this (randomly selected)
 		 * reciprocal lattice point */
@@ -844,6 +851,8 @@ static int check_twinning(UnitCell *c1, UnitCell *c2, int verbose)
 		STATUS("%i duplicates.\n", n_dup);
 	}
 
+	gsl_rng_free(rng);
+
 	if ( n_dup > 10 ) return 1;
 	return 0;
 }
diff --git a/libcrystfel/src/utils.c b/libcrystfel/src/utils.c
index 2f5107e3..4c85bf9d 100644
--- a/libcrystfel/src/utils.c
+++ b/libcrystfel/src/utils.c
@@ -3,11 +3,11 @@
  *
  * Utility stuff
  *
- * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY,
- *                  a research centre of the Helmholtz Association.
+ * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY,
+ *                       a research centre of the Helmholtz Association.
  *
  * Authors:
- *   2009-2012 Thomas White <taw@physics.org>
+ *   2009-2014 Thomas White <taw@physics.org>
  *
  * This file is part of CrystFEL.
  *
@@ -165,26 +165,26 @@ void progress_bar(int val, int total, const char *text)
 }
 
 
-double random_flat(double max)
+double random_flat(gsl_rng *rng, double max)
 {
-	return max * (double)random()/RAND_MAX;
+	return max * gsl_rng_uniform(rng);
 }
 
 
-double flat_noise(double expected, double width)
+double flat_noise(gsl_rng *rng, double expected, double width)
 {
-	double noise = random_flat(2.0*width);
+	double noise = random_flat(rng, 2.0*width);
 	return expected+noise-width;
 }
 
 
-double gaussian_noise(double expected, double stddev)
+double gaussian_noise(gsl_rng *rng, double expected, double stddev)
 {
 	double x1, x2, noise;
 
 	/* A uniformly distributed random number between 0 and 1 */
-	x1 = ((double)random()/RAND_MAX);
-	x2 = ((double)random()/RAND_MAX);
+	x1 = gsl_rng_uniform(rng);
+	x2 = gsl_rng_uniform(rng);
 
 	noise = sqrt(-2.0*log(x1)) * cos(2.0*M_PI*x2);
 
@@ -192,14 +192,14 @@ double gaussian_noise(double expected, double stddev)
 }
 
 
-static int fake_poisson_noise(double expected)
+static int fake_poisson_noise(gsl_rng *rng, double expected)
 {
-	double rf = gaussian_noise(expected, sqrt(expected));
+	double rf = gaussian_noise(rng, expected, sqrt(expected));
 	return (int)rf;
 }
 
 
-int poisson_noise(double expected)
+int poisson_noise(gsl_rng *rng, double expected)
 {
 	double L;
 	int k = 0;
@@ -207,7 +207,7 @@ int poisson_noise(double expected)
 
 	/* For large values of the mean, we get big problems with arithmetic.
 	 * In such cases, fall back on a Gaussian with the right variance. */
-	if ( expected > 100.0 ) return fake_poisson_noise(expected);
+	if ( expected > 100.0 ) return fake_poisson_noise(rng, expected);
 
 	L = exp(-expected);
 
@@ -216,7 +216,7 @@ int poisson_noise(double expected)
 		double r;
 
 		k++;
-		r = (double)random()/(double)RAND_MAX;
+		r = gsl_rng_uniform(rng);
 		p *= r;
 
 	} while ( p > L );
@@ -282,14 +282,14 @@ struct quaternion normalise_quaternion(struct quaternion q)
  *
  * Returns: a randomly generated, normalised, quaternion.
  **/
-struct quaternion random_quaternion()
+struct quaternion random_quaternion(gsl_rng *rng)
 {
 	struct quaternion q;
 
-	q.w = 2.0*(double)random()/RAND_MAX - 1.0;
-	q.x = 2.0*(double)random()/RAND_MAX - 1.0;
-	q.y = 2.0*(double)random()/RAND_MAX - 1.0;
-	q.z = 2.0*(double)random()/RAND_MAX - 1.0;
+	q.w = 2.0*gsl_rng_uniform(rng) - 1.0;
+	q.x = 2.0*gsl_rng_uniform(rng) - 1.0;
+	q.y = 2.0*gsl_rng_uniform(rng) - 1.0;
+	q.z = 2.0*gsl_rng_uniform(rng) - 1.0;
 	q = normalise_quaternion(q);
 
 	return q;
diff --git a/libcrystfel/src/utils.h b/libcrystfel/src/utils.h
index f4477717..81270916 100644
--- a/libcrystfel/src/utils.h
+++ b/libcrystfel/src/utils.h
@@ -3,11 +3,11 @@
  *
  * Utility stuff
  *
- * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY,
- *                  a research centre of the Helmholtz Association.
+ * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY,
+ *                       a research centre of the Helmholtz Association.
  *
  * Authors:
- *   2009-2012 Thomas White <taw@physics.org>
+ *   2009-2014 Thomas White <taw@physics.org>
  *
  * This file is part of CrystFEL.
  *
@@ -40,6 +40,7 @@
 #include <pthread.h>
 #include <gsl/gsl_matrix.h>
 #include <gsl/gsl_vector.h>
+#include <gsl/gsl_rng.h>
 
 
 #include "thread-pool.h"
@@ -89,7 +90,7 @@ struct quaternion {
 
 extern struct quaternion normalise_quaternion(struct quaternion q);
 extern double quaternion_modulus(struct quaternion q);
-extern struct quaternion random_quaternion(void);
+extern struct quaternion random_quaternion(gsl_rng *rng);
 extern int quaternion_valid(struct quaternion q);
 extern struct rvec quat_rot(struct rvec q, struct quaternion z);
 
@@ -114,10 +115,10 @@ extern int assplode(const char *a, const char *delims, char ***pbits,
                     AssplodeFlag flags);
 
 extern void progress_bar(int val, int total, const char *text);
-extern double random_flat(double max);
-extern double flat_noise(double expected, double width);
-extern double gaussian_noise(double expected, double stddev);
-extern int poisson_noise(double expected);
+extern double random_flat(gsl_rng *rng, double max);
+extern double flat_noise(gsl_rng *rng, double expected, double width);
+extern double gaussian_noise(gsl_rng *rng, double expected, double stddev);
+extern int poisson_noise(gsl_rng *rng, double expected);
 
 /* Keep these ones inline, to avoid function call overhead */
 static inline struct quaternion invalid_quaternion(void)