From 4eb0015369115791fc721ca97736d30fc38d8ae1 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Tue, 27 Feb 2018 15:58:00 +0100 Subject: Switch to simpler scaling algorithm --- src/partialator.c | 69 +++------- src/post-refinement.c | 8 ++ src/scaling.c | 346 +++++--------------------------------------------- src/scaling.h | 5 +- 4 files changed, 61 insertions(+), 367 deletions(-) diff --git a/src/partialator.c b/src/partialator.c index b5f24b1d..7c3a481c 100644 --- a/src/partialator.c +++ b/src/partialator.c @@ -1221,22 +1221,21 @@ int main(int argc, char *argv[]) STATUS("Checking patterns.\n"); //early_rejection(crystals, n_crystals); - /* Initial rejection, figures of merit etc */ + /* Create reference data set if we don't already have one */ if ( reference == NULL ) { - if ( !no_scale ) { - scale_all(crystals, n_crystals, nthreads, pmodel); - } full = merge_intensities(crystals, n_crystals, nthreads, pmodel, min_measurements, push_res, 1); } else { full = reference; } - check_rejection(crystals, n_crystals, full, max_B); + /* Scale everything to the reference */ if ( !no_scale ) { - scale_all_to_reference(crystals, n_crystals, full); + scale_all_to_reference(crystals, n_crystals, full, nthreads); } + /* Check rejection and write figures of merit */ + check_rejection(crystals, n_crystals, full, max_B); show_all_residuals(crystals, n_crystals, full); write_pgraph(full, crystals, n_crystals, 0); write_specgraph(full, crystals[0], 0); @@ -1248,41 +1247,19 @@ int main(int argc, char *argv[]) if ( !no_pr ) { refine_all(crystals, n_crystals, full, nthreads, pmodel); + } else if ( !no_scale ) { + scale_all_to_reference(crystals, n_crystals, full, nthreads); } - if ( !no_scale ) { - - if ( reference == NULL ) { - - /* No reference -> XSCALE-like algorithm */ - full = merge_intensities(crystals, n_crystals, - nthreads, pmodel, - min_measurements, - push_res, 1); - - } else { - - /* Reference -> Ginn-style linear algorithm */ - if ( !no_scale ) { - scale_all_to_reference(crystals, n_crystals, - reference); - } - - } - } - - check_rejection(crystals, n_crystals, full, max_B); - + /* Create new reference if needed */ if ( reference == NULL ) { - if ( !no_scale ) { - scale_all(crystals, n_crystals, nthreads, pmodel); - } reflist_free(full); full = merge_intensities(crystals, n_crystals, nthreads, pmodel, min_measurements, push_res, 1); } /* else full still equals reference */ + check_rejection(crystals, n_crystals, full, max_B); show_all_residuals(crystals, n_crystals, full); write_pgraph(full, crystals, n_crystals, i+1); write_specgraph(full, crystals[0], i+1); @@ -1319,34 +1296,24 @@ int main(int argc, char *argv[]) } } - + /* Final merge */ + STATUS("Final merge...\n"); if ( reference == NULL ) { - if ( !no_scale ) { - scale_all(crystals, n_crystals, nthreads, pmodel); - } reflist_free(full); full = merge_intensities(crystals, n_crystals, nthreads, pmodel, min_measurements, push_res, 1); - } /* else full still equals reference */ - if ( !no_scale ) { - scale_all_to_reference(crystals, n_crystals, full); + } else { + scale_all_to_reference(crystals, n_crystals, reference, nthreads); + full = merge_intensities(crystals, n_crystals, nthreads, + pmodel, min_measurements, push_res, 1); } + + /* Write final figures of merit (no rejection any more) */ show_all_residuals(crystals, n_crystals, full); write_pgraph(full, crystals, n_crystals, n_iter+1); write_specgraph(full, crystals[0], n_iter+1); - //STATUS("Final profile radius: %e\n", crystal_get_profile_radius(crystals[0])); - - /* If we've been using a reference up to now, it's time to actually - * scale and merge the final version */ - if ( reference != NULL ) { - STATUS("Starting final reference-free merge\n"); - scale_all_to_reference(crystals, n_crystals, reference); - STATUS("Scaling done\n"); - full = merge_intensities(crystals, n_crystals, nthreads, - pmodel, min_measurements, push_res, 1); - STATUS("Done\n"); - } /* else we just did it */ + STATUS("Final profile radius: %e\n", crystal_get_profile_radius(crystals[0])); /* Output results */ STATUS("Writing overall results to %s\n", outfile); diff --git a/src/post-refinement.c b/src/post-refinement.c index af701575..7402b91c 100644 --- a/src/post-refinement.c +++ b/src/post-refinement.c @@ -381,6 +381,8 @@ static void do_pr_refine(Crystal *cr, const RefList *full, int n_params = 0; int n_iter = 0; int status; + int r; + double G; double residual_start, residual_free_start; residual_start = residual(cr, full, 0, NULL, NULL); @@ -456,6 +458,12 @@ static void do_pr_refine(Crystal *cr, const RefList *full, apply_parameters(min->x, initial, rv, cr); update_predictions(cr); calculate_partialities(cr, PMODEL_XSPHERE); + r = linear_scale(full, crystal_get_reflections(cr), &G); + if ( r == 0 ) { + crystal_set_osf(cr, G); + } else { + fprintf(stderr, "Scaling failure after refinement.\n"); + } if ( verbose ) { STATUS("PR final: dev = %10.5e, free dev = %10.5e\n", diff --git a/src/scaling.c b/src/scaling.c index b2f063ec..bf8857ac 100644 --- a/src/scaling.c +++ b/src/scaling.c @@ -48,184 +48,6 @@ #include "scaling.h" -/* Maximum number of iterations of NLSq to do for each image per macrocycle. */ -#define MAX_CYCLES (10) - - -/* Apply the given shift to the 'k'th parameter of 'image'. */ -static void apply_shift(Crystal *cr, int k, double shift) -{ - double t; - - switch ( k ) { - - case GPARAM_BFAC : - t = crystal_get_Bfac(cr); - t += shift; - crystal_set_Bfac(cr, t); - break; - - case GPARAM_OSF : - t = -log(crystal_get_osf(cr)); - t += shift; - crystal_set_osf(cr, exp(-t)); - break; - - default : - ERROR("No shift defined for parameter %i\n", k); - abort(); - - } -} - - -/* Perform one cycle of scaling of 'cr' against 'full' */ -static double scale_iterate(Crystal *cr, const RefList *full, - PartialityModel pmodel, int *nr) -{ - gsl_matrix *M; - gsl_vector *v; - gsl_vector *shifts; - int param; - Reflection *refl; - RefListIterator *iter; - RefList *reflections; - double max_shift; - int nref = 0; - int num_params = 0; - enum gparam rv[32]; - double G, B; - - *nr = 0; - - rv[num_params++] = GPARAM_OSF; - rv[num_params++] = GPARAM_BFAC; - - M = gsl_matrix_calloc(num_params, num_params); - v = gsl_vector_calloc(num_params); - - reflections = crystal_get_reflections(cr); - G = crystal_get_osf(cr); - B = crystal_get_Bfac(cr); - - /* Scaling terms */ - for ( refl = first_refl(reflections, &iter); - refl != NULL; - refl = next_refl(refl, iter) ) - { - signed int ha, ka, la; - double I_full, delta_I, esd; - double w; - double I_partial; - int k; - double p, L, s; - double fx; - Reflection *match; - double gradients[num_params]; - - /* If reflection is free-flagged, don't use it here */ - if ( get_flag(refl) ) continue; - - /* Find the full version */ - get_indices(refl, &ha, &ka, &la); - match = find_refl(full, ha, ka, la); - if ( match == NULL ) continue; - - /* Merged intensitty */ - I_full = get_intensity(match); - - /* Actual measurement of this reflection from this pattern */ - I_partial = get_intensity(refl); - esd = get_esd_intensity(refl); - p = get_partiality(refl); - - /* Scale only using strong reflections */ - if ( I_partial <= 3.0*esd ) continue; /* Also because of log */ - if ( get_redundancy(match) < 2 ) continue; - if ( I_full <= 0 ) continue; /* Because log */ - if ( p <= 0.0 ) continue; /* Because of log */ - - L = get_lorentz(refl); - s = resolution(crystal_get_cell(cr), ha, ka, la); - - /* Calculate the weight for this reflection */ - w = 1.0; - - /* Calculate all gradients for this reflection */ - for ( k=0; k k ) continue; - - 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); - gsl_matrix_set(M, g, k, M_curr + M_c); - - } - - fx = -log(G) + log(p) - log(L) - B*s*s + log(I_full); - delta_I = log(I_partial) - fx; - v_c = w * delta_I * gradients[k]; - v_curr = gsl_vector_get(v, k); - gsl_vector_set(v, k, v_curr + v_c); - - } - - nref++; - } - - *nr = nref; - - if ( nref < num_params ) { - crystal_set_user_flag(cr, PRFLAG_FEWREFL); - gsl_matrix_free(M); - gsl_vector_free(v); - return 0.0; - } - - max_shift = 0.0; - shifts = solve_svd(v, M, NULL, 0); - if ( shifts != NULL ) { - - for ( param=0; param max_shift ) max_shift = fabs(shift); - } - - } else { - crystal_set_user_flag(cr, PRFLAG_SOLVEFAIL); - } - - gsl_matrix_free(M); - gsl_vector_free(v); - gsl_vector_free(shifts); - - return max_shift; -} - - double log_residual(Crystal *cr, const RefList *full, int free, int *pn_used, const char *filename) { @@ -292,32 +114,6 @@ double log_residual(Crystal *cr, const RefList *full, int free, } -static void do_scale_refine(Crystal *cr, const RefList *full, - PartialityModel pmodel, int *nr) -{ - int i, done; - double old_dev; - - old_dev = log_residual(cr, full, 0, NULL, NULL); - - i = 0; - done = 0; - do { - - double dev; - - scale_iterate(cr, full, pmodel, nr); - - dev = log_residual(cr, full, 0, 0, NULL); - if ( fabs(dev - old_dev) < dev*0.01 ) done = 1; - - i++; - old_dev = dev; - - } while ( i < MAX_CYCLES && !done ); -} - - struct scale_args { RefList *full; @@ -341,8 +137,14 @@ struct queue_args static void scale_crystal(void *task, int id) { struct scale_args *pargs = task; - do_scale_refine(pargs->crystal, pargs->full, pargs->pmodel, - &pargs->n_reflections); + int r; + double G; + + /* Simple iterative algorithm */ + r = linear_scale(pargs->full, crystal_get_reflections(pargs->crystal), &G); + if ( r == 0 ) { + crystal_set_osf(pargs->crystal, G); + } /* else don't change it */ } @@ -375,91 +177,6 @@ static void done_crystal(void *vqargs, void *task) } -static double total_log_r(Crystal **crystals, int n_crystals, RefList *full, - int *ninc) -{ - int i; - double total = 0.0; - int n = 0; - - for ( i=0; i= 0.01*old_res) && (niter < 10) ); - - if ( niter == 10 ) { - ERROR("Too many iterations - giving up!\n"); - } -} - - /* Calculates G, by which list2 should be multiplied to fit list1 */ int linear_scale(const RefList *list1, const RefList *list2, double *G) { @@ -490,6 +207,7 @@ int linear_scale(const RefList *list1, const RefList *list2, double *G) { signed int h, k, l; double Ih1, Ih2; + double esd1, esd2; nb++; get_indices(refl1, &h, &k, &l); @@ -500,8 +218,13 @@ int linear_scale(const RefList *list1, const RefList *list2, double *G) Ih1 = get_intensity(refl1); Ih2 = get_intensity(refl2); + esd1 = get_esd_intensity(refl1); + esd2 = get_esd_intensity(refl2); + if ( Ih1 <= 3.0*esd1 ) continue; + if ( Ih2 <= 3.0*esd2 ) continue; if ( (Ih1 <= 0.0) || (Ih2 <= 0.0) ) continue; + if ( get_redundancy(refl1) < 2 ) continue; if ( isnan(Ih1) || isinf(Ih1) ) continue; if ( isnan(Ih2) || isinf(Ih2) ) continue; @@ -537,7 +260,7 @@ int linear_scale(const RefList *list1, const RefList *list2, double *G) if ( isnan(*G) ) { ERROR("Scaling gave NaN (%i pairs)\n", n); - abort(); + *G = 1.0; return 1; } @@ -550,25 +273,24 @@ int linear_scale(const RefList *list1, const RefList *list2, double *G) void scale_all_to_reference(Crystal **crystals, int n_crystals, - RefList *reference) + RefList *reference, int nthreads) { - int i; - - for ( i=0; i