aboutsummaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
Diffstat (limited to 'src')
-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
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);
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-18 06:57:19 +0000