process_hkl: Restructure and remove all intensity analysis stuff

3 files changed, 130 insertions, 284 deletions

diff --git a/src/likelihood.c b/src/likelihood.c
index 417ebea3..b6a37994 100644
--- a/src/likelihood.c
+++ b/src/likelihood.c
@@ -17,12 +17,6 @@
 #include "statistics.h"
 #include "utils.h"
 
-void detwin_intensities(const double *model, double *new_pattern,
-                        const unsigned int *model_counts,
-                        ReflItemList *items)
-{
-	/* Placeholder... */
-}
 
 void scale_intensities(const double *model, double *new_pattern,
                        const unsigned int *model_counts,
diff --git a/src/likelihood.h b/src/likelihood.h
index ad6585eb..08956fe3 100644
--- a/src/likelihood.h
+++ b/src/likelihood.h
@@ -20,10 +20,6 @@
 #include "utils.h"
 
 
-extern void detwin_intensities(const double *model, double *new_pattern,
-	                       const unsigned int *model_counts,
-                               ReflItemList *items);
-
 extern void scale_intensities(const double *model, double *new_pattern,
 	                      const unsigned int *model_counts,
                               ReflItemList *items, double f0,
diff --git a/src/process_hkl.c b/src/process_hkl.c
index 29288154..b5b6df9a 100644
--- a/src/process_hkl.c
+++ b/src/process_hkl.c
@@ -56,16 +56,7 @@ static void show_help(const char *s)
 "      --stop-after=<n>      Stop after processing n patterns.  Zero means\n"
 "                             keep going until the end of the input, and is\n"
 "                             the default.\n"
-"  -c, --compare-with=<file> Compare with reflection intensities in this file\n"
 "\n"
-"  -e, --output-every=<n>    Analyse figures of merit after every n patterns\n"
-"                             Default: 1000.  A value of zero means to do the\n"
-"                             analysis only after reading all the patterns.\n"
-"  -r, --rvsq                Output lists of R vs |q| (\"Luzzatti plots\")\n"
-"                             when analysing figures of merit.\n"
-"      --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"
 "  -y, --symmetry=<sym>      Merge according to point group <sym>.\n"
@@ -73,113 +64,6 @@ static void show_help(const char *s)
 }
 
 
-static void write_RvsQ(const char *name, double *ref, double *trueref,
-                       unsigned int *counts, double scale, UnitCell *cell)
-{
-	FILE *fh;
-	double smax, sbracket;
-	signed int h, k, l;
-
-	fh = fopen(name, "w");
-
-	smax = 0.0;
-	for ( h=-INDMAX; h<INDMAX; h++ ) {
-	for ( k=-INDMAX; k<INDMAX; k++ ) {
-	for ( l=-INDMAX; l<INDMAX; l++ ) {
-		double s = 2.0*resolution(cell, h, k, l);
-		if ( (lookup_count(counts, h, k, l) > 0) && (s > smax) ) {
-			smax = s;
-		}
-	}
-	}
-	}
-
-	for ( sbracket=0.0; sbracket<smax; sbracket+=smax/RVQDV ) {
-
-		double top = 0.0;
-		double bot = 0.0;
-		int nhits = 0;
-		int nrefl = 0;
-		double R;
-		double hits_per_refl;
-
-		for ( h=-INDMAX; h<INDMAX; h++ ) {
-		for ( k=-INDMAX; k<INDMAX; k++ ) {
-		for ( l=-INDMAX; l<INDMAX; l++ ) {
-
-			double s;
-			int c;
-
-			if ( (h==0) && (k==0) && (l==0) ) continue;
-
-			c = lookup_count(counts, h, k, l);
-			s = 2.0*resolution(cell, h, k, l);
-			if ((s>=sbracket) && (s<sbracket+smax/RVQDV) && (c>0)) {
-
-				double obs, calc, obsi;
-
-				obs = lookup_intensity(ref, h, k, l);
-				calc = lookup_intensity(trueref, h, k, l);
-
-				obsi = obs / (double)c;
-				top += pow(obsi - scale*calc, 2.0);
-				bot += pow(obsi, 2.0);
-				nhits += c;
-				nrefl++;
-
-			}
-
-		}
-		}
-		}
-
-		R = sqrt(top/bot);
-		hits_per_refl = nrefl ? (double)nhits/nrefl : 0;
-
-		fprintf(fh, "%8.5f %8.5f %5.2f\n", sbracket+smax/(2.0*RVQDV),
-		                                R, hits_per_refl);
-
-	}
-	fclose(fh);
-}
-
-
-static void process_reflections(double *ref, unsigned int *counts,
-                                double *trueref, unsigned int *truecounts,
-                                unsigned int n_patterns,
-                                UnitCell *cell, int do_rvsq)
-{
-	int j;
-	double mean_counts;
-	int ctot = 0;
-	int nmeas = 0;
-	double R, scale;
-	FILE *fh;
-
-	for ( j=0; j<LIST_SIZE; j++ ) {
-		ctot += counts[j];
-		if ( counts[j] > 0 ) nmeas++;
-	}
-	mean_counts = (double)ctot/nmeas;
-
-	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);
-
-	if ( do_rvsq ) {
-		/* Record graph of R against q for this N */
-		char name[64];
-		snprintf(name, 63, "R_vs_q-%u.dat", n_patterns);
-		write_RvsQ(name, ref, trueref, counts, scale, cell);
-	}
-
-	fh = fopen("results/convergence.dat", "a");
-	fprintf(fh, "%u %5.2f %5.2f\n", n_patterns, R*100.0, mean_counts);
-	fclose(fh);
-}
-
-
 /* Note "holo" needn't actually be a holohedral point group, if you want to try
  * something strange like resolving from a low-symmetry group into an even
  * lower symmetry one.
@@ -321,34 +205,132 @@ static void merge_pattern(double *model, const double *new,
 }
 
 
+static void merge_all(FILE *fh, double *model, unsigned int *model_counts,
+                      int config_maxonly, int config_scale, int config_sum,
+                      int config_stopafter,
+                      ReflItemList *twins, const char *holo, const char *mero,
+                      int n_total_patterns)
+{
+	char *rval;
+	float f0;
+	int n_nof0 = 0;
+	int f0_valid = 0;
+	int n_patterns = 0;
+	double *new_pattern = new_list_intensity();
+	ReflItemList *items = new_items();
+
+	do {
+
+		char line[1024];
+		signed int h, k, l;
+		float intensity;
+		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;
+			}
+
+			/* Assume a default I0 if we don't have one by now */
+			if ( config_scale && !f0_valid ) {
+				n_nof0++;
+				f0 = 1.0;
+			}
+
+			/* Scale if requested */
+			if ( config_scale ) {
+				scale_intensities(model, new_pattern,
+				                  model_counts, items, f0,
+				                  f0_valid);
+			}
+
+			/* Start of second or later pattern */
+			merge_pattern(model, new_pattern, model_counts,
+			              items, config_maxonly, config_sum, twins,
+			              holo, mero);
+
+			if ( n_patterns == config_stopafter ) break;
+
+			n_patterns++;
+			clear_items(items);
+
+			progress_bar(n_patterns, n_total_patterns, "Merging");
+
+			f0_valid = 0;
+
+		}
+
+		if ( strncmp(line, "f0 = ", 5) == 0 ) {
+			r = sscanf(line, "f0 = %f", &f0);
+			if ( r != 1 ) {
+				f0 = 1.0;
+				f0_valid = 0;
+				continue;
+			}
+			f0_valid = 1;
+		}
+
+		r = sscanf(line, "%i %i %i %f", &h, &k, &l, &intensity);
+		if ( r != 4 ) continue;
+
+		if ( (h==0) && (k==0) && (l==0) ) continue;
+
+		if ( find_item(items, h, k, l) != 0 ) {
+			ERROR("More than one measurement for %i %i %i in"
+			      " pattern number %i\n", h, k, l, n_patterns);
+		}
+		set_intensity(new_pattern, h, k, l, intensity);
+		add_item(items, h, k, l);
+
+	} while ( rval != NULL );
+
+	delete_items(items);
+	free(new_pattern);
+
+	STATUS("%i patterns had no f0 valid value.\n", n_nof0);
+}
+
+
+static int count_patterns(FILE *fh)
+{
+	char *rval;
+
+	int n_total_patterns = 0;
+	do {
+		char line[1024];
+
+		rval = fgets(line, 1023, fh);
+		if ( (strncmp(line, "Reflections from indexing", 25) == 0)
+		    || (strncmp(line, "New pattern", 11) == 0) ) {
+		    n_total_patterns++;
+		}
+	} while ( rval != NULL );
+
+	return n_total_patterns;
+}
+
+
 int main(int argc, char *argv[])
 {
 	int c;
 	char *filename = NULL;
 	char *output = NULL;
 	FILE *fh;
-	unsigned int n_patterns;
-	double *model, *trueref = NULL;
+	double *model;
 	unsigned int *model_counts;
-	char *rval;
 	UnitCell *cell;
 	int config_maxonly = 0;
-	int config_every = 1000;
-	int config_rvsq = 0;
 	int config_stopafter = 0;
 	int config_sum = 0;
-	int config_detwin = 0;
 	int config_scale = 0;
-	char *intfile = NULL;
-	double *new_pattern = NULL;
 	unsigned int n_total_patterns;
-	unsigned int *truecounts = NULL;
 	char *sym = NULL;
 	char *pdb = NULL;
-	float f0;
-	int f0_valid;
-	int n_nof0 = 0;
-	ReflItemList *items;
 	ReflItemList *twins;
 	int i;
 
@@ -359,11 +341,8 @@ int main(int argc, char *argv[])
 		{"output",             1, NULL,               'o'},
 		{"max-only",           0, &config_maxonly,     1},
 		{"output-every",       1, NULL,               'e'},
-		{"rvsq",               0, NULL,               'r'},
 		{"stop-after",         1, NULL,               's'},
-		{"compare-with",       0, NULL,               'c'},
 		{"sum",                0, &config_sum,         1},
-		{"detwin",             0, &config_detwin,      1},
 		{"scale",              0, &config_scale,       1},
 		{"symmetry",           1, NULL,               'y'},
 		{"pdb",                1, NULL,               'p'},
@@ -387,23 +366,10 @@ int main(int argc, char *argv[])
 			output = strdup(optarg);
 			break;
 
-
-		case 'r' :
-			config_rvsq = 1;
-			break;
-
-		case 'e' :
-			config_every = atoi(optarg);
-			break;
-
 		case 's' :
 			config_stopafter = atoi(optarg);
 			break;
 
-		case 'c' :
-			intfile = strdup(optarg);
-			break;
-
 		case 'p' :
 			pdb = strdup(optarg);
 			break;
@@ -426,15 +392,6 @@ int main(int argc, char *argv[])
 		return 1;
 	}
 
-	if ( intfile != NULL ) {
-		truecounts = new_list_count();
-		STATUS("Comparing against '%s'\n", intfile);
-		trueref = read_reflections(intfile, truecounts, NULL);
-		free(intfile);
-	} else {
-		trueref = NULL;
-	}
-
 	if ( pdb == NULL ) {
 		pdb = strdup("molecule.pdb");
 	}
@@ -445,37 +402,10 @@ int main(int argc, char *argv[])
 
 	model = new_list_intensity();
 	model_counts = new_list_count();
-	new_pattern = new_list_intensity();
-	items = new_items();
 
 	cell = load_cell_from_pdb(pdb);
 	free(pdb);
 
-	if ( strcmp(filename, "-") == 0 ) {
-		fh = stdin;
-	} else {
-		fh = fopen(filename, "r");
-	}
-	free(filename);
-	if ( fh == NULL ) {
-		ERROR("Failed to open input file\n");
-		return 1;
-	}
-
-	/* Count the number of patterns in the file */
-	n_total_patterns = 0;
-	do {
-		char line[1024];
-
-		rval = fgets(line, 1023, fh);
-		if ( (strncmp(line, "Reflections from indexing", 25) == 0)
-		    || (strncmp(line, "New pattern", 11) == 0) ) {
-		    n_total_patterns++;
-		}
-	} while ( rval != NULL );
-	rewind(fh);
-	STATUS("There are %i patterns to process\n", n_total_patterns);
-
 	/* Show useful symmetry information */
 	const char *holo = get_holohedral(sym);
 	int np = num_general_equivs(holo) / num_general_equivs(sym);
@@ -496,112 +426,38 @@ int main(int argc, char *argv[])
 		twins = NULL;
 	}
 
-	n_patterns = 0;
-	f0_valid = 0;
-	do {
-
-		char line[1024];
-		signed int h, k, l;
-		float intensity;
-		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;
-			}
-
-			/* Detwin before scaling */
-			if ( config_detwin ) {
-				detwin_intensities(model, new_pattern,
-				                   model_counts, items);
-			}
-
-			/* Assume a default I0 if we don't have one by now */
-			if ( config_scale && !f0_valid ) {
-				n_nof0++;
-				f0 = 1.0;
-			}
-
-			/* Scale if requested */
-			if ( config_scale ) {
-				scale_intensities(model, new_pattern,
-				                  model_counts, items, f0,
-				                  f0_valid);
-			}
-
-			/* Start of second or later pattern */
-			merge_pattern(model, new_pattern, model_counts,
-			              items, config_maxonly, config_sum, twins,
-			              holo, sym);
-
-			if ( (trueref != NULL) && config_every
-			    && (n_patterns % config_every == 0) ) {
-				process_reflections(model, model_counts,
-				                    trueref, truecounts,
-				                    n_patterns, cell,
-				                    config_rvsq);
-			}
-
-			if ( n_patterns == config_stopafter ) break;
-
-			n_patterns++;
-			clear_items(items);
-
-			progress_bar(n_patterns, n_total_patterns, "Merging");
-
-			f0_valid = 0;
-
-		}
-
-		if ( strncmp(line, "f0 = ", 5) == 0 ) {
-			r = sscanf(line, "f0 = %f", &f0);
-			if ( r != 1 ) {
-				f0 = 1.0;
-				f0_valid = 0;
-				continue;
-			}
-			f0_valid = 1;
-		}
-
-		r = sscanf(line, "%i %i %i %f", &h, &k, &l, &intensity);
-		if ( r != 4 ) continue;
-
-		if ( (h==0) && (k==0) && (l==0) ) continue;
+	/* Open the data stream */
+	if ( strcmp(filename, "-") == 0 ) {
+		fh = stdin;
+	} else {
+		fh = fopen(filename, "r");
+	}
+	free(filename);
+	if ( fh == NULL ) {
+		ERROR("Failed to open input file\n");
+		return 1;
+	}
 
-		if ( find_item(items, h, k, l) != 0 ) {
-			ERROR("More than one measurement for %i %i %i in"
-			      " pattern number %i\n", h, k, l, n_patterns);
-		}
-		set_intensity(new_pattern, h, k, l, intensity);
-		add_item(items, h, k, l);
+	/* 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);
+	rewind(fh);
 
-	} while ( rval != NULL );
+	merge_all(fh, model, model_counts,
+	          config_maxonly, config_scale, config_sum, config_stopafter,
+                  twins, holo, sym, n_total_patterns);
+	rewind(fh);
 
 	fclose(fh);
 
-	if ( trueref != NULL ) {
-		process_reflections(model, model_counts, trueref, truecounts,
-		                    n_patterns, cell, config_rvsq);
-	}
-
 	if ( output != NULL ) {
 		write_reflections(output, model_counts, model, NULL,
 		                  0, cell, 1);
 	}
 
-	STATUS("There were %u patterns.\n", n_patterns);
-	STATUS("%i had no f0 valid value.\n", n_nof0);
-
-	delete_items(items);
 	free(sym);
 	free(model);
 	free(model_counts);
-	free(new_pattern);
 	free(output);
 	free(cell);