From c8a51805efd3deb068b17b5f79c1afcfd5e8d73d Mon Sep 17 00:00:00 2001 From: Thomas White Date: Fri, 2 Mar 2018 16:38:43 +0100 Subject: Restore old scaling code --- src/scaling.c | 310 +++++++++++++++++++++++++++++++++++++++++++++++++++++++--- src/scaling.h | 3 + 2 files changed, 301 insertions(+), 12 deletions(-) (limited to 'src') diff --git a/src/scaling.c b/src/scaling.c index 4f92f501..0ad50634 100644 --- a/src/scaling.c +++ b/src/scaling.c @@ -48,6 +48,179 @@ #include "scaling.h" +/* 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, 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) { @@ -114,11 +287,35 @@ double log_residual(Crystal *cr, const RefList *full, int free, } +static void do_scale_refine(Crystal *cr, const RefList *full, 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, 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 < 10 && !done ); +} + + struct scale_args { RefList *full; Crystal *crystal; - PartialityModel pmodel; int n_reflections; }; @@ -137,16 +334,7 @@ struct queue_args static void scale_crystal(void *task, int id) { struct scale_args *pargs = task; - int r; - double G; - - /* Simple iterative algorithm */ - r = linear_scale(pargs->full, crystal_get_reflections(pargs->crystal), &G, 0); - if ( r == 0 ) { - crystal_set_osf(pargs->crystal, G); - } else { - crystal_set_user_flag(pargs->crystal, PRFLAG_SCALEBAD); - } + do_scale_refine(pargs->crystal, pargs->full, &pargs->n_reflections); } @@ -179,6 +367,104 @@ 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"); + } +} + + +static void scale_crystal_linear(void *task, int id) +{ + struct scale_args *pargs = task; + int r; + double G; + + /* Simple iterative algorithm */ + r = linear_scale(pargs->full, crystal_get_reflections(pargs->crystal), &G, 0); + if ( r == 0 ) { + crystal_set_osf(pargs->crystal, G); + } else { + crystal_set_user_flag(pargs->crystal, PRFLAG_SCALEBAD); + } +} + + /* Calculates G, by which list2 should be multiplied to fit list1 */ int linear_scale(const RefList *list1, const RefList *list2, double *G, int complain_loudly) @@ -333,6 +619,6 @@ void scale_all_to_reference(Crystal **crystals, int n_crystals, qargs.n_started = 0; qargs.n_done = 0; qargs.n_reflections = 0; - run_threads(nthreads, scale_crystal, get_crystal, done_crystal, + run_threads(nthreads, scale_crystal_linear, get_crystal, done_crystal, &qargs, n_crystals, 0, 0, 0); } diff --git a/src/scaling.h b/src/scaling.h index b345029a..620518ef 100644 --- a/src/scaling.h +++ b/src/scaling.h @@ -44,6 +44,9 @@ extern double log_residual(Crystal *cr, const RefList *full, int free, extern int linear_scale(const RefList *list1, const RefList *list2, double *G, int complain_loudly); +extern void scale_all(Crystal **crystals, int n_crystals, int nthreads, + PartialityModel pmodel); + extern void scale_all_to_reference(Crystal **crystals, int n_crystals, RefList *reference, int nthreads); -- cgit v1.2.3