aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--Makefile.am5
-rw-r--r--src/hrs-scaling.c170
-rw-r--r--src/partialator.c127
-rw-r--r--src/rejection.c83
-rw-r--r--src/rejection.h42
5 files changed, 368 insertions, 59 deletions
diff --git a/Makefile.am b/Makefile.am
index c1d28436..db3d3aac 100644
--- a/Makefile.am
+++ b/Makefile.am
@@ -92,7 +92,7 @@ src_render_hkl_SOURCES = src/render_hkl.c
endif
src_partialator_SOURCES = src/partialator.c src/post-refinement.c \
- src/hrs-scaling.c
+ src/hrs-scaling.c src/rejection.c
src_ambigator_SOURCES = src/ambigator.c
@@ -135,7 +135,8 @@ EXTRA_DIST += src/dw-hdfsee.h src/hdfsee.h src/render_hkl.h \
src/post-refinement.h src/hrs-scaling.h src/scaling-report.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/im-sandbox.h src/process_image.h src/multihistogram.h \
+ src/rejection.h
crystfeldir = $(datadir)/crystfel
crystfel_DATA = data/diffraction.cl data/hdfsee.ui
diff --git a/src/hrs-scaling.c b/src/hrs-scaling.c
index 7f15f54d..0225d625 100644
--- a/src/hrs-scaling.c
+++ b/src/hrs-scaling.c
@@ -3,11 +3,11 @@
*
* Intensity scaling using generalised HRS target function
*
- * 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.
*
@@ -37,6 +37,7 @@
#include <gsl/gsl_vector.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_eigen.h>
+#include <gsl/gsl_fit.h>
#include "image.h"
#include "peaks.h"
@@ -45,6 +46,7 @@
#include "cell.h"
#include "utils.h"
#include "reflist.h"
+#include "cell-utils.h"
/* Minimum partiality of a reflection for it to be used for scaling */
@@ -63,6 +65,7 @@ struct scale_queue_args
Crystal **crystals;
int n_started;
PartialityModel pmodel;
+ UnitCell *cell;
};
@@ -71,6 +74,7 @@ struct scale_worker_args
Crystal *crystal;
RefList *reference;
PartialityModel pmodel;
+ int crystal_number;
};
@@ -83,6 +87,7 @@ static void *create_scale_job(void *vqargs)
wargs->reference = qargs->reference;
wargs->pmodel = qargs->pmodel;
+ wargs->crystal_number = qargs->n_started;
wargs->crystal = qargs->crystals[qargs->n_started++];
return wargs;
@@ -96,9 +101,35 @@ static void run_scale_job(void *vwargs, int cookie)
RefList *reference = wargs->reference;
Reflection *refl;
RefListIterator *iter;
- double num = 0.0;
- double den = 0.0;
- double g;
+ int n = 0;
+ double *x;
+ double *y;
+ double *w;
+ int max_n = 256;
+ double c0, c1, cov00, cov01, cov11, chisq;
+ double G, B;
+ int r;
+ int debug = 0;
+ FILE *fh;
+
+ /* If this crystal's scaling was dodgy, it doesn't contribute to the
+ * merged intensities */
+ if ( crystal_get_user_flag(cr) != 0 ) return;
+
+ x = malloc(max_n*sizeof(double));
+ y = malloc(max_n*sizeof(double));
+ w = malloc(max_n*sizeof(double));
+ if ( (x==NULL) || (y==NULL) || (w==NULL) ) {
+ ERROR("Failed to allocate memory for scaling.\n");
+ return;
+ }
+
+ if ( wargs->crystal_number == 0 ) {
+ char fn[256];
+ snprintf(fn, 255, "scale-debug.%i", wargs->crystal_number);
+ debug = 1;
+ fh = fopen(fn, "w");
+ }
for ( refl = first_refl(crystal_get_reflections(cr), &iter);
refl != NULL;
@@ -107,30 +138,80 @@ static void run_scale_job(void *vwargs, int cookie)
signed int h, k, l;
double Ih, Ihl, corr;
Reflection *r;
+ double res;
+
+ get_indices(refl, &h, &k, &l);
+ res = resolution(crystal_get_cell(cr), h, k, l);
if ( (get_partiality(refl) < MIN_PART_SCALE)
- || (get_intensity(refl) < 3.0*get_esd_intensity(refl)) ) {
+ || (get_intensity(refl) < 5.0*get_esd_intensity(refl)) ) {
continue;
}
/* Look up by asymmetric indices */
- get_indices(refl, &h, &k, &l);
r = find_refl(reference, h, k, l);
- if ( r == NULL ) continue;
+ if ( r == NULL ) {
+ continue;
+ }
Ih = get_intensity(r);
-
corr = get_partiality(refl) * get_lorentz(refl);
-
Ihl = get_intensity(refl) / corr;
- num += Ih * Ihl;
- den += Ih * Ih;
+ if ( Ihl <= 0.0 ) continue;
+ if ( Ih <= 0.0 ) continue;
+ if ( isnan(Ih) || isinf(Ihl) ) continue;
+ if ( isnan(Ih) || isinf(Ih) ) 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;
+ }
+ }
+
+ x[n] = res*res;
+ y[n] = log(Ihl/Ih);
+ w[n] = 1.0;//log(get_esd_intensity(refl)*corr);
+ if ( debug ) {
+ fprintf(fh, "%i %10.2e %10.2f %10.2f %.2f %.2f %.2f\n",
+ n, x[n], y[n], w[n], Ihl, Ih, res/1e9);
+ }
+ n++;
}
- g = num / den;
- crystal_set_osf(cr, g); /* If it's NaN, it'll get rejected later */
+ if ( n < 2 ) {
+ crystal_set_user_flag(cr, 1);
+ return;
+ }
+
+ r = gsl_fit_wlinear(x, 1, w, 1, y, 1, n, &c0, &c1,
+ &cov00, &cov01, &cov11, &chisq);
+
+ G = 1.0/exp(c0);
+ B = -c1/2.0;
+ if ( debug ) {
+ printf("intercept = %e, gradient = %e\n", c0, c1);
+ printf("Scale factor = %f\n", G);
+ printf("Relative B = %.2f A^2\n", B * 1e20);
+ fclose(fh);
+ }
+
+ free(x);
+ free(y);
+ free(w);
+ if ( r || isnan(c0) ) {
+ crystal_set_user_flag(cr, 1);
+ return;
+ }
+
+ crystal_set_osf(cr, G);
+ crystal_set_Bfac(cr, B);
}
@@ -174,6 +255,7 @@ struct merge_worker_args
RefList *full;
pthread_rwlock_t *full_lock;
PartialityModel pmodel;
+ int crystal_number;
};
@@ -187,6 +269,7 @@ static void *create_merge_job(void *vqargs)
wargs->full_lock = &qargs->full_lock;
wargs->pmodel = qargs->pmodel;
+ wargs->crystal_number = qargs->n_started;
wargs->crystal = qargs->crystals[qargs->n_started++];
return wargs;
@@ -200,13 +283,23 @@ static void run_merge_job(void *vwargs, int cookie)
RefList *full = wargs->full;
Reflection *refl;
RefListIterator *iter;
- double G;
+ double G, B;
+ FILE *fh;
+ int debug = 0;
/* 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);
+ B = crystal_get_Bfac(cr);
+
+ if ( wargs->crystal_number == 0 ) {
+ char fn[256];
+ snprintf(fn, 255, "merge-debug.%i", wargs->crystal_number);
+ debug = 1;
+ fh = fopen(fn, "w");
+ }
for ( refl = first_refl(crystal_get_reflections(cr), &iter);
refl != NULL;
@@ -216,7 +309,7 @@ static void run_merge_job(void *vwargs, int cookie)
signed int h, k, l;
double num, den;
int red;
- double Ihl, corr;
+ double Ihl, corr, res;
if ( get_partiality(refl) < MIN_PART_MERGE ) continue;
@@ -262,19 +355,27 @@ static void run_merge_job(void *vwargs, int cookie)
}
+ res = resolution(crystal_get_cell(cr), h, k, l);
corr = get_partiality(refl) * get_lorentz(refl);
Ihl = get_intensity(refl) / corr;
- num += Ihl / G;
+ num += Ihl * G * exp(2.0*B*res);
den += 1.0;
red++;
+ if ( debug ) {
+ fprintf(fh, "%.2f %.2f %f %e\n", res/1e9,
+ Ihl*G*exp(2.0*B*res),
+ G, B);
+ }
+
set_temp1(f, num);
set_temp2(f, den);
set_redundancy(f, red);
unlock_reflection(f);
}
+ if ( debug ) fclose(fh);
}
@@ -449,10 +550,33 @@ static void reject_outliers(double *old_osfs, int n, Crystal **crystals)
int i;
for ( i=0; i<n; i++ ) {
+
double osf = crystal_get_osf(crystals[i]);
- if ( isnan(osf) || (osf < 0.0) || (osf > 3.0) ) {
+ double Bfac = crystal_get_Bfac(crystals[i]);
+ int bad = 0;
+
+ if ( isnan(osf) || (osf < 0.0) || (osf > 10.0) ) bad = 1;
+ if ( isnan(Bfac) || (Bfac<-40e-20) || (Bfac>40e-20) ) bad = 1;
+
+ if ( bad ) {
crystal_set_user_flag(crystals[i], 1);
+ //STATUS("crystal %i bad: osf=%f B=%f A^2\n",
+ // i, osf, Bfac*1e20);
+ } else {
+ //STATUS("OK %i: osf=%f B=%f A^2\n",
+ // i, osf, Bfac*1e20);
}
+
+ }
+}
+
+
+static void reset_scaling_flag(Crystal *crystal)
+{
+ /* Every pattern gets another chance at being scaled, but stays flagged
+ * if it was flagged for any other reason */
+ if ( crystal_get_user_flag(crystal) == 1 ) {
+ crystal_set_user_flag(crystal, 0);
}
}
@@ -490,8 +614,9 @@ RefList *scale_intensities(Crystal **crystals, int n,
int done;
for ( i=0; i<n; i++ ) {
- crystal_set_user_flag(crystals[i], 0);
+ reset_scaling_flag(crystals[i]);
crystal_set_osf(crystals[i], 1.0);
+ crystal_set_Bfac(crystals[i], 0.0);
}
if ( noscale ) {
@@ -518,26 +643,27 @@ RefList *scale_intensities(Crystal **crystals, int n,
for ( j=0; j<n; j++ ) {
old_osfs[j] = crystal_get_osf(crystals[j]);
- crystal_set_user_flag(crystals[j], 0);
+ reset_scaling_flag(crystals[j]);
}
iterate_scale(crystals, n, full, n_threads, pmodel);
+ reject_outliers(old_osfs, n, crystals);
/* Normalise the scale factors */
for ( j=0; j<n; j++ ) {
double osf = crystal_get_osf(crystals[j]);
- if ( !isnan(osf) ) {
+ if ( crystal_get_user_flag(crystals[j]) == 0 ) {
total_sf += osf;
n_sf++;
}
}
norm_sf = total_sf / n_sf;
+ STATUS("norm = %f\n", norm_sf);
for ( j=0; j<n; j++ ) {
crystal_set_osf(crystals[j],
crystal_get_osf(crystals[j])/norm_sf);
}
- reject_outliers(old_osfs, n, crystals);
done = test_convergence(old_osfs, n, crystals);
/* Generate list for next iteration */
diff --git a/src/partialator.c b/src/partialator.c
index e5f965ed..0f0c2631 100644
--- a/src/partialator.c
+++ b/src/partialator.c
@@ -3,11 +3,11 @@
*
* Scaling and post refinement for coherent nanocrystallography
*
- * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY,
* a research centre of the Helmholtz Association.
*
* Authors:
- * 2010-2013 Thomas White <taw@physics.org>
+ * 2010-2015 Thomas White <taw@physics.org>
*
* This file is part of CrystFEL.
*
@@ -55,6 +55,7 @@
#include "post-refinement.h"
#include "hrs-scaling.h"
#include "scaling-report.h"
+#include "rejection.h"
static void show_help(const char *s)
@@ -192,11 +193,29 @@ static void display_progress(int n_images, int n_crystals)
static const char *str_flags(Crystal *cr)
{
- if ( crystal_get_user_flag(cr) ) {
- return "N";
- }
+ switch ( crystal_get_user_flag(cr) ) {
+
+ case 0 :
+ return "OK";
+
+ case 1 :
+ return "bad scaling";
+
+ case 2 :
+ return "not enough reflections";
+
+ case 3 :
+ return "PR solve failed";
- return "-";
+ case 4 :
+ return "PR lost too many reflections";
+
+ case 5 :
+ return "Early rejection";
+
+ default :
+ return "Unknown flag";
+ }
}
@@ -225,6 +244,62 @@ static RefList *apply_max_adu(RefList *list, double max_adu)
}
+static void show_duds(Crystal **crystals, int n_crystals)
+{
+ int j;
+ int n_dud = 0;
+ int n_noscale = 0;
+ int n_noref = 0;
+ int n_solve = 0;
+ int n_lost = 0;
+ int n_early = 0;
+
+ for ( j=0; j<n_crystals; j++ ) {
+ int flag;
+ flag = crystal_get_user_flag(crystals[j]);
+ if ( flag != 0 ) n_dud++;
+ switch ( flag ) {
+
+ case 0:
+ break;
+
+ case 1:
+ n_noscale++;
+ break;
+
+ case 2:
+ n_noref++;
+ break;
+
+ case 3:
+ n_solve++;
+ break;
+
+ case 4:
+ n_lost++;
+ break;
+
+ case 5:
+ n_early++;
+ break;
+
+ default:
+ STATUS("Unknown flag %i\n", flag);
+ break;
+ }
+ }
+
+ if ( n_dud ) {
+ STATUS("%i bad crystals:\n", n_dud);
+ STATUS(" %i scaling failed.\n", n_noscale);
+ STATUS(" %i not enough reflections.\n", n_noref);
+ STATUS(" %i solve failed.\n", n_solve);
+ STATUS(" %i lost too many reflections.\n", n_lost);
+ STATUS(" %i early rejection.\n", n_early);
+ }
+}
+
+
int main(int argc, char *argv[])
{
int c;
@@ -462,6 +537,7 @@ int main(int argc, char *argv[])
as = asymmetric_indices(cr_refl, sym);
crystal_set_reflections(cr, as);
+ crystal_set_user_flag(cr, 0);
reflist_free(cr_refl);
n_crystals++;
@@ -499,6 +575,10 @@ int main(int argc, char *argv[])
}
}
+ /* Make a first pass at cutting out crap */
+ STATUS("Checking patterns.\n");
+ early_rejection(crystals, n_crystals);
+
/* Make initial estimates */
STATUS("Performing initial scaling.\n");
if ( noscale ) STATUS("Scale factors fixed at 1.\n");
@@ -515,14 +595,10 @@ int main(int argc, char *argv[])
infile, cmdline);
sr_iteration(sr, 0, &srdata);
+ show_duds(crystals, n_crystals);
+
/* Iterate */
for ( i=0; i<n_iter; i++ ) {
- int n_noscale = 0;
- int n_noref = 0;
- int n_solve = 0;
- int n_lost = 0;
- int n_dud = 0;
- int j;
STATUS("Post refinement cycle %i of %i\n", i+1, n_iter);
@@ -532,29 +608,8 @@ int main(int argc, char *argv[])
refine_all(crystals, n_crystals, full, nthreads, pmodel,
&srdata);
- for ( j=0; j<n_crystals; j++ ) {
- int flag;
- flag = crystal_get_user_flag(crystals[j]);
- if ( flag != 0 ) n_dud++;
- if ( flag == 1 ) {
- n_noscale++;
- } else if ( flag == 2 ) {
- n_noref++;
- } else if ( flag == 3 ) {
- n_solve++;
- } else if ( flag == 4 ) {
- n_lost++;
- }
- }
+ show_duds(crystals, n_crystals);
- if ( n_dud ) {
- STATUS("%i crystals could not be refined this cycle.\n",
- n_dud);
- STATUS("%i scaling failed.\n", n_noscale);
- STATUS("%i not enough reflections.\n", n_noref);
- STATUS("%i solve failed.\n", n_solve);
- STATUS("%i lost too many reflections.\n", n_lost);
- }
/* Re-estimate all the full intensities */
reflist_free(full);
full = scale_intensities(crystals, n_crystals, nthreads,
@@ -577,9 +632,11 @@ int main(int argc, char *argv[])
if ( fh == NULL ) {
ERROR("Couldn't open partialator.params!\n");
} else {
+ fprintf(fh, " cr OSF relB div flag\n");
for ( i=0; i<n_crystals; i++ ) {
- fprintf(fh, "%4i %5.2f %8.5e %s\n", i,
+ fprintf(fh, "%4i %10.5f %10.2f %8.5e %s\n", i,
crystal_get_osf(crystals[i]),
+ crystal_get_Bfac(crystals[i])*1e20,
crystal_get_image(crystals[i])->div,
str_flags(crystals[i]));
}
diff --git a/src/rejection.c b/src/rejection.c
new file mode 100644
index 00000000..d413666f
--- /dev/null
+++ b/src/rejection.c
@@ -0,0 +1,83 @@
+/*
+ * rejection.c
+ *
+ * Crystal rejection for 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 "crystal.h"
+#include "reflist.h"
+
+
+static double mean_intensity(RefList *list)
+{
+ Reflection *refl;
+ RefListIterator *iter;
+ double total = 0.0;
+ int n = 0;
+
+ for ( refl = first_refl(list, &iter);
+ refl != NULL;
+ refl = next_refl(refl, iter) )
+ {
+ total += get_intensity(refl);
+ n++;
+ }
+
+ return total/n;
+}
+
+
+/* Reject really obvious outliers */
+void early_rejection(Crystal **crystals, int n)
+{
+ int i;
+ double m = 0.0;
+ FILE *fh = fopen("reject.dat", "w");
+ int n_flag = 0;
+
+ for ( i=0; i<n; i++ ) {
+ double u;
+ RefList *list = crystal_get_reflections(crystals[i]);
+ u = mean_intensity(list);
+ m += u;
+ fprintf(fh, "%i %f\n", i, u);
+ if ( u < 100.0 ) {
+ crystal_set_user_flag(crystals[i], 5);
+ n_flag++;
+ }
+ }
+ fclose(fh);
+
+ STATUS("Mean intensity/peak = %f ADU\n", m/n);
+ STATUS("%i crystals flagged\n", n_flag);
+}
diff --git a/src/rejection.h b/src/rejection.h
new file mode 100644
index 00000000..891d6216
--- /dev/null
+++ b/src/rejection.h
@@ -0,0 +1,42 @@
+/*
+ * rejection.h
+ *
+ * Crystal rejection for 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 REJECTION_H
+#define REJECTION_H
+
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+
+#include "crystal.h"
+
+extern void early_rejection(Crystal **crystals, int n);
+
+#endif /* REJECTION_H */
generated by cgit v1.2.3 (git 2.25.1) at 2024-06-28 05:51:08 +0000