aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--src/hrs-scaling.c139
-rw-r--r--src/post-refinement.c1
2 files changed, 105 insertions, 35 deletions
diff --git a/src/hrs-scaling.c b/src/hrs-scaling.c
index ac9091ed..912993a5 100644
--- a/src/hrs-scaling.c
+++ b/src/hrs-scaling.c
@@ -3,11 +3,11 @@
*
* Intensity scaling using generalised HRS target function
*
- * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY,
* a research centre of the Helmholtz Association.
*
* Authors:
- * 2010-2013 Thomas White <taw@physics.org>
+ * 2010-2014 Thomas White <taw@physics.org>
*
* This file is part of CrystFEL.
*
@@ -212,6 +212,10 @@ static void run_merge_job(void *vwargs, int cookie)
RefListIterator *iter;
double G;
+ /* If this crystal's scaling was dodgy, it doesn't contribute to the
+ * merged intensities */
+ if ( crystal_get_user_flag(cr) != 0 ) return;
+
G = crystal_get_osf(cr);
for ( refl = first_refl(crystal_get_reflections(cr), &iter);
@@ -366,6 +370,10 @@ static void run_esd_job(void *vwargs, int cookie)
RefListIterator *iter;
double G;
+ /* If this crystal's scaling was dodgy, it doesn't contribute to the
+ * merged intensities */
+ if ( crystal_get_user_flag(cr) != 0 ) return;
+
G = crystal_get_osf(cr);
for ( refl = first_refl(crystal_get_reflections(cr), &iter);
@@ -446,29 +454,67 @@ static void calculate_esds(Crystal **crystals, int n, RefList *full,
}
-static double max_outlier_change(Crystal **crystals, int n, double *old_osfs)
+struct osfcheck
+{
+ int n; /* Crystal number */
+ double old_osf;
+ double new_osf;
+ double change;
+};
+
+
+
+static int compare_osf_change(const void *av, const void *bv)
+{
+ struct osfcheck *a = (struct osfcheck *)av;
+ struct osfcheck *b = (struct osfcheck *)bv;
+ if ( a->change < b->change ) return 1;
+ return -1;
+}
+
+
+static void reject_outliers(struct osfcheck *och, int n, Crystal **crystals)
+{
+ int i;
+
+ for ( i=0; i<n; i++ ) {
+ if ( (och[i].new_osf < 0.0) || (och[i].new_osf > 3.0) ) {
+ crystal_set_user_flag(crystals[och[i].n], 1);
+ crystal_set_osf(crystals[och[i].n],
+ och[i].old_osf);
+ } else {
+ crystal_set_user_flag(crystals[och[i].n], 0);
+ }
+ }
+}
+
+
+static int test_convergence(struct osfcheck *och, int n, Crystal **crystals)
{
- double rdtot = 0.0;
- double rdmax = 0.0;
- double rdmean;
- int j;
-
- for ( j=0; j<n; j++ ) {
- double r;
- r = crystal_get_osf(crystals[j]) / old_osfs[j];
- rdtot += r;
+ int i;
+ double total_change = 0.0;
+ double mean_change;
+ int n_change = 0;
+
+ for ( i=0; i<n; i++ ) {
+ och[i].change = fabs(och[i].old_osf - och[i].new_osf);
}
- rdmean = rdtot / n;
- for ( j=0; j<n; j++ ) {
- double r;
- r = crystal_get_osf(crystals[j]) / old_osfs[j];
- old_osfs[j] = crystal_get_osf(crystals[j]);
- if ( fabs(r-rdmean) > rdmax ) {
- rdmax = fabs(r-rdmean);
+
+ /* Sort the crystals according to how much they changed */
+ qsort(och, n, sizeof(struct osfcheck), compare_osf_change);
+
+ for ( i=0; i<n; i++ ) {
+ if ( crystal_get_user_flag(crystals[och[i].n]) == 0 ) {
+ total_change += och[i].change;
+ n_change++;
}
+
}
+ mean_change = total_change / n_change;
- return rdmax;
+ STATUS("Mean OSF change = %f\n", mean_change);
+
+ return mean_change < 0.01;
}
@@ -480,11 +526,11 @@ RefList *scale_intensities(Crystal **crystals, int n, RefList *gref,
int i;
double max_corr;
RefList *full = NULL;
- double *old_osfs;
- double rdval;
+ struct osfcheck *och;
+ int done;
- old_osfs = malloc(n*sizeof(double));
- if ( old_osfs == NULL ) return NULL;
+ och = malloc(n*sizeof(struct osfcheck));
+ if ( och == NULL ) return NULL;
for ( i=0; i<n; i++ ) crystal_set_osf(crystals[i], 1.0);
@@ -492,7 +538,7 @@ RefList *scale_intensities(Crystal **crystals, int n, RefList *gref,
full = lsq_intensities(crystals, n, n_threads, pmodel);
calculate_esds(crystals, n, full, n_threads, min_redundancy,
pmodel);
- free(old_osfs);
+ free(och);
return full;
}
@@ -501,15 +547,20 @@ RefList *scale_intensities(Crystal **crystals, int n, RefList *gref,
full = lsq_intensities(crystals, n, n_threads, pmodel);
}
- for ( i=0; i<n; i++ ) {
- old_osfs[i] = crystal_get_osf(crystals[i]);
- }
-
/* Iterate */
i = 0;
do {
RefList *reference;
+ double total_sf = 0.0;
+ double norm_sf;
+ int j;
+
+ for ( j=0; j<n; j++ ) {
+ och[j].n = j;
+ och[j].old_osf = crystal_get_osf(crystals[j]);
+ crystal_set_user_flag(crystals[j], 0);
+ }
/* Refine against reference or current "full" estimates */
if ( gref != NULL ) {
@@ -521,11 +572,29 @@ RefList *scale_intensities(Crystal **crystals, int n, RefList *gref,
max_corr = iterate_scale(crystals, n, reference, n_threads,
pmodel);
- rdval = max_outlier_change(crystals, n, old_osfs);
+ /* Normalise the scale factors */
+ for ( j=0; j<n; j++ ) {
+ total_sf += crystal_get_osf(crystals[j]);
+ }
+ norm_sf = total_sf / n;
+ for ( j=0; j<n; j++ ) {
+ crystal_set_osf(crystals[j],
+ crystal_get_osf(crystals[j])/norm_sf);
+ }
+
+ for ( j=0; j<n; j++ ) {
+ och[j].new_osf = crystal_get_osf(crystals[j]);
+ }
+
+ reject_outliers(och, n, crystals);
- STATUS("Scaling iteration %2i: max correction = %5.2f "
- "dev of worst outlier = %5.2f\n",
- i+1, max_corr, rdval);
+ done = test_convergence(och, n, crystals);
+
+ FILE *fh = fopen("scale-factors.log", "a");
+ for ( j=0; j<n; j++ ) {
+ fprintf(fh, "%5.2f\n", crystal_get_osf(crystals[j]));
+ }
+ fclose(fh);
/* No reference -> generate list for next iteration */
if ( gref == NULL ) {
@@ -535,7 +604,7 @@ RefList *scale_intensities(Crystal **crystals, int n, RefList *gref,
i++;
- } while ( (rdval > 0.05) && (i < MAX_CYCLES) );
+ } while ( !done && (i < MAX_CYCLES) );
if ( i == MAX_CYCLES ) {
ERROR("Warning: Scaling did not converge.\n");
@@ -547,6 +616,6 @@ RefList *scale_intensities(Crystal **crystals, int n, RefList *gref,
calculate_esds(crystals, n, full, n_threads, min_redundancy, pmodel);
- free(old_osfs);
+ free(och);
return full;
}
diff --git a/src/post-refinement.c b/src/post-refinement.c
index 0f001a21..26299f2c 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -581,6 +581,7 @@ static double pr_iterate(Crystal *cr, const RefList *full,
//STATUS("Shift %i: %e\n", param, shift);
if ( fabs(shift) > max_shift ) max_shift = fabs(shift);
}
+ crystal_set_user_flag(cr, 0);
} else {
crystal_set_user_flag(cr, 2);
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-16 09:16:21 +0000