process_hkl: Generate intensity histograms

2 files changed, 114 insertions, 155 deletions

diff --git a/scripts/frequency b/scripts/frequency
index 26e29bb4..d6caf72c 100755
--- a/scripts/frequency
+++ b/scripts/frequency
@@ -2,168 +2,37 @@
 
 use strict;
 
-my $bins = 200;
 my $nplots = 32;
+my $sym = "6/mmm";
 
-open(FH, $ARGV[0]);
+my $stream = $ARGV[0];
 
-my $line;
-my %counts;
-
-# --- Step 1: Count reflections and find the most popular ones
-
-while ( $line = <FH> ) {
-
-	my $h;
-	my $k;
-	my $l;
-
-	chomp($line);
-
-	if ( $line =~ /^\s+(\d+)\s+(\d+)\s+(\d+)/ ) {
-
-		my $key;
-
-		$h = $1;
-		$k = $2;
-		$l = $3;
-
-		$key = $h.",".$k.",".$l;
-		$counts{$key}++;
-
-	}
-
-}
-
-my $ref;
-my $i;
-$i = 0;
-foreach $ref ( sort {-($counts{$a} <=> $counts{$b})} keys %counts )
-{
-	$i++;
-	if ( $i > $nplots ) {
-		delete($counts{$ref});
-	} else {
-		printf("%s %i\n", $ref, $counts{$ref});
-	}
-}
-
-
-# --- Step 2: Work out histogram bins for each reflection
-
-seek(FH, 0, 0);
-
-my %mins;
-my %maxs;
+#system("~/crystfel/src/process_hkl -i ".$stream." -o test.hkl -y ".$sym);
 
+open(FH, "test.hkl");
+my $max = 0;
+my $h;
+my $k;
+my $l;
+my $line;
 while ( $line = <FH> ) {
 
-	my $key;
-	my $intensity;
-
-	chomp($line);
-
-	if ( $line =~ /^\s+(\d+)\s+(\d+)\s+(\d+)\s+([\d\.]+)/ ) {
-
-		my $h;
-		my $k;
-		my $l;
+	chomp $line;
 
-		$h = $1;
-		$k = $2;
-		$l = $3;
-		$intensity = $4;
-
-		$key = $h.",".$k.",".$l;
-
-		next unless exists($counts{$key});
-
-		unless ( exists($mins{$key}) ) {
-			$mins{$key} = $intensity;
-		}
-		unless ( exists($maxs{$key}) ) {
-			$maxs{$key} = $intensity;
-		}
-		if ( $intensity < $mins{$key} ) {
-			$mins{$key} = $intensity;
-		}
-		if ( $intensity > $maxs{$key} ) {
-			$maxs{$key} = $intensity;
-		}
-
-	}
-
-}
-
-# --- Step 3: Bin the counts and plot
-
-open(GP, "| gnuplot");
-printf(GP "set term postscript enhanced font \"Helvetica,20\"\n");
-printf(GP "set output \"histo.ps\"\n");
-printf(GP "set xtics nomirror out rotate by -60\n");
-#printf(GP "set logscale x\n");
-foreach $ref ( sort {-($counts{$a} <=> $counts{$b})} keys %counts )
-{
-	my $max = $maxs{$ref};
-	my $min = $mins{$ref};
-	my $binsize = ($max+1-$min)/$bins;
-	my @hist;
-
-	printf("%s %f -> %f\n", $ref, $mins{$ref}, $maxs{$ref});
-
-	for ( $i=0; $i<$bins; $i++ ) {
-		$hist[$i] = 0;
-	}
-
-	seek(FH, 0, 0);
-	while ( $line = <FH> ) {
-
-		my $bin;
-		my $h;
-		my $k;
-		my $l;
-		my $intensity;
-
-		chomp($line);
-
-		if ( $line =~ /^\s+(\d+)\s+(\d+)\s+(\d+)\s+([\d\.]+)/ ) {
-
-			my $key;
+	if ( $line =~ /\s+(\d+)$/ ) {
 
+		my $n = $1;
+		if ( $n > $max ) {
+			$line =~ /^\s+(\d+)\s+(\d+)\s+(\d+)\s+/;
 			$h = $1;
 			$k = $2;
 			$l = $3;
-			$intensity = $4;
-
-			$key = $h.",".$k.",".$l;
-
-			next unless ( $ref eq $key );
-
-			$bin = int(($intensity-$min)/$binsize);
-			$hist[$bin]++;
-
+			$max = $n;
 		}
 
 	}
-
-	open(OFH, "> ".$ref.".dat");
-	for ( $i=0; $i<$bins; $i++ ) {
-		printf(OFH "%f\t%i\n", $min+(($i+0.5)*$binsize), $hist[$i]);
-	}
-	close(OFH);
-
-	my $nref = $ref;
-	$nref =~ s/\,/\ /g;
-	my $title = "t \"".$nref." reflection\"";
-	printf(GP "plot [] [0:20]\"".$ref.".dat\" u 1:2 w histeps ".$title."\n");
-
 }
 
-close(FH);
-close(GP);
+printf("%i %i %i = %i\n", $h, $k, $l, $max);
 
-foreach $ref ( sort {-($counts{$a} <=> $counts{$b})} keys %counts ) {
-	unlink($ref.".dat");
-}
-system("ps2pdf histo.ps");
-unlink("histo.ps");
+exit 0;
diff --git a/src/process_hkl.c b/src/process_hkl.c
index ec23e1cd..05c71985 100644
--- a/src/process_hkl.c
+++ b/src/process_hkl.c
@@ -29,8 +29,8 @@
 #include "symmetry.h"
 
 
-/* Number of divisions for R vs |q| graphs */
-#define RVQDV (20)
+/* Number of divisions for intensity histograms */
+#define NBINS (200)
 
 
 static void show_help(const char *s)
@@ -56,6 +56,8 @@ 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"
+"  -g, --histogram=<h,k,l>   Calculate the histogram of measurements for this\n"
+"                             reflection.\n"
 "\n"
 "      --scale               Scale each pattern for best fit with the current\n"
 "                             model.\n"
@@ -64,6 +66,47 @@ static void show_help(const char *s)
 }
 
 
+static void plot_histogram(double *vals, int n)
+{
+	int i;
+	double max = -INFINITY;
+	double min = +INFINITY;
+	double step;
+	int histo[NBINS];
+	FILE *fh;
+
+	fh = fopen("histogram.dat", "w");
+	if ( fh == NULL ) {
+		ERROR("Couldn't open 'histogram.dat'\n");
+		return;
+	}
+
+	for ( i=0; i<n; i++ ) {
+		if ( vals[i] > max ) max = vals[i];
+		if ( vals[i] < min ) min = vals[i];
+	}
+	STATUS("%f %f\n", min, max);
+
+	for ( i=0; i<NBINS; i++ ) {
+		histo[i] = 0;
+	}
+
+	step = (max-min)/NBINS;
+
+	for ( i=0; i<n; i++ ) {
+		int bin;
+		bin = (vals[i]-min)/step;
+		histo[bin]++;
+	}
+
+	for ( i=0; i<NBINS; i++ ) {
+		fprintf(fh, "%f %i\n", min+step*i, histo[i]);
+	}
+
+	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.
@@ -160,7 +203,9 @@ 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)
+                          const char *holo, const char *mero, double *hist_vals,
+                          signed int hist_h, signed int hist_k,
+                          signed int hist_l, int *hist_n)
 {
 	int i;
 	int twin;
@@ -216,6 +261,14 @@ static void merge_pattern(double *model, ReflItemList *observed,
 		/* Increase count count */
 		integrate_count(model_counts, h, k, l, 1);
 
+		if ( hist_vals != NULL ) {
+			int p = *hist_n;
+			if ( (h==hist_h) && (k==hist_k) && (l==hist_l) ) {
+				hist_vals[p] = intensity;
+				*hist_n = p+1;
+			}
+		}
+
 	}
 
 	/* Dump the reflections in this pattern into the overall list */
@@ -230,7 +283,9 @@ static void merge_all(FILE *fh, double **pmodel, ReflItemList **pobserved,
                       int config_maxonly, int config_scale, int config_sum,
                       int config_stopafter,
                       ReflItemList *twins, const char *holo, const char *mero,
-                      int n_total_patterns)
+                      int n_total_patterns, double *hist_vals,
+                      signed int hist_h, signed int hist_k, signed int hist_l,
+                      int *hist_i)
 {
 	char *rval;
 	float f0;
@@ -277,7 +332,9 @@ static void merge_all(FILE *fh, double **pmodel, ReflItemList **pobserved,
 			/* Start of second or later pattern */
 			merge_pattern(model, observed, new_pattern, items,
 			              counts, config_maxonly,
-			              twins, holo, mero);
+			              twins, holo, mero,
+			              hist_vals, hist_h, hist_k, hist_l,
+			              hist_i);
 
 			if ( n_patterns == config_stopafter ) break;
 
@@ -380,6 +437,10 @@ int main(int argc, char *argv[])
 	ReflItemList *observed;
 	int i;
 	const char *holo = NULL;
+	char *histo = NULL;
+	signed int hist_h, hist_k, hist_l;
+	double *hist_vals = NULL;
+	int hist_i;
 
 	/* Long options */
 	const struct option longopts[] = {
@@ -393,11 +454,12 @@ int main(int argc, char *argv[])
 		{"scale",              0, &config_scale,       1},
 		{"symmetry",           1, NULL,               'y'},
 		{"pdb",                1, NULL,               'p'},
+		{"histogram",          1, NULL,               'g'},
 		{0, 0, NULL, 0}
 	};
 
 	/* Short options */
-	while ((c = getopt_long(argc, argv, "hi:e:ro:p:y:",
+	while ((c = getopt_long(argc, argv, "hi:e:ro:p:y:g:",
 	                        longopts, NULL)) != -1) {
 
 		switch (c) {
@@ -425,6 +487,10 @@ int main(int argc, char *argv[])
 			sym = strdup(optarg);
 			break;
 
+		case 'g' :
+			histo = strdup(optarg);
+			break;
+
 		case 0 :
 			break;
 
@@ -475,6 +541,22 @@ int main(int argc, char *argv[])
 		holo = strdup("1");
 	}
 
+	if ( histo != NULL ) {
+		int r;
+		r = sscanf(histo, "%i,%i,%i", &hist_h, &hist_k, &hist_l);
+		if ( r != 3 ) {
+			ERROR("Invalid indices for '--histogram'\n");
+			return 1;
+		}
+		hist_vals = malloc(10*1024*sizeof(double));
+		free(histo);
+		STATUS("Histogramming %i %i %i -> ", hist_h, hist_k, hist_l);
+		/* Put into the asymmetric cell for the target group */
+		get_asymm(hist_h, hist_k, hist_l,
+		          &hist_h, &hist_k, &hist_l, sym);
+		STATUS("%i %i %i\n", hist_h, hist_k, hist_l);
+	}
+
 	/* Open the data stream */
 	if ( strcmp(filename, "-") == 0 ) {
 		fh = stdin;
@@ -492,13 +574,21 @@ int main(int argc, char *argv[])
 	STATUS("There are %i patterns to process\n", n_total_patterns);
 	rewind(fh);
 
+	hist_i = 0;
 	merge_all(fh, &model, &observed, &counts,
 	          config_maxonly, config_scale, config_sum, config_stopafter,
-                  twins, holo, sym, n_total_patterns);
+                  twins, holo, sym, n_total_patterns,
+                  hist_vals, hist_h, hist_k, hist_l, &hist_i);
 	rewind(fh);
 
 	fclose(fh);
 
+	if ( hist_vals != NULL ) {
+		STATUS("%i %i %i was seen %i times.\n", hist_h, hist_k, hist_l,
+		                                        hist_i);
+		plot_histogram(hist_vals, hist_i);
+	}
+
 	if ( output != NULL ) {
 		write_reflections(output, observed, model, NULL, counts, cell);
 	}