From 90ee3c269580104f2d16d28aeaa565063f6fc1f2 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Fri, 17 Jan 2014 16:52:57 +0100 Subject: RNG overhaul Previously, we were using random(), which is really really bad. --- libcrystfel/src/detector.c | 4 ++-- libcrystfel/src/detector.h | 2 +- libcrystfel/src/reax.c | 15 ++++++++++++--- libcrystfel/src/utils.c | 40 ++++++++++++++++++++-------------------- libcrystfel/src/utils.h | 17 +++++++++-------- 5 files changed, 44 insertions(+), 34 deletions(-) (limited to 'libcrystfel') 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 + * 2009-2014 Thomas White * * 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 + * 2009-2014 Thomas White * * This file is part of CrystFEL. * @@ -40,6 +40,7 @@ #include #include #include +#include #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) -- cgit v1.2.3