diff options
| -rw-r--r-- | libcrystfel/src/geometry.h | 4 | ||||
| -rw-r--r-- | src/partialator.c | 33 | ||||
| -rw-r--r-- | src/post-refinement.c | 74 | ||||
| -rw-r--r-- | src/post-refinement.h | 10 | ||||
| -rw-r--r-- | tests/pr_p_gradient_check.c | 2 |
5 files changed, 70 insertions, 53 deletions
diff --git a/libcrystfel/src/geometry.h b/libcrystfel/src/geometry.h index cecdc4f1..152c0e47 100644 --- a/libcrystfel/src/geometry.h +++ b/libcrystfel/src/geometry.h @@ -77,7 +77,9 @@ enum gparam { GPARAM_DIV, GPARAM_DETX, GPARAM_DETY, - GPARAM_CLEN + GPARAM_CLEN, + GPARAM_OSF, + GPARAM_BFAC }; diff --git a/src/partialator.c b/src/partialator.c index 9bdb684b..19449590 100644 --- a/src/partialator.c +++ b/src/partialator.c @@ -85,6 +85,7 @@ struct refine_args RefList *full; Crystal *crystal; PartialityModel pmodel; + int no_scale; struct prdata prdata; }; @@ -105,7 +106,8 @@ static void refine_image(void *task, int id) struct refine_args *pargs = task; Crystal *cr = pargs->crystal; - pargs->prdata = pr_refine(cr, pargs->full, pargs->pmodel); + pargs->prdata = pr_refine(cr, pargs->full, pargs->pmodel, + pargs->no_scale); } @@ -148,10 +150,6 @@ static void refine_all(Crystal **crystals, int n_crystals, struct refine_args task_defaults; struct queue_args qargs; - /* If the partiality model is "p=1", this refinement is really, really - * easy... */ - if ( pmodel == PMODEL_UNITY ) return; - task_defaults.full = full; task_defaults.crystal = NULL; task_defaults.pmodel = pmodel; @@ -641,15 +639,8 @@ int main(int argc, char *argv[]) // early_rejection(crystals, n_crystals); /* Make initial estimates */ - STATUS("Performing initial scaling.\n"); - if ( noscale ) { - STATUS("Skipping scaling step (--no-scale).\n"); - full = lsq_intensities(crystals, n_crystals, nthreads, pmodel, - min_measurements); - } else { - full = scale_intensities(crystals, n_crystals, nthreads, pmodel, - min_measurements); - } + full = lsq_intensities(crystals, n_crystals, nthreads, pmodel, + min_measurements); check_rejection(crystals, n_crystals); @@ -668,11 +659,11 @@ int main(int argc, char *argv[]) /* Iterate */ for ( i=0; i<n_iter; i++ ) { - STATUS("Post refinement cycle %i of %i\n", i+1, n_iter); + STATUS("Refinement cycle %i of %i\n", i+1, n_iter); srdata.n_filtered = 0; - /* Refine the geometry of all patterns to get the best fit */ + /* Refine all crystals to get the best fit */ refine_all(crystals, n_crystals, full, nthreads, pmodel, &srdata); @@ -681,14 +672,8 @@ int main(int argc, char *argv[]) /* Re-estimate all the full intensities */ reflist_free(full); - if ( noscale ) { - STATUS("Skipping scaling step (--no-scale).\n"); - full = lsq_intensities(crystals, n_crystals, nthreads, - pmodel, min_measurements); - } else { - full = scale_intensities(crystals, n_crystals, nthreads, - pmodel, min_measurements); - } + full = lsq_intensities(crystals, n_crystals, nthreads, + pmodel, min_measurements); check_rejection(crystals, n_crystals); diff --git a/src/post-refinement.c b/src/post-refinement.c index 8e6f6108..5a1053b8 100644 --- a/src/post-refinement.c +++ b/src/post-refinement.c @@ -174,14 +174,16 @@ static double volume_fraction(double rlow, double rhigh, double pr, } -/* Return the gradient of partiality wrt parameter 'k' given the current status - * of 'image'. */ -double p_gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel) +/* Return the gradient of p/G wrt parameter 'k' given the current status + * of the crystal. */ +double gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel) { double glow, ghigh; double rlow, rhigh, p; struct image *image = crystal_get_image(cr); double R = crystal_get_profile_radius(cr); + double osf = crystal_get_osf(cr); + double gr; get_partial(refl, &rlow, &rhigh, &p); @@ -196,7 +198,8 @@ double p_gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel) Rglow = volume_fraction_rgradient(rlow, R, pmodel); Rghigh = volume_fraction_rgradient(rhigh, R, pmodel); - return 4.0*psph/(3.0*D) + (4.0*R/(3.0*D))*(Rglow - Rghigh); + gr = 4.0*psph/(3.0*D) + (4.0*R/(3.0*D))*(Rglow - Rghigh); + return gr / osf; } @@ -214,11 +217,17 @@ double p_gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel) get_symmetric_indices(refl, &hs, &ks, &ls); ds = 2.0 * resolution(crystal_get_cell(cr), hs, ks, ls); - return (ds/2.0)*(glow+ghigh) - 4.0*R*psph*ds/(3.0*D*D); + gr = (ds/2.0)*(glow+ghigh) - 4.0*R*psph*ds/(3.0*D*D); + return gr / osf; } - return r_gradient(crystal_get_cell(cr), k, refl, image) * (glow-ghigh); + if ( k == GPARAM_OSF ) { + return -p/(osf*osf); + } + + gr = r_gradient(crystal_get_cell(cr), k, refl, image) * (glow-ghigh); + return gr / osf; } @@ -274,6 +283,18 @@ static void apply_shift(Crystal *cr, int k, double shift) crystal_set_profile_radius(cr, t); break; + case GPARAM_OSF : + if ( isnan(shift) ) { + ERROR("Refusing nan shift of OSF\n"); + } else if (shift > 100.0 ) { + ERROR("Refusing large (%e) OSF shift.\n", shift); + } else { + t = crystal_get_osf(cr); + t += shift; + crystal_set_osf(cr, t); + } + break; + case GPARAM_ASX : case GPARAM_ASY : case GPARAM_ASZ : @@ -296,7 +317,7 @@ static void apply_shift(Crystal *cr, int k, double shift) /* Perform one cycle of post refinement on 'image' against 'full' */ static double pr_iterate(Crystal *cr, const RefList *full, - PartialityModel pmodel, int *n_filtered) + PartialityModel pmodel, int no_scale, int *n_filtered) { gsl_matrix *M; gsl_vector *v; @@ -310,6 +331,7 @@ static double pr_iterate(Crystal *cr, const RefList *full, const int verbose = 0; int num_params = 0; enum gparam rv[32]; + double G; *n_filtered = 0; @@ -327,13 +349,19 @@ static double pr_iterate(Crystal *cr, const RefList *full, rv[num_params++] = GPARAM_CSZ; } - STATUS("Refining %i parameters.\n", num_params); + /* If we are scaling, refine scale factors (duh) */ + if ( !no_scale ) { + rv[num_params++] = GPARAM_OSF; + //rv[num_params++] = GPARAM_BFAC; FIXME: Next + } reflections = crystal_get_reflections(cr); M = gsl_matrix_calloc(num_params, num_params); v = gsl_vector_calloc(num_params); + G = crystal_get_osf(cr); + /* Construct the equations, one per reflection in this image */ for ( refl = first_refl(reflections, &iter); refl != NULL; @@ -344,7 +372,7 @@ static double pr_iterate(Crystal *cr, const RefList *full, double w; double I_partial; int k; - double p, l; + double p, L; Reflection *match; double gradients[num_params]; @@ -359,19 +387,22 @@ static double pr_iterate(Crystal *cr, const RefList *full, I_full = get_intensity(match); - /* Actual measurement of this reflection from this pattern? */ - I_partial = get_intensity(refl) / crystal_get_osf(cr); p = get_partiality(refl); - l = get_lorentz(refl); + L = get_lorentz(refl); + + /* Actual measurement of this reflection from this pattern? */ + I_partial = get_intensity(refl); /* Calculate the weight for this reflection */ - w = pow(get_esd_intensity(refl), 2.0); - w += l * p * I_full * pow(get_esd_intensity(match), 2.0); - w = pow(w, -1.0); + //w = pow(get_esd_intensity(refl), 2.0); + //w += (p/L) * I_full * pow(get_esd_intensity(match), 2.0); + //w = pow(w, -1.0); + w = 1.0; /* Calculate all gradients for this reflection */ for ( k=0; k<num_params; k++ ) { - gradients[k] = p_gradient(cr, rv[k], refl, pmodel) * l; + gradients[k] = gradient(cr, rv[k], refl, pmodel); + gradients[k] *= I_full / L; } for ( k=0; k<num_params; k++ ) { @@ -386,8 +417,7 @@ static double pr_iterate(Crystal *cr, const RefList *full, /* Matrix is symmetric */ if ( g > k ) continue; - M_c = gradients[g] * gradients[k]; - M_c *= w * pow(I_full, 2.0); + M_c = w * gradients[g] * gradients[k]; M_curr = gsl_matrix_get(M, k, g); gsl_matrix_set(M, k, g, M_curr + M_c); @@ -395,8 +425,8 @@ static double pr_iterate(Crystal *cr, const RefList *full, } - delta_I = I_partial - (l * p * I_full); - v_c = w * delta_I * I_full * gradients[k]; + delta_I = I_partial - (p * I_full / (L*G)); + v_c = w * delta_I * gradients[k]; v_curr = gsl_vector_get(v, k); gsl_vector_set(v, k, v_curr + v_c); @@ -486,7 +516,7 @@ static double guide_dev(Crystal *cr, const RefList *full) struct prdata pr_refine(Crystal *cr, const RefList *full, - PartialityModel pmodel) + PartialityModel pmodel, int no_scale) { double dev; int i; @@ -518,7 +548,7 @@ struct prdata pr_refine(Crystal *cr, const RefList *full, cell_get_reciprocal(crystal_get_cell(cr), &asx, &asy, &asz, &bsx, &bsy, &bsz, &csx, &csy, &csz); - pr_iterate(cr, full, pmodel, &prdata.n_filtered); + pr_iterate(cr, full, pmodel, no_scale, &prdata.n_filtered); update_partialities(cr, pmodel); diff --git a/src/post-refinement.h b/src/post-refinement.h index caf5f9af..844afbd1 100644 --- a/src/post-refinement.h +++ b/src/post-refinement.h @@ -3,11 +3,11 @@ * * Post refinement * - * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * * Authors: - * 2010-2014 Thomas White <taw@physics.org> + * 2010-2015 Thomas White <taw@physics.org> * * This file is part of CrystFEL. * @@ -51,10 +51,10 @@ struct prdata extern struct prdata pr_refine(Crystal *cr, const RefList *full, - PartialityModel pmodel); + PartialityModel pmodel, int no_scale); /* Exported so it can be poked by tests/pr_p_gradient_check */ -extern double p_gradient(Crystal *cr, int k, Reflection *refl, - PartialityModel pmodel); +extern double gradient(Crystal *cr, int k, Reflection *refl, + PartialityModel pmodel); #endif /* POST_REFINEMENT_H */ diff --git a/tests/pr_p_gradient_check.c b/tests/pr_p_gradient_check.c index affddeab..1ae7b049 100644 --- a/tests/pr_p_gradient_check.c +++ b/tests/pr_p_gradient_check.c @@ -303,7 +303,7 @@ static double test_gradients(Crystal *cr, double incr_val, int refine, grad = (grad1 + grad2) / 2.0; i++; - cgrad = p_gradient(cr, refine, refl, pmodel); + cgrad = gradient(cr, refine, refl, pmodel); get_partial(refl, &r1, &r2, &p); |