aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2018-02-27 13:44:29 +0100
committerThomas White <taw@physics.org>2018-02-27 17:12:42 +0100
commitbc24110f035c62b3e7e873a2b84c47bcf7414f19 (patch)
tree9208ba54ae6178437665aafa5178bb3868e118bf
parentde9fbc4335a5d2e3896381403f4b1b5b19dcf444 (diff)
Scaling fixes
-rw-r--r--src/merge.c2
-rw-r--r--src/partialator.c91
-rw-r--r--src/post-refinement.c2
-rw-r--r--src/rejection.c2
-rw-r--r--src/scaling.c10
5 files changed, 71 insertions, 36 deletions
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<n_iter; i++ ) {
@@ -1237,10 +1245,7 @@ int main(int argc, char *argv[])
STATUS("Scaling and refinement cycle %i of %i\n", i+1, n_iter);
if ( !no_pr ) {
-
- r = (reference != NULL) ? reference : full;
- refine_all(crystals, n_crystals, r, nthreads, pmodel);
-
+ refine_all(crystals, n_crystals, full, nthreads, pmodel);
}
if ( !no_scale ) {
@@ -1248,8 +1253,6 @@ int main(int argc, char *argv[])
if ( reference == NULL ) {
/* No reference -> 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);