From d813960912ddc76d5cff0501b4f7783c3eeca95d Mon Sep 17 00:00:00 2001 From: Thomas White Date: Tue, 8 Nov 2016 16:20:58 +0100 Subject: Separate partiality calculation from prediction --- libcrystfel/src/geometry.c | 130 ++++++++++++++++++++++++-------------- libcrystfel/src/geometry.h | 14 ++-- libcrystfel/src/integration.c | 3 +- libcrystfel/src/predict-refine.c | 8 +-- src/partial_sim.c | 3 +- src/partialator.c | 3 +- src/post-refinement.c | 9 ++- tests/pr_p_gradient_check.c | 20 ++++-- tests/prediction_gradient_check.c | 7 +- tests/prof2d_check.c | 6 +- 10 files changed, 126 insertions(+), 77 deletions(-) diff --git a/libcrystfel/src/geometry.c b/libcrystfel/src/geometry.c index b968dc06..576655de 100644 --- a/libcrystfel/src/geometry.c +++ b/libcrystfel/src/geometry.c @@ -266,15 +266,12 @@ static double partiality(PartialityModel pmodel, static Reflection *check_reflection(struct image *image, Crystal *cryst, - PartialityModel pmodel, signed int h, signed int k, signed int l, double xl, double yl, double zl, Reflection *updateme) { - const int output = 0; double tl; - double rlow, rhigh; /* "Excitation error" */ - double part; /* Partiality */ + double rlow, rhigh; /* "Excitation error" */ double klow, khigh; /* Wavenumber */ Reflection *refl; double cet, cez; /* Centre of Ewald sphere */ @@ -312,20 +309,6 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst, && (fabs(rlow) > pr) && (fabs(rhigh) > pr) ) return NULL; - /* Calculate partiality */ - part = partiality(pmodel, h, k, l, image->serial, rlow, rhigh, pr); - - if ( isnan(part) && ((h!=0) || (k!=0) || (l!=0)) ) { - ERROR("Assigning NAN partiality!\n"); - ERROR("%3i %3i %3i rlow = %e, rhigh = %e\n", - h, k, l, rlow, rhigh); - ERROR("div + m = %e, R = %e, bw = %e\n", del, pr, image->bw); - ERROR("D = %e\n", rlow - rhigh); - ERROR("klow = %e khigh = %e\n", klow, khigh); - ERROR("distance(..) = %e\n", distance(cet, cez, tl, zl)); - abort(); - } - if ( updateme == NULL ) { refl = reflection_new(h, k, l); } else { @@ -373,15 +356,12 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst, } } - set_partial(refl, rlow, rhigh, part); + set_partial(refl, rlow, rhigh, 1.0); /* Actual partiality set by + * calculate_partialities() */ set_lorentz(refl, 1.0); set_symmetric_indices(refl, h, k, l); set_redundancy(refl, 1); - if ( output ) { - printf("%3i %3i %3i %6f %5.2f\n", h, k, l, 0.0, part); - } - return refl; } @@ -472,15 +452,7 @@ double r_gradient(UnitCell *cell, int k, Reflection *refl, struct image *image) } -RefList *find_intersections(struct image *image, Crystal *cryst, - PartialityModel pmodel) -{ - return find_intersections_to_res(image, cryst, pmodel, INFINITY); -} - - -RefList *find_intersections_to_res(struct image *image, Crystal *cryst, - PartialityModel pmodel, double max_res) +RefList *predict_to_res(Crystal *cryst, double max_res) { double ax, ay, az; double bx, by, bz; @@ -509,7 +481,7 @@ RefList *find_intersections_to_res(struct image *image, Crystal *cryst, cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz); - mres = largest_q(image); + mres = largest_q(crystal_get_image(cryst)); if ( mres > max_res ) mres = max_res; hmax = mres * modulus(ax, ay, az); @@ -544,7 +516,7 @@ RefList *find_intersections_to_res(struct image *image, Crystal *cryst, yl = h*asy + k*bsy + l*csy; zl = h*asz + k*bsz + l*csz; - refl = check_reflection(image, cryst, pmodel, + refl = check_reflection(crystal_get_image(cryst), cryst, h, k, l, xl, yl, zl, NULL); if ( refl != NULL ) { @@ -559,12 +531,12 @@ RefList *find_intersections_to_res(struct image *image, Crystal *cryst, } -static void set_unity_partialities(Crystal *cryst) +static void set_unity_partialities(RefList *list) { Reflection *refl; RefListIterator *iter; - for ( refl = first_refl(crystal_get_reflections(cryst), &iter); + for ( refl = first_refl(list, &iter); refl != NULL; refl = next_refl(refl, iter) ) { @@ -575,14 +547,85 @@ static void set_unity_partialities(Crystal *cryst) /** - * update_partialities: + * calculate_partialities: * @cryst: A %Crystal * @pmodel: A %PartialityModel * - * Updates the partialities of @cryst's reflections according to @pmodel and + * Calculates the partialities for the reflections in @cryst, given the current + * crystal and image parameters. The crystal's image and reflection lists + * must be set. The specified %PartialityModel will be used. + * + * You must not have changed the crystal or image parameters since you last + * called predict_to_res() or update_predictions(), because this function + * relies on the limiting wavelength values calculated by those functions. + */ +void calculate_partialities(Crystal *cryst, PartialityModel pmodel) +{ + RefList *list; + Reflection *refl; + RefListIterator *iter; + double pr; + struct image *image; + + list = crystal_get_reflections(cryst); + if ( list == NULL ) { + ERROR("No reflections for partiality calculation!\n"); + return; + } + + if ( pmodel == PMODEL_UNITY ) { + set_unity_partialities(list); + return; + } + + image = crystal_get_image(cryst); + if ( image == NULL ) { + ERROR("No image for partiality calculation!\n"); + return; + } + + pr = crystal_get_profile_radius(cryst); + + for ( refl = first_refl(crystal_get_reflections(cryst), &iter); + refl != NULL; + refl = next_refl(refl, iter) ) + { + double rlow, rhigh, part; + signed int h, k, l; + + get_symmetric_indices(refl, &h, &k, &l); + get_partial(refl, &rlow, &rhigh, &part); + + part = partiality(pmodel, h, k, l, image->serial, + rlow, rhigh, pr); + + if ( isnan(part) && ((h!=0) || (k!=0) || (l!=0)) ) { + ERROR("Assigning NAN partiality!\n"); + ERROR("%3i %3i %3i rlow = %e, rhigh = %e\n", + h, k, l, rlow, rhigh); + ERROR("R = %e, bw = %e\n", pr, image->bw); + ERROR("D = %e\n", rlow - rhigh); + abort(); + } + + set_partiality(refl, part); + + } +} + + +/** + * update_predictions: + * @cryst: A %Crystal + * + * Updates the predicted reflections (positions and excitation errors, but not + * the actual partialities) of @cryst's reflections according to * the current state of the crystal (e.g. its unit cell parameters). + * + * If you need to update the partialities as well, call calculate_partialities() + * afterwards. */ -void update_partialities(Crystal *cryst, PartialityModel pmodel) +void update_predictions(Crystal *cryst) { Reflection *refl; RefListIterator *iter; @@ -591,11 +634,6 @@ void update_partialities(Crystal *cryst, PartialityModel pmodel) double csx, csy, csz; struct image *image = crystal_get_image(cryst); - if ( pmodel == PMODEL_UNITY ) { - set_unity_partialities(cryst); - return; - } - cell_get_reciprocal(crystal_get_cell(cryst), &asx, &asy, &asz, &bsx, &bsy, &bsz, &csx, &csy, &csz); @@ -613,8 +651,8 @@ void update_partialities(Crystal *cryst, PartialityModel pmodel) yl = h*asy + k*bsy + l*csy; zl = h*asz + k*bsz + l*csz; - check_reflection(image, cryst, pmodel, - h, k, l, xl, yl, zl, refl); + check_reflection(image, cryst, h, k, l, + xl, yl, zl, refl); } } diff --git a/libcrystfel/src/geometry.h b/libcrystfel/src/geometry.h index 6b035175..4fbd0bca 100644 --- a/libcrystfel/src/geometry.h +++ b/libcrystfel/src/geometry.h @@ -3,12 +3,12 @@ * * Geometry of diffraction * - * Copyright © 2013-2015 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2013-2016 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * Copyright © 2012 Richard Kirian * * Authors: - * 2010-2015 Thomas White + * 2010-2016 Thomas White * 2012 Richard Kirian * * This file is part of CrystFEL. @@ -85,15 +85,13 @@ enum gparam { }; -extern RefList *find_intersections(struct image *image, Crystal *cryst, - PartialityModel pmodel); -extern RefList *find_intersections_to_res(struct image *image, Crystal *cryst, - PartialityModel pmodel, - double max_res); +extern RefList *predict_to_res(Crystal *cryst, double max_res); + +extern void calculate_partialities(Crystal *cryst, PartialityModel pmodel); extern double r_gradient(UnitCell *cell, int k, Reflection *refl, struct image *image); -extern void update_partialities(Crystal *cryst, PartialityModel pmodel); +extern void update_predictions(Crystal *cryst); extern void polarisation_correction(RefList *list, UnitCell *cell, struct image *image); diff --git a/libcrystfel/src/integration.c b/libcrystfel/src/integration.c index d50f8d6b..96eec92a 100644 --- a/libcrystfel/src/integration.c +++ b/libcrystfel/src/integration.c @@ -1673,8 +1673,7 @@ void integrate_all_4(struct image *image, IntegrationMethod meth, image->features); crystal_set_resolution_limit(image->crystals[i], res); - list = find_intersections_to_res(image, image->crystals[i], - pmodel, res+push_res); + list = predict_to_res(image->crystals[i], res+push_res); crystal_set_reflections(image->crystals[i], list); } diff --git a/libcrystfel/src/predict-refine.c b/libcrystfel/src/predict-refine.c index e702cf33..3246dbec 100644 --- a/libcrystfel/src/predict-refine.c +++ b/libcrystfel/src/predict-refine.c @@ -262,7 +262,7 @@ static int pair_peaks(struct image *image, Crystal *cr, * doesn't fall on this panel. We're only interested * in how far away it is from the peak location. * The predicted position and excitation errors will be - * filled in by update_partialities(). */ + * filled in by update_predictions(). */ set_panel(refl, f->p); rps[n].refl = refl; @@ -275,7 +275,7 @@ static int pair_peaks(struct image *image, Crystal *cr, /* Get the excitation errors and detector positions for the candidate * reflections */ crystal_set_reflections(cr, all_reflist); - update_partialities(cr, PMODEL_SCSPHERE); + update_predictions(cr); /* Pass over the peaks again, keeping only the ones which look like * good pairings */ @@ -339,7 +339,7 @@ void refine_radius(Crystal *cr, struct image *image) return; } crystal_set_reflections(cr, reflist); - update_partialities(cr, PMODEL_SCSPHERE); + update_predictions(cr); crystal_set_reflections(cr, NULL); qsort(rps, n_acc, sizeof(struct reflpeak), cmpd2); @@ -635,7 +635,7 @@ int refine_prediction(struct image *image, Crystal *cr) /* Refine */ for ( i=0; iimage.filename, 255, "dummy.h5"); } - reflections = find_intersections(&wargs->image, cr, PMODEL_SCSPHERE); + reflections = predict_to_res(cr, largest_q(&wargs->image)); crystal_set_reflections(cr, reflections); + calculate_partialities(cr, PMODEL_SCSPHERE); for ( i=0; in_ref[i] = 0; diff --git a/src/partialator.c b/src/partialator.c index aded835d..428b0d68 100644 --- a/src/partialator.c +++ b/src/partialator.c @@ -1135,7 +1135,8 @@ int main(int argc, char *argv[]) crystal_set_image(cryst, &images[i]); /* Now it's safe to do the following */ - update_partialities(cryst, pmodel); + update_predictions(cryst); + calculate_partialities(cryst, pmodel); } } diff --git a/src/post-refinement.c b/src/post-refinement.c index 487ce872..bf30d299 100644 --- a/src/post-refinement.c +++ b/src/post-refinement.c @@ -606,7 +606,8 @@ static void write_residual_graph(Crystal *cr, const RefList *full) cell_set_reciprocal(cell, asx, asy, asz, bsx, bsy, bsz, csx, csy, csz); - update_partialities(cr, PMODEL_SCSPHERE); + update_predictions(cr); + calculate_partialities(cr, PMODEL_SCSPHERE); res = residual(cr, full, 0, &n, NULL); fprintf(fh, "%i %e %e %i\n", i, asx, res, n); } @@ -614,7 +615,8 @@ static void write_residual_graph(Crystal *cr, const RefList *full) cell_set_reciprocal(cell, orig_asx, asy, asz, bsx, bsy, bsz, csx, csy, csz); - update_partialities(cr, PMODEL_SCSPHERE); + update_predictions(cr); + calculate_partialities(cr, PMODEL_SCSPHERE); fclose(fh); } @@ -648,7 +650,8 @@ static void do_pr_refine(Crystal *cr, const RefList *full, pr_iterate(cr, full, pmodel, NULL, verbose); - update_partialities(cr, pmodel); + update_predictions(cr); + calculate_partialities(cr, pmodel); dev = residual(cr, full, 0, 0, NULL); if ( fabs(dev - old_dev) < dev*0.0001 ) done = 1; diff --git a/tests/pr_p_gradient_check.c b/tests/pr_p_gradient_check.c index 9c2a6633..5322fcca 100644 --- a/tests/pr_p_gradient_check.c +++ b/tests/pr_p_gradient_check.c @@ -178,7 +178,9 @@ static void calc_either_side(Crystal *cr, double incr_val, /* Crystal properties */ cr_new = new_shifted_crystal(cr, refine, -incr_val); - compare = find_intersections(image, cr_new, pmodel); + compare = predict_to_res(cr_new, largest_q(image)); + crystal_set_reflections(cr_new, compare); + calculate_partialities(cr_new, pmodel); scan_partialities(crystal_get_reflections(cr), compare, valid, vals, 0, pmodel); cell_free(crystal_get_cell(cr_new)); @@ -186,7 +188,9 @@ static void calc_either_side(Crystal *cr, double incr_val, reflist_free(compare); cr_new = new_shifted_crystal(cr, refine, +incr_val); - compare = find_intersections(image, cr_new, pmodel); + compare = predict_to_res(cr_new, largest_q(image)); + crystal_set_reflections(cr_new, compare); + calculate_partialities(cr_new, pmodel); scan_partialities(crystal_get_reflections(cr), compare, valid, vals, 2, pmodel); cell_free(crystal_get_cell(cr_new)); @@ -196,18 +200,24 @@ static void calc_either_side(Crystal *cr, double incr_val, } else { struct image im_moved; + Crystal *cr_new = crystal_copy(cr); + crystal_set_image(cr_new, &im_moved); /* "Image" properties */ im_moved = *image; shift_parameter(&im_moved, refine, -incr_val); - compare = find_intersections(&im_moved, cr, pmodel); + compare = predict_to_res(cr_new, largest_q(&im_moved)); + crystal_set_reflections(cr_new, compare); + calculate_partialities(cr_new, pmodel); scan_partialities(crystal_get_reflections(cr), compare, valid, vals, 0, pmodel); reflist_free(compare); im_moved = *image; shift_parameter(&im_moved, refine, +incr_val); - compare = find_intersections(&im_moved, cr, pmodel); + compare = predict_to_res(cr_new, largest_q(&im_moved)); + crystal_set_reflections(cr_new, compare); + calculate_partialities(cr_new, pmodel); scan_partialities(crystal_get_reflections(cr), compare, valid, vals, 2, pmodel); reflist_free(compare); @@ -237,7 +247,7 @@ static double test_gradients(Crystal *cr, double incr_val, int refine, int n_line; double cc; - reflections = find_intersections(crystal_get_image(cr), cr, pmodel); + reflections = predict_to_res(cr, largest_q(crystal_get_image(cr))); crystal_set_reflections(cr, reflections); nref = num_reflections(reflections); diff --git a/tests/prediction_gradient_check.c b/tests/prediction_gradient_check.c index b9f23217..0ffc07ca 100644 --- a/tests/prediction_gradient_check.c +++ b/tests/prediction_gradient_check.c @@ -185,14 +185,14 @@ static void calc_either_side(Crystal *cr, double incr_val, Crystal *cr_new; cr_new = new_shifted_crystal(cr, refine, -incr_val); - compare = find_intersections(image, cr_new, PMODEL_SCSPHERE); + compare = predict_to_res(cr_new, largest_q(image)); scan(crystal_get_reflections(cr), compare, valid, vals, 0); cell_free(crystal_get_cell(cr_new)); crystal_free(cr_new); reflist_free(compare); cr_new = new_shifted_crystal(cr, refine, +incr_val); - compare = find_intersections(image, cr_new, PMODEL_SCSPHERE); + compare = predict_to_res(cr_new, largest_q(image)); scan(crystal_get_reflections(cr), compare, valid, vals, 2); cell_free(crystal_get_cell(cr_new)); crystal_free(cr_new); @@ -221,8 +221,7 @@ static double test_gradients(Crystal *cr, double incr_val, int refine, int n_line; double cc; - reflections = find_intersections(crystal_get_image(cr), cr, - PMODEL_SCSPHERE); + reflections = predict_to_res(cr, largest_q(crystal_get_image(cr))); crystal_set_reflections(cr, reflections); nref = num_reflections(reflections); diff --git a/tests/prof2d_check.c b/tests/prof2d_check.c index 9fe8b356..8b9de871 100644 --- a/tests/prof2d_check.c +++ b/tests/prof2d_check.c @@ -3,11 +3,11 @@ * * Check 2D profile fitting * - * Copyright © 2013-2014 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2013-2016 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * * Authors: - * 2013-2014 Thomas White + * 2013-2014,2016 Thomas White * * This file is part of CrystFEL. * @@ -135,7 +135,7 @@ int main(int argc, char *argv[]) image.n_crystals = 1; image.crystals = &cr; - list = find_intersections(&image, cr, PMODEL_SCSPHERE); + list = predict_to_res(cr, largest_q(&image)); crystal_set_reflections(cr, list); for ( fs=0; fs