diff options
| -rw-r--r-- | src/partial_sim.c | 158 |
1 files changed, 121 insertions, 37 deletions
diff --git a/src/partial_sim.c b/src/partial_sim.c index 8fc2fb21..86ad5536 100644 --- a/src/partial_sim.c +++ b/src/partial_sim.c @@ -21,6 +21,7 @@ #include <unistd.h> #include <getopt.h> #include <assert.h> +#include <pthread.h> #include "utils.h" #include "reflist-utils.h" @@ -29,6 +30,7 @@ #include "detector.h" #include "geometry.h" #include "stream.h" +#include "thread-pool.h" static void mess_up_cell(UnitCell *cell) @@ -64,7 +66,8 @@ static void mess_up_cell(UnitCell *cell) * according to "full" */ static void calculate_partials(RefList *partial, double osf, RefList *full, const SymOpList *sym, - int random_intensities) + int random_intensities, + pthread_mutex_t *full_lock) { Reflection *refl; RefListIterator *iter; @@ -81,13 +84,24 @@ static void calculate_partials(RefList *partial, double osf, get_asymm(sym, h, k, l, &h, &k, &l); p = get_partiality(refl); + pthread_mutex_lock(full_lock); rfull = find_refl(full, h, k, l); + pthread_mutex_unlock(full_lock); + if ( rfull == NULL ) { if ( random_intensities ) { + + /* The full reflection is immutable (in this + * program) once created, but creating it must + * be an atomic operation. So do the whole + * thing under lock. */ + pthread_mutex_lock(full_lock); rfull = add_refl(full, h, k, l); If = fabs(gaussian_noise(0.0, 1000.0)); set_int(rfull, If); set_redundancy(rfull, 1); + pthread_mutex_unlock(full_lock); + } else { set_redundancy(refl, 0); If = 0.0; @@ -134,6 +148,87 @@ static void show_help(const char *s) } + +struct queue_args +{ + RefList *full; + pthread_mutex_t full_lock; + + int n_done; + int n_to_do; + + SymOpList *sym; + int random_intensities; + UnitCell *cell; + + struct image *template_image; + + FILE *stream; +}; + + +struct worker_args +{ + struct queue_args *qargs; + struct image image; +}; + + +static void *create_job(void *vqargs) +{ + struct worker_args *wargs; + struct queue_args *qargs = vqargs; + + wargs = malloc(sizeof(struct worker_args)); + + wargs->qargs = qargs; + wargs->image = *qargs->template_image; + + return wargs; +} + + +static void run_job(void *vwargs, int cookie) +{ + double osf; + struct quaternion orientation; + struct worker_args *wargs = vwargs; + struct queue_args *qargs = wargs->qargs; + + osf = gaussian_noise(1.0, 0.3); + + /* Set up a random orientation */ + orientation = random_quaternion(); + wargs->image.indexed_cell = cell_rotate(qargs->cell, orientation); + + snprintf(wargs->image.filename, 255, "dummy.h5"); + wargs->image.reflections = find_intersections(&wargs->image, + wargs->image.indexed_cell); + calculate_partials(wargs->image.reflections, osf, qargs->full, + qargs->sym, qargs->random_intensities, + &qargs->full_lock); + + /* Give a slightly incorrect cell in the stream */ + mess_up_cell(wargs->image.indexed_cell); +} + + +static void finalise_job(void *vqargs, void *vwargs) +{ + struct worker_args *wargs = vwargs; + struct queue_args *qargs = vqargs; + + write_chunk(qargs->stream, &wargs->image, STREAM_INTEGRATED); + + reflist_free(wargs->image.reflections); + cell_free(wargs->image.indexed_cell); + free(wargs); + + qargs->n_done++; + progress_bar(qargs->n_done, qargs->n_to_do, "Simulating"); +} + + int main(int argc, char *argv[]) { int c; @@ -148,13 +243,13 @@ int main(int argc, char *argv[]) char *sym_str = NULL; SymOpList *sym; UnitCell *cell = NULL; - struct quaternion orientation; - struct image image; FILE *ofh; int n = 2; - int i; int random_intensities = 0; char *save_file = NULL; + struct queue_args qargs; + struct image image; + int n_threads = 1; /* Long options */ const struct option longopts[] = { @@ -170,7 +265,7 @@ int main(int argc, char *argv[]) }; /* Short options */ - while ((c = getopt_long(argc, argv, "hi:o:b:p:g:y:n:r:", + while ((c = getopt_long(argc, argv, "hi:o:b:p:g:y:n:r:j:", longopts, NULL)) != -1) { switch (c) { @@ -210,6 +305,10 @@ int main(int argc, char *argv[]) save_file = strdup(optarg); break; + case 'j' : + n_threads = atoi(optarg); + break; + case 0 : break; @@ -218,6 +317,11 @@ int main(int argc, char *argv[]) } } + if ( n_threads < 1 ) { + ERROR("Invalid number of threads.\n"); + return 1; + } + /* Load beam */ if ( beamfile == NULL ) { ERROR("You need to provide a beam parameters file.\n"); @@ -318,38 +422,18 @@ int main(int argc, char *argv[]) full = reflist_new(); } - for ( i=0; i<n; i++ ) { - - double osf; - - //if ( random() > RAND_MAX/2 ) { - // osf = 1.0; - //} else { - // osf = 2.0; - //} - //STATUS("Image %i scale factor %f\n", i, osf); - osf = gaussian_noise(1.0, 0.3); - - /* Set up a random orientation */ - orientation = random_quaternion(); - image.indexed_cell = cell_rotate(cell, orientation); - - snprintf(image.filename, 255, "dummy.h5"); - image.reflections = find_intersections(&image, - image.indexed_cell); - calculate_partials(image.reflections, osf, full, sym, - random_intensities); - - /* Give a slightly incorrect cell in the stream */ - mess_up_cell(image.indexed_cell); - write_chunk(ofh, &image, STREAM_INTEGRATED); - - reflist_free(image.reflections); - cell_free(image.indexed_cell); - - progress_bar(i+1, n, "Simulating"); - - } + qargs.full = full; + pthread_mutex_init(&qargs.full_lock, NULL); + qargs.n_to_do = n; + qargs.n_done = 0; + qargs.sym = sym; + qargs.random_intensities = random_intensities; + qargs.cell = cell; + qargs.template_image = ℑ + qargs.stream = ofh; + + run_threads(n_threads, run_job, create_job, finalise_job, + &qargs, n, n, 1, 0); if ( random_intensities ) { STATUS("Writing full intensities to %s\n", save_file); |