From b57deca8061316d2fc09dedcf8a91fa7523f081f Mon Sep 17 00:00:00 2001
From: Thomas White <taw@physics.org>
Date: Mon, 1 Sep 2014 14:54:09 +0200
Subject: WIP on prediction

---
 libcrystfel/src/geometry.c |  5 +++-
 src/process_image.c        | 73 +++++++++++++++++++++++++++++++++++++++++++++-
 2 files changed, 76 insertions(+), 2 deletions(-)

diff --git a/libcrystfel/src/geometry.c b/libcrystfel/src/geometry.c
index 8952486d..250f889c 100644
--- a/libcrystfel/src/geometry.c
+++ b/libcrystfel/src/geometry.c
@@ -151,6 +151,7 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst,
 	double pr;
 	double D;
 	double del;
+	double rlowuc, rhighuc;  /* Values before clamping */
 
 	/* Don't predict 000 */
 	if ( abs(h)+abs(k)+abs(l) == 0 ) return NULL;
@@ -197,6 +198,8 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst,
 
 	D = rlow - rhigh;
 
+	rlowuc = rlow;  rhighuc = rhigh;
+
 	/* If the "lower" Ewald sphere is a long way away, use the
 	 * position at which the Ewald sphere would just touch the
 	 * reflection.
@@ -243,7 +246,7 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst,
 		set_detector_pos(refl, 0.0, xda, yda);
 	}
 
-	set_partial(refl, rlow, rhigh, part, clamp_low, clamp_high);
+	set_partial(refl, rlowuc, rhighuc, part, clamp_low, clamp_high);
 	set_lorentz(refl, 1.0);
 	set_symmetric_indices(refl, h, k, l);
 	set_redundancy(refl, 1);
diff --git a/src/process_image.c b/src/process_image.c
index 1d1dac82..f4392efe 100644
--- a/src/process_image.c
+++ b/src/process_image.c
@@ -34,6 +34,8 @@
 #include <stdlib.h>
 #include <hdf5.h>
 #include <gsl/gsl_errno.h>
+#include <gsl/gsl_statistics_double.h>
+#include <gsl/gsl_sort.h>
 #include <unistd.h>
 
 #include "utils.h"
@@ -50,6 +52,63 @@
 #include "integration.h"
 
 
+static int cmpd2(const void *av, const void *bv)
+{
+	double *ap, *bp;
+	double a, b;
+
+	ap = (double *)av;
+	bp = (double *)bv;
+
+	a = ap[1];
+	b = bp[1];
+
+	if ( a > b ) return -1;
+	return 1;
+}
+
+
+static void refine_radius(Crystal *cr)
+{
+	Reflection *refl;
+	RefListIterator *iter;
+	FILE *fh;
+	double vals[num_reflections(crystal_get_reflections(cr))*2];
+	int n = 0;
+	int i;
+
+	fh = fopen("graph.dat", "w");
+
+	for ( refl = first_refl(crystal_get_reflections(cr), &iter);
+	      refl != NULL;
+	      refl = next_refl(refl, iter) )
+	{
+		double i = get_intensity(refl);
+		double rlow, rhigh, p;
+		int cl, ch;
+
+		get_partial(refl, &rlow, &rhigh, &p, &cl, &ch);
+		fprintf(fh, "%e %10f %e %e %f %i %i\n", (rhigh+rlow)/2.0, i,
+		                           rlow, rhigh, p, cl, ch);
+
+		vals[(2*n)+0] = i;
+		vals[(2*n)+1] = (rhigh+rlow)/2.0;
+		n++;
+
+	}
+
+	qsort(vals, n, sizeof(double)*2, cmpd2);
+
+	for ( i=0; i<n-1; i++ ) {
+		double mean = gsl_stats_mean(vals, 2, i);
+		double var = gsl_stats_variance_m(vals, 2, i, mean);
+		printf("%.2f  %i\n", mean/var, i);
+	}
+
+	fclose(fh);
+}
+
+
 void process_image(const struct index_args *iargs, struct pattern_args *pargs,
                    Stream *st, int cookie, const char *tmpdir, int results_pipe)
 {
@@ -169,7 +228,7 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs,
 
 	/* Default parameters */
 	image.div = 0.0;
-	image.bw = 0.001;
+	image.bw = 0.00000001;
 	STATUS("Warning: div, bw and pr are hardcoded in this version.\n");
 	for ( i=0; i<image.n_crystals; i++ ) {
 		crystal_set_profile_radius(image.crystals[i], 0.01e9);
@@ -183,6 +242,18 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs,
 	                iargs->int_diag, iargs->int_diag_h,
 	                iargs->int_diag_k, iargs->int_diag_l, results_pipe);
 
+	for ( i=0; i<image.n_crystals; i++ ) {
+		refine_radius(image.crystals[i]);
+		reflist_free(crystal_get_reflections(image.crystals[i]));
+	}
+
+	integrate_all_4(&image, iargs->int_meth, PMODEL_SCSPHERE,
+		                iargs->push_res,
+			        iargs->ir_inn, iargs->ir_mid, iargs->ir_out,
+			        iargs->int_diag, iargs->int_diag_h,
+			        iargs->int_diag_k, iargs->int_diag_l,
+			        results_pipe);
+
 	write_chunk(st, &image, hdfile,
 	            iargs->stream_peaks, iargs->stream_refls,
 	            pargs->filename_p_e->ev);
-- 
cgit v1.2.3