From bc24110f035c62b3e7e873a2b84c47bcf7414f19 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Tue, 27 Feb 2018 13:44:29 +0100 Subject: Scaling fixes --- src/merge.c | 2 +- src/partialator.c | 91 +++++++++++++++++++++++++++++++++------------------ src/post-refinement.c | 2 +- src/rejection.c | 2 +- src/scaling.c | 10 ++++-- 5 files changed, 71 insertions(+), 36 deletions(-) (limited to 'src') diff --git a/src/merge.c b/src/merge.c index 33221bda..099bad2d 100644 --- a/src/merge.c +++ b/src/merge.c @@ -189,7 +189,7 @@ static void run_merge_job(void *vwargs, int cookie) } /* Total (multiplicative) correction factor */ - corr = (1.0/G) * exp(B*res*res) * get_lorentz(refl) + corr = G * exp(B*res*res) * get_lorentz(refl) / get_partiality(refl); if ( isnan(corr) ) { ERROR("NaN corr:\n"); diff --git a/src/partialator.c b/src/partialator.c index 46a5a25f..124ce2e4 100644 --- a/src/partialator.c +++ b/src/partialator.c @@ -637,7 +637,7 @@ static void write_to_pgraph(FILE *fh, RefList *list, RefList *full, Crystal *cr, pcalc = get_partiality(refl); /* Observed partiality */ - corr = (1.0/G) * exp(B*res*res) * get_lorentz(refl); + corr = G * exp(B*res*res) * get_lorentz(refl); Ipart = get_intensity(refl) * corr; pobs = Ipart / get_intensity(match); @@ -689,7 +689,7 @@ static void write_specgraph(RefList *full, Crystal *crystal, int in) match = find_refl(full, h, k, l); if ( match == NULL ) continue; - corr = (1.0/G) * exp(B*res*res) * get_lorentz(refl); + corr = G * exp(B*res*res) * get_lorentz(refl); Ipart = get_intensity(refl) * corr; Ifull = get_intensity(match); esd = get_esd_intensity(match); @@ -847,7 +847,6 @@ int main(int argc, char *argv[]) double max_B = 1e-18; char *rfile = NULL; RefList *reference = NULL; - RefList *r; /* Long options */ const struct option longopts[] = { @@ -1223,13 +1222,22 @@ int main(int argc, char *argv[]) //early_rejection(crystals, n_crystals); /* Initial rejection, figures of merit etc */ - full = merge_intensities(crystals, n_crystals, nthreads, pmodel, - min_measurements, push_res, 1); + if ( reference == NULL ) { + 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); - r = (reference != NULL) ? reference : full; - show_all_residuals(crystals, n_crystals, r); - write_pgraph(r, crystals, n_crystals, 0); - write_specgraph(r, crystals[0], 0); + + if ( !no_scale ) { + scale_all_to_reference(crystals, n_crystals, full); + } + + show_all_residuals(crystals, n_crystals, full); + write_pgraph(full, crystals, n_crystals, 0); + write_specgraph(full, crystals[0], 0); /* Iterate */ for ( i=0; i XSCALE-like algorithm */ - scale_all(crystals, n_crystals, nthreads, pmodel); - reflist_free(full); full = merge_intensities(crystals, n_crystals, nthreads, pmodel, min_measurements, @@ -1258,21 +1261,27 @@ int main(int argc, char *argv[]) } else { /* Reference -> Ginn-style linear algorithm */ - scale_all_to_reference(crystals, n_crystals, - reference); + if ( !no_scale ) { + scale_all_to_reference(crystals, n_crystals, + reference); + } } } check_rejection(crystals, n_crystals, full, max_B); - reflist_free(full); - full = merge_intensities(crystals, n_crystals, nthreads, - pmodel, min_measurements, - push_res, 1); - r = (reference != NULL) ? reference : full; - show_all_residuals(crystals, n_crystals, r); - write_pgraph(r, crystals, n_crystals, i+1); - write_specgraph(r, crystals[0], i+1); + + if ( reference == NULL ) { + 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 */ + + show_all_residuals(crystals, n_crystals, full); + write_pgraph(full, crystals, n_crystals, i+1); + write_specgraph(full, crystals[0], i+1); if ( output_everycycle ) { @@ -1306,12 +1315,32 @@ int main(int argc, char *argv[]) } } - full = merge_intensities(crystals, n_crystals, nthreads, pmodel, - min_measurements, push_res, 1); - r = (reference != NULL) ? reference : full; - show_all_residuals(crystals, n_crystals, r); - write_pgraph(r, crystals, n_crystals, n_iter+1); - write_specgraph(r, crystals[0], n_iter+1); + + if ( reference == NULL ) { + 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); + } + 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 */ /* Output results */ STATUS("Writing overall results to %s\n", outfile); diff --git a/src/post-refinement.c b/src/post-refinement.c index e1a2e05a..8e2241f6 100644 --- a/src/post-refinement.c +++ b/src/post-refinement.c @@ -92,7 +92,7 @@ double residual(Crystal *cr, const RefList *full, int free, double B = crystal_get_Bfac(cr); UnitCell *cell = crystal_get_cell(cr); - if ( linear_scale(crystal_get_reflections(cr), full, &G) ) { + if ( linear_scale(full, crystal_get_reflections(cr), &G) ) { return INFINITY; } diff --git a/src/rejection.c b/src/rejection.c index a44b697f..9d41f9b4 100644 --- a/src/rejection.c +++ b/src/rejection.c @@ -137,7 +137,7 @@ static void check_cc(Crystal *cr, RefList *full) pcalc = get_partiality(refl); /* Observed partiality */ - corr = (1.0/G) * exp(B*res*res) * get_lorentz(refl); + corr = G * exp(B*res*res) * get_lorentz(refl); Ipart = get_intensity(refl) * corr; pobs = Ipart / get_intensity(match); diff --git a/src/scaling.c b/src/scaling.c index 97503b89..ab913229 100644 --- a/src/scaling.c +++ b/src/scaling.c @@ -515,9 +515,9 @@ int linear_scale(const RefList *list1, const RefList *list2, double *G) } } - x[n] = Ih2; + x[n] = Ih2 / get_partiality(refl2); y[n] = Ih1; - w[n] = get_partiality(refl1); + w[n] = get_partiality(refl2); n++; } @@ -534,6 +534,12 @@ int linear_scale(const RefList *list1, const RefList *list2, double *G) return 1; } + if ( isnan(*G) ) { + ERROR("Scaling gave NaN (%i pairs)\n", n); + abort(); + return 1; + } + free(x); free(y); free(w); -- cgit v1.2.3