diff options
| -rw-r--r-- | Makefile.am | 4 | ||||
| -rw-r--r-- | src/partialator.c | 79 | ||||
| -rw-r--r-- | src/post-refinement.c | 334 | ||||
| -rw-r--r-- | src/post-refinement.h | 8 | ||||
| -rw-r--r-- | src/scaling.c | 439 | ||||
| -rw-r--r-- | src/scaling.h | 47 |
6 files changed, 548 insertions, 363 deletions
diff --git a/Makefile.am b/Makefile.am index e8782455..1940affc 100644 --- a/Makefile.am +++ b/Makefile.am @@ -95,7 +95,7 @@ src_render_hkl_SOURCES = src/render_hkl.c endif src_partialator_SOURCES = src/partialator.c src/post-refinement.c \ - src/merge.c src/rejection.c + src/merge.c src/rejection.c src/scaling.c src_ambigator_SOURCES = src/ambigator.c @@ -133,7 +133,7 @@ EXTRA_DIST += src/dw-hdfsee.h src/hdfsee.h src/render_hkl.h \ src/cl-utils.h src/hdfsee-render.h src/diffraction.h \ src/diffraction-gpu.h src/pattern_sim.h src/list_tmp.h \ src/im-sandbox.h src/process_image.h src/multihistogram.h \ - src/rejection.h + src/rejection.h src/scaling.h crystfeldir = $(datadir)/crystfel crystfel_DATA = data/diffraction.cl data/hdfsee.ui diff --git a/src/partialator.c b/src/partialator.c index 26c72696..6de69cd2 100644 --- a/src/partialator.c +++ b/src/partialator.c @@ -54,6 +54,7 @@ #include <cell-utils.h> #include "version.h" +#include "scaling.h" #include "post-refinement.h" #include "merge.h" #include "rejection.h" @@ -564,23 +565,35 @@ static void write_pgraph(RefList *full, Crystal **crystals, int n_crystals, } -static void normalise_scales(Crystal **crystals, int n_crystals) +static void all_residuals(Crystal **crystals, int n_crystals, RefList *full, + double *presidual, double *pfree_residual, + double *plog_residual, double *pfree_log_residual) { int i; - int n = 0; - double sum_B = 0.0; - double mean_B; + + *presidual = 0.0; + *pfree_residual = 0.0; + *plog_residual = 0.0; + *pfree_log_residual = 0.0; for ( i=0; i<n_crystals; i++ ) { - if ( crystal_get_user_flag(crystals[i]) != 0 ) continue; - sum_B += crystal_get_Bfac(crystals[i]); - n++; - } - mean_B = sum_B / n; - STATUS("Mean B = %e A^2\n", mean_B*1e20); - for ( i=0; i<n_crystals; i++ ) { - double B = crystal_get_Bfac(crystals[i]); - crystal_set_Bfac(crystals[i], B-mean_B); + + double r, free_r, log_r, free_log_r; + + if ( crystal_get_user_flag(crystals[i]) ) continue; + + r = residual(crystals[i], full, 0, NULL, NULL); + free_r = residual(crystals[i], full, 1, NULL, NULL); + log_r = log_residual(crystals[i], full, 0, NULL, NULL); + free_log_r = log_residual(crystals[i], full, 1, NULL, NULL); + + if ( isnan(r) || isnan(free_r) + || isnan(log_r) || isnan(free_log_r) ) continue; + + *presidual += r; + *pfree_residual += free_r; + *plog_residual += log_r; + *pfree_log_residual += free_log_r; } } @@ -993,40 +1006,50 @@ int main(int argc, char *argv[]) STATUS("Checking patterns.\n"); //early_rejection(crystals, n_crystals); - /* Make initial estimates */ + /* Initial rejection, figures of merit etc */ full = merge_intensities(crystals, n_crystals, nthreads, pmodel, min_measurements, push_res, 1); - check_rejection(crystals, n_crystals, full); - + show_all_residuals(crystals, n_crystals, full); write_pgraph(full, crystals, n_crystals, 0); /* Iterate */ for ( i=0; i<n_iter; i++ ) { - show_all_residuals(crystals, n_crystals, full); - - STATUS("Refinement cycle %i of %i\n", i+1, n_iter); + STATUS("Scaling and refinement cycle %i of %i\n", i+1, n_iter); - /* Refine all crystals to get the best fit */ - refine_all(crystals, n_crystals, full, nthreads, pmodel, - no_scale, no_pr, max_B); + if ( !no_scale ) { + scale_all(crystals, n_crystals, nthreads, pmodel, + max_B); + reflist_free(full); + full = merge_intensities(crystals, n_crystals, nthreads, + pmodel, min_measurements, + push_res, 1); + } - show_all_residuals(crystals, n_crystals, full); + if ( !no_pr ) { + refine_all(crystals, n_crystals, full, nthreads, + pmodel); + reflist_free(full); + full = merge_intensities(crystals, n_crystals, nthreads, + pmodel, min_measurements, + push_res, 1); + } check_rejection(crystals, n_crystals, full); - normalise_scales(crystals, n_crystals); - - /* Re-estimate all the full intensities */ reflist_free(full); full = merge_intensities(crystals, n_crystals, nthreads, - pmodel, min_measurements, push_res, 1); - + pmodel, min_measurements, + push_res, 1); + show_all_residuals(crystals, n_crystals, full); write_pgraph(full, crystals, n_crystals, i+1); } + full = merge_intensities(crystals, n_crystals, nthreads, pmodel, + min_measurements, push_res, 1); show_all_residuals(crystals, n_crystals, full); + write_pgraph(full, crystals, n_crystals, n_iter+1); /* Output results */ STATUS("Writing overall results to %s\n", outfile); diff --git a/src/post-refinement.c b/src/post-refinement.c index cf3108db..487ce872 100644 --- a/src/post-refinement.c +++ b/src/post-refinement.c @@ -216,21 +216,9 @@ double gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel) struct image *image = crystal_get_image(cr); double R = crystal_get_profile_radius(cr); double gr; - signed int hi, ki, li; - double s; - get_indices(refl, &hi, &ki, &li); - s = resolution(crystal_get_cell(cr), hi, ki, li); get_partial(refl, &rlow, &rhigh, &p); - if ( k == GPARAM_OSF ) { - return 1.0; - } - - if ( k == GPARAM_BFAC ) { - return -s*s; - } - if ( k == GPARAM_R ) { double Rglow, Rghigh; @@ -347,18 +335,6 @@ static void apply_shift(Crystal *cr, int k, double shift) crystal_set_profile_radius(cr, t); break; - case GPARAM_BFAC : - t = crystal_get_Bfac(cr); - t += shift; - crystal_set_Bfac(cr, t); - break; - - case GPARAM_OSF : - t = crystal_get_osf(cr); - t += shift; - crystal_set_osf(cr, t); - break; - case GPARAM_ASX : case GPARAM_ASY : case GPARAM_ASZ : @@ -379,141 +355,6 @@ static void apply_shift(Crystal *cr, int k, double shift) } -/* Perform one cycle of scaling of 'cr' against 'full' */ -static double scale_iterate(Crystal *cr, const RefList *full, - PartialityModel pmodel) -{ - 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; - - 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<num_params; k++ ) { - gradients[k] = gradient(cr, rv[k], refl, pmodel); - } - - for ( k=0; k<num_params; k++ ) { - - int g; - double v_c, v_curr; - - for ( g=0; g<num_params; g++ ) { - - double M_c, M_curr; - - /* Matrix is symmetric */ - if ( g > 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 = 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++; - } - - 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<num_params; param++ ) { - double shift = gsl_vector_get(shifts, param); - apply_shift(cr, rv[param], shift); - if ( fabs(shift) > 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; -} - - /* Perform one cycle of post refinement on 'image' against 'full' */ static double pr_iterate(Crystal *cr, const RefList *full, PartialityModel pmodel, @@ -667,74 +508,8 @@ static double pr_iterate(Crystal *cr, const RefList *full, } -static double log_residual(Crystal *cr, const RefList *full, int free, - int *pn_used, const char *filename) -{ - double dev = 0.0; - double G, B; - Reflection *refl; - RefListIterator *iter; - int n_used = 0; - FILE *fh = NULL; - - G = crystal_get_osf(cr); - B = crystal_get_Bfac(cr); - if ( filename != NULL ) { - fh = fopen(filename, "a"); - if ( fh == NULL ) { - ERROR("Failed to open '%s'\n", filename); - } - } - - for ( refl = first_refl(crystal_get_reflections(cr), &iter); - refl != NULL; - refl = next_refl(refl, iter) ) - { - double p, L, s, w; - signed int h, k, l; - Reflection *match; - double esd, I_full, I_partial; - double fx, dc; - - if ( free && !get_flag(refl) ) continue; - - get_indices(refl, &h, &k, &l); - match = find_refl(full, h, k, l); - if ( match == NULL ) continue; - - p = get_partiality(refl); - L = get_lorentz(refl); - I_partial = get_intensity(refl); - I_full = get_intensity(match); - esd = get_esd_intensity(refl); - s = resolution(crystal_get_cell(cr), h, k, l); - - 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 */ - - fx = G + log(p) - log(L) - B*s*s + log(I_full); - dc = log(I_partial) - fx; - w = 1.0; - dev += w*dc*dc; - - if ( fh != NULL ) { - fprintf(fh, "%4i %4i %4i %e %e\n", - h, k, l, s, dev); - } - - } - - if ( fh != NULL ) fclose(fh); - - if ( pn_used != NULL ) *pn_used = n_used; - return dev; -} - - -static double residual(Crystal *cr, const RefList *full, int free, - int *pn_used, const char *filename) +double residual(Crystal *cr, const RefList *full, int free, + int *pn_used, const char *filename) { double dev = 0.0; double G, B; @@ -799,39 +574,6 @@ static double residual(Crystal *cr, const RefList *full, int free, } -void all_residuals(Crystal **crystals, int n_crystals, RefList *full, - double *presidual, double *pfree_residual, - double *plog_residual, double *pfree_log_residual) -{ - int i; - - *presidual = 0.0; - *pfree_residual = 0.0; - *plog_residual = 0.0; - *pfree_log_residual = 0.0; - - for ( i=0; i<n_crystals; i++ ) { - - double r, free_r, log_r, free_log_r; - - if ( crystal_get_user_flag(crystals[i]) ) continue; - - r = residual(crystals[i], full, 0, NULL, NULL); - free_r = residual(crystals[i], full, 1, NULL, NULL); - log_r = log_residual(crystals[i], full, 0, NULL, NULL); - free_log_r = log_residual(crystals[i], full, 1, NULL, NULL); - - if ( isnan(r) || isnan(free_r) - || isnan(log_r) || isnan(free_log_r) ) continue; - - *presidual += r; - *pfree_residual += free_r; - *plog_residual += log_r; - *pfree_log_residual += free_log_r; - } -} - - static void write_residual_graph(Crystal *cr, const RefList *full) { int i; @@ -877,49 +619,6 @@ static void write_residual_graph(Crystal *cr, const RefList *full) } -static void do_scale_refine(Crystal *cr, const RefList *full, - PartialityModel pmodel, int verbose) -{ - int i, done; - double old_dev; - - old_dev = log_residual(cr, full, 0, NULL, NULL); - - if ( verbose ) { - STATUS("Initial G=%.2f, B=%e\n", - crystal_get_osf(cr), crystal_get_Bfac(cr)); - STATUS("Scaling initial dev = %10.5e, free dev = %10.5e\n", - old_dev, log_residual(cr, full, 1, NULL, NULL)); - } - - i = 0; - done = 0; - do { - - double dev; - - scale_iterate(cr, full, pmodel); - - dev = log_residual(cr, full, 0, 0, NULL); - if ( fabs(dev - old_dev) < dev*0.01 ) done = 1; - - if ( verbose ) { - STATUS("Scaling %2i: dev = %10.5e, free dev = %10.5e\n", - i+1, dev, log_residual(cr, full, 0, 0, NULL)); - } - - i++; - old_dev = dev; - - } while ( i < MAX_CYCLES && !done ); - - if ( verbose ) { - STATUS("Final G=%.2f, B=%e\n", crystal_get_osf(cr), - crystal_get_Bfac(cr)); - } -} - - static void do_pr_refine(Crystal *cr, const RefList *full, PartialityModel pmodel, int verbose) { @@ -977,8 +676,7 @@ static void do_pr_refine(Crystal *cr, const RefList *full, static struct prdata pr_refine(Crystal *cr, const RefList *full, - PartialityModel pmodel, int no_scale, int no_pr, - double max_B) + PartialityModel pmodel) { int verbose = 0; struct prdata prdata; @@ -986,23 +684,11 @@ static struct prdata pr_refine(Crystal *cr, const RefList *full, prdata.refined = 0; prdata.n_filtered = 0; - if ( !no_scale ) { - do_scale_refine(cr, full, pmodel, verbose); - } - - /* Reject if B factor modulus is very large */ - if ( fabs(crystal_get_Bfac(cr)) > max_B ) { - crystal_set_user_flag(cr, PRFLAG_BIGB); - return prdata; - } - if ( verbose ) { write_residual_graph(cr, full); } - if ( !no_pr ) { - do_pr_refine(cr, full, pmodel, verbose); - } + do_pr_refine(cr, full, pmodel, verbose); if ( crystal_get_user_flag(cr) == 0 ) { prdata.refined = 1; @@ -1017,9 +703,6 @@ struct refine_args RefList *full; Crystal *crystal; PartialityModel pmodel; - int no_scale; - int no_pr; - double max_B; struct prdata prdata; }; @@ -1039,8 +722,7 @@ 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->no_scale, pargs->no_pr, pargs->max_B); + pargs->prdata = pr_refine(cr, pargs->full, pargs->pmodel); } @@ -1072,8 +754,7 @@ static void done_image(void *vqargs, void *task) void refine_all(Crystal **crystals, int n_crystals, - RefList *full, int nthreads, PartialityModel pmodel, - int no_scale, int no_pr, double max_B) + RefList *full, int nthreads, PartialityModel pmodel) { struct refine_args task_defaults; struct queue_args qargs; @@ -1083,9 +764,6 @@ void refine_all(Crystal **crystals, int n_crystals, task_defaults.pmodel = pmodel; task_defaults.prdata.refined = 0; task_defaults.prdata.n_filtered = 0; - task_defaults.no_scale = no_scale; - task_defaults.no_pr = no_pr; - task_defaults.max_B = max_B; qargs.task_defaults = task_defaults; qargs.n_started = 0; diff --git a/src/post-refinement.h b/src/post-refinement.h index b3509501..7f395013 100644 --- a/src/post-refinement.h +++ b/src/post-refinement.h @@ -57,15 +57,13 @@ enum prflag extern const char *str_prflag(enum prflag flag); extern void refine_all(Crystal **crystals, int n_crystals, - RefList *full, int nthreads, PartialityModel pmodel, - int no_scale, int no_pr, double max_B); + RefList *full, int nthreads, PartialityModel pmodel); /* Exported so it can be poked by tests/pr_p_gradient_check */ extern double gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel); -extern void all_residuals(Crystal **crystals, int n_crystals, RefList *full, - double *residual, double *free_residual, - double *log_residual, double *free_log_residual); +extern double residual(Crystal *cr, const RefList *full, int free, + int *pn_used, const char *filename); #endif /* POST_REFINEMENT_H */ diff --git a/src/scaling.c b/src/scaling.c new file mode 100644 index 00000000..b3b9137e --- /dev/null +++ b/src/scaling.c @@ -0,0 +1,439 @@ +/* + * scaling.c + * + * Scaling + * + * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. + * + * Authors: + * 2010-2015 Thomas White <taw@physics.org> + * + * This file is part of CrystFEL. + * + * CrystFEL is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * CrystFEL is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with CrystFEL. If not, see <http://www.gnu.org/licenses/>. + * + */ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + + +#include <stdlib.h> +#include <assert.h> +#include <gsl/gsl_matrix.h> +#include <gsl/gsl_vector.h> +#include <gsl/gsl_linalg.h> +#include <gsl/gsl_eigen.h> +#include <gsl/gsl_blas.h> + +#include "merge.h" +#include "post-refinement.h" +#include "symmetry.h" +#include "cell.h" +#include "cell-utils.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 = crystal_get_osf(cr); + t += shift; + crystal_set_osf(cr, 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) +{ + 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; + + 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<num_params; k++ ) { + + if ( rv[k] == GPARAM_OSF ) { + gradients[k] = 1.0; + } else if ( rv[k] == GPARAM_BFAC ) { + gradients[k] = -s*s; + } else { + ERROR("Unrecognised scaling gradient.\n"); + abort(); + } + } + + for ( k=0; k<num_params; k++ ) { + + int g; + double v_c, v_curr; + + for ( g=0; g<num_params; g++ ) { + + double M_c, M_curr; + + /* Matrix is symmetric */ + if ( g > 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 = 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++; + } + + 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<num_params; param++ ) { + double shift = gsl_vector_get(shifts, param); + apply_shift(cr, rv[param], shift); + if ( fabs(shift) > 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) +{ + double dev = 0.0; + double G, B; + Reflection *refl; + RefListIterator *iter; + int n_used = 0; + FILE *fh = NULL; + + G = crystal_get_osf(cr); + B = crystal_get_Bfac(cr); + if ( filename != NULL ) { + fh = fopen(filename, "a"); + if ( fh == NULL ) { + ERROR("Failed to open '%s'\n", filename); + } + } + + for ( refl = first_refl(crystal_get_reflections(cr), &iter); + refl != NULL; + refl = next_refl(refl, iter) ) + { + double p, L, s, w; + signed int h, k, l; + Reflection *match; + double esd, I_full, I_partial; + double fx, dc; + + if ( free && !get_flag(refl) ) continue; + + get_indices(refl, &h, &k, &l); + match = find_refl(full, h, k, l); + if ( match == NULL ) continue; + + p = get_partiality(refl); + L = get_lorentz(refl); + I_partial = get_intensity(refl); + I_full = get_intensity(match); + esd = get_esd_intensity(refl); + s = resolution(crystal_get_cell(cr), h, k, l); + + 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 */ + + fx = G + log(p) - log(L) - B*s*s + log(I_full); + dc = log(I_partial) - fx; + w = 1.0; + dev += w*dc*dc; + + if ( fh != NULL ) { + fprintf(fh, "%4i %4i %4i %e %e\n", + h, k, l, s, dev); + } + + } + + if ( fh != NULL ) fclose(fh); + + if ( pn_used != NULL ) *pn_used = n_used; + return dev; +} + + +static void do_scale_refine(Crystal *cr, const RefList *full, + PartialityModel pmodel) +{ + 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); + + 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; + Crystal *crystal; + PartialityModel pmodel; + double max_B; +}; + + +struct queue_args +{ + int n_started; + int n_done; + Crystal **crystals; + int n_crystals; + 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->pmodel); + + /* Reject if B factor modulus is very large */ + if ( fabs(crystal_get_Bfac(pargs->crystal)) > pargs->max_B ) { + crystal_set_user_flag(pargs->crystal, PRFLAG_BIGB); + } +} + + +static void *get_crystal(void *vqargs) +{ + struct scale_args *task; + struct queue_args *qargs = vqargs; + + task = malloc(sizeof(struct scale_args)); + memcpy(task, &qargs->task_defaults, sizeof(struct scale_args)); + + task->crystal = qargs->crystals[qargs->n_started]; + + qargs->n_started++; + + return task; +} + + +static void done_crystal(void *vqargs, void *task) +{ + struct queue_args *qa = vqargs; + + qa->n_done++; + + progress_bar(qa->n_done, qa->n_crystals, "Scaling"); + free(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<n_crystals; i++ ) { + + double r; + if ( crystal_get_user_flag(crystals[i]) ) continue; + r = log_residual(crystals[i], full, 0, NULL, NULL); + if ( isnan(r) ) continue; + total += r; + n++; + + } + + if ( ninc != NULL ) *ninc = n; + return total; +} + + +/* Perform iterative scaling, all the way to convergence */ +void scale_all(Crystal **crystals, int n_crystals, int nthreads, + PartialityModel pmodel, double max_B) +{ + struct scale_args task_defaults; + struct queue_args qargs; + double old_res, new_res; + + task_defaults.crystal = NULL; + task_defaults.pmodel = pmodel; + task_defaults.max_B = max_B; + + qargs.task_defaults = task_defaults; + qargs.n_crystals = n_crystals; + qargs.crystals = crystals; + + /* Don't have threads which are doing nothing */ + if ( n_crystals < nthreads ) nthreads = n_crystals; + + new_res = INFINITY; + do { + RefList *full; + int ninc; + double bef_res; + + full = merge_intensities(crystals, n_crystals, nthreads, + pmodel, 2, INFINITY, 0); + old_res = new_res; + bef_res = total_log_r(crystals, n_crystals, full, NULL); + + qargs.task_defaults.full = full; + qargs.n_started = 0; + qargs.n_done = 0; + run_threads(nthreads, scale_crystal, get_crystal, done_crystal, + &qargs, n_crystals, 0, 0, 0); + + new_res = total_log_r(crystals, n_crystals, full, &ninc); + STATUS("Log residual went from %e to %e, %i crystals\n", + bef_res, new_res, ninc); + + reflist_free(full); + + } while ( fabs(new_res-old_res) >= 0.01*old_res ); +} diff --git a/src/scaling.h b/src/scaling.h new file mode 100644 index 00000000..f8eb4e9e --- /dev/null +++ b/src/scaling.h @@ -0,0 +1,47 @@ +/* + * scaling.h + * + * Scaling + * + * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. + * + * Authors: + * 2010-2015 Thomas White <taw@physics.org> + * + * This file is part of CrystFEL. + * + * CrystFEL is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * CrystFEL is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with CrystFEL. If not, see <http://www.gnu.org/licenses/>. + * + */ + +#ifndef SCALING_H +#define SCALING_H + + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + + +#include "crystal.h" +#include "geometry.h" + +extern double log_residual(Crystal *cr, const RefList *full, int free, + int *pn_used, const char *filename); + +extern void scale_all(Crystal **crystals, int n_crystals, int nthreads, + PartialityModel pmodel, double max_B); + +#endif /* SCALING_H */ |