From 294965d42b309e98c8952d3a5dea753af21713a6 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Wed, 2 May 2018 17:37:12 +0200 Subject: Preparation for adjusting B factor during post-refinement Add --no-Bscale option to partialator, and pass down as far as needed residual() no longer does scaling: call scale_one_crystal() first if necessary scale_one() replaces old linear_scale() function to scale a pair of RefLists (but so far does the same as the old function) --- src/partialator.c | 35 +++++++++++++++++++++-------- src/post-refinement.c | 55 ++++++++++++++++++++++++++-------------------- src/post-refinement.h | 4 ++-- src/scaling.c | 27 ++++++++++++++++++----- src/scaling.h | 18 ++++++++++----- tests/linear_scale_check.c | 10 ++++----- 6 files changed, 99 insertions(+), 50 deletions(-) diff --git a/src/partialator.c b/src/partialator.c index 0f3bd104..e775caef 100644 --- a/src/partialator.c +++ b/src/partialator.c @@ -317,6 +317,7 @@ static void show_help(const char *s) " --stop-after= Stop after merging crystals.\n" " -n, --iterations= Run cycles of scaling and post-refinement.\n" " --no-scale Disable scale factor (G, B) refinement.\n" +" --no-Bscale Disable B factor scaling.\n" " --no-pr Disable orientation/physics refinement.\n" " -m, --model= Specify partiality model.\n" " --min-measurements= Minimum number of measurements to require.\n" @@ -717,6 +718,7 @@ static void all_residuals(Crystal **crystals, int n_crystals, RefList *full, if ( crystal_get_user_flag(crystals[i]) ) continue; + /* Scaling should have been done right before calling this */ 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); @@ -754,6 +756,7 @@ struct log_qargs Crystal **crystals; int n_crystals; RefList *full; + int scaleflags; }; @@ -761,6 +764,7 @@ struct log_args { Crystal *cr; RefList *full; + int scaleflags; int iter; int cnum; }; @@ -780,6 +784,7 @@ static void *get_log_task(void *vp) task->full = qargs->full; task->iter = qargs->iter; task->cnum = qargs->next; + task->scaleflags = qargs->scaleflags; qargs->next += 20; return task; @@ -790,7 +795,8 @@ static void write_logs(void *vp, int cookie) { struct log_args *args = vp; write_specgraph(args->cr, args->full, args->iter, args->cnum); - write_gridscan(args->cr, args->full, args->iter, args->cnum); + write_gridscan(args->cr, args->full, args->iter, args->cnum, + args->scaleflags); write_test_logs(args->cr, args->full, args->iter, args->cnum); } @@ -803,16 +809,17 @@ static void done_log(void *qargs, void *vp) static void write_logs_parallel(Crystal **crystals, int n_crystals, - RefList *full, int iter, int n_threads) + RefList *full, int iter, int n_threads, + int scaleflags) { struct log_qargs qargs; - qargs.iter = iter; qargs.next = 0; qargs.full = full; qargs.crystals = crystals; qargs.n_crystals = n_crystals; + qargs.scaleflags = scaleflags; STATUS("Writing logs...\n"); run_threads(n_threads, write_logs, get_log_task, done_log, &qargs, @@ -839,6 +846,7 @@ int main(int argc, char *argv[]) int n_crystals_seen = 0; char cmdline[1024]; int no_scale = 0; + int no_Bscale = 0; int no_pr = 0; Stream *st; Crystal **crystals; @@ -867,6 +875,7 @@ int main(int argc, char *argv[]) double force_bandwidth = -1.0; double force_radius = -1.0; char *audit_info; + int scaleflags = 0; /* Long options */ const struct option longopts[] = { @@ -895,6 +904,7 @@ int main(int argc, char *argv[]) {"force-radius", 1, NULL, 11}, {"no-scale", 0, &no_scale, 1}, + {"no-Bscale", 0, &no_Bscale, 1}, {"no-pr", 0, &no_pr, 1}, {"no-polarisation", 0, &polarisation, 0}, {"no-polarization", 0, &polarisation, 0}, /* compat */ @@ -1146,6 +1156,10 @@ int main(int argc, char *argv[]) "partialities (--model=unity).\n"); } + if ( no_Bscale ) { + scaleflags |= SCALE_NO_B; + } + if ( !no_logs ) { int r = mkdir("pr-logs", S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH); if ( r ) { @@ -1327,7 +1341,7 @@ int main(int argc, char *argv[]) /* Create reference data set if we don't already have one */ if ( reference == NULL ) { STATUS("Initial scaling...\n"); - scale_all(crystals, n_crystals, nthreads); + scale_all(crystals, n_crystals, nthreads, scaleflags); full = merge_intensities(crystals, n_crystals, nthreads, min_measurements, push_res, 1); } else { @@ -1339,7 +1353,8 @@ int main(int argc, char *argv[]) show_all_residuals(crystals, n_crystals, full); if ( !no_pr && !no_logs ) { write_pgraph(full, crystals, n_crystals, 0, ""); - write_logs_parallel(crystals, n_crystals, full, 0, nthreads); + write_logs_parallel(crystals, n_crystals, full, 0, nthreads, + scaleflags); } /* Iterate */ @@ -1349,14 +1364,15 @@ int main(int argc, char *argv[]) if ( !no_pr ) { refine_all(crystals, n_crystals, full, nthreads, pmodel, - 0, i+1, no_logs, sym, amb); + 0, i+1, no_logs, sym, amb, scaleflags); } /* Create new reference if needed */ if ( reference == NULL ) { reflist_free(full); if ( !no_scale ) { - scale_all(crystals, n_crystals, nthreads); + scale_all(crystals, n_crystals, nthreads, + scaleflags); } full = merge_intensities(crystals, n_crystals, nthreads, min_measurements, @@ -1403,7 +1419,7 @@ int main(int argc, char *argv[]) if ( reference == NULL ) { reflist_free(full); if ( !no_scale ) { - scale_all(crystals, n_crystals, nthreads); + scale_all(crystals, n_crystals, nthreads, scaleflags); } full = merge_intensities(crystals, n_crystals, nthreads, min_measurements, @@ -1417,7 +1433,8 @@ int main(int argc, char *argv[]) show_all_residuals(crystals, n_crystals, full); if ( !no_pr && !no_logs ) { write_pgraph(full, crystals, n_crystals, -1, ""); - write_logs_parallel(crystals, n_crystals, full, -1, nthreads); + write_logs_parallel(crystals, n_crystals, full, -1, nthreads, + scaleflags); } /* Output results */ diff --git a/src/post-refinement.c b/src/post-refinement.c index f260da7f..293fe332 100644 --- a/src/post-refinement.c +++ b/src/post-refinement.c @@ -95,10 +95,6 @@ double residual(Crystal *cr, const RefList *full, int free, double B = crystal_get_Bfac(cr); UnitCell *cell = crystal_get_cell(cr); - if ( linear_scale(full, crystal_get_reflections(cr), &G, complain) ) { - return GSL_NAN; - } - for ( refl = first_refl(crystal_get_reflections(cr), &iter); refl != NULL; refl = next_refl(refl, iter) ) @@ -232,6 +228,7 @@ struct rf_priv { int verbose; int serial; gsl_vector *initial; + int scaleflags; /* For freeing later */ gsl_vector *vals; @@ -349,6 +346,11 @@ static double residual_f(const gsl_vector *v, void *pp) update_predictions(cr); calculate_partialities(cr, PMODEL_XSPHERE); + if ( scale_one_crystal(cr, pv->full, pv->scaleflags) ) { + crystal_free(cr); + if ( pv->verbose ) STATUS("Bad scaling\n"); + return GSL_NAN; + } res = residual(cr, pv->full, 0, NULL, NULL); cell_free(crystal_get_cell(cr)); @@ -463,8 +465,8 @@ static void reindex_cell(UnitCell *cell, SymOpList *amb, int idx) } -void try_reindex(Crystal *crin, const RefList *full, - SymOpList *sym, SymOpList *amb) +static void try_reindex(Crystal *crin, const RefList *full, + SymOpList *sym, SymOpList *amb, int scaleflags) { RefList *list; Crystal *cr; @@ -474,6 +476,7 @@ void try_reindex(Crystal *crin, const RefList *full, if ( sym == NULL || amb == NULL ) return; + if ( scale_one_crystal(crin, full, scaleflags) ) return; residual_original = residual(crin, full, 0, NULL, NULL); cr = crystal_copy(crin); @@ -494,6 +497,7 @@ void try_reindex(Crystal *crin, const RefList *full, update_predictions(cr); calculate_partialities(cr, PMODEL_XSPHERE); + if ( scale_one_crystal(cr, full, scaleflags) ) return; residual_flipped = residual(cr, full, 0, NULL, NULL); if ( residual_flipped < residual_original ) { @@ -637,6 +641,7 @@ void write_specgraph(Crystal *crystal, const RefList *full, static gsl_multimin_fminimizer *setup_minimiser(Crystal *cr, const RefList *full, int verbose, int serial, + int scaleflags, struct rf_priv *priv) { gsl_multimin_fminimizer *min; @@ -654,6 +659,7 @@ static gsl_multimin_fminimizer *setup_minimiser(Crystal *cr, const RefList *full priv->full = full; priv->verbose = verbose; priv->serial = serial; + priv->scaleflags = scaleflags; priv->f.f = residual_f; priv->f.n = n_params; @@ -678,7 +684,7 @@ static gsl_multimin_fminimizer *setup_minimiser(Crystal *cr, const RefList *full static void write_grid(Crystal *cr, const RefList *full, - signed int cycle, int serial, + signed int cycle, int serial, int scaleflags, const enum gparam par1, const enum gparam par2, const char *name) { @@ -690,7 +696,7 @@ static void write_grid(Crystal *cr, const RefList *full, int idx1, idx2; int i; - min = setup_minimiser(cr, full, 0, serial, &priv); + min = setup_minimiser(cr, full, 0, serial, scaleflags, &priv); idx1 = 99; idx2 = 99; @@ -747,13 +753,13 @@ static void write_grid(Crystal *cr, const RefList *full, void write_gridscan(Crystal *cr, const RefList *full, - signed int cycle, int serial) + signed int cycle, int serial, int scaleflags) { - write_grid(cr, full, cycle, serial, + write_grid(cr, full, cycle, serial, scaleflags, GPARAM_ANG1, GPARAM_ANG2, "ang1-ang2"); - write_grid(cr, full, cycle, serial, + write_grid(cr, full, cycle, serial, scaleflags, GPARAM_ANG1, GPARAM_WAVELENGTH, "ang1-wave"); - write_grid(cr, full, cycle, serial, + write_grid(cr, full, cycle, serial, scaleflags, GPARAM_R, GPARAM_WAVELENGTH, "R-wave"); } @@ -761,19 +767,21 @@ void write_gridscan(Crystal *cr, const RefList *full, static void do_pr_refine(Crystal *cr, const RefList *full, PartialityModel pmodel, int verbose, int serial, int cycle, int write_logs, - SymOpList *sym, SymOpList *amb) + SymOpList *sym, SymOpList *amb, int scaleflags) { gsl_multimin_fminimizer *min; struct rf_priv priv; int n_iter = 0; int status; - int r; - double G; double residual_start, residual_free_start; FILE *fh = NULL; - try_reindex(cr, full, sym, amb); + try_reindex(cr, full, sym, amb, scaleflags); + if ( scale_one_crystal(cr, full, scaleflags) ) { + ERROR("Bad scaling at start of refinement.\n"); + return; + } residual_start = residual(cr, full, 0, NULL, NULL); residual_free_start = residual(cr, full, 1, NULL, NULL); @@ -782,7 +790,7 @@ static void do_pr_refine(Crystal *cr, const RefList *full, residual_start, residual_free_start); } - min = setup_minimiser(cr, full, verbose, serial, &priv); + min = setup_minimiser(cr, full, verbose, serial, scaleflags, &priv); if ( verbose ) { double res = residual_f(min->x, &priv); @@ -902,10 +910,7 @@ static void do_pr_refine(Crystal *cr, const RefList *full, apply_parameters(min->x, priv.initial, priv.rv, cr); update_predictions(cr); calculate_partialities(cr, PMODEL_XSPHERE); - r = linear_scale(full, crystal_get_reflections(cr), &G, 0); - if ( r == 0 ) { - crystal_set_osf(cr, G); - } + scale_one_crystal(cr, full, scaleflags); if ( verbose ) { @@ -922,7 +927,7 @@ static void do_pr_refine(Crystal *cr, const RefList *full, } if ( write_logs ) { - write_gridscan(cr, full, cycle, serial); + write_gridscan(cr, full, cycle, serial, scaleflags); write_specgraph(cr, full, cycle, serial); write_test_logs(cr, full, cycle, serial); } @@ -950,6 +955,7 @@ struct refine_args int no_logs; SymOpList *sym; SymOpList *amb; + int scaleflags; }; @@ -974,7 +980,7 @@ static void refine_image(void *task, int id) do_pr_refine(cr, pargs->full, pargs->pmodel, pargs->verbose, pargs->serial, pargs->cycle, write_logs, - pargs->sym, pargs->amb); + pargs->sym, pargs->amb, pargs->scaleflags); if ( crystal_get_user_flag(cr) == 0 ) { pargs->prdata.refined = 1; @@ -1013,7 +1019,7 @@ static void done_image(void *vqargs, void *task) void refine_all(Crystal **crystals, int n_crystals, RefList *full, int nthreads, PartialityModel pmodel, int verbose, int cycle, int no_logs, - SymOpList *sym, SymOpList *amb) + SymOpList *sym, SymOpList *amb, int scaleflags) { struct refine_args task_defaults; struct queue_args qargs; @@ -1027,6 +1033,7 @@ void refine_all(Crystal **crystals, int n_crystals, task_defaults.no_logs = no_logs; task_defaults.sym = sym; task_defaults.amb = amb; + task_defaults.scaleflags = scaleflags; qargs.task_defaults = task_defaults; qargs.n_started = 0; diff --git a/src/post-refinement.h b/src/post-refinement.h index 8c345729..b8923d2c 100644 --- a/src/post-refinement.h +++ b/src/post-refinement.h @@ -61,10 +61,10 @@ extern const char *str_prflag(enum prflag flag); extern void refine_all(Crystal **crystals, int n_crystals, RefList *full, int nthreads, PartialityModel pmodel, int verbose, int cycle, int no_logs, - SymOpList *sym, SymOpList *amb); + SymOpList *sym, SymOpList *amb, int scaleflags); extern void write_gridscan(Crystal *cr, const RefList *full, - int cycle, int serial); + int cycle, int serial, int scaleflags); extern void write_specgraph(Crystal *crystal, const RefList *full, signed int cycle, int serial); diff --git a/src/scaling.c b/src/scaling.c index 39541926..5f896e65 100644 --- a/src/scaling.c +++ b/src/scaling.c @@ -389,7 +389,7 @@ static double total_log_r(Crystal **crystals, int n_crystals, RefList *full, /* Perform iterative scaling, all the way to convergence */ -void scale_all(Crystal **crystals, int n_crystals, int nthreads) +void scale_all(Crystal **crystals, int n_crystals, int nthreads, int no_Bscale) { struct scale_args task_defaults; struct queue_args qargs; @@ -448,10 +448,11 @@ void scale_all(Crystal **crystals, int n_crystals, int nthreads) } -/* 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) +/* 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) { + int complain_loudly = 0; const Reflection *refl1; const Reflection *refl2; RefListIterator *iter; @@ -474,6 +475,8 @@ int linear_scale(const RefList *list1, const RefList *list2, double *G, int n_part = 0; int n_nom = 0; + *B = 0.0; /* FIXME */ + x = malloc(max_n*sizeof(double)); w = malloc(max_n*sizeof(double)); y = malloc(max_n*sizeof(double)); @@ -556,6 +559,7 @@ int linear_scale(const RefList *list1, const RefList *list2, double *G, if ( r ) { ERROR("Scaling failed.\n"); + *G = 1.0; return 1; } @@ -567,7 +571,7 @@ int linear_scale(const RefList *list1, const RefList *list2, double *G, int i; for ( i=0; i + * 2010-2018 Thomas White * * This file is part of CrystFEL. * @@ -38,12 +38,20 @@ #include "crystal.h" #include "geometry.h" +enum ScaleFlags +{ + SCALE_NO_B, /* Don't use Debye-Waller part */ +}; + extern double log_residual(Crystal *cr, const RefList *full, int free, int *pn_used, const char *filename); -extern int linear_scale(const RefList *list1, const RefList *list2, double *G, - int complain_loudly); +extern int scale_one(const RefList *list1, const RefList *list2, int flags, + double *G, double *B); + +extern int scale_one_crystal(Crystal *cr, const RefList *list2, int flags); -extern void scale_all(Crystal **crystals, int n_crystals, int nthreads); +extern void scale_all(Crystal **crystals, int n_crystals, int nthreads, + int flags); #endif /* SCALING_H */ diff --git a/tests/linear_scale_check.c b/tests/linear_scale_check.c index 4cb6a1ba..5c723849 100644 --- a/tests/linear_scale_check.c +++ b/tests/linear_scale_check.c @@ -3,11 +3,11 @@ * * Check that linear scaling works * - * Copyright © 2017 Deutsches Elektronen-Synchrotron DESY, - * a research centre of the Helmholtz Association. + * Copyright © 2017-2018 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. * * Authors: - * 2017 Thomas White + * 2017-2018 Thomas White * * This file is part of CrystFEL. * @@ -45,7 +45,7 @@ int main(int argc, char *argv[]) gsl_rng *rng; RefList *list1; RefList *list2; - double G; + double G, B; int r; list1 = reflist_new(); @@ -70,7 +70,7 @@ int main(int argc, char *argv[]) set_partiality(refl2, 1.0); } - r = linear_scale(list1, list2, &G, 1); + r = scale_one(list1, list2, SCALE_NO_B, &G, &B); STATUS("Scaling result: %i, G = %f\n", r, G); return fail; -- cgit v1.2.3