RNG overhaul

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

6 files changed, 82 insertions, 42 deletions

diff --git a/src/diffraction.c b/src/diffraction.c
index 826aac77..92e18521 100644
--- a/src/diffraction.c
+++ b/src/diffraction.c
@@ -456,7 +456,7 @@ static struct sample *get_gaussian_spectrum(double eV_cen, double eV_step,
 }
 
 
-static int add_sase_noise(struct sample *spectrum, int nsteps)
+static int add_sase_noise(struct sample *spectrum, int nsteps, gsl_rng *rng)
 {
 	struct sample *noise;
 	int i, j;
@@ -483,7 +483,7 @@ static int add_sase_noise(struct sample *spectrum, int nsteps)
 		noise[i].weight = 0.0;
 
 		/* Gaussian noise with mean = 0, std = 1 */
-		gaussianNoise[i] = gaussian_noise(noise_mean, noise_sigma);
+		gaussianNoise[i] = gaussian_noise(rng, noise_mean, noise_sigma);
 		gaussianNoise[i+nsteps] = gaussianNoise[i];
 		gaussianNoise[i+2*nsteps] = gaussianNoise[i];
 	}
@@ -546,7 +546,7 @@ struct sample *generate_tophat(struct image *image)
 }
 
 
-struct sample *generate_SASE(struct image *image)
+struct sample *generate_SASE(struct image *image, gsl_rng *rng)
 {
 	struct sample *spectrum;
 	int i;
@@ -555,7 +555,7 @@ struct sample *generate_SASE(struct image *image)
 	const double jitter_sigma_eV = 8.0;
 
 	/* Central wavelength jitters with Gaussian distribution */
-	eV_cen = gaussian_noise(ph_lambda_to_eV(image->lambda),
+	eV_cen = gaussian_noise(rng, ph_lambda_to_eV(image->lambda),
 	                        jitter_sigma_eV);
 
 	/* Convert FWHM to standard deviation.  Note that bandwidth is taken to
@@ -571,7 +571,7 @@ struct sample *generate_SASE(struct image *image)
 	spectrum = get_gaussian_spectrum(eV_cen, eV_step, sigma, spec_size);
 
 	/* Add SASE-type noise to Gaussian spectrum */
-	add_sase_noise(spectrum, spec_size);
+	add_sase_noise(spectrum, spec_size, rng);
 
 	/* Normalise intensity (before taking restricted number of samples) */
 	double total_weight = 0.0;
diff --git a/src/diffraction.h b/src/diffraction.h
index a903346e..5b67c198 100644
--- a/src/diffraction.h
+++ b/src/diffraction.h
@@ -35,10 +35,13 @@
 #ifndef DIFFRACTION_H
 #define DIFFRACTION_H
 
+#include <gsl/gsl_rng.h>
+
 #include "image.h"
 #include "cell.h"
 #include "symmetry.h"
 
+
 typedef enum {
 	GRADIENT_MOSAIC,
 	GRADIENT_INTERPOLATE,
@@ -52,6 +55,6 @@ extern void get_diffraction(struct image *image, int na, int nb, int nc,
 
 extern struct sample *generate_tophat(struct image *image);
 
-extern struct sample *generate_SASE(struct image *image);
+extern struct sample *generate_SASE(struct image *image, gsl_rng *rng);
 
 #endif	/* DIFFRACTION_H */
diff --git a/src/get_hkl.c b/src/get_hkl.c
index 454737b0..8ba77ebf 100644
--- a/src/get_hkl.c
+++ b/src/get_hkl.c
@@ -97,6 +97,9 @@ static void poisson_reflections(RefList *list, double adu_per_photon)
 {
 	Reflection *refl;
 	RefListIterator *iter;
+	gsl_rng *rng;
+
+	rng = gsl_rng_alloc(gsl_rng_mt19937);
 
 	for ( refl = first_refl(list, &iter);
 	      refl != NULL;
@@ -106,10 +109,12 @@ static void poisson_reflections(RefList *list, double adu_per_photon)
 
 		val = get_intensity(refl);
 
-		c = adu_per_photon * poisson_noise(val/adu_per_photon);
+		c = adu_per_photon * poisson_noise(rng, val/adu_per_photon);
 		set_intensity(refl, c);
 
 	}
+
+	gsl_rng_free(rng);
 }
 
 
diff --git a/src/partial_sim.c b/src/partial_sim.c
index 7a3517a3..d9c9ae9f 100644
--- a/src/partial_sim.c
+++ b/src/partial_sim.c
@@ -3,11 +3,11 @@
  *
  * Generate partials for testing scaling
  *
- * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY,
  *                       a research centre of the Helmholtz Association.
  *
  * Authors:
- *   2011-2013 Thomas White <taw@physics.org>
+ *   2011-2014 Thomas White <taw@physics.org>
  *
  * This file is part of CrystFEL.
  *
@@ -39,6 +39,7 @@
 #include <getopt.h>
 #include <assert.h>
 #include <pthread.h>
+#include <gsl/gsl_rng.h>
 
 #include "image.h"
 #include "utils.h"
@@ -54,7 +55,7 @@
 #define NBINS 50
 
 
-static void mess_up_cell(Crystal *cr, double cnoise)
+static void mess_up_cell(Crystal *cr, double cnoise, gsl_rng *rng)
 {
 	double ax, ay, az;
 	double bx, by, bz;
@@ -65,15 +66,15 @@ static void mess_up_cell(Crystal *cr, double cnoise)
 	//cell_print(cell);
 
 	cell_get_reciprocal(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
-	ax = flat_noise(ax, cnoise*fabs(ax)/100.0);
-	ay = flat_noise(ay, cnoise*fabs(ay)/100.0);
-	az = flat_noise(az, cnoise*fabs(az)/100.0);
-	bx = flat_noise(bx, cnoise*fabs(bx)/100.0);
-	by = flat_noise(by, cnoise*fabs(by)/100.0);
-	bz = flat_noise(bz, cnoise*fabs(bz)/100.0);
-	cx = flat_noise(cx, cnoise*fabs(cx)/100.0);
-	cy = flat_noise(cy, cnoise*fabs(cy)/100.0);
-	cz = flat_noise(cz, cnoise*fabs(cz)/100.0);
+	ax = flat_noise(rng, ax, cnoise*fabs(ax)/100.0);
+	ay = flat_noise(rng, ay, cnoise*fabs(ay)/100.0);
+	az = flat_noise(rng, az, cnoise*fabs(az)/100.0);
+	bx = flat_noise(rng, bx, cnoise*fabs(bx)/100.0);
+	by = flat_noise(rng, by, cnoise*fabs(by)/100.0);
+	bz = flat_noise(rng, bz, cnoise*fabs(bz)/100.0);
+	cx = flat_noise(rng, cx, cnoise*fabs(cx)/100.0);
+	cy = flat_noise(rng, cy, cnoise*fabs(cy)/100.0);
+	cz = flat_noise(rng, cz, cnoise*fabs(cz)/100.0);
 	cell_set_reciprocal(cell, ax, ay, az, bx, by, bz, cx, cy, cz);
 
 	//STATUS("Changed:\n");
@@ -89,7 +90,7 @@ static void calculate_partials(Crystal *cr,
                                pthread_rwlock_t *full_lock,
                                unsigned long int *n_ref, double *p_hist,
                                double *p_max, double max_q, double full_stddev,
-                               double noise_stddev)
+                               double noise_stddev, gsl_rng *rng)
 {
 	Reflection *refl;
 	RefListIterator *iter;
@@ -124,7 +125,7 @@ static void calculate_partials(Crystal *cr,
 				rfull = find_refl(full, h, k, l);
 				if ( rfull == NULL ) {
 					rfull = add_refl(full, h, k, l);
-					If = fabs(gaussian_noise(0.0,
+					If = fabs(gaussian_noise(rng, 0.0,
 					                         full_stddev));
 					set_intensity(rfull, If);
 					set_redundancy(rfull, 1);
@@ -160,7 +161,7 @@ static void calculate_partials(Crystal *cr,
 			       res, bin,  p);
 		}
 
-		Ip = gaussian_noise(Ip, noise_stddev);
+		Ip = gaussian_noise(rng, Ip, noise_stddev);
 
 		set_intensity(refl, Ip);
 		set_esd_intensity(refl, noise_stddev);
@@ -226,6 +227,7 @@ struct queue_args
 	double p_max[NBINS];
 
 	Stream *stream;
+	gsl_rng **rngs;
 };
 
 
@@ -280,14 +282,15 @@ static void run_job(void *vwargs, int cookie)
 	crystal_set_image(cr, &wargs->image);
 
 	do {
-		osf = gaussian_noise(1.0, qargs->osf_stddev);
+		osf = gaussian_noise(qargs->rngs[cookie], 1.0,
+		                     qargs->osf_stddev);
 	} while ( osf <= 0.0 );
 	crystal_set_osf(cr, osf);
 	crystal_set_mosaicity(cr, 0.0);
 	crystal_set_profile_radius(cr, wargs->image.beam->profile_radius);
 
 	/* Set up a random orientation */
-	orientation = random_quaternion();
+	orientation = random_quaternion(qargs->rngs[cookie]);
 	crystal_set_cell(cr, cell_rotate(qargs->cell, orientation));
 
 	snprintf(wargs->image.filename, 255, "dummy.h5");
@@ -305,10 +308,10 @@ static void run_job(void *vwargs, int cookie)
 	                   &qargs->full_lock,
 	                   wargs->n_ref, wargs->p_hist, wargs->p_max,
 	                   qargs->max_q, qargs->full_stddev,
-	                   qargs->noise_stddev);
+	                   qargs->noise_stddev, qargs->rngs[cookie]);
 
 	/* Give a slightly incorrect cell in the stream */
-	mess_up_cell(cr, qargs->cnoise);
+	mess_up_cell(cr, qargs->cnoise, qargs->rngs[cookie]);
 
 	image_add_crystal(&wargs->image, cr);
 }
@@ -367,6 +370,7 @@ int main(int argc, char *argv[])
 	double osf_stddev = 2.0;
 	double full_stddev = 1000.0;
 	double noise_stddev = 20.0;
+	gsl_rng *rng_for_seeds;
 
 	/* Long options */
 	const struct option longopts[] = {
@@ -623,6 +627,18 @@ int main(int argc, char *argv[])
 	qargs.noise_stddev = noise_stddev;
 	qargs.max_q = largest_q(&image);
 
+	qargs.rngs = malloc(n_threads * sizeof(gsl_rng *));
+	if ( qargs.rngs == NULL ) {
+		ERROR("Failed to allocate RNGs\n");
+		return 1;
+	}
+	rng_for_seeds = gsl_rng_alloc(gsl_rng_mt19937);
+	for ( i=0; i<n_threads; i++ ) {
+		qargs.rngs[i] = gsl_rng_alloc(gsl_rng_mt19937);
+		gsl_rng_set(qargs.rngs[i], gsl_rng_get(rng_for_seeds));
+	}
+	gsl_rng_free(rng_for_seeds);
+
 	for ( i=0; i<NBINS; i++ ) {
 		qargs.n_ref[i] = 0;
 		qargs.p_hist[i] = 0.0;
@@ -683,6 +699,9 @@ int main(int argc, char *argv[])
 
 	}
 
+	for ( i=0; i<n_threads; i++ ) {
+		gsl_rng_free(qargs.rngs[i]);
+	}
 	pthread_rwlock_destroy(&qargs.full_lock);
 	close_stream(stream);
 	cell_free(cell);
diff --git a/src/pattern_sim.c b/src/pattern_sim.c
index 81a94c37..5ea328bf 100644
--- a/src/pattern_sim.c
+++ b/src/pattern_sim.c
@@ -3,11 +3,13 @@
  *
  * Simulate diffraction patterns from small crystals
  *
- * 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>
+ *   2013-2014 Chun Hong Yoon <chun.hong.yoon@desy.de>
+ *   2013      Alexandra Tolstikova
  *
  * This file is part of CrystFEL.
  *
@@ -278,6 +280,7 @@ int main(int argc, char *argv[])
 	char *sym_str = NULL;
 	SymOpList *sym;
 	int nsamples = 3;
+	gsl_rng *rng;
 
 	/* Long options */
 	const struct option longopts[] = {
@@ -407,15 +410,6 @@ int main(int argc, char *argv[])
 
 	}
 
-	if ( config_random ) {
-		FILE *fh;
-		unsigned int seed;
-		fh = fopen("/dev/urandom", "r");
-		fread(&seed, sizeof(seed), 1, fh);
-		fclose(fh);
-		srand(seed);
-	}
-
 	if ( random_size == 1 ) {
 		ERROR("You must specify both --min-size and --max-size.\n");
 		return 1;
@@ -565,6 +559,16 @@ int main(int argc, char *argv[])
 	image.features = NULL;
 	image.flags = NULL;
 
+	rng = gsl_rng_alloc(gsl_rng_mt19937);
+	if ( config_random ) {
+		FILE *fh;
+		unsigned long int seed;
+		fh = fopen("/dev/urandom", "r");
+		fread(&seed, sizeof(seed), 1, fh);
+		fclose(fh);
+		gsl_rng_set(rng, seed);
+	}
+
 	powder = calloc(image.width*image.height, sizeof(*powder));
 
 	/* Splurge a few useful numbers */
@@ -597,7 +601,7 @@ int main(int argc, char *argv[])
 
 		/* Read quaternion from stdin */
 		if ( config_randomquat ) {
-			orientation = random_quaternion();
+			orientation = random_quaternion(rng);
 		} else {
 			orientation = read_quaternion();
 		}
@@ -620,7 +624,7 @@ int main(int argc, char *argv[])
 			break;
 
 			case SPECTRUM_SASE :
-			image.spectrum = generate_SASE(&image);
+			image.spectrum = generate_SASE(&image, rng);
 			break;
 
 		}
@@ -650,7 +654,7 @@ int main(int argc, char *argv[])
 			goto skip;
 		}
 
-		record_image(&image, !config_nonoise);
+		record_image(&image, !config_nonoise, rng);
 
 		if ( powder_fn != NULL ) {
 
@@ -723,5 +727,7 @@ skip:
 	free(sym_str);
 	free_symoplist(sym);
 
+	gsl_rng_free(rng);
+
 	return 0;
 }
diff --git a/src/scaling-report.c b/src/scaling-report.c
index 89b2f3f2..85ed0ad2 100644
--- a/src/scaling-report.c
+++ b/src/scaling-report.c
@@ -196,6 +196,7 @@ static void partiality_graph(cairo_t *cr, Crystal **crystals, int n,
 	double pcalcmin[nbins];
 	double pcalcmax[nbins];
 	int num_nondud;
+	gsl_rng *rng;
 
 	show_text_simple(cr, "Observed partiality", -20.0, g_height/2.0,
 	                      NULL, -M_PI_2, J_CENTER);
@@ -223,6 +224,10 @@ static void partiality_graph(cairo_t *cr, Crystal **crystals, int n,
 		num_nondud++;
 	}
 
+	/* The reflections chosen for the graph will be the same every time
+	 * (given the same sequence of input reflections, scalabilities etc) */
+	rng = gsl_rng_alloc(gsl_rng_mt19937);
+
 	cairo_set_source_rgb(cr, 0.0, 0.7, 0.0);
 	prob = 1.0 / num_nondud;
 	for ( i=0; i<n; i++ ) {
@@ -279,13 +284,15 @@ static void partiality_graph(cairo_t *cr, Crystal **crystals, int n,
 
 			bin = nbins * pcalc;
 
-			if ( random_flat(1.0) < prob ) {
+			if ( random_flat(rng, 1.0) < prob ) {
 				plot_point(cr, g_width, g_height, pcalc, pobs);
 			}
 		}
 
 	}
 
+	gsl_rng_free(rng);
+
 	cairo_new_path(cr);
 	cairo_rectangle(cr, 0.0, 0.0, g_width, g_height);
 	cairo_clip(cr);