diff options
author | Thomas White <taw@physics.org> | 2013-03-07 10:36:26 +0100 |
---|---|---|
committer | Thomas White <taw@physics.org> | 2013-03-07 10:36:26 +0100 |
commit | 33260b9a7bf27a87afa1dde789f0b0d9ce28eaab (patch) | |
tree | 686fcdfa791e2f4a14eb2f5575971021cde0d86f /src/partialator.c | |
parent | 22cd1d08e6eeb6c8e43a709021746c057f6661d7 (diff) | |
parent | 1e03ed982741fdc576ec5a915da120450df20499 (diff) |
Merge branch 'tom/multicrystal'
Conflicts:
libcrystfel/src/mosflm.c
Diffstat (limited to 'src/partialator.c')
-rw-r--r-- | src/partialator.c | 235 |
1 files changed, 139 insertions, 96 deletions
diff --git a/src/partialator.c b/src/partialator.c index c2f6c299..a4af3c18 100644 --- a/src/partialator.c +++ b/src/partialator.c @@ -3,11 +3,11 @@ * * Scaling and post refinement for coherent nanocrystallography * - * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY, - * a research centre of the Helmholtz Association. + * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. * * Authors: - * 2010-2012 Thomas White <taw@physics.org> + * 2010-2013 Thomas White <taw@physics.org> * * This file is part of CrystFEL. * @@ -87,16 +87,16 @@ static void show_help(const char *s) struct refine_args { RefList *full; - struct image *image; + Crystal *crystal; }; struct queue_args { - int n; + int n_started; int n_done; - int n_total_patterns; - struct image *images; + Crystal **crystals; + int n_crystals; struct refine_args task_defaults; }; @@ -104,10 +104,9 @@ struct queue_args static void refine_image(void *task, int id) { struct refine_args *pargs = task; - struct image *image = pargs->image; - image->id = id; + Crystal *cr = pargs->crystal; - pr_refine(image, pargs->full); + pr_refine(cr, pargs->full); } @@ -119,9 +118,9 @@ static void *get_image(void *vqargs) task = malloc(sizeof(struct refine_args)); memcpy(task, &qargs->task_defaults, sizeof(struct refine_args)); - task->image = &qargs->images[qargs->n]; + task->crystal = qargs->crystals[qargs->n_started]; - qargs->n++; + qargs->n_started++; return task; } @@ -133,12 +132,12 @@ static void done_image(void *vqargs, void *task) qargs->n_done++; - progress_bar(qargs->n_done, qargs->n_total_patterns, "Refining"); + progress_bar(qargs->n_done, qargs->n_crystals, "Refining"); free(task); } -static void refine_all(struct image *images, int n_total_patterns, +static void refine_all(Crystal **crystals, int n_crystals, struct detector *det, RefList *full, int nthreads) { @@ -146,19 +145,19 @@ static void refine_all(struct image *images, int n_total_patterns, struct queue_args qargs; task_defaults.full = full; - task_defaults.image = NULL; + task_defaults.crystal = NULL; qargs.task_defaults = task_defaults; - qargs.n = 0; + qargs.n_started = 0; qargs.n_done = 0; - qargs.n_total_patterns = n_total_patterns; - qargs.images = images; + qargs.n_crystals = n_crystals; + qargs.crystals = crystals; /* Don't have threads which are doing nothing */ - if ( n_total_patterns < nthreads ) nthreads = n_total_patterns; + if ( n_crystals < nthreads ) nthreads = n_crystals; run_threads(nthreads, refine_image, get_image, done_image, - &qargs, n_total_patterns, 0, 0, 0); + &qargs, n_crystals, 0, 0, 0); } @@ -201,13 +200,14 @@ static int select_scalable_reflections(RefList *list, RefList *reference) } -static void select_reflections_for_refinement(struct image *images, int n, +static void select_reflections_for_refinement(Crystal **crystals, int n, RefList *full, int have_reference) { int i; for ( i=0; i<n; i++ ) { + RefList *reflist; Reflection *refl; RefListIterator *iter; int n_acc = 0; @@ -216,9 +216,10 @@ static void select_reflections_for_refinement(struct image *images, int n, int n_fewmatch = 0; int n_ref = 0; - if ( images[i].pr_dud ) continue; + if ( crystal_get_user_flag(crystals[i]) ) continue; - for ( refl = first_refl(images[i].reflections, &iter); + reflist = crystal_get_reflections(crystals[i]); + for ( refl = first_refl(reflist, &iter); refl != NULL; refl = next_refl(refl, iter) ) { @@ -272,6 +273,20 @@ static void select_reflections_for_refinement(struct image *images, int n, } +static void display_progress(int n_images, int n_crystals) +{ + if ( !isatty(STDERR_FILENO) ) return; + if ( tcgetpgrp(STDERR_FILENO) != getpgrp() ) return; + + pthread_mutex_lock(&stderr_lock); + fprintf(stderr, "\r%i images loaded, %i crystals.", + n_images, n_crystals); + pthread_mutex_unlock(&stderr_lock); + + fflush(stdout); +} + + int main(int argc, char *argv[]) { int c; @@ -280,16 +295,15 @@ int main(int argc, char *argv[]) char *geomfile = NULL; char *sym_str = NULL; SymOpList *sym; - FILE *fh; int nthreads = 1; struct detector *det; int i; - int n_total_patterns; struct image *images; int n_iter = 10; struct beam_params *beam = NULL; RefList *full; - int n_usable_patterns = 0; + int n_images = 0; + int n_crystals = 0; int nobs; char *reference_file = NULL; RefList *reference = NULL; @@ -298,6 +312,8 @@ int main(int argc, char *argv[]) char cmdline[1024]; SRContext *sr; int noscale = 0; + Stream *st; + Crystal **crystals; /* Long options */ const struct option longopts[] = { @@ -386,19 +402,15 @@ int main(int argc, char *argv[]) return 1; } - /* Sanitise input filename and open */ if ( infile == NULL ) { infile = strdup("-"); } - if ( strcmp(infile, "-") == 0 ) { - fh = stdin; - } else { - fh = fopen(infile, "r"); - } - if ( fh == NULL ) { - ERROR("Failed to open input file '%s'\n", infile); + st = open_stream_for_read(infile); + if ( st == NULL ) { + ERROR("Failed to open input stream '%s'\n", infile); return 1; } + /* Don't free "infile", because it's needed for the scaling report */ /* Sanitise output filename */ if ( outfile == NULL ) { @@ -438,86 +450,116 @@ int main(int argc, char *argv[]) } - n_total_patterns = count_patterns(fh); - if ( n_total_patterns == 0 ) { - ERROR("No patterns to process.\n"); - return 1; - } - STATUS("There are %i patterns to process\n", n_total_patterns); - gsl_set_error_handler_off(); - images = malloc(n_total_patterns * sizeof(struct image)); - if ( images == NULL ) { - ERROR("Couldn't allocate memory for images.\n"); - return 1; - } - /* Fill in what we know about the images so far */ - rewind(fh); - nobs = 0; - for ( i=0; i<n_total_patterns; i++ ) { + n_images = 0; + n_crystals = 0; + images = NULL; + crystals = NULL; + do { RefList *as; + int i; + struct image *images_new; struct image *cur; - cur = &images[n_usable_patterns]; + images_new = realloc(images, (n_images+1)*sizeof(struct image)); + if ( images_new == NULL ) { + ERROR("Failed to allocate memory for image list.\n"); + return 1; + } + images = images_new; + cur = &images[n_images]; cur->det = det; - - if ( read_chunk(fh, cur) != 0 ) { - /* Should not happen, because we counted the patterns - * earlier. */ - ERROR("Failed to read chunk from the input stream.\n"); - return 1; + if ( read_chunk(st, cur) != 0 ) { + break; } /* Won't be needing this, if it exists */ image_feature_list_free(cur->features); cur->features = NULL; - - /* "n_usable_patterns" will not be incremented in this case */ - if ( cur->indexed_cell == NULL ) continue; - - /* Fill in initial estimates of stuff */ cur->div = beam->divergence; cur->bw = beam->bandwidth; cur->width = det->max_fs; cur->height = det->max_ss; - cur->osf = 1.0; - cur->profile_radius = beam->profile_radius; - cur->pr_dud = 0; - - /* Muppet proofing */ cur->data = NULL; cur->flags = NULL; cur->beam = NULL; - /* This is the raw list of reflections */ - as = asymmetric_indices(cur->reflections, sym); - reflist_free(cur->reflections); - cur->reflections = as; + n_images++; + + for ( i=0; i<cur->n_crystals; i++ ) { + + Crystal *cr; + Crystal **crystals_new; + RefList *cr_refl; + + crystals_new = realloc(crystals, + (n_crystals+1)*sizeof(Crystal *)); + if ( crystals_new == NULL ) { + ERROR("Failed to allocate memory for crystal " + "list.\n"); + return 1; + } + crystals = crystals_new; + crystals[n_crystals] = cur->crystals[i]; + cr = crystals[n_crystals]; + + /* Image pointer will change due to later reallocs */ + crystal_set_image(cr, NULL); + + /* Fill in initial estimates of stuff */ + crystal_set_osf(cr, 1.0); + crystal_set_profile_radius(cr, beam->profile_radius); + crystal_set_user_flag(cr, 0); + + /* This is the raw list of reflections */ + cr_refl = crystal_get_reflections(cr); + as = asymmetric_indices(cr_refl, sym); + crystal_set_reflections(cr, as); + reflist_free(cr_refl); - update_partialities(cur); + n_crystals++; - nobs += select_scalable_reflections(cur->reflections, - reference); + } + + display_progress(n_images, n_crystals); + + } while ( 1 ); + + close_stream(st); + + /* Fill in image pointers */ + nobs = 0; + for ( i=0; i<n_images; i++ ) { + int j; + for ( j=0; j<images[i].n_crystals; j++ ) { + + Crystal *cryst; + RefList *as; - progress_bar(i, n_total_patterns-1, "Loading pattern data"); - n_usable_patterns++; + cryst = images[i].crystals[j]; + crystal_set_image(cryst, &images[i]); + /* Now it's safe to do the following */ + update_partialities(cryst); + as = crystal_get_reflections(cryst); + nobs += select_scalable_reflections(as, reference); + + } } - fclose(fh); /* Make initial estimates */ - STATUS("Performing initial scaling.\n"); + STATUS("\nPerforming initial scaling.\n"); if ( noscale ) STATUS("Scale factors fixed at 1.\n"); - full = scale_intensities(images, n_usable_patterns, reference, + full = scale_intensities(crystals, n_crystals, reference, nthreads, noscale); - sr = sr_titlepage(images, n_usable_patterns, "scaling-report.pdf", + sr = sr_titlepage(crystals, n_crystals, "scaling-report.pdf", infile, cmdline); - sr_iteration(sr, 0, images, n_usable_patterns, full); + sr_iteration(sr, 0, crystals, n_crystals, full); /* Iterate */ for ( i=0; i<n_iter; i++ ) { @@ -534,54 +576,55 @@ int main(int argc, char *argv[]) } /* Refine the geometry of all patterns to get the best fit */ - select_reflections_for_refinement(images, n_usable_patterns, + select_reflections_for_refinement(crystals, n_crystals, comp, have_reference); - refine_all(images, n_usable_patterns, det, comp, nthreads); + refine_all(crystals, n_crystals, det, comp, nthreads); nobs = 0; - for ( j=0; j<n_usable_patterns; j++ ) { + for ( j=0; j<n_crystals; j++ ) { - struct image *cur = &images[j]; + Crystal *cr = crystals[j]; + RefList *rf = crystal_get_reflections(cr); - nobs += select_scalable_reflections(cur->reflections, - reference); + nobs += select_scalable_reflections(rf, reference); } /* Re-estimate all the full intensities */ reflist_free(full); - full = scale_intensities(images, n_usable_patterns, + full = scale_intensities(crystals, n_crystals, reference, nthreads, noscale); - sr_iteration(sr, i+1, images, n_usable_patterns, full); + sr_iteration(sr, i+1, crystals, n_crystals, full); } sr_finish(sr); n_dud = 0; - for ( i=0; i<n_usable_patterns; i++ ) { - if ( images[i].pr_dud ) n_dud++; + for ( i=0; i<n_crystals; i++ ) { + if ( crystal_get_user_flag(crystals[i]) ) n_dud++; } - STATUS("%i images could not be refined on the last cycle.\n", n_dud); + STATUS("%i crystals could not be refined on the last cycle.\n", n_dud); /* Output results */ write_reflist(outfile, full); /* Clean up */ - for ( i=0; i<n_usable_patterns; i++ ) { - reflist_free(images[i].reflections); + for ( i=0; i<n_crystals; i++ ) { + reflist_free(crystal_get_reflections(crystals[i])); + crystal_free(crystals[i]); } reflist_free(full); free(sym); free(outfile); free_detector_geometry(det); free(beam); + free(crystals); if ( reference != NULL ) { reflist_free(reference); } - for ( i=0; i<n_usable_patterns; i++ ) { - cell_free(images[i].indexed_cell); + for ( i=0; i<n_images; i++ ) { free(images[i].filename); } free(images); |