aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2015-11-11 17:31:34 +0100
committerThomas White <taw@physics.org>2015-11-18 11:17:57 +0100
commita2b11362e440d3487520a39284ae72fcb4cb37f5 (patch)
tree8e9aaf9c965468b8e160d62c8ba09ca7aa4c95b9
parentaf7c721a1c75692ecf9ba9840a5847af90a48b53 (diff)
partialator: Scale (strictly) using strong reflections only
-rw-r--r--Makefile.am4
-rw-r--r--src/partialator.c79
-rw-r--r--src/post-refinement.c334
-rw-r--r--src/post-refinement.h8
-rw-r--r--src/scaling.c439
-rw-r--r--src/scaling.h47
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 */