aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--src/partial_sim.c158
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 = &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);