diff options
author | Thomas White <taw@bitwiz.org.uk> | 2010-09-19 23:05:25 +0200 |
---|---|---|
committer | Thomas White <taw@physics.org> | 2012-02-22 15:26:58 +0100 |
commit | 1757b74cf13ceb57ef056a0140c3e97358e640f0 (patch) | |
tree | 6653155ba53137b7a35122c7829b2731b0f9d172 | |
parent | 09e36df6c83b23b434238e9aa2bef29cd04712f0 (diff) |
facetron: Add threads
-rw-r--r-- | src/facetron.c | 341 |
1 files changed, 275 insertions, 66 deletions
diff --git a/src/facetron.c b/src/facetron.c index 722f8195..7da96153 100644 --- a/src/facetron.c +++ b/src/facetron.c @@ -20,20 +20,39 @@ #include <string.h> #include <unistd.h> #include <getopt.h> -#include <errno.h> +#include <pthread.h> +#include <sys/time.h> +#include <assert.h> -#include "image.h" -#include "cell.h" +#include "utils.h" #include "hdf5-file.h" -#include "detector.h" -#include "geometry.h" +#include "symmetry.h" + + +#define MAX_THREADS (256) + +struct process_args +{ + char *filename; + int id; + + /* Thread control */ + pthread_mutex_t control_mutex; /* Protects the scary stuff below */ + int start; + int finish; + int done; + + UnitCell *cell; + struct detector *det; + char *sym; +}; static void show_help(const char *s) { printf("Syntax: %s [options]\n\n", s); printf( -"Profile fitting for coherent nanocrystallography.\n" +"Post-refinement and profile fitting for coherent nanocrystallography.\n" "\n" " -h, --help Display this help message.\n" "\n" @@ -43,7 +62,83 @@ static void show_help(const char *s) " -x, --prefix=<p> Prefix filenames from input file with <p>.\n" " --basename Remove the directory parts of the filenames.\n" " --no-check-prefix Don't attempt to correct the --prefix.\n" -); +" -j <n> Run <n> analyses in parallel.\n"); +} + + +static void process_image(struct process_args *pargs) +{ + struct hdfile *hdfile; + struct image image; + + image.features = NULL; + image.data = NULL; + image.flags = NULL; + image.indexed_cell = NULL; + image.id = pargs->id; + image.filename = pargs->filename; + image.hits = NULL; + image.n_hits = 0; + image.det = pargs->det; + + /* View head-on (unit cell is tilted) */ + image.orientation.w = 1.0; + image.orientation.x = 0.0; + image.orientation.y = 0.0; + image.orientation.z = 0.0; + + STATUS("Processing '%s'\n", pargs->filename); + + hdfile = hdfile_open(pargs->filename); + if ( hdfile == NULL ) { + return; + } else if ( hdfile_set_first_image(hdfile, "/") ) { + ERROR("Couldn't select path\n"); + hdfile_close(hdfile); + return; + } + + hdf5_read(hdfile, &image, 1); + + + free(image.data); + cell_free(pargs->cell); + if ( image.flags != NULL ) free(image.flags); + hdfile_close(hdfile); +} + + +static void *worker_thread(void *pargsv) +{ + struct process_args *pargs = pargsv; + int finish; + + do { + + int wakeup; + + process_image(pargs); + + pthread_mutex_lock(&pargs->control_mutex); + pargs->done = 1; + pargs->start = 0; + pthread_mutex_unlock(&pargs->control_mutex); + + /* Go to sleep until told to exit or process next image */ + do { + + pthread_mutex_lock(&pargs->control_mutex); + /* Either of these can result in the thread waking up */ + wakeup = pargs->start || pargs->finish; + finish = pargs->finish; + pthread_mutex_unlock(&pargs->control_mutex); + usleep(20000); + + } while ( !wakeup ); + + } while ( !pargs->finish ); + + return NULL; } @@ -116,28 +211,21 @@ static int find_chunk(FILE *fh, UnitCell **cell, char **filename) } -static char *twiddle_filename(char *filename, int config_basename, - const char *prefix) +static void add_to_mean(UnitCell *cell, double *ast, double *bst, double *cst, + double *alst, double *best, double *gast) { - char *f3; - size_t len; - - if ( config_basename ) { - char *f2; - /* Happy with either the GNU or non-GNU versions of basename() - * here. Because _GNU_SOURCE is defined in configure.ac, we - * should have the GNU version. */ - f2 = strdup(basename(filename)); - free(filename); - filename = f2; - } - - len = strlen(prefix)+strlen(filename)+1; - f3 = malloc(len); - snprintf(f3, len, "%s%s", prefix, filename); - - free(filename); - return f3; + double asx, asy, asz; + double bsx, bsy, bsz; + double csx, csy, csz; + + cell_get_reciprocal(cell, &asx, &asy, &asz, &bsx, &bsy, &bsz, + &csx, &csy, &csz); + *ast += modulus(asx, asy, asz); + *bst += modulus(bsx, bsy, bsz); + *cst += modulus(csx, csy, csz); + *alst += angle_between(bsx, bsy, bsz, csx, csy, csz); + *best += angle_between(asx, asy, asz, csx, csy, csz); + *gast += angle_between(asx, asy, asz, bsx, bsy, bsz); } @@ -145,14 +233,21 @@ int main(int argc, char *argv[]) { int c; char *infile = NULL; + char *geomfile = NULL; FILE *fh; - UnitCell *cell; - char *filename; - struct detector *det; - char *geometry = NULL; + int rval; + int n_images; char *prefix = NULL; + int nthreads = 1; + pthread_t workers[MAX_THREADS]; + struct process_args *worker_args[MAX_THREADS]; + int worker_active[MAX_THREADS]; int config_basename = 0; int config_checkprefix = 1; + struct detector *det; + int i; + char *sym = NULL; + double as, bs, cs, als, bes, gas; /* Long options */ const struct option longopts[] = { @@ -166,7 +261,8 @@ int main(int argc, char *argv[]) }; /* Short options */ - while ((c = getopt_long(argc, argv, "hi:g:x:", longopts, NULL)) != -1) { + while ((c = getopt_long(argc, argv, "hi:g:x:j:", + longopts, NULL)) != -1) { switch (c) { case 'h' : @@ -178,13 +274,17 @@ int main(int argc, char *argv[]) break; case 'g' : - geometry = strdup(optarg); + geomfile = strdup(optarg); break; case 'x' : prefix = strdup(optarg); break; + case 'j' : + nthreads = atoi(optarg); + break; + case 0 : break; @@ -194,11 +294,18 @@ int main(int argc, char *argv[]) } - if ( infile == NULL ) infile = strdup("-"); - fh = fopen(infile, "r"); + + if ( infile == NULL ) { + infile = strdup("-"); + } + if ( strcmp(infile, "-") == 0 ) { + fh = stdin; + } else { + fh = fopen(infile, "r"); + } if ( fh == NULL ) { - ERROR("Couldn't open input stream '%s'\n", infile); - return ENOENT; + ERROR("Failed to open input file '%s'\n", infile); + return 1; } free(infile); @@ -210,51 +317,153 @@ int main(int argc, char *argv[]) } } - det = get_detector_geometry(geometry); + det = get_detector_geometry(geomfile); if ( det == NULL ) { - ERROR("Failed to read detector geometry from '%s'\n", geometry); + ERROR("Failed to read detector geometry from '%s'\n", geomfile); return 1; } - free(geometry); + free(geomfile); - /* Loop over all successfully indexed patterns */ - while ( find_chunk(fh, &cell, &filename) == 0 ) { + sym = strdup("6/mmm"); /* FIXME: Should be on command line */ - struct image image; - struct hdfile *hdfile; - struct reflhit *hits; - int np; + as = 0.0; bs = 0.0; cs = 0.0; als = 0.0; bes = 0.0; gas = 0.0; - filename = twiddle_filename(filename, config_basename, prefix); + /* Initialise worker arguments */ + for ( i=0; i<nthreads; i++ ) { - STATUS("Integrating intensities from '%s'\n", filename); + worker_args[i] = malloc(sizeof(struct process_args)); + worker_args[i]->filename = malloc(1024); + worker_active[i] = 0; + worker_args[i]->det = det; + pthread_mutex_init(&worker_args[i]->control_mutex, NULL); + worker_args[i]->sym = sym; - image.det = det; + } - hdfile = hdfile_open(filename); - if ( hdfile == NULL ) { - return ENOENT; - } else if ( hdfile_set_image(hdfile, "/data/data0") ) { - ERROR("Couldn't select path\n"); - return ENOENT; - } + n_images = 0; - hdf5_read(hdfile, &image, 0); + /* Start threads off */ + for ( i=0; i<nthreads; i++ ) { - hits = find_intersections(&image, cell, 5.0e-3, 0.0001, &np, 1); + struct process_args *pargs; + int r; + int rval; + char *filename; + UnitCell *cell; - hdfile_close(hdfile); - cell_free(cell); + pargs = worker_args[i]; + + /* Get the next filename */ + rval = find_chunk(fh, &cell, &filename); + if ( rval == 1 ) break; + add_to_mean(cell, &as, &bs, &cs, &als, &bes, &gas); + if ( config_basename ) { + char *tmp; + tmp = basename(filename); + free(filename); + filename = tmp; + } + snprintf(pargs->filename, 1023, "%s%s", + prefix, filename); + pargs->cell = cell; free(filename); - free(image.data); - free(image.flags); + + n_images++; + + pthread_mutex_lock(&pargs->control_mutex); + pargs->done = 0; + pargs->start = 1; + pargs->finish = 0; + pthread_mutex_unlock(&pargs->control_mutex); + + worker_active[i] = 1; + r = pthread_create(&workers[i], NULL, worker_thread, pargs); + if ( r != 0 ) { + worker_active[i] = 0; + ERROR("Couldn't start thread %i\n", i); + } + + } + + /* Keep threads busy until the end of the data */ + do { + + int i; + rval = 0; + + for ( i=0; i<nthreads; i++ ) { + + struct process_args *pargs; + int done; + char *filename; + UnitCell *cell; + + /* Spend time working, not managing threads */ + usleep(100000); + + /* Are we using this thread record at all? */ + if ( !worker_active[i] ) continue; + + /* Has the thread finished yet? */ + pargs = worker_args[i]; + pthread_mutex_lock(&pargs->control_mutex); + done = pargs->done; + pthread_mutex_unlock(&pargs->control_mutex); + if ( !done ) continue; + + /* Get the next filename */ + rval = find_chunk(fh, &cell, &filename); + if ( rval == 1 ) break; + add_to_mean(cell, &as, &bs, &cs, &als, &bes, &gas); + if ( config_basename ) { + char *tmp; + tmp = basename(filename); + free(filename); + filename = tmp; + } + snprintf(pargs->filename, 1023, "%s%s", + prefix, filename); + pargs->cell = cell; + free(filename); + + n_images++; + + STATUS("Done %i images\n", n_images); + + /* Wake the thread up ... */ + pthread_mutex_lock(&pargs->control_mutex); + pargs->done = 0; + pargs->start = 1; + pthread_mutex_unlock(&pargs->control_mutex); + + } + + } while ( rval == 0 ); + + /* Join threads */ + for ( i=0; i<nthreads; i++ ) { + + if ( !worker_active[i] ) goto free; + + /* Tell the thread to exit */ + struct process_args *pargs = worker_args[i]; + pthread_mutex_lock(&pargs->control_mutex); + pargs->finish = 1; + pthread_mutex_unlock(&pargs->control_mutex); + + /* Wait for it to join */ + pthread_join(workers[i], NULL); + + free: + if ( worker_args[i]->filename != NULL ) { + free(worker_args[i]->filename); + } } - free(det->panels); - free(det); fclose(fh); - free(prefix); + + STATUS("There were %i images.\n", n_images); return 0; } |