aboutsummaryrefslogtreecommitdiff
path: root/src/post-refinement.c
diff options
context:
space:
mode:
Diffstat (limited to 'src/post-refinement.c')
-rw-r--r--src/post-refinement.c128
1 files changed, 97 insertions, 31 deletions
diff --git a/src/post-refinement.c b/src/post-refinement.c
index 9369cfe9..165080b4 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -34,6 +34,7 @@
#include <stdlib.h>
#include <assert.h>
#include <gsl/gsl_multimin.h>
+#include <gsl/gsl_fit.h>
#include "image.h"
#include "post-refinement.h"
@@ -78,35 +79,108 @@ 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 dev = 0.0;
- double G, B;
+ double G;
Reflection *refl;
RefListIterator *iter;
- FILE *fh = NULL;
int n_used = 0;
+ double num = 0.0;
+ double den = 0.0;
- if ( filename != NULL ) {
- fh = fopen(filename, "w");
- if ( fh == NULL ) {
- ERROR("Failed to open '%s'\n", filename);
- }
+ if ( linear_scale(crystal_get_reflections(cr), full, &G) ) {
+ return INFINITY;
}
- G = crystal_get_osf(cr);
- B = crystal_get_Bfac(cr);
-
for ( refl = first_refl(crystal_get_reflections(cr), &iter);
refl != NULL;
refl = next_refl(refl, iter) )
{
- double p, L, s, w;
+ double p, w;
signed int h, k, l;
Reflection *match;
- double esd, I_full, I_partial;
- double fx, dc;
+ double I_full;
+ double int1, int2;
if ( free && !get_flag(refl) ) continue;
@@ -115,34 +189,26 @@ double residual(Crystal *cr, const RefList *full, int free,
if ( match == NULL ) continue;
I_full = get_intensity(match);
- if ( get_redundancy(match) < 2 ) continue;
+ if ( get_redundancy(match) < 1 ) continue; /* FIXME: 2 */
p = get_partiality(refl);
//if ( p < 0.2 ) continue;
- I_partial = get_intensity(refl);
- esd = get_esd_intensity(refl);
- //if ( I_partial < 3.0*esd ) continue;
+ int1 = get_intensity(refl) / G;
+ int2 = p*I_full;
+ w = 1.0;
- L = get_lorentz(refl);
- s = resolution(crystal_get_cell(cr), h, k, l);
+ if ( int1 + int2 < 0.0 ) continue;
+
+ num += fabs(int1 - int2) * w;
+ den += (int1 + int2) * w/2.0;
- fx = exp(G)*p*exp(-B*s*s)*I_full/L;
- dc = I_partial - fx;
- w = (s/1e9)*(s/1e9)/(esd*esd);
- dev += w*dc*dc;
n_used++;
- 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;
+ return num/(den*sqrt(2));
}