process_hkl: Add --reference and --outstream options

7 files changed, 113 insertions, 99 deletions

diff --git a/Makefile.am b/Makefile.am
index bf20389d..6dedd3c3 100644
--- a/Makefile.am
+++ b/Makefile.am
@@ -5,7 +5,7 @@ EXTRA_DIST = configure src/cell.h src/hdf5-file.h src/image.h \
              data/displaywindow.ui src/dirax.h src/peaks.h src/index.h \
              src/filters.h src/diffraction-gpu.h src/cl-utils.h \
              data/defs.h src/parameters-lcls.tmp \
-             data/diffraction.cl src/likelihood.h src/symmetry.h \
+             data/diffraction.cl src/symmetry.h \
              src/povray.h src/index-priv.h src/geometry.h src/templates.h \
              data/sfac/Ca.nff data/sfac/C.nff data/sfac/Fe.nff data/sfac/H.nff \
              data/sfac/Mg.nff data/sfac/N.nff data/sfac/O.nff data/sfac/P.nff \
diff --git a/Makefile.in b/Makefile.in
index 19434a30..ec902cc9 100644
--- a/Makefile.in
+++ b/Makefile.in
@@ -202,7 +202,7 @@ EXTRA_DIST = configure src/cell.h src/hdf5-file.h src/image.h \
              data/displaywindow.ui src/dirax.h src/peaks.h src/index.h \
              src/filters.h src/diffraction-gpu.h src/cl-utils.h \
              data/defs.h src/parameters-lcls.tmp \
-             data/diffraction.cl src/likelihood.h src/symmetry.h \
+             data/diffraction.cl src/symmetry.h \
              src/povray.h src/index-priv.h src/geometry.h src/templates.h \
              data/sfac/Ca.nff data/sfac/C.nff data/sfac/Fe.nff data/sfac/H.nff \
              data/sfac/Mg.nff data/sfac/N.nff data/sfac/O.nff data/sfac/P.nff \
diff --git a/src/Makefile.am b/src/Makefile.am
index 7fdf6dc9..66871249 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -17,7 +17,7 @@ endif
 pattern_sim_LDADD = @LIBS@
 
 process_hkl_SOURCES = process_hkl.c sfac.c statistics.c cell.c utils.c \
-                      reflections.c likelihood.c symmetry.c
+                      reflections.c symmetry.c
 process_hkl_LDADD = @LIBS@
 
 indexamajig_SOURCES = indexamajig.c hdf5-file.c utils.c cell.c image.c \
diff --git a/src/Makefile.in b/src/Makefile.in
index 4b4da1f7..2fa6d54c 100644
--- a/src/Makefile.in
+++ b/src/Makefile.in
@@ -116,7 +116,7 @@ powder_plot_OBJECTS = $(am_powder_plot_OBJECTS)
 powder_plot_DEPENDENCIES =
 am_process_hkl_OBJECTS = process_hkl.$(OBJEXT) sfac.$(OBJEXT) \
 	statistics.$(OBJEXT) cell.$(OBJEXT) utils.$(OBJEXT) \
-	reflections.$(OBJEXT) likelihood.$(OBJEXT) symmetry.$(OBJEXT)
+	reflections.$(OBJEXT) symmetry.$(OBJEXT)
 process_hkl_OBJECTS = $(am_process_hkl_OBJECTS)
 process_hkl_DEPENDENCIES =
 am_render_hkl_OBJECTS = render_hkl.$(OBJEXT) cell.$(OBJEXT) \
@@ -255,7 +255,7 @@ pattern_sim_SOURCES = pattern_sim.c diffraction.c utils.c image.c \
 	$(am__append_2)
 pattern_sim_LDADD = @LIBS@
 process_hkl_SOURCES = process_hkl.c sfac.c statistics.c cell.c utils.c \
-                      reflections.c likelihood.c symmetry.c
+                      reflections.c symmetry.c
 
 process_hkl_LDADD = @LIBS@
 indexamajig_SOURCES = indexamajig.c hdf5-file.c utils.c cell.c image.c \
@@ -423,7 +423,6 @@ distclean-compile:
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/image.Po@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/index.Po@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/indexamajig.Po@am__quote@
-@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/likelihood.Po@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/pattern_sim.Po@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/peaks.Po@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/povray.Po@am__quote@
diff --git a/src/likelihood.c b/src/likelihood.c
deleted file mode 100644
index cb5595f8..00000000
--- a/src/likelihood.c
+++ /dev/null
@@ -1,41 +0,0 @@
-/*
- * 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 scale_intensities(const double *model, ReflItemList *model_items,
-                       double *new_pattern, ReflItemList *new_items,
-                       double f0, int f0_valid)
-{
-	double s;
-	unsigned int i;
-	ReflItemList *items;
-
-	items = intersection_items(model_items, new_items);
-	s = stat_scale_intensity(model, new_pattern, items);
-	delete_items(items);
-	if ( f0_valid ) printf("%f %f\n", s, f0);
-
-	/* NaN -> abort */
-	if ( isnan(s) ) return;
-
-	/* Multiply the new pattern up by "s" */
-	for ( i=0; i<LIST_SIZE; i++ ) {
-		new_pattern[i] *= s;
-	}
-}
diff --git a/src/likelihood.h b/src/likelihood.h
deleted file mode 100644
index 70aef098..00000000
--- a/src/likelihood.h
+++ /dev/null
@@ -1,28 +0,0 @@
-/*
- * 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
-
-
-#include "utils.h"
-
-
-extern void scale_intensities(const double *model, ReflItemList *model_items,
-                              double *new_pattern, ReflItemList *new_items,
-                              double f0, int f0_valid);
-
-
-#endif	/* LIKELIHOOD_H */
diff --git a/src/process_hkl.c b/src/process_hkl.c
index 67e17212..d78cd2fc 100644
--- a/src/process_hkl.c
+++ b/src/process_hkl.c
@@ -25,7 +25,6 @@
 #include "statistics.h"
 #include "sfac.h"
 #include "reflections.h"
-#include "likelihood.h"
 #include "symmetry.h"
 
 
@@ -64,6 +63,10 @@ static void show_help(const char *s)
 "      --scale               Scale each pattern for best fit with the current\n"
 "                             model.\n"
 "  -y, --symmetry=<sym>      Merge according to point group <sym>.\n"
+"      --reference=<file>    Compare against intensities from <file> when\n"
+"                             scaling or resolving ambiguities.\n"
+"      --outstream=<file>    Write an annotated version of the input stream\n"
+"                             to <file>.\n"
 );
 }
 
@@ -205,22 +208,33 @@ static void merge_pattern(double *model, ReflItemList *observed,
                           const double *new,  ReflItemList *items,
                           unsigned int *model_counts,  int mo,
                           ReflItemList *twins,
-                          const char *holo, const char *mero, double *hist_vals,
+                          const char *holo, const char *mero,
+                          ReflItemList *reference, const double *reference_i,
+                          double *hist_vals,
                           signed int hist_h, signed int hist_k,
                           signed int hist_l, int *hist_n, double *devs,
-                          double *tots, double *means)
+                          double *tots, double *means, FILE *outfh)
 {
 	int i;
 	int twin;
 	ReflItemList *sym_items = new_items();
 
 	if ( twins != NULL ) {
-		twin = resolve_twin(model, observed, new, items,
-		                    twins, holo, mero);
+		if ( reference != NULL ) {
+			twin = resolve_twin(reference_i, reference, new, items,
+			                    twins, holo, mero);
+		} else {
+			twin = resolve_twin(model, observed, new, items,
+			                    twins, holo, mero);
+		}
 	} else {
 		twin = 0;
 	}
 
+	if ( outfh != NULL ) {
+		fprintf(outfh, "twin=%i\n", twin);
+	}
+
 	for ( i=0; i<num_items(items); i++ ) {
 
 		double intensity;
@@ -288,20 +302,46 @@ static void merge_pattern(double *model, ReflItemList *observed,
 }
 
 
+static void scale_intensities(const double *model, ReflItemList *model_items,
+                              double *new_pattern, ReflItemList *new_items,
+                              double f0, int f0_valid)
+{
+	double s;
+	unsigned int i;
+	ReflItemList *items;
+
+	items = intersection_items(model_items, new_items);
+	s = stat_scale_intensity(model, new_pattern, items);
+	delete_items(items);
+	if ( f0_valid ) printf("%f %f\n", s, f0);
+
+	/* NaN -> abort */
+	if ( isnan(s) ) return;
+
+	/* Multiply the new pattern up by "s" */
+	for ( i=0; i<LIST_SIZE; i++ ) {
+		new_pattern[i] *= s;
+	}
+}
+
+
 static void merge_all(FILE *fh, double **pmodel, ReflItemList **pobserved,
                       unsigned int **pcounts,
                       int config_maxonly, int config_scale, int config_sum,
                       int config_startafter, int config_stopafter,
                       ReflItemList *twins, const char *holo, const char *mero,
-                      int n_total_patterns, double *hist_vals,
+                      int n_total_patterns,
+                      ReflItemList *reference, double *reference_i,
+                      double *hist_vals,
                       signed int hist_h, signed int hist_k, signed int hist_l,
-                      int *hist_i, double *devs, double *tots, double *means)
+                      int *hist_i, double *devs, double *tots, double *means,
+                      FILE *outfh)
 {
 	char *rval;
 	float f0;
 	int n_nof0 = 0;
 	int f0_valid = 0;
-	int n_patterns = 0;
+	int n_patterns = 1;
 	double *new_pattern = new_list_intensity();
 	ReflItemList *items = new_items();
 	ReflItemList *observed = new_items();
@@ -336,14 +376,7 @@ static void merge_all(FILE *fh, double **pmodel, ReflItemList **pobserved,
 		int r;
 
 		rval = fgets(line, 1023, fh);
-		if ( (strncmp(line, "Reflections from indexing", 25) == 0)
-		    || (strncmp(line, "New pattern", 11) == 0) ) {
-
-			/* Start of first pattern? */
-			if ( n_patterns == 0 ) {
-				n_patterns++;
-				continue;
-			}
+		if ( strcmp(line, "\n") == 0 ) {
 
 			/* Assume a default I0 if we don't have one by now */
 			if ( config_scale && !f0_valid ) {
@@ -353,17 +386,25 @@ static void merge_all(FILE *fh, double **pmodel, ReflItemList **pobserved,
 
 			/* Scale if requested */
 			if ( config_scale ) {
-				scale_intensities(model, observed,
-				                  new_pattern, items,
-				                  f0, f0_valid);
+				if ( reference == NULL ) {
+					scale_intensities(model, observed,
+					                  new_pattern, items,
+					                  f0, f0_valid);
+				} else {
+					scale_intensities(reference_i,
+					                  reference,
+					                  new_pattern, items,
+					                  f0, f0_valid);
+				}
 			}
 
 			/* Start of second or later pattern */
 			merge_pattern(model, observed, new_pattern, items,
 			              counts, config_maxonly,
 			              twins, holo, mero,
+			              reference, reference_i,
 			              hist_vals, hist_h, hist_k, hist_l,
-			              hist_i, devs, tots, means);
+			              hist_i, devs, tots, means, outfh);
 
 			if ( n_patterns == config_stopafter ) break;
 
@@ -377,6 +418,10 @@ static void merge_all(FILE *fh, double **pmodel, ReflItemList **pobserved,
 
 		}
 
+		if ( outfh != NULL ) {
+			fprintf(outfh, "%s", line);
+		}
+
 		if ( strncmp(line, "f0 = ", 5) == 0 ) {
 			r = sscanf(line, "f0 = %f", &f0);
 			if ( r != 1 ) {
@@ -472,6 +517,11 @@ int main(int argc, char *argv[])
 	signed int hist_h, hist_k, hist_l;
 	double *hist_vals = NULL;
 	int hist_i;
+	char *outstream = NULL;
+	char *reference = NULL;
+	ReflItemList *reference_items;
+	double *reference_i;
+	FILE *outfh = NULL;
 
 	/* Long options */
 	const struct option longopts[] = {
@@ -488,11 +538,13 @@ int main(int argc, char *argv[])
 		{"pdb",                1, NULL,               'p'},
 		{"histogram",          1, NULL,               'g'},
 		{"rmerge",             0, &config_rmerge,      1},
+		{"outstream",          1, NULL,               'a'},
+		{"reference",          1, NULL,               'r'},
 		{0, 0, NULL, 0}
 	};
 
 	/* Short options */
-	while ((c = getopt_long(argc, argv, "hi:e:ro:p:y:g:f:",
+	while ((c = getopt_long(argc, argv, "hi:e:ro:p:y:g:f:a:r:",
 	                        longopts, NULL)) != -1) {
 
 		switch (c) {
@@ -528,6 +580,14 @@ int main(int argc, char *argv[])
 			histo = strdup(optarg);
 			break;
 
+		case 'r' :
+			reference = strdup(optarg);
+			break;
+
+		case 'a' :
+			outstream = strdup(optarg);
+			break;
+
 		case 0 :
 			break;
 
@@ -612,6 +672,28 @@ int main(int argc, char *argv[])
 		return 1;
 	}
 
+	/* Read the reference reflections */
+	if ( reference != NULL ) {
+		reference_i = new_list_intensity();
+		reference_items = read_reflections(reference, reference_i,
+		                                   NULL, NULL);
+		if ( reference_items == NULL ) {
+			ERROR("Couldn't read '%s'\n", reference);
+			return 1;
+		}
+	} else {
+		reference_items = NULL;
+		reference_i = NULL;
+	}
+
+	if ( outstream != NULL ) {
+		outfh = fopen(outstream, "w");
+		if ( outfh == NULL ) {
+			ERROR("Couldn't open '%s'\n", outstream);
+			return 1;
+		}
+	}
+
 	/* Count the number of patterns in the file */
 	n_total_patterns = count_patterns(fh);
 	STATUS("There are %i patterns to process\n", n_total_patterns);
@@ -622,7 +704,9 @@ int main(int argc, char *argv[])
 	          config_maxonly, config_scale, config_sum,
 	          config_startafter, config_stopafter,
                   twins, holo, sym, n_total_patterns,
-                  hist_vals, hist_h, hist_k, hist_l, &hist_i, NULL, NULL, NULL);
+                  reference_items, reference_i,
+                  hist_vals, hist_h, hist_k, hist_l, &hist_i, NULL, NULL, NULL,
+                  outfh);
 	rewind(fh);
 
 	if ( hist_vals != NULL ) {
@@ -651,8 +735,8 @@ int main(int argc, char *argv[])
 		merge_all(fh, &model, &observed, &counts,
 		          config_maxonly, config_scale, 0,
 		          config_startafter, config_stopafter, twins, holo, sym,
-		          n_total_patterns,
-		          NULL, 0, 0, 0, NULL, devs, tots, model);
+		          n_total_patterns, reference_items, reference_i,
+		          NULL, 0, 0, 0, NULL, devs, tots, model, NULL);
 
 		for ( i=0; i<num_items(observed); i++ ) {