From d1a2c71235611ad878fa9bc705c6eff1a04b3600 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Thu, 3 May 2018 14:35:32 +0200 Subject: Simplify scaling --- src/scaling.c | 350 ++++++++++----------------------------------- src/scaling.h | 8 +- tests/linear_scale_check.c | 17 ++- 3 files changed, 92 insertions(+), 283 deletions(-) diff --git a/src/scaling.c b/src/scaling.c index 92f618ae..29b0c2b4 100644 --- a/src/scaling.c +++ b/src/scaling.c @@ -48,195 +48,38 @@ #include "scaling.h" -/* Apply the given shift to the 'k'th parameter of 'image'. */ -static void apply_shift(Crystal *cr, int k, double shift) +struct scale_args { - 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(); - - } -} + RefList *full; + Crystal *crystal; +}; -/* Perform one cycle of scaling of 'cr' against 'full' */ -static double scale_iterate(Crystal *cr, const RefList *full, int *nr) +struct queue_args { - 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; -} + int n_started; + int n_done; + Crystal **crystals; + int n_crystals; + struct scale_args task_defaults; +}; -static void do_scale_refine(Crystal *cr, const RefList *full, int *nr) +static void scale_crystal(void *task, int id) { + struct scale_args *pargs = task; int i, done; double old_dev; - old_dev = log_residual(cr, full, 0, NULL, NULL); + old_dev = log_residual(pargs->crystal, pargs->full, 0, NULL, NULL); i = 0; done = 0; do { double dev; - - scale_iterate(cr, full, nr); - - dev = log_residual(cr, full, 0, 0, NULL); + scale_one_crystal(pargs->crystal, pargs->full, 0); + dev = log_residual(pargs->crystal, pargs->full, 0, 0, NULL); if ( fabs(dev - old_dev) < dev*0.01 ) done = 1; i++; @@ -246,32 +89,6 @@ static void do_scale_refine(Crystal *cr, const RefList *full, int *nr) } -struct scale_args -{ - RefList *full; - Crystal *crystal; - int n_reflections; -}; - - -struct queue_args -{ - int n_started; - int n_done; - Crystal **crystals; - int n_crystals; - long long int n_reflections; - struct scale_args task_defaults; -}; - - -static void scale_crystal(void *task, int id) -{ - struct scale_args *pargs = task; - do_scale_refine(pargs->crystal, pargs->full, &pargs->n_reflections); -} - - static void *get_crystal(void *vqargs) { struct scale_args *task; @@ -291,11 +108,7 @@ static void *get_crystal(void *vqargs) static void done_crystal(void *vqargs, void *task) { struct queue_args *qa = vqargs; - struct scale_args *wargs = task; - qa->n_done++; - qa->n_reflections += wargs->n_reflections; - progress_bar(qa->n_done, qa->n_crystals, "Scaling"); free(task); } @@ -353,11 +166,8 @@ void scale_all(Crystal **crystals, int n_crystals, int nthreads, int no_Bscale) qargs.task_defaults.full = full; qargs.n_started = 0; qargs.n_done = 0; - qargs.n_reflections = 0; run_threads(nthreads, scale_crystal, get_crystal, done_crystal, &qargs, n_crystals, 0, 0, 0); - STATUS("%lli reflections went into the scaling.\n", - qargs.n_reflections); new_res = total_log_r(crystals, n_crystals, full, &ninc); STATUS("Log residual went from %e to %e, %i crystals\n", @@ -382,13 +192,11 @@ void scale_all(Crystal **crystals, int n_crystals, int nthreads, int no_Bscale) } -/* Calculates G and B, by which list2 should be multiplied to fit list1 */ -int scale_one(const RefList *list1, const RefList *list2, int flags, - double *G, double *B) +/* Calculates G and B, by which cr's reflections should be multiplied to fit reference */ +int scale_one_crystal(Crystal *cr, const RefList *listR, int flags) { int complain_loudly = 0; - const Reflection *refl1; - const Reflection *refl2; + const Reflection *reflS; RefListIterator *iter; int max_n = 256; int n = 0; @@ -396,20 +204,24 @@ int scale_one(const RefList *list1, const RefList *list2, int flags, double *y; double *w; int r; - double cov11; - double sumsq; - int n_esd1 = 0; - int n_esd2 = 0; - int n_ih1 = 0; - int n_ih2 = 0; - int n_nan1 = 0; - int n_nan2 = 0; - int n_inf1 = 0; - int n_inf2 = 0; + double cov00, cov01, cov11, chisq; + int n_esdS = 0; + int n_esdR = 0; + int n_ihS = 0; + int n_ihR = 0; + int n_nanS = 0; + int n_nanR = 0; + int n_infS = 0; + int n_infR = 0; int n_part = 0; int n_nom = 0; + RefList *listS = crystal_get_reflections(cr); + UnitCell *cell = crystal_get_cell(cr); + double G, B; - *B = 0.0; /* FIXME */ + assert(cell != NULL); + assert(listR != NULL); + assert(listS != NULL); x = malloc(max_n*sizeof(double)); w = malloc(max_n*sizeof(double)); @@ -420,38 +232,40 @@ int scale_one(const RefList *list1, const RefList *list2, int flags, } int nb = 0; - for ( refl2 = first_refl_const(list2, &iter); - refl2 != NULL; - refl2 = next_refl_const(refl2, iter) ) + for ( reflS = first_refl_const(listS, &iter); + reflS != NULL; + reflS = next_refl_const(reflS, iter) ) { signed int h, k, l; - double Ih1, Ih2; - double esd1, esd2; + const Reflection *reflR; + double IhR, IhS, esdS, pS, LS; + double s; + nb++; - get_indices(refl2, &h, &k, &l); - refl1 = find_refl(list1, h, k, l); - if ( refl1 == NULL ) { + get_indices(reflS, &h, &k, &l); + reflR = find_refl(listR, h, k, l); + if ( reflR == NULL ) { n_nom++; continue; } - Ih1 = get_intensity(refl1); - Ih2 = get_intensity(refl2); - esd1 = get_esd_intensity(refl1); - esd2 = get_esd_intensity(refl2); + s = resolution(cell, h, k, l); + + IhR = get_intensity(reflR); + IhS = get_intensity(reflS); + esdS = get_esd_intensity(reflS); + pS = get_partiality(reflS); + LS = get_lorentz(reflS); /* Problem cases in approximate descending order of severity */ - if ( isnan(Ih1) ) { n_nan1++; continue; } - if ( isinf(Ih1) ) { n_inf1++; continue; } - if ( isnan(Ih2) ) { n_nan2++; continue; } - if ( isinf(Ih2) ) { n_inf2++; continue; } - if ( get_partiality(refl2) < 0.3 ) { n_part++; continue; } - //if ( Ih1 <= 0.0 ) { n_ih1++; continue; } - if ( Ih2 <= 0.0 ) { n_ih2++; continue; } - //if ( Ih1 <= 3.0*esd1 ) { n_esd1++; continue; } - if ( Ih2 <= 3.0*esd2 ) { n_esd2++; continue; } - //if ( get_redundancy(refl1) < 2 ) continue; + if ( isnan(IhR) ) { n_nanR++; continue; } + if ( isinf(IhR) ) { n_infR++; continue; } + if ( isnan(IhS) ) { n_nanS++; continue; } + if ( isinf(IhS) ) { n_infS++; continue; } + if ( pS < 0.3 ) { n_part++; continue; } + if ( IhS <= 0.0 ) { n_ihS++; continue; } + if ( IhS <= 3.0*esdS ) { n_esdS++; continue; } if ( n == max_n ) { max_n *= 2; @@ -464,9 +278,9 @@ int scale_one(const RefList *list1, const RefList *list2, int flags, } } - x[n] = Ih2 / get_partiality(refl2); - y[n] = Ih1; - w[n] = get_partiality(refl2); + x[n] = s*s; + y[n] = log(IhR) - log(pS) + log(LS) - log(IhS); + w[n] = pS; n++; } @@ -474,62 +288,48 @@ int scale_one(const RefList *list1, const RefList *list2, int flags, if ( n < 2 ) { if ( complain_loudly ) { ERROR("Not enough reflections for scaling (had %i, but %i remain)\n", nb, n); - if ( n_esd1 ) ERROR("%i reference reflection esd\n", n_esd1); - if ( n_esd2 ) ERROR("%i subject reflection esd\n", n_esd2); - if ( n_ih1 ) ERROR("%i reference reflection intensity\n", n_ih1); - if ( n_ih2 ) ERROR("%i subject reflection intensity\n", n_ih2); - if ( n_nan1 ) ERROR("%i reference reflection nan\n", n_nan1); - if ( n_nan2 ) ERROR("%i subject reflection nan\n", n_nan2); - if ( n_inf1 ) ERROR("%i reference reflection inf\n", n_inf1); - if ( n_inf2 ) ERROR("%i subject reflection inf\n", n_inf2); + if ( n_esdR ) ERROR("%i reference reflection esd\n", n_esdR); + if ( n_esdS ) ERROR("%i subject reflection esd\n", n_esdS); + if ( n_ihR ) ERROR("%i reference reflection intensity\n", n_ihR); + if ( n_ihS ) ERROR("%i subject reflection intensity\n", n_ihS); + if ( n_nanR ) ERROR("%i reference reflection nan\n", n_nanR); + if ( n_nanS ) ERROR("%i subject reflection nan\n", n_nanS); + if ( n_infR ) ERROR("%i reference reflection inf\n", n_infR); + if ( n_infS ) ERROR("%i subject reflection inf\n", n_infS); if ( n_part ) ERROR("%i subject reflection partiality\n", n_part); if ( n_nom ) ERROR("%i no match in reference list\n", n_nom); } - *G = 1.0; return 1; } - r = gsl_fit_wmul(x, 1, w, 1, y, 1, n, G, &cov11, &sumsq); + r = gsl_fit_wlinear(x, 1, w, 1, y, 1, n, &G, &B, &cov00, &cov01, &cov11, &chisq); if ( r ) { ERROR("Scaling failed.\n"); - *G = 1.0; return 1; } - if ( isnan(*G) ) { + if ( isnan(G) ) { if ( complain_loudly ) { ERROR("Scaling gave NaN (%i pairs)\n", n); if ( n < 10 ) { int i; for ( i=0; i