process_hkl: Implement --scale option

7 files changed, 129 insertions, 38 deletions

diff --git a/src/compare_hkl.c b/src/compare_hkl.c
index bb07d17e..e7752078 100644
--- a/src/compare_hkl.c
+++ b/src/compare_hkl.c
@@ -53,7 +53,6 @@ int main(int argc, char *argv[])
 	double scale, R;
 	unsigned int *c1;
 	unsigned int *c2;
-	unsigned int *cjoint;
 	int i;
 	int nc1, nc2, ncom;
 
@@ -134,18 +133,16 @@ int main(int argc, char *argv[])
 	}
 	}
 
-	cjoint = new_list_count();
 	nc1 = 0;
 	nc2 = 0;
 	ncom = 0;
 	for ( i=0; i<IDIM*IDIM*IDIM; i++ ) {
-		cjoint[i] = c1[i] && c2[i];
 		nc1 += c1[i];
 		nc2 += c2[i];
-		ncom += cjoint[i];
+		ncom += c1[i] && c2[i];
 	}
 	STATUS("%i,%i reflections: %i in common\n", nc1, nc2, ncom);
-	R = stat_r2(ref1, ref2, cjoint, IDIM*IDIM*IDIM, &scale);
+	R = stat_r2(ref1, c1, ref2, c2, &scale);
 	STATUS("R2 = %5.4f %% (scale=%5.2f)\n", R*100.0, scale);
 
 	if ( outfile != NULL ) {
diff --git a/src/likelihood.c b/src/likelihood.c
new file mode 100644
index 00000000..e8c62010
--- /dev/null
+++ b/src/likelihood.c
@@ -0,0 +1,44 @@
+/*
+ * likelihood.c
+ *
+ * Likelihood maximisation
+ *
+ * (c) 2006-2010 Thomas White <taw@physics.org>
+ *
+ * Part of CrystFEL - crystallography with a FEL
+ *
+ */
+
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include "statistics.h"
+#include "utils.h"
+
+void detwin_intensities(const double *model, double *new_pattern,
+                        const unsigned int *model_counts,
+                        unsigned int *new_counts)
+{
+	/* Placeholder... */
+}
+
+void scale_intensities(const double *model, double *new_pattern,
+                       const unsigned int *model_counts,
+                       unsigned int *new_counts)
+{
+	double s;
+	unsigned int i;
+
+	s = stat_scale_intensity(model, model_counts, new_pattern, new_counts);
+	printf("%f\n", s);
+
+	/* NaN -> abort */
+	if ( isnan(s) ) return;
+
+	/* Multiply the new pattern up by "s" */
+	for ( i=0; i<LIST_SIZE; i++ ) {
+		new_counts[i] *= s;
+	}
+}
diff --git a/src/likelihood.h b/src/likelihood.h
new file mode 100644
index 00000000..9a773829
--- /dev/null
+++ b/src/likelihood.h
@@ -0,0 +1,29 @@
+/*
+ * likelihood.h
+ *
+ * Likelihood maximisation
+ *
+ * (c) 2006-2010 Thomas White <taw@physics.org>
+ *
+ * Part of CrystFEL - crystallography with a FEL
+ *
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#ifndef LIKELIHOOD_H
+#define LIKELIHOOD_H
+
+
+extern void detwin_intensities(const double *model, double *new_pattern,
+	                       const unsigned int *model_counts,
+                               unsigned int *new_counts);
+
+extern void scale_intensities(const double *model, double *new_pattern,
+	                      const unsigned int *model_counts,
+                              unsigned int *new_counts);
+
+
+#endif	/* LIKELIHOOD_H */
diff --git a/src/list_tmp.h b/src/list_tmp.h
index 2b8f703d..749d7e7a 100644
--- a/src/list_tmp.h
+++ b/src/list_tmp.h
@@ -18,6 +18,8 @@
  *
  */
 
+#include <stdio.h>
+
 #define ERROR_T(...) fprintf(stderr, __VA_ARGS__)
 
 static inline void LABEL(integrate)(TYPE *ref, signed int h,
diff --git a/src/process_hkl.c b/src/process_hkl.c
index dfb84ee1..d09fca22 100644
--- a/src/process_hkl.c
+++ b/src/process_hkl.c
@@ -66,6 +66,8 @@ static void show_help(const char *s)
 "      --detwin              Correlate each new pattern with the current\n"
 "                             model and choose the best fitting out of the\n"
 "                             allowable twins.\n"
+"      --scale               Scale each pattern for best fit with the current\n"
+"                             model.\n"
 );
 }
 
@@ -141,8 +143,9 @@ static void write_RvsQ(const char *name, double *ref, double *trueref,
 }
 
 
-static void process_reflections(double *ref, double *trueref,
-                                unsigned int *counts, unsigned int n_patterns,
+static void process_reflections(double *ref, unsigned int *counts,
+                                double *trueref, unsigned int *truecounts,
+                                unsigned int n_patterns,
                                 UnitCell *cell, int do_rvsq, int do_zoneaxis)
 {
 	int j;
@@ -158,7 +161,7 @@ static void process_reflections(double *ref, double *trueref,
 	}
 	mean_counts = (double)ctot/nmeas;
 
-	R = stat_r2(ref, trueref, counts, LIST_SIZE, &scale);
+	R = stat_r2(ref, counts, trueref, truecounts, &scale);
 	STATUS("%8u: R=%5.2f%% (sf=%7.4e)  mean meas/refl=%5.2f,"
 	       " %i reflections measured\n",
 	       n_patterns, R*100.0, scale, mean_counts, nmeas);
@@ -236,10 +239,12 @@ int main(int argc, char *argv[])
 	int config_zoneaxis = 0;
 	int config_sum = 0;
 	int config_detwin = 0;
+	int config_scale = 0;
 	char *intfile = NULL;
 	double *new_pattern = NULL;
 	unsigned int *new_counts = NULL;
 	unsigned int n_total_patterns;
+	unsigned int *truecounts = NULL;
 
 	/* Long options */
 	const struct option longopts[] = {
@@ -254,6 +259,7 @@ int main(int argc, char *argv[])
 		{"compare-with",       0, NULL,               'c'},
 		{"sum",                0, &config_sum,         1},
 		{"detwin",             0, &config_detwin,      1},
+		{"scale",              0, &config_scale,      1},
 		{0, 0, NULL, 0}
 	};
 
@@ -313,8 +319,9 @@ int main(int argc, char *argv[])
 	}
 
 	if ( intfile != NULL ) {
+		truecounts = new_list_count();
 		STATUS("Comparing against '%s'\n", intfile);
-		trueref = read_reflections(intfile, NULL);
+		trueref = read_reflections(intfile, truecounts);
 		free(intfile);
 	} else {
 		trueref = NULL;
@@ -369,19 +376,27 @@ int main(int argc, char *argv[])
 				continue;
 			}
 
+			/* Detwin before scaling */
 			if ( config_detwin ) {
 				detwin_intensities(model, new_pattern,
 				                   model_counts, new_counts);
 			}
 
+			/* Scale if requested */
+			if ( config_scale ) {
+				scale_intensities(model, new_pattern,
+				                  model_counts, new_counts);
+			}
+
 			/* Start of second or later pattern */
 			merge_pattern(model, new_pattern, model_counts,
 			              new_counts, config_maxonly, config_sum);
 
 			if (config_every && (n_patterns % config_every == 0)) {
-				process_reflections(model, trueref,
-				                    model_counts, n_patterns,
-				                    cell, config_rvsq,
+				process_reflections(model, model_counts,
+				                    trueref, truecounts,
+				                    n_patterns, cell,
+				                    config_rvsq,
 				                    config_zoneaxis);
 			}
 
@@ -412,8 +427,9 @@ int main(int argc, char *argv[])
 	fclose(fh);
 
 	if ( trueref != NULL ) {
-		process_reflections(model, trueref, model_counts, n_patterns,
-		                    cell, config_rvsq, config_zoneaxis);
+		process_reflections(model, model_counts, trueref, truecounts,
+		                    n_patterns, cell, config_rvsq,
+		                    config_zoneaxis);
 	}
 
 	if ( output != NULL ) {
diff --git a/src/statistics.c b/src/statistics.c
index be499814..406f1f56 100644
--- a/src/statistics.c
+++ b/src/statistics.c
@@ -18,26 +18,26 @@
 #include <stdlib.h>
 
 #include "statistics.h"
+#include "utils.h"
 
 
-/* By what (best-fitted) factor must the list "list" be multiplied by,
- * if it were to be merged with "target"? */
-static double stat_scale_intensity(double *obs, double *calc, unsigned int *c,
-                                   int size)
+double stat_scale_intensity(const double *ref1, const unsigned int *c1,
+                            const double *ref2, const unsigned int *c2)
 {
 	double top = 0.0;
 	double bot = 0.0;
 	int i;
 
 	/* Start from i=1 -> skip central beam */
-	for ( i=1; i<size; i++ ) {
+	for ( i=1; i<LIST_SIZE; i++ ) {
 
-		if ( c[i] > 0 ) {
-			double obsi;
-			obsi = obs[i] / (double)c[i];
-			top += obsi * calc[i];
-			bot += calc[i] * calc[i];
-		} /* else reflection not measured so don't worry about it */
+		if ( c1[i] && c2[i] ) {
+			double i1, i2;
+			i1 = ref1[i] / (double)c1[i];
+			i2 = ref2[i] / (double)c2[i];
+			top += i1 * i2;
+			bot += i2 * i2;
+		} /* else reflection not common so don't worry about it */
 
 	}
 
@@ -45,28 +45,27 @@ static double stat_scale_intensity(double *obs, double *calc, unsigned int *c,
 }
 
 
-double stat_r2(double *obs, double *calc, unsigned int *c, int size,
-               double *scalep)
+double stat_r2(const double *ref1, const unsigned int *c1,
+               const double *ref2, const unsigned int *c2, double *scalep)
 {
 	double top = 0.0;
 	double bot = 0.0;
 	double scale;
 	int i;
-	scale = stat_scale_intensity(obs, calc, c, size);
+	scale = stat_scale_intensity(ref1, c1, ref2, c2);
 	*scalep = scale;
 
 	/* Start from i=1 -> skip central beam */
-	for ( i=1; i<size; i++ ) {
-
-		if ( c[i] > 0 ) {
+	for ( i=1; i<LIST_SIZE; i++ ) {
 
-			double obsi;
+		if ( c1[i] && c2[i] ) {
 
-			obsi = obs[i] / (double)c[i];
-			obsi = obsi / scale;
+			double i1, i2;
+			i1 = ref1[i] / (scale*(double)c1[i]);
+			i2 = ref2[i] / (scale*(double)c2[i]);
 
-			top += pow(fabs(obsi - calc[i]), 2.0);
-			bot += pow(obsi, 2.0);
+			top += pow(fabs(i1 - i2), 2.0);
+			bot += pow(i1, 2.0);
 
 		} /* else reflection not measured so don't worry about it */
 
diff --git a/src/statistics.h b/src/statistics.h
index 5db7b045..fe7f57e1 100644
--- a/src/statistics.h
+++ b/src/statistics.h
@@ -16,7 +16,11 @@
 #ifndef STATISTICS_H
 #define STATISTICS_H
 
-double stat_r2(double *obs, double *calc, unsigned int *c, int size,
-               double *scalep);
+extern double stat_scale_intensity(const double *ref1, const unsigned int *c1,
+                                   const double *ref2, const unsigned int *c2);
+
+extern double stat_r2(const double *ref1, const unsigned int *c1,
+                      const double *ref2, const unsigned int *c2,
+                      double *scalep);
 
 #endif	/* STATISTICS_H */