aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--src/partialator.c51
-rw-r--r--src/post-refinement.c94
-rw-r--r--src/scaling.c101
-rw-r--r--src/scaling.h7
4 files changed, 147 insertions, 106 deletions
diff --git a/src/partialator.c b/src/partialator.c
index 49605bf9..af04d935 100644
--- a/src/partialator.c
+++ b/src/partialator.c
@@ -654,7 +654,8 @@ static void write_specgraph(RefList *full, Crystal *crystal, int in)
char tmp[256];
Reflection *refl;
RefListIterator *iter;
- double G, B;
+ double G = crystal_get_osf(crystal);
+ double B = crystal_get_Bfac(crystal);
UnitCell *cell;
struct image *image = crystal_get_image(crystal);
@@ -671,8 +672,6 @@ static void write_specgraph(RefList *full, Crystal *crystal, int in)
fprintf(fh, "khalf/m pcalc pobs\n");
- G = crystal_get_osf(crystal);
- B = crystal_get_Bfac(crystal);
cell = crystal_get_cell(crystal);
for ( refl = first_refl(crystal_get_reflections(crystal), &iter);
@@ -848,6 +847,7 @@ 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[] = {
@@ -1226,8 +1226,10 @@ int main(int argc, char *argv[])
full = merge_intensities(crystals, n_crystals, nthreads, pmodel,
min_measurements, push_res, 1);
check_rejection(crystals, n_crystals, full, max_B);
- show_all_residuals(crystals, n_crystals, full);
- write_pgraph(full, crystals, n_crystals, 0);
+ 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);
/* Iterate */
for ( i=0; i<n_iter; i++ ) {
@@ -1235,22 +1237,36 @@ int main(int argc, char *argv[])
STATUS("Scaling and refinement cycle %i of %i\n", i+1, n_iter);
if ( !no_scale ) {
- scale_all(crystals, n_crystals, nthreads, pmodel);
- reflist_free(full);
- full = merge_intensities(crystals, n_crystals, nthreads,
- pmodel, min_measurements,
- push_res, 1);
+
+ 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,
+ push_res, 1);
+
+ } else {
+
+ /* Reference -> Ginn-style linear algorithm */
+ scale_all_to_reference(crystals, n_crystals,
+ reference);
+
+ }
}
if ( !no_pr ) {
- RefList *r = (reference != NULL) ? reference : full;
+ r = (reference != NULL) ? reference : full;
refine_all(crystals, n_crystals, r, nthreads, pmodel);
reflist_free(full);
full = merge_intensities(crystals, n_crystals, nthreads,
pmodel, min_measurements,
push_res, 1);
+
}
check_rejection(crystals, n_crystals, full, max_B);
@@ -1258,9 +1274,10 @@ int main(int argc, char *argv[])
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, i+1);
- write_specgraph(full, crystals[0], i+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 ( output_everycycle ) {
@@ -1296,8 +1313,10 @@ int main(int argc, char *argv[])
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);
+ 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);
/* Output results */
STATUS("Writing overall results to %s\n", outfile);
diff --git a/src/post-refinement.c b/src/post-refinement.c
index 165080b4..5624b15d 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -44,6 +44,7 @@
#include "cell.h"
#include "cell-utils.h"
#include "reflist-utils.h"
+#include "scaling.h"
struct prdata
@@ -79,104 +80,23 @@ const char *str_prflag(enum prflag flag)
}
-static int linear_scale(RefList *list1, const RefList *list2, double *G)
-{
- Reflection *refl1;
- Reflection *refl2;
- RefListIterator *iter;
- int max_n = 256;
- int n = 0;
- double *x;
- double *y;
- double *w;
- int r;
- double cov11;
- double sumsq;
-
- x = malloc(max_n*sizeof(double));
- w = malloc(max_n*sizeof(double));
- y = malloc(max_n*sizeof(double));
- if ( (x==NULL) || (y==NULL) || (w==NULL) ) {
- ERROR("Failed to allocate memory for scaling.\n");
- return 1;
- }
-
- for ( refl1 = first_refl(list1, &iter);
- refl1 != NULL;
- refl1 = next_refl(refl1, iter) )
- {
- signed int h, k, l;
- double Ih1, Ih2;
-
- get_indices(refl1, &h, &k, &l);
- refl2 = find_refl(list2, h, k, l);
- if ( refl2 == NULL ) continue;
-
- Ih1 = get_intensity(refl1);
- Ih2 = get_intensity(refl2);
-
- if ( (Ih1 <= 0.0) || (Ih2 <= 0.0) ) continue;
- if ( isnan(Ih1) || isinf(Ih1) ) continue;
- if ( isnan(Ih2) || isinf(Ih2) ) continue;
-
- if ( n == max_n ) {
- max_n *= 2;
- x = realloc(x, max_n*sizeof(double));
- y = realloc(y, max_n*sizeof(double));
- w = realloc(w, max_n*sizeof(double));
- if ( (x==NULL) || (y==NULL) || (w==NULL) ) {
- ERROR("Failed to allocate memory for scaling.\n");
- return 1;
- }
- }
-
- x[n] = Ih2;
- y[n] = Ih1;
- w[n] = get_partiality(refl1);
- n++;
-
- }
-
- if ( n < 2 ) {
- ERROR("Not enough reflections for scaling\n");
- return 1;
- }
-
- r = gsl_fit_wmul(x, 1, w, 1, y, 1, n, G, &cov11, &sumsq);
-
- if ( r ) {
- ERROR("Scaling failed.\n");
- return 1;
- }
-
- free(x);
- free(y);
- free(w);
-
- return 0;
-}
-
-
-
double residual(Crystal *cr, const RefList *full, int free,
int *pn_used, const char *filename)
{
- double G;
Reflection *refl;
RefListIterator *iter;
int n_used = 0;
double num = 0.0;
double den = 0.0;
-
- if ( linear_scale(crystal_get_reflections(cr), full, &G) ) {
- return INFINITY;
- }
+ double G = crystal_get_osf(cr);
+ double B = crystal_get_Bfac(cr);
+ UnitCell *cell = crystal_get_cell(cr);
for ( refl = first_refl(crystal_get_reflections(cr), &iter);
refl != NULL;
refl = next_refl(refl, iter) )
{
- double p, w;
+ double p, w, corr, res;
signed int h, k, l;
Reflection *match;
double I_full;
@@ -185,6 +105,7 @@ double residual(Crystal *cr, const RefList *full, int free,
if ( free && !get_flag(refl) ) continue;
get_indices(refl, &h, &k, &l);
+ res = resolution(cell, h, k, l);
match = find_refl(full, h, k, l);
if ( match == NULL ) continue;
I_full = get_intensity(match);
@@ -194,7 +115,8 @@ double residual(Crystal *cr, const RefList *full, int free,
p = get_partiality(refl);
//if ( p < 0.2 ) continue;
- int1 = get_intensity(refl) / G;
+ corr = exp(-G) * exp(B*res*res) * get_lorentz(refl);
+ int1 = get_intensity(refl) * corr;
int2 = p*I_full;
w = 1.0;
diff --git a/src/scaling.c b/src/scaling.c
index 77a64384..0908cfb3 100644
--- a/src/scaling.c
+++ b/src/scaling.c
@@ -3,11 +3,11 @@
*
* Scaling
*
- * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY,
* a research centre of the Helmholtz Association.
*
* Authors:
- * 2010-2015 Thomas White <taw@physics.org>
+ * 2010-2017 Thomas White <taw@physics.org>
*
* This file is part of CrystFEL.
*
@@ -38,6 +38,7 @@
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_blas.h>
+#include <gsl/gsl_fit.h>
#include "merge.h"
#include "post-refinement.h"
@@ -456,3 +457,99 @@ void scale_all(Crystal **crystals, int n_crystals, int nthreads,
ERROR("Too many iterations - giving up!\n");
}
}
+
+
+int linear_scale(RefList *list1, const RefList *list2, double *G)
+{
+ Reflection *refl1;
+ Reflection *refl2;
+ RefListIterator *iter;
+ int max_n = 256;
+ int n = 0;
+ double *x;
+ double *y;
+ double *w;
+ int r;
+ double cov11;
+ double sumsq;
+
+ x = malloc(max_n*sizeof(double));
+ w = malloc(max_n*sizeof(double));
+ y = malloc(max_n*sizeof(double));
+ if ( (x==NULL) || (y==NULL) || (w==NULL) ) {
+ ERROR("Failed to allocate memory for scaling.\n");
+ return 1;
+ }
+
+ int nb = 0;
+ for ( refl1 = first_refl(list1, &iter);
+ refl1 != NULL;
+ refl1 = next_refl(refl1, iter) )
+ {
+ signed int h, k, l;
+ double Ih1, Ih2;
+ nb++;
+
+ get_indices(refl1, &h, &k, &l);
+ refl2 = find_refl(list2, h, k, l);
+ if ( refl2 == NULL ) {
+ continue;
+ }
+
+ Ih1 = get_intensity(refl1);
+ Ih2 = get_intensity(refl2);
+
+ if ( (Ih1 <= 0.0) || (Ih2 <= 0.0) ) continue;
+ if ( isnan(Ih1) || isinf(Ih1) ) continue;
+ if ( isnan(Ih2) || isinf(Ih2) ) continue;
+
+ if ( n == max_n ) {
+ max_n *= 2;
+ x = realloc(x, max_n*sizeof(double));
+ y = realloc(y, max_n*sizeof(double));
+ w = realloc(w, max_n*sizeof(double));
+ if ( (x==NULL) || (y==NULL) || (w==NULL) ) {
+ ERROR("Failed to allocate memory for scaling.\n");
+ return 1;
+ }
+ }
+
+ x[n] = Ih2;
+ y[n] = Ih1;
+ w[n] = get_partiality(refl1);
+ n++;
+
+ }
+
+ if ( n < 2 ) {
+ ERROR("Not enough reflections for scaling (had %i, but %i remain)\n", nb, n);
+ return 1;
+ }
+
+ r = gsl_fit_wmul(x, 1, w, 1, y, 1, n, G, &cov11, &sumsq);
+
+ if ( r ) {
+ ERROR("Scaling failed.\n");
+ return 1;
+ }
+
+ free(x);
+ free(y);
+ free(w);
+
+ return 0;
+}
+
+
+void scale_all_to_reference(Crystal **crystals, int n_crystals,
+ RefList *reference)
+{
+ int i;
+
+ for ( i=0; i<n_crystals; i++ ) {
+ double G;
+ linear_scale(crystal_get_reflections(crystals[i]), reference, &G);
+ crystal_set_osf(crystals[i], -log(G));
+ crystal_set_Bfac(crystals[i], 0.0);
+ }
+}
diff --git a/src/scaling.h b/src/scaling.h
index 534043d9..7b055003 100644
--- a/src/scaling.h
+++ b/src/scaling.h
@@ -3,11 +3,11 @@
*
* Scaling
*
- * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY,
* a research centre of the Helmholtz Association.
*
* Authors:
- * 2010-2015 Thomas White <taw@physics.org>
+ * 2010-2017 Thomas White <taw@physics.org>
*
* This file is part of CrystFEL.
*
@@ -44,4 +44,7 @@ extern double log_residual(Crystal *cr, const RefList *full, int free,
extern void scale_all(Crystal **crystals, int n_crystals, int nthreads,
PartialityModel pmodel);
+extern void scale_all_to_reference(Crystal **crystals, int n_crystals,
+ RefList *reference);
+
#endif /* SCALING_H */
generated by cgit v1.2.3 (git 2.25.1) at 2024-07-02 14:03:57 +0000