aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2015-05-07 17:03:01 +0200
committerThomas White <taw@physics.org>2015-05-13 13:48:33 +0200
commitcf68fc854e148d410bba1540f1738290e6f5313e (patch)
tree17f72cacca80da30bb6c0ca91fcbc4b43a1fbdbe
parentb7540780d46b92f9ce96de74eb4c2e05b2369983 (diff)
partialator: Move scaling calculation into PR proper
-rw-r--r--libcrystfel/src/geometry.h4
-rw-r--r--src/partialator.c33
-rw-r--r--src/post-refinement.c74
-rw-r--r--src/post-refinement.h10
-rw-r--r--tests/pr_p_gradient_check.c2
5 files changed, 70 insertions, 53 deletions
diff --git a/libcrystfel/src/geometry.h b/libcrystfel/src/geometry.h
index cecdc4f1..152c0e47 100644
--- a/libcrystfel/src/geometry.h
+++ b/libcrystfel/src/geometry.h
@@ -77,7 +77,9 @@ enum gparam {
GPARAM_DIV,
GPARAM_DETX,
GPARAM_DETY,
- GPARAM_CLEN
+ GPARAM_CLEN,
+ GPARAM_OSF,
+ GPARAM_BFAC
};
diff --git a/src/partialator.c b/src/partialator.c
index 9bdb684b..19449590 100644
--- a/src/partialator.c
+++ b/src/partialator.c
@@ -85,6 +85,7 @@ struct refine_args
RefList *full;
Crystal *crystal;
PartialityModel pmodel;
+ int no_scale;
struct prdata prdata;
};
@@ -105,7 +106,8 @@ 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->prdata = pr_refine(cr, pargs->full, pargs->pmodel,
+ pargs->no_scale);
}
@@ -148,10 +150,6 @@ static void refine_all(Crystal **crystals, int n_crystals,
struct refine_args task_defaults;
struct queue_args qargs;
- /* If the partiality model is "p=1", this refinement is really, really
- * easy... */
- if ( pmodel == PMODEL_UNITY ) return;
-
task_defaults.full = full;
task_defaults.crystal = NULL;
task_defaults.pmodel = pmodel;
@@ -641,15 +639,8 @@ int main(int argc, char *argv[])
// early_rejection(crystals, n_crystals);
/* Make initial estimates */
- STATUS("Performing initial scaling.\n");
- if ( noscale ) {
- STATUS("Skipping scaling step (--no-scale).\n");
- full = lsq_intensities(crystals, n_crystals, nthreads, pmodel,
- min_measurements);
- } else {
- full = scale_intensities(crystals, n_crystals, nthreads, pmodel,
- min_measurements);
- }
+ full = lsq_intensities(crystals, n_crystals, nthreads, pmodel,
+ min_measurements);
check_rejection(crystals, n_crystals);
@@ -668,11 +659,11 @@ int main(int argc, char *argv[])
/* Iterate */
for ( i=0; i<n_iter; i++ ) {
- STATUS("Post refinement cycle %i of %i\n", i+1, n_iter);
+ STATUS("Refinement cycle %i of %i\n", i+1, n_iter);
srdata.n_filtered = 0;
- /* Refine the geometry of all patterns to get the best fit */
+ /* Refine all crystals to get the best fit */
refine_all(crystals, n_crystals, full, nthreads, pmodel,
&srdata);
@@ -681,14 +672,8 @@ int main(int argc, char *argv[])
/* Re-estimate all the full intensities */
reflist_free(full);
- if ( noscale ) {
- STATUS("Skipping scaling step (--no-scale).\n");
- full = lsq_intensities(crystals, n_crystals, nthreads,
- pmodel, min_measurements);
- } else {
- full = scale_intensities(crystals, n_crystals, nthreads,
- pmodel, min_measurements);
- }
+ full = lsq_intensities(crystals, n_crystals, nthreads,
+ pmodel, min_measurements);
check_rejection(crystals, n_crystals);
diff --git a/src/post-refinement.c b/src/post-refinement.c
index 8e6f6108..5a1053b8 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -174,14 +174,16 @@ static double volume_fraction(double rlow, double rhigh, double pr,
}
-/* Return the gradient of partiality wrt parameter 'k' given the current status
- * of 'image'. */
-double p_gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel)
+/* Return the gradient of p/G wrt parameter 'k' given the current status
+ * of the crystal. */
+double gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel)
{
double glow, ghigh;
double rlow, rhigh, p;
struct image *image = crystal_get_image(cr);
double R = crystal_get_profile_radius(cr);
+ double osf = crystal_get_osf(cr);
+ double gr;
get_partial(refl, &rlow, &rhigh, &p);
@@ -196,7 +198,8 @@ double p_gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel)
Rglow = volume_fraction_rgradient(rlow, R, pmodel);
Rghigh = volume_fraction_rgradient(rhigh, R, pmodel);
- return 4.0*psph/(3.0*D) + (4.0*R/(3.0*D))*(Rglow - Rghigh);
+ gr = 4.0*psph/(3.0*D) + (4.0*R/(3.0*D))*(Rglow - Rghigh);
+ return gr / osf;
}
@@ -214,11 +217,17 @@ double p_gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel)
get_symmetric_indices(refl, &hs, &ks, &ls);
ds = 2.0 * resolution(crystal_get_cell(cr), hs, ks, ls);
- return (ds/2.0)*(glow+ghigh) - 4.0*R*psph*ds/(3.0*D*D);
+ gr = (ds/2.0)*(glow+ghigh) - 4.0*R*psph*ds/(3.0*D*D);
+ return gr / osf;
}
- return r_gradient(crystal_get_cell(cr), k, refl, image) * (glow-ghigh);
+ if ( k == GPARAM_OSF ) {
+ return -p/(osf*osf);
+ }
+
+ gr = r_gradient(crystal_get_cell(cr), k, refl, image) * (glow-ghigh);
+ return gr / osf;
}
@@ -274,6 +283,18 @@ static void apply_shift(Crystal *cr, int k, double shift)
crystal_set_profile_radius(cr, t);
break;
+ case GPARAM_OSF :
+ if ( isnan(shift) ) {
+ ERROR("Refusing nan shift of OSF\n");
+ } else if (shift > 100.0 ) {
+ ERROR("Refusing large (%e) OSF shift.\n", shift);
+ } else {
+ t = crystal_get_osf(cr);
+ t += shift;
+ crystal_set_osf(cr, t);
+ }
+ break;
+
case GPARAM_ASX :
case GPARAM_ASY :
case GPARAM_ASZ :
@@ -296,7 +317,7 @@ static void apply_shift(Crystal *cr, int k, double shift)
/* Perform one cycle of post refinement on 'image' against 'full' */
static double pr_iterate(Crystal *cr, const RefList *full,
- PartialityModel pmodel, int *n_filtered)
+ PartialityModel pmodel, int no_scale, int *n_filtered)
{
gsl_matrix *M;
gsl_vector *v;
@@ -310,6 +331,7 @@ static double pr_iterate(Crystal *cr, const RefList *full,
const int verbose = 0;
int num_params = 0;
enum gparam rv[32];
+ double G;
*n_filtered = 0;
@@ -327,13 +349,19 @@ static double pr_iterate(Crystal *cr, const RefList *full,
rv[num_params++] = GPARAM_CSZ;
}
- STATUS("Refining %i parameters.\n", num_params);
+ /* If we are scaling, refine scale factors (duh) */
+ if ( !no_scale ) {
+ rv[num_params++] = GPARAM_OSF;
+ //rv[num_params++] = GPARAM_BFAC; FIXME: Next
+ }
reflections = crystal_get_reflections(cr);
M = gsl_matrix_calloc(num_params, num_params);
v = gsl_vector_calloc(num_params);
+ G = crystal_get_osf(cr);
+
/* Construct the equations, one per reflection in this image */
for ( refl = first_refl(reflections, &iter);
refl != NULL;
@@ -344,7 +372,7 @@ static double pr_iterate(Crystal *cr, const RefList *full,
double w;
double I_partial;
int k;
- double p, l;
+ double p, L;
Reflection *match;
double gradients[num_params];
@@ -359,19 +387,22 @@ static double pr_iterate(Crystal *cr, const RefList *full,
I_full = get_intensity(match);
- /* Actual measurement of this reflection from this pattern? */
- I_partial = get_intensity(refl) / crystal_get_osf(cr);
p = get_partiality(refl);
- l = get_lorentz(refl);
+ L = get_lorentz(refl);
+
+ /* Actual measurement of this reflection from this pattern? */
+ I_partial = get_intensity(refl);
/* Calculate the weight for this reflection */
- w = pow(get_esd_intensity(refl), 2.0);
- w += l * p * I_full * pow(get_esd_intensity(match), 2.0);
- w = pow(w, -1.0);
+ //w = pow(get_esd_intensity(refl), 2.0);
+ //w += (p/L) * I_full * pow(get_esd_intensity(match), 2.0);
+ //w = pow(w, -1.0);
+ w = 1.0;
/* Calculate all gradients for this reflection */
for ( k=0; k<num_params; k++ ) {
- gradients[k] = p_gradient(cr, rv[k], refl, pmodel) * l;
+ gradients[k] = gradient(cr, rv[k], refl, pmodel);
+ gradients[k] *= I_full / L;
}
for ( k=0; k<num_params; k++ ) {
@@ -386,8 +417,7 @@ static double pr_iterate(Crystal *cr, const RefList *full,
/* Matrix is symmetric */
if ( g > k ) continue;
- M_c = gradients[g] * gradients[k];
- M_c *= w * pow(I_full, 2.0);
+ 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);
@@ -395,8 +425,8 @@ static double pr_iterate(Crystal *cr, const RefList *full,
}
- delta_I = I_partial - (l * p * I_full);
- v_c = w * delta_I * I_full * gradients[k];
+ delta_I = I_partial - (p * I_full / (L*G));
+ v_c = w * delta_I * gradients[k];
v_curr = gsl_vector_get(v, k);
gsl_vector_set(v, k, v_curr + v_c);
@@ -486,7 +516,7 @@ static double guide_dev(Crystal *cr, const RefList *full)
struct prdata pr_refine(Crystal *cr, const RefList *full,
- PartialityModel pmodel)
+ PartialityModel pmodel, int no_scale)
{
double dev;
int i;
@@ -518,7 +548,7 @@ struct prdata pr_refine(Crystal *cr, const RefList *full,
cell_get_reciprocal(crystal_get_cell(cr), &asx, &asy, &asz,
&bsx, &bsy, &bsz, &csx, &csy, &csz);
- pr_iterate(cr, full, pmodel, &prdata.n_filtered);
+ pr_iterate(cr, full, pmodel, no_scale, &prdata.n_filtered);
update_partialities(cr, pmodel);
diff --git a/src/post-refinement.h b/src/post-refinement.h
index caf5f9af..844afbd1 100644
--- a/src/post-refinement.h
+++ b/src/post-refinement.h
@@ -3,11 +3,11 @@
*
* Post refinement
*
- * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY,
* a research centre of the Helmholtz Association.
*
* Authors:
- * 2010-2014 Thomas White <taw@physics.org>
+ * 2010-2015 Thomas White <taw@physics.org>
*
* This file is part of CrystFEL.
*
@@ -51,10 +51,10 @@ struct prdata
extern struct prdata pr_refine(Crystal *cr, const RefList *full,
- PartialityModel pmodel);
+ PartialityModel pmodel, int no_scale);
/* Exported so it can be poked by tests/pr_p_gradient_check */
-extern double p_gradient(Crystal *cr, int k, Reflection *refl,
- PartialityModel pmodel);
+extern double gradient(Crystal *cr, int k, Reflection *refl,
+ PartialityModel pmodel);
#endif /* POST_REFINEMENT_H */
diff --git a/tests/pr_p_gradient_check.c b/tests/pr_p_gradient_check.c
index affddeab..1ae7b049 100644
--- a/tests/pr_p_gradient_check.c
+++ b/tests/pr_p_gradient_check.c
@@ -303,7 +303,7 @@ static double test_gradients(Crystal *cr, double incr_val, int refine,
grad = (grad1 + grad2) / 2.0;
i++;
- cgrad = p_gradient(cr, refine, refl, pmodel);
+ cgrad = gradient(cr, refine, refl, pmodel);
get_partial(refl, &r1, &r2, &p);