From 2f5b44635ac37f988b3b281c5aa0aada5a9d9d03 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Tue, 8 Feb 2011 14:24:03 +0100 Subject: Don't constantly re-integrate peaks --- src/partialator.c | 242 ++++++++++++++------------------------------------ src/post-refinement.c | 105 ++++++++++------------ src/post-refinement.h | 33 +------ src/reflist.c | 27 ------ src/reflist.h | 2 - 5 files changed, 119 insertions(+), 290 deletions(-) diff --git a/src/partialator.c b/src/partialator.c index a4f3381c..37c8fbd6 100644 --- a/src/partialator.c +++ b/src/partialator.c @@ -37,10 +37,6 @@ #include "reflist.h" -/* Maximum number of iterations of NLSq to do for each image per macrocycle. */ -#define MAX_CYCLES (100) - - static void show_help(const char *s) { printf("Syntax: %s [options]\n\n", s); @@ -57,9 +53,6 @@ static void show_help(const char *s) " initial values for parameters, and nominal\n" " wavelengths if no per-shot value is found in \n" " an HDF5 file.\n" -" -x, --prefix=

Prefix filenames from input file with

.\n" -" --basename Remove the directory parts of the filenames.\n" -" --no-check-prefix Don't attempt to correct the --prefix.\n" " -y, --symmetry= Merge according to symmetry .\n" " -n, --iterations= Run cycles of scaling and post-refinement.\n" "\n" @@ -83,58 +76,8 @@ static void refine_image(int mytask, void *tasks) struct refine_args *all_args = tasks; struct refine_args *pargs = &all_args[mytask]; struct image *image = pargs->image; - double nominal_photon_energy = pargs->image->beam->photon_energy; - struct hdfile *hdfile; - int i; - double dev, last_dev; - RefList *reflections; - - hdfile = hdfile_open(image->filename); - if ( hdfile == NULL ) { - ERROR("Couldn't open '%s'\n", image->filename); - return; - } else if ( hdfile_set_image(hdfile, "/data/data0") ) { - ERROR("Couldn't select path\n"); - hdfile_close(hdfile); - return; - } - - if ( hdf5_read(hdfile, pargs->image, 0, nominal_photon_energy) ) { - ERROR("Couldn't read '%s'\n", image->filename); - hdfile_close(hdfile); - return; - } - - double a, b, c, al, be, ga; - cell_get_parameters(image->indexed_cell, &a, &b, &c, &al, &be, &ga); - STATUS("Initial cell: %5.2f %5.2f %5.2f nm %5.2f %5.2f %5.2f deg\n", - a/1.0e-9, b/1.0e-9, c/1.0e-9, - rad2deg(al), rad2deg(be), rad2deg(ga)); - - /* FIXME: Don't do this each time */ - reflections = find_intersections(image, image->indexed_cell, 0); - dev = +INFINITY; - i = 0; - do { - last_dev = dev; - dev = pr_iterate(image, pargs->i_full, pargs->sym, reflections); - STATUS("Iteration %2i: mean dev = %5.2f\n", i, dev); - i++; - } while ( (fabs(last_dev - dev) > 1.0) && (i < MAX_CYCLES) ); - mean_partial_dev(image, reflections, pargs->sym, - pargs->i_full, pargs->graph); - if ( pargs->pgraph ) { - fprintf(pargs->pgraph, "%5i %5.2f\n", mytask, dev); - } - free(image->data); - if ( image->flags != NULL ) free(image->flags); - hdfile_close(hdfile); - reflist_free(reflections); - - /* Muppet proofing */ - image->data = NULL; - image->flags = NULL; + pr_refine(image, pargs->i_full, pargs->sym); } @@ -165,87 +108,6 @@ static void refine_all(struct image *images, int n_total_patterns, } -static void uniquify(Reflection *refl, const char *sym) -{ - signed int h, k, l; - signed int ha, ka, la; - - get_indices(refl, &h, &k, &l); - get_asymm(h, k, l, &ha, &ka, &la, sym); - set_indices(refl, h, k, l); -} - - -/* FIXME: Get rid of this */ -static void integrate_image(struct image *image, ReflItemList *obs, - const char *sym) -{ - RefList *reflections; - Reflection *refl; - struct hdfile *hdfile; - double nominal_photon_energy = image->beam->photon_energy; - - hdfile = hdfile_open(image->filename); - if ( hdfile == NULL ) { - ERROR("Couldn't open '%s'\n", image->filename); - return; - } else if ( hdfile_set_image(hdfile, "/data/data0") ) { - ERROR("Couldn't select path\n"); - hdfile_close(hdfile); - return; - } - - if ( hdf5_read(hdfile, image, 0, nominal_photon_energy) ) { - ERROR("Couldn't read '%s'\n", image->filename); - hdfile_close(hdfile); - return; - } - - /* Figure out which spots should appear in this pattern */ - reflections = find_intersections(image, image->indexed_cell, 0); - - /* For each reflection, estimate the partiality */ - for ( refl = first_refl(reflections); - refl != NULL; - refl = next_refl(refl) ) { - - signed int h, k, l; - float i_partial; - float xc, yc; - double x, y; - - uniquify(refl, sym); - get_indices(refl, &h, &k, &l); - - /* Don't attempt to use spots with very small - * partialities, since it won't be accurate. */ - if ( get_partiality(refl) < 0.1 ) continue; - - /* Actual measurement of this reflection from this pattern? */ - get_detector_pos(refl, &x, &y); - if ( integrate_peak(image, x, y, - &xc, &yc, &i_partial, NULL, NULL, 1, 0) ) { - delete_refl(refl); - continue; - } - - set_int(refl, i_partial); - - if ( !find_item(obs, h, k, l) ) add_item(obs, h, k, l); - - } - image->reflections = reflections; - - free(image->data); - if ( image->flags != NULL ) free(image->flags); - hdfile_close(hdfile); - - /* Muppet proofing */ - image->data = NULL; - image->flags = NULL; -} - - /* Decide which reflections can be scaled */ static void select_scalable_reflections(struct image *images, int n) { @@ -280,12 +142,9 @@ int main(int argc, char *argv[]) char *infile = NULL; char *outfile = NULL; char *geomfile = NULL; - char *prefix = NULL; char *sym = NULL; FILE *fh; int nthreads = 1; - int config_basename = 0; - int config_checkprefix = 1; struct detector *det; unsigned int *cts; ReflItemList *obs; @@ -303,9 +162,6 @@ int main(int argc, char *argv[]) {"output", 1, NULL, 'o'}, {"geometry", 1, NULL, 'g'}, {"beam", 1, NULL, 'b'}, - {"prefix", 1, NULL, 'x'}, - {"basename", 0, &config_basename, 1}, - {"no-check-prefix", 0, &config_checkprefix, 0}, {"symmetry", 1, NULL, 'y'}, {"iterations", 1, NULL, 'n'}, {0, 0, NULL, 0} @@ -329,10 +185,6 @@ int main(int argc, char *argv[]) geomfile = strdup(optarg); break; - case 'x' : - prefix = strdup(optarg); - break; - case 'j' : nthreads = atoi(optarg); break; @@ -387,15 +239,6 @@ int main(int argc, char *argv[]) outfile = strdup("facetron.hkl"); } - /* Sanitise prefix */ - if ( prefix == NULL ) { - prefix = strdup(""); - } else { - if ( config_checkprefix ) { - prefix = check_prefix(prefix); - } - } - if ( sym == NULL ) sym = strdup("1"); /* Get detector geometry */ @@ -427,7 +270,11 @@ int main(int argc, char *argv[]) UnitCell *cell; char *filename; - char *fnamereal; + char *rval; + char line[1024]; + RefList *peaks; + RefList *transfer; + Reflection *refl; if ( find_chunk(fh, &cell, &filename) == 1 ) { ERROR("Couldn't get all of the filenames and cells" @@ -437,38 +284,83 @@ int main(int argc, char *argv[]) images[i].indexed_cell = cell; - /* Mangle the filename now */ - if ( config_basename ) { - char *tmp; - tmp = safe_basename(filename); - free(filename); - filename = tmp; - } - fnamereal = malloc(1024); - snprintf(fnamereal, 1023, "%s%s", prefix, filename); - free(filename); - - images[i].filename = fnamereal; + images[i].filename = filename; images[i].div = beam->divergence; images[i].bw = beam->bandwidth; images[i].det = det; images[i].beam = beam; images[i].osf = 1.0; images[i].profile_radius = 0.005e9; + images[i].reflections = reflist_new(); /* Muppet proofing */ images[i].data = NULL; images[i].flags = NULL; - /* Get reflections from this image. - * FIXME: Use the ones from the stream */ - integrate_image(&images[i], obs, sym); + /* Read integrated intensities from pattern */ + peaks = reflist_new(); + do { + + Reflection *refl; + signed int h, k, l; + float intensity; + int r; + + rval = fgets(line, 1023, fh); + chomp(line); + + if ( (strlen(line) == 0) || (rval == NULL) ) break; + + r = sscanf(line, "%i %i %i %f", &h, &k, &l, &intensity); + if ( r != 4 ) continue; + + /* Not interested in the central beam */ + if ( (h==0) && (k==0) && (l==0) ) continue; + + refl = add_refl(peaks, h, k, l); + set_int(refl, intensity); + + } while ( (strlen(line) != 0) && (rval != NULL) ); + + /* Calculate initial partialities and fill in intensities from + * the stream */ + transfer = find_intersections(&images[i], cell, 0); + images[i].reflections = reflist_new(); + + for ( refl = first_refl(transfer); + refl != NULL; + refl = next_refl(refl) ) { + + Reflection *peak; + Reflection *new; + signed int h, k, l, ha, ka, la; + double r1, r2, p, x, y; + int clamp1, clamp2; + + get_indices(refl, &h, &k, &l); + peak = find_refl(peaks, h, k, l); + if ( peak == NULL ) { + if ( (h==0) && (k==0) && (l==0) ) continue; + STATUS("%3i %3i %3i not found\n", h, k, l); + continue; + } + + get_asymm(h, k, l, &ha, &ka, &la, sym); + add_item(obs, ha, ka, la); + new = add_refl(images[i].reflections, ha, ka, la); + get_partial(refl, &r1, &r2, &p, &clamp1, &clamp2); + get_detector_pos(refl, &x, &y); + set_partial(new, r1, r2, p, clamp1, clamp2); + set_detector_pos(new, 0.0, x, y); + + } + reflist_free(peaks); + reflist_free(transfer); progress_bar(i, n_total_patterns-1, "Loading pattern data"); } fclose(fh); - free(prefix); cts = new_list_count(); diff --git a/src/post-refinement.c b/src/post-refinement.c index 25664e69..93b65810 100644 --- a/src/post-refinement.c +++ b/src/post-refinement.c @@ -28,6 +28,28 @@ #include "cell.h" +/* Maximum number of iterations of NLSq to do for each image per macrocycle. */ +#define MAX_CYCLES (100) + + +/* Refineable parameters */ +enum { + REF_SCALE, + REF_DIV, + REF_R, + REF_ASX, + REF_BSX, + REF_CSX, + REF_ASY, + REF_BSY, + REF_CSY, + REF_ASZ, + REF_BSZ, + REF_CSZ, + NUM_PARAMS +}; + + /* Returns dp/dr at "r" */ static double partiality_gradient(double r, double profile_radius) { @@ -65,7 +87,7 @@ static double partiality_rgradient(double r, double profile_radius) /* Return the gradient of parameter 'k' given the current status of 'image'. */ -double gradient(struct image *image, int k, Reflection *refl, double r) +static double gradient(struct image *image, int k, Reflection *refl, double r) { double ds, tt, azi; double nom, den; @@ -187,7 +209,7 @@ static void apply_cell_shift(UnitCell *cell, int k, double shift) /* Apply the given shift to the 'k'th parameter of 'image'. */ -void apply_shift(struct image *image, int k, double shift) +static void apply_shift(struct image *image, int k, double shift) { switch ( k ) { @@ -226,67 +248,19 @@ void apply_shift(struct image *image, int k, double shift) } -double mean_partial_dev(struct image *image, RefList *reflections, - const char *sym, double *i_full, FILE *graph) -{ - int n_used; - double delta_I = 0.0; - Reflection *refl; - - n_used = 0; - for ( refl = first_refl(reflections); - refl != NULL; - refl = next_refl(refl) ) { - - signed int hind, kind, lind; - signed int ha, ka, la; - double I_full; - float I_partial; - float xc, yc; - double x, y; - double p; - - get_indices(refl, &hind, &kind, &lind); - - if ( !get_scalable(refl) ) continue; - - /* Actual measurement of this reflection from this pattern? */ - get_detector_pos(refl, &x, &y); - /* FIXME: Get rid of this */ - if ( integrate_peak(image, x, y, - &xc, &yc, &I_partial, NULL, NULL, - 1, 0) ) { - continue; - } - - get_asymm(hind, kind, lind, &ha, &ka, &la, sym); - I_full = lookup_intensity(i_full, ha, ka, la); - p = get_partiality(refl); - delta_I += fabs(I_partial - (p * I_full / image->osf)); - n_used++; - - if ( graph != NULL ) { - fprintf(graph, "%3i %3i %3i %5.2f (at %5.2f,%5.2f)" - " %5.2f %5.2f\n", - hind, kind, lind, I_partial/p, xc, yc, - p, I_partial / I_full); - } - - } - - return delta_I / (double)n_used; -} - - /* Perform one cycle of post refinement on 'image' against 'i_full' */ -double pr_iterate(struct image *image, double *i_full, const char *sym, - RefList *reflections) +static double pr_iterate(struct image *image, const double *i_full, + const char *sym) { gsl_matrix *M; gsl_vector *v; gsl_vector *shifts; int param; Reflection *refl; + RefList *reflections; + double max_shift; + + reflections = image->reflections; M = gsl_matrix_calloc(NUM_PARAMS, NUM_PARAMS); v = gsl_vector_calloc(NUM_PARAMS); @@ -355,14 +329,31 @@ double pr_iterate(struct image *image, double *i_full, const char *sym, shifts = gsl_vector_alloc(NUM_PARAMS); gsl_linalg_HH_solve(M, v, shifts); + + max_shift = 0.0; for ( param=0; param max_shift ) max_shift = fabs(shift); } gsl_matrix_free(M); gsl_vector_free(v); gsl_vector_free(shifts); - return mean_partial_dev(image, reflections, sym, i_full, NULL); + return max_shift; +} + + +void pr_refine(struct image *image, const double *i_full, const char *sym) +{ + double max_shift; + int i; + + i = 0; + do { + max_shift = pr_iterate(image, i_full, sym); + STATUS("Iteration %2i: max shift = %5.2f\n", i, max_shift); + i++; + } while ( (max_shift > 0.01) && (i < MAX_CYCLES) ); } diff --git a/src/post-refinement.h b/src/post-refinement.h index b3c6c074..a10fafe4 100644 --- a/src/post-refinement.h +++ b/src/post-refinement.h @@ -3,7 +3,7 @@ * * Post refinement * - * (c) 2006-2010 Thomas White + * (c) 2006-2011 Thomas White * * Part of CrystFEL - crystallography with a FEL * @@ -21,36 +21,11 @@ #include #include "image.h" -#include "reflist.h" +#include "utils.h" -/* Refineable parameters */ -enum { - REF_SCALE, - REF_DIV, - REF_R, - REF_ASX, - REF_BSX, - REF_CSX, - REF_ASY, - REF_BSY, - REF_CSY, - REF_ASZ, - REF_BSZ, - REF_CSZ, - NUM_PARAMS -}; - -/* Apply the given shift to the 'k'th parameter of 'image'. */ -void apply_shift(struct image *image, int k, double shift); - - -double mean_partial_dev(struct image *image, RefList *reflections, - const char *sym, double *i_full, FILE *graph); - - -double pr_iterate(struct image *image, double *i_full, const char *sym, - RefList *reflections); +extern void pr_refine(struct image *image, const double *i_full, + const char *sym); #endif /* POST_REFINEMENT_H */ diff --git a/src/reflist.c b/src/reflist.c index d297c4d0..2bb4e5f6 100644 --- a/src/reflist.c +++ b/src/reflist.c @@ -19,7 +19,6 @@ struct _reflection { /* Listy stuff */ - RefList *list; unsigned int serial; /* Unique serial number, key */ struct _reflection *child[2]; /* Child nodes */ struct _reflection *parent; /* Parent node */ @@ -224,32 +223,6 @@ void set_partial(Reflection *refl, double r1, double r2, double p, } -void set_indices(Reflection *refl, signed int h, signed int k, signed int l) -{ - /* Tricky, because the indices determine the position in the tree */ - Reflection copy; - Reflection *new; - Reflection transfer; - - /* Copy all data */ - copy = *refl; - - /* Delete and re-add with new indices */ - delete_refl(refl); - new = add_refl(copy.list, h, k, l); - - /* Transfer data back */ - transfer = *new; - *new = copy; - new->list = transfer.list; - new->parent = transfer.parent; - new->child[0] = transfer.child[0]; - new->child[1] = transfer.child[1]; - new->h = transfer.h; new->k = transfer.k; new->l = transfer.l; - new->serial = transfer.serial; -} - - void set_int(Reflection *refl, double intensity) { refl->intensity = intensity; diff --git a/src/reflist.h b/src/reflist.h index c7b2ad42..78aadf2a 100644 --- a/src/reflist.h +++ b/src/reflist.h @@ -46,8 +46,6 @@ extern void set_detector_pos(Reflection *refl, double exerr, double x, double y); extern void set_partial(Reflection *refl, double r1, double r2, double p, double clamp_low, double clamp_high); -extern void set_indices(Reflection *refl, - signed int h, signed int k, signed int l); extern void set_int(Reflection *refl, double intensity); extern void set_scalable(Reflection *refl, int scalable); -- cgit v1.2.3