aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--libcrystfel/src/detector.c4
-rw-r--r--libcrystfel/src/detector.h2
-rw-r--r--libcrystfel/src/reax.c15
-rw-r--r--libcrystfel/src/utils.c40
-rw-r--r--libcrystfel/src/utils.h17
-rw-r--r--src/diffraction.c10
-rw-r--r--src/diffraction.h5
-rw-r--r--src/get_hkl.c7
-rw-r--r--src/partial_sim.c57
-rw-r--r--src/pattern_sim.c36
-rw-r--r--src/scaling-report.c9
-rw-r--r--tests/cell_check.c12
-rw-r--r--tests/centering_check.c93
-rw-r--r--tests/gpu_sim_check.c13
-rw-r--r--tests/integration_check.c11
-rw-r--r--tests/pr_l_gradient_check.c11
-rw-r--r--tests/pr_p_gradient_check.c11
-rw-r--r--tests/pr_pl_gradient_check.c11
-rw-r--r--tests/prof2d_check.c28
-rw-r--r--tests/ring_check.c31
-rw-r--r--tests/transformation_check.c12
21 files changed, 277 insertions, 158 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)
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);
diff --git a/tests/cell_check.c b/tests/cell_check.c
index f62cd7be..009b4c1f 100644
--- a/tests/cell_check.c
+++ b/tests/cell_check.c
@@ -3,7 +3,11 @@
*
* Check that unit cells work correctly
*
- * Copyright © 2012 Thomas White <taw@physics.org>
+ * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY,
+ * a research centre of the Helmholtz Association.
+ *
+ * Authors:
+ * 2012-2014 Thomas White <taw@physics.org>
*
* This file is part of CrystFEL.
*
@@ -49,6 +53,9 @@ int main(int argc, char *argv[])
double ax, ay, az;
double bx, by, bz;
double cx, cy, cz;
+ gsl_rng *rng;
+
+ rng = gsl_rng_alloc(gsl_rng_mt19937);
cell = cell_new_from_parameters(27.155e-9, 28.155e-9, 10.987e-9,
deg2rad(90.0),
@@ -56,7 +63,7 @@ int main(int argc, char *argv[])
deg2rad(120.0));
if ( cell == NULL ) return 1;
- orientation = random_quaternion();
+ orientation = random_quaternion(rng);
cell = cell_rotate(cell, orientation);
cell_get_reciprocal(cell, &asx, &asy, &asz,
@@ -112,6 +119,7 @@ int main(int argc, char *argv[])
cx, cy, cz);
cell_print(cell);
+ gsl_rng_free(rng);
return fail;
}
diff --git a/tests/centering_check.c b/tests/centering_check.c
index f4073072..887311dd 100644
--- a/tests/centering_check.c
+++ b/tests/centering_check.c
@@ -3,7 +3,11 @@
*
* Check that centering of cells works
*
- * Copyright © 2012 Thomas White <taw@physics.org>
+ * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY,
+ * a research centre of the Helmholtz Association.
+ *
+ * Authors:
+ * 2012-2014 Thomas White <taw@physics.org>
*
* This file is part of CrystFEL.
*
@@ -53,7 +57,7 @@ static int check_cell(UnitCell *cell, const char *text)
static int check_centering(double a, double b, double c,
double al, double be, double ga,
- LatticeType latt, char cen, char ua)
+ LatticeType latt, char cen, char ua, gsl_rng *rng)
{
UnitCell *cell, *cref;
UnitCell *n;
@@ -77,7 +81,7 @@ static int check_centering(double a, double b, double c,
cell_set_centering(cref, cen);
cell_set_unique_axis(cref, ua);
- cell = cell_rotate(cref, random_quaternion());
+ cell = cell_rotate(cref, random_quaternion(rng));
if ( cell == NULL ) return 1;
cell_free(cref);
@@ -110,9 +114,9 @@ static int check_centering(double a, double b, double c,
do {
- h = flat_noise(0, 30);
- k = flat_noise(0, 30);
- l = flat_noise(0, 30);
+ h = gsl_rng_uniform_int(rng, 30);
+ k = gsl_rng_uniform_int(rng, 30);
+ l = gsl_rng_uniform_int(rng, 30);
} while ( forbidden_reflection(cell, h, k, l) );
@@ -156,9 +160,9 @@ static int check_centering(double a, double b, double c,
int f = 0;
long int ih, ik, il;
- h = flat_noise(0, 30);
- k = flat_noise(0, 30);
- l = flat_noise(0, 30);
+ h = gsl_rng_uniform_int(rng, 30);
+ k = gsl_rng_uniform_int(rng, 30);
+ l = gsl_rng_uniform_int(rng, 30);
x = h*asx + k*bsx + l*csx;
y = h*asy + k*bsy + l*csy;
@@ -199,110 +203,115 @@ static int check_centering(double a, double b, double c,
int main(int argc, char *argv[])
{
int fail = 0;
+ gsl_rng *rng;
+
+ rng = gsl_rng_alloc(gsl_rng_mt19937);
/* Triclinic P */
fail += check_centering(50e-10, 55e-10, 70e-10, 67.0, 70.0, 77.0,
- L_TRICLINIC, 'P', '*');
+ L_TRICLINIC, 'P', '*', rng);
/* Monoclinic P */
fail += check_centering(10e-10, 20e-10, 30e-10, 100.0, 90.0, 90.0,
- L_MONOCLINIC, 'P', 'a');
+ L_MONOCLINIC, 'P', 'a', rng);
fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 100.0, 90.0,
- L_MONOCLINIC, 'P', 'b');
+ L_MONOCLINIC, 'P', 'b', rng);
fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 90.0, 100.0,
- L_MONOCLINIC, 'P', 'c');
+ L_MONOCLINIC, 'P', 'c', rng);
/* Monoclinic "C"-centered, unique axis a, three cell choices */
fail += check_centering(10e-10, 20e-10, 30e-10, 100.0, 90.0, 90.0,
- L_MONOCLINIC, 'B', 'a');
+ L_MONOCLINIC, 'B', 'a', rng);
fail += check_centering(10e-10, 20e-10, 30e-10, 100.0, 90.0, 90.0,
- L_MONOCLINIC, 'C', 'a');
+ L_MONOCLINIC, 'C', 'a', rng);
fail += check_centering(10e-10, 20e-10, 30e-10, 100.0, 90.0, 90.0,
- L_MONOCLINIC, 'I', 'a');
+ L_MONOCLINIC, 'I', 'a', rng);
/* Monoclinic "C"-centered, unique axis b, three cell choices */
fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 100.0, 90.0,
- L_MONOCLINIC, 'C', 'b');
+ L_MONOCLINIC, 'C', 'b', rng);
fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 100.0, 90.0,
- L_MONOCLINIC, 'A', 'b');
+ L_MONOCLINIC, 'A', 'b', rng);
fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 100.0, 90.0,
- L_MONOCLINIC, 'I', 'b');
+ L_MONOCLINIC, 'I', 'b', rng);
/* Monoclinic "C"-centered, unique axis c, three cell choices */
fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 90.0, 100.0,
- L_MONOCLINIC, 'A', 'c');
+ L_MONOCLINIC, 'A', 'c', rng);
fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 90.0, 100.0,
- L_MONOCLINIC, 'B', 'c');
+ L_MONOCLINIC, 'B', 'c', rng);
fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 90.0, 100.0,
- L_MONOCLINIC, 'I', 'c');
+ L_MONOCLINIC, 'I', 'c', rng);
/* Orthorhombic P */
fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 90.0, 90.0,
- L_ORTHORHOMBIC, 'P', '*');
+ L_ORTHORHOMBIC, 'P', '*', rng);
/* Orthorhombic A */
fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 90.0, 90.0,
- L_ORTHORHOMBIC, 'A', '*');
+ L_ORTHORHOMBIC, 'A', '*', rng);
/* Orthorhombic B */
fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 90.0, 90.0,
- L_ORTHORHOMBIC, 'B', '*');
+ L_ORTHORHOMBIC, 'B', '*', rng);
/* Orthorhombic C */
fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 90.0, 90.0,
- L_ORTHORHOMBIC, 'C', '*');
+ L_ORTHORHOMBIC, 'C', '*', rng);
/* Orthorhombic I */
fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 90.0, 90.0,
- L_ORTHORHOMBIC, 'I', '*');
+ L_ORTHORHOMBIC, 'I', '*', rng);
/* Orthorhombic F */
fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 90.0, 90.0,
- L_ORTHORHOMBIC, 'F', '*');
+ L_ORTHORHOMBIC, 'F', '*', rng);
/* Tetragonal P */
fail += check_centering(10e-10, 30e-10, 30e-10, 90.0, 90.0, 90.0,
- L_TETRAGONAL, 'P', 'a');
+ L_TETRAGONAL, 'P', 'a', rng);
fail += check_centering(30e-10, 10e-10, 30e-10, 90.0, 90.0, 90.0,
- L_TETRAGONAL, 'P', 'b');
+ L_TETRAGONAL, 'P', 'b', rng);
fail += check_centering(30e-10, 30e-10, 10e-10, 90.0, 90.0, 90.0,
- L_TETRAGONAL, 'P', 'c');
+ L_TETRAGONAL, 'P', 'c', rng);
/* Tetragonal I */
fail += check_centering(10e-10, 30e-10, 30e-10, 90.0, 90.0, 90.0,
- L_TETRAGONAL, 'I', 'a');
+ L_TETRAGONAL, 'I', 'a', rng);
fail += check_centering(30e-10, 10e-10, 30e-10, 90.0, 90.0, 90.0,
- L_TETRAGONAL, 'I', 'b');
+ L_TETRAGONAL, 'I', 'b', rng);
fail += check_centering(30e-10, 30e-10, 10e-10, 90.0, 90.0, 90.0,
- L_TETRAGONAL, 'I', 'c');
+ L_TETRAGONAL, 'I', 'c', rng);
/* Rhombohedral R */
fail += check_centering(10e-10, 10e-10, 10e-10, 60.0, 60.0, 60.0,
- L_RHOMBOHEDRAL, 'R', '*');
+ L_RHOMBOHEDRAL, 'R', '*', rng);
/* Hexagonal P */
fail += check_centering(30e-10, 10e-10, 10e-10, 120.0, 90.0, 90.0,
- L_HEXAGONAL, 'P', 'a');
+ L_HEXAGONAL, 'P', 'a', rng);
fail += check_centering(10e-10, 30e-10, 10e-10, 90.0, 120.0, 90.0,
- L_HEXAGONAL, 'P', 'b');
+ L_HEXAGONAL, 'P', 'b', rng);
fail += check_centering(10e-10, 10e-10, 30e-10, 90.0, 90.0, 120.0,
- L_HEXAGONAL, 'P', 'c');
+ L_HEXAGONAL, 'P', 'c', rng);
/* Hexagonal H (PDB-speak for rhombohedral) */
fail += check_centering(20e-10, 20e-10, 40e-10, 90.0, 90.0, 120.0,
- L_HEXAGONAL, 'H', 'c');
+ L_HEXAGONAL, 'H', 'c', rng);
/* Cubic P */
fail += check_centering(30e-10, 30e-10, 30e-10, 90.0, 90.0, 90.0,
- L_CUBIC, 'P', '*');
+ L_CUBIC, 'P', '*', rng);
/* Cubic I */
fail += check_centering(30e-10, 30e-10, 30e-10, 90.0, 90.0, 90.0,
- L_CUBIC, 'I', '*');
+ L_CUBIC, 'I', '*', rng);
/* Cubic F */
fail += check_centering(30e-10, 30e-10, 30e-10, 90.0, 90.0, 90.0,
- L_CUBIC, 'F', '*');
+ L_CUBIC, 'F', '*', rng);
+
+ gsl_rng_free(rng);
return fail;
}
diff --git a/tests/gpu_sim_check.c b/tests/gpu_sim_check.c
index c183a2a2..d2475213 100644
--- a/tests/gpu_sim_check.c
+++ b/tests/gpu_sim_check.c
@@ -3,7 +3,11 @@
*
* Check that GPU simulation agrees with CPU version
*
- * Copyright © 2012 Thomas White <taw@physics.org>
+ * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY,
+ * a research centre of the Helmholtz Association.
+ *
+ * Authors:
+ * 2012-2014 Thomas White <taw@physics.org>
*
* This file is part of CrystFEL.
*
@@ -80,6 +84,9 @@ int main(int argc, char *argv[])
double start, end;
double gpu_time, cpu_time;
SymOpList *sym;
+ gsl_rng *rng;
+
+ rng = gsl_rng_alloc(gsl_rng_mt19937);
gctx = setup_gpu(1, NULL, NULL, NULL, 0);
if ( gctx == NULL ) {
@@ -90,7 +97,7 @@ int main(int argc, char *argv[])
cell_raw = cell_new_from_parameters(28.1e-9, 28.1e-9, 16.5e-9,
deg2rad(90.0), deg2rad(90.0), deg2rad(120.0));
- cell = cell_rotate(cell_raw, random_quaternion());
+ cell = cell_rotate(cell_raw, random_quaternion(rng));
gpu_image.width = 1024;
gpu_image.height = 1024;
@@ -217,5 +224,7 @@ int main(int argc, char *argv[])
}
+ gsl_rng_free(rng);
+
return 0;
}
diff --git a/tests/integration_check.c b/tests/integration_check.c
index f1b5384b..0196742f 100644
--- a/tests/integration_check.c
+++ b/tests/integration_check.c
@@ -43,7 +43,7 @@ int main(int argc, char *argv[])
struct image image;
int fs, ss;
FILE *fh;
- unsigned int seed;
+ unsigned long int seed;
int fail = 0;
const int w = 128;
const int h = 128;
@@ -57,11 +57,14 @@ int main(int argc, char *argv[])
int i;
Histogram *hi;
double esd_sum = 0.0;
+ gsl_rng *rng;
+
+ rng = gsl_rng_alloc(gsl_rng_mt19937);
fh = fopen("/dev/urandom", "r");
fread(&seed, sizeof(seed), 1, fh);
fclose(fh);
- srand(seed);
+ gsl_rng_set(rng, seed);
image.flags = NULL;
image.beam = NULL;
@@ -110,7 +113,7 @@ int main(int argc, char *argv[])
for ( fs=0; fs<w; fs++ ) {
for ( ss=0; ss<h; ss++ ) {
- image.dp[0][fs+w*ss] = 10.0*poisson_noise(40);
+ image.dp[0][fs+w*ss] = 10.0*poisson_noise(rng, 40);
if ( (fs-64)*(fs-64) + (ss-64)*(ss-64) > 2 ) continue;
//image.dp[0][fs+w*ss] += 10.0*poisson_noise(10);
}
@@ -124,7 +127,7 @@ int main(int argc, char *argv[])
cell_set_centering(cell, 'P');
cell_set_parameters(cell, 800.0e-10, 800.0e-10, 800.0e-10,
deg2rad(90.0), deg2rad(90.0), deg2rad(90.0));
- cell = cell_rotate(cell, random_quaternion());
+ cell = cell_rotate(cell, random_quaternion(rng));
ic.halfw = ir_out;
ic.image = &image;
diff --git a/tests/pr_l_gradient_check.c b/tests/pr_l_gradient_check.c
index a4fd4526..4b3894ba 100644
--- a/tests/pr_l_gradient_check.c
+++ b/tests/pr_l_gradient_check.c
@@ -3,11 +3,11 @@
*
* Check Lorentz factor gradients for post refinement
*
- * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY,
* a research centre of the Helmholtz Association.
*
* Authors:
- * 2012-2013 Thomas White <taw@physics.org>
+ * 2012-2014 Thomas White <taw@physics.org>
*
* This file is part of CrystFEL.
*
@@ -284,6 +284,7 @@ int main(int argc, char *argv[])
int quiet = 0;
int plot = 0;
int c;
+ gsl_rng *rng;
const struct option longopts[] = {
{"quiet", 0, &quiet, 1},
@@ -334,12 +335,14 @@ int main(int argc, char *argv[])
deg2rad(90.0),
deg2rad(90.0));
+ rng = gsl_rng_alloc(gsl_rng_mt19937);
+
for ( i=0; i<1; i++ ) {
UnitCell *rot;
double val;
- orientation = random_quaternion();
+ orientation = random_quaternion(rng);
rot = cell_rotate(cell, orientation);
crystal_set_cell(cr, rot);
@@ -354,5 +357,7 @@ int main(int argc, char *argv[])
}
+ gsl_rng_free(rng);
+
return fail;
}
diff --git a/tests/pr_p_gradient_check.c b/tests/pr_p_gradient_check.c
index 9adf4780..c0582036 100644
--- a/tests/pr_p_gradient_check.c
+++ b/tests/pr_p_gradient_check.c
@@ -3,11 +3,11 @@
*
* Check partiality gradients for post refinement
*
- * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY,
* a research centre of the Helmholtz Association.
*
* Authors:
- * 2012-2013 Thomas White <taw@physics.org>
+ * 2012-2014 Thomas White <taw@physics.org>
*
* This file is part of CrystFEL.
*
@@ -397,6 +397,7 @@ int main(int argc, char *argv[])
int quiet = 0;
int plot = 0;
int c;
+ gsl_rng *rng;
const struct option longopts[] = {
{"quiet", 0, &quiet, 1},
@@ -447,12 +448,14 @@ int main(int argc, char *argv[])
deg2rad(90.0),
deg2rad(90.0));
+ rng = gsl_rng_alloc(gsl_rng_mt19937);
+
for ( i=0; i<1; i++ ) {
UnitCell *rot;
double val;
- orientation = random_quaternion();
+ orientation = random_quaternion(rng);
rot = cell_rotate(cell, orientation);
crystal_set_cell(cr, rot);
@@ -511,5 +514,7 @@ int main(int argc, char *argv[])
}
+ gsl_rng_free(rng);
+
return fail;
}
diff --git a/tests/pr_pl_gradient_check.c b/tests/pr_pl_gradient_check.c
index f603fff0..5ef01bf5 100644
--- a/tests/pr_pl_gradient_check.c
+++ b/tests/pr_pl_gradient_check.c
@@ -3,11 +3,11 @@
*
* Check partiality gradients for post refinement
*
- * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY,
* a research centre of the Helmholtz Association.
*
* Authors:
- * 2012-2013 Thomas White <taw@physics.org>
+ * 2012-2014 Thomas White <taw@physics.org>
*
* This file is part of CrystFEL.
*
@@ -399,6 +399,7 @@ int main(int argc, char *argv[])
int quiet = 0;
int plot = 0;
int c;
+ gsl_rng *rng;
const struct option longopts[] = {
{"quiet", 0, &quiet, 1},
@@ -449,12 +450,14 @@ int main(int argc, char *argv[])
deg2rad(90.0),
deg2rad(90.0));
+ rng = gsl_rng_alloc(gsl_rng_mt19937);
+
for ( i=0; i<1; i++ ) {
UnitCell *rot;
double val;
- orientation = random_quaternion();
+ orientation = random_quaternion(rng);
rot = cell_rotate(cell, orientation);
crystal_set_cell(cr, rot);
@@ -514,5 +517,7 @@ int main(int argc, char *argv[])
}
+ gsl_rng_free(rng);
+
return fail;
}
diff --git a/tests/prof2d_check.c b/tests/prof2d_check.c
index c1718a92..f1665a6c 100644
--- a/tests/prof2d_check.c
+++ b/tests/prof2d_check.c
@@ -3,7 +3,11 @@
*
* Check 2D profile fitting
*
- * Copyright © 2013 Thomas White <taw@physics.org>
+ * Copyright © 2013-2014 Deutsches Elektronen-Synchrotron DESY,
+ * a research centre of the Helmholtz Association.
+ *
+ * Authors:
+ * 2013-2014 Thomas White <taw@physics.org>
*
* This file is part of CrystFEL.
*
@@ -48,7 +52,7 @@ int main(int argc, char *argv[])
struct image image;
int fs, ss;
FILE *fh;
- unsigned int seed;
+ unsigned long int seed;
int fail = 0;
const int w = 1024;
const int h = 1024;
@@ -65,11 +69,14 @@ int main(int argc, char *argv[])
int n = 0;
int n_strong = 0;
int n_weak = 0;
+ gsl_rng *rng;
+
+ rng = gsl_rng_alloc(gsl_rng_mt19937);
fh = fopen("/dev/urandom", "r");
fread(&seed, sizeof(seed), 1, fh);
fclose(fh);
- srand(seed);
+ gsl_rng_set(rng, seed);
image.flags = NULL;
image.beam = NULL;
@@ -121,7 +128,7 @@ int main(int argc, char *argv[])
cell_set_centering(cell, 'P');
cell_set_parameters(cell, 800.0e-10, 800.0e-10, 800.0e-10,
deg2rad(90.0), deg2rad(90.0), deg2rad(90.0));
- cell = cell_rotate(cell, random_quaternion());
+ cell = cell_rotate(cell, random_quaternion(rng));
cr = crystal_new();
crystal_set_profile_radius(cr, 0.001e9);
@@ -136,7 +143,7 @@ int main(int argc, char *argv[])
for ( fs=0; fs<w; fs++ ) {
for ( ss=0; ss<h; ss++ ) {
- image.dp[0][fs+w*ss] = 10.0*poisson_noise(40);
+ image.dp[0][fs+w*ss] = 10.0*poisson_noise(rng, 40);
}
}
@@ -154,11 +161,11 @@ int main(int argc, char *argv[])
get_indices(refl, &hj, &kj, &lj);
if ( lj % 2 ) {
const int pk_ph = 1000;
- ADD_PX(fs, ss, 10.0*poisson_noise(pk_ph));
- ADD_PX(fs-1, ss, 10.0*poisson_noise(pk_ph));
- ADD_PX(fs+1, ss, 10.0*poisson_noise(pk_ph));
- ADD_PX(fs, ss-1, 10.0*poisson_noise(pk_ph));
- ADD_PX(fs, ss+1, 10.0*poisson_noise(pk_ph));
+ ADD_PX(fs, ss, 10.0*poisson_noise(rng, pk_ph));
+ ADD_PX(fs-1, ss, 10.0*poisson_noise(rng, pk_ph));
+ ADD_PX(fs+1, ss, 10.0*poisson_noise(rng, pk_ph));
+ ADD_PX(fs, ss-1, 10.0*poisson_noise(rng, pk_ph));
+ ADD_PX(fs, ss+1, 10.0*poisson_noise(rng, pk_ph));
n_strong++;
} else {
/* Absent peak */
@@ -226,6 +233,7 @@ int main(int argc, char *argv[])
free(image.det);
free(image.dp[0]);
free(image.dp);
+ gsl_rng_free(rng);
if ( fail ) return 1;
diff --git a/tests/ring_check.c b/tests/ring_check.c
index ac871581..f7e8c473 100644
--- a/tests/ring_check.c
+++ b/tests/ring_check.c
@@ -1,10 +1,14 @@
/*
- * integration_check.c
+ * ring_check.c
*
* Check peak integration
*
- * Copyright © 2012 Thomas White <taw@physics.org>
- * Copyright © 2012 Andrew Martin <andrew.martin@desy.de>
+ * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY,
+ * a research centre of the Helmholtz Association.
+ *
+ * Authors:
+ * 2011-2014 Thomas White <taw@physics.org>
+ * 2012 Andrew Martin <andrew.martin@desy.de>
*
* This file is part of CrystFEL.
*
@@ -41,7 +45,7 @@
/* The third integration check draws a Poisson background and checks that, on
* average, it gets subtracted by the background subtraction. */
static void third_integration_check(struct image *image, int n_trials,
- int *fail)
+ int *fail, gsl_rng *rng)
{
double mean_intensity = 0.0;
double mean_bg = 0.0;
@@ -59,7 +63,8 @@ static void third_integration_check(struct image *image, int n_trials,
for ( fs=0; fs<image->width; fs++ ) {
for ( ss=0; ss<image->height; ss++ ) {
- image->data[fs+image->width*ss] = poisson_noise(1000.0);
+ image->data[fs+image->width*ss]
+ = poisson_noise(rng, 1000.0);
}
}
@@ -100,7 +105,7 @@ static void third_integration_check(struct image *image, int n_trials,
* top of it, then checks that the intensity of the peak is correctly recovered
* accounting for the background. */
static void fourth_integration_check(struct image *image, int n_trials,
- int *fail)
+ int *fail, gsl_rng *rng)
{
double mean_intensity = 0.0;
double mean_sigma = 0.0;
@@ -118,7 +123,7 @@ static void fourth_integration_check(struct image *image, int n_trials,
for ( fs=0; fs<image->width; fs++ ) {
for ( ss=0; ss<image->height; ss++ ) {
int idx = fs+image->width*ss;
- image->data[idx] = poisson_noise(1000.0);
+ image->data[idx] = poisson_noise(rng, 1000.0);
if ( (fs-64)*(fs-64) + (ss-64)*(ss-64) > 9*9 ) continue;
image->data[idx] += 1000.0;
pcount++;
@@ -162,16 +167,19 @@ int main(int argc, char *argv[])
double fsp, ssp, intensity, sigma;
int fs, ss;
FILE *fh;
- unsigned int seed;
+ unsigned long int seed;
int fail = 0;
const int n_trials = 100;
int r, npx;
double ex;
+ gsl_rng *rng;
+
+ rng = gsl_rng_alloc(gsl_rng_mt19937);
fh = fopen("/dev/urandom", "r");
fread(&seed, sizeof(seed), 1, fh);
fclose(fh);
- srand(seed);
+ gsl_rng_set(rng, seed);
image.data = malloc(128*128*sizeof(float));
image.flags = NULL;
@@ -261,10 +269,10 @@ int main(int argc, char *argv[])
}
/* Third check: Poisson background should get mostly subtracted */
- third_integration_check(&image, n_trials, &fail);
+ third_integration_check(&image, n_trials, &fail, rng);
/* Fourth check: peak on Poisson background */
- fourth_integration_check(&image, n_trials, &fail);
+ fourth_integration_check(&image, n_trials, &fail, rng);
/* Fifth check: uniform peak on uniform background */
npx = 0;
@@ -307,6 +315,7 @@ int main(int argc, char *argv[])
free(image.det->panels);
free(image.det);
free(image.data);
+ gsl_rng_free(rng);
if ( fail ) return 1;
diff --git a/tests/transformation_check.c b/tests/transformation_check.c
index 7d25aa04..3affb652 100644
--- a/tests/transformation_check.c
+++ b/tests/transformation_check.c
@@ -3,7 +3,11 @@
*
* Check that unit cell transformations work
*
- * Copyright © 2012 Thomas White <taw@physics.org>
+ * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY,
+ * a research centre of the Helmholtz Association.
+ *
+ * Authors:
+ * 2012-2014 Thomas White <taw@physics.org>
*
* This file is part of CrystFEL.
*
@@ -111,6 +115,9 @@ int main(int argc, char *argv[])
int fail = 0;
UnitCell *cell, *cref;
UnitCellTransformation *tfn;
+ gsl_rng *rng;
+
+ rng = gsl_rng_alloc(gsl_rng_mt19937);
cref = cell_new_from_parameters(50e-10, 55e-10, 70e-10,
deg2rad(67.0),
@@ -118,7 +125,7 @@ int main(int argc, char *argv[])
deg2rad(77.0));
if ( cref == NULL ) return 1;
- cell = cell_rotate(cref, random_quaternion());
+ cell = cell_rotate(cref, random_quaternion(rng));
if ( cell == NULL ) return 1;
cell_free(cref);
@@ -171,6 +178,7 @@ int main(int argc, char *argv[])
tfn_free(tfn);
cell_free(cell);
+ gsl_rng_free(rng);
return fail;
}
generated by cgit v1.2.3 (git 2.25.1) at 2024-06-26 18:53:16 +0000