compare_hkl: Take symmetry into account

6 files changed, 83 insertions, 75 deletions

diff --git a/src/Makefile.am b/src/Makefile.am
index b9b70c22..a61ebb6b 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -39,7 +39,7 @@ get_hkl_SOURCES = get_hkl.c sfac.c cell.c utils.c reflections.c symmetry.c
 get_hkl_LDADD = @LIBS@
 
 compare_hkl_SOURCES = compare_hkl.c sfac.c cell.c utils.c reflections.c \
-                      statistics.c
+                      statistics.c symmetry.c
 compare_hkl_LDADD = @LIBS@
 
 powder_plot_SOURCES = powder_plot.c cell.c utils.c image.c hdf5-file.c \
diff --git a/src/Makefile.in b/src/Makefile.in
index 7c8c1028..7cd5226c 100644
--- a/src/Makefile.in
+++ b/src/Makefile.in
@@ -60,7 +60,7 @@ calibrate_detector_OBJECTS = $(am_calibrate_detector_OBJECTS)
 calibrate_detector_DEPENDENCIES =
 am_compare_hkl_OBJECTS = compare_hkl.$(OBJEXT) sfac.$(OBJEXT) \
 	cell.$(OBJEXT) utils.$(OBJEXT) reflections.$(OBJEXT) \
-	statistics.$(OBJEXT)
+	statistics.$(OBJEXT) symmetry.$(OBJEXT)
 compare_hkl_OBJECTS = $(am_compare_hkl_OBJECTS)
 compare_hkl_DEPENDENCIES =
 am_get_hkl_OBJECTS = get_hkl.$(OBJEXT) sfac.$(OBJEXT) cell.$(OBJEXT) \
@@ -256,7 +256,7 @@ indexamajig_LDADD = @LIBS@
 get_hkl_SOURCES = get_hkl.c sfac.c cell.c utils.c reflections.c symmetry.c
 get_hkl_LDADD = @LIBS@
 compare_hkl_SOURCES = compare_hkl.c sfac.c cell.c utils.c reflections.c \
-                      statistics.c
+                      statistics.c symmetry.c
 
 compare_hkl_LDADD = @LIBS@
 powder_plot_SOURCES = powder_plot.c cell.c utils.c image.c hdf5-file.c \
diff --git a/src/compare_hkl.c b/src/compare_hkl.c
index 445798d6..27ca438a 100644
--- a/src/compare_hkl.c
+++ b/src/compare_hkl.c
@@ -25,6 +25,7 @@
 #include "sfac.h"
 #include "reflections.h"
 #include "statistics.h"
+#include "symmetry.h"
 
 
 static void show_help(const char *s)
@@ -45,6 +46,7 @@ int main(int argc, char *argv[])
 	int c;
 	double *ref1;
 	double *ref2;
+	double *ref2_transformed;
 	double *out;
 	UnitCell *cell;
 	char *outfile = NULL;
@@ -64,7 +66,7 @@ int main(int argc, char *argv[])
 	};
 
 	/* Short options */
-	while ((c = getopt_long(argc, argv, "ho:", longopts, NULL)) != -1) {
+	while ((c = getopt_long(argc, argv, "ho:y:", longopts, NULL)) != -1) {
 
 		switch (c) {
 		case 'h' :
@@ -110,36 +112,44 @@ int main(int argc, char *argv[])
 		return 1;
 	}
 
-	/* Find common reflections */
-	icommon = intersection_items(i1, i2);
-	ncom = num_items(icommon);
 
 	/* List for output scale factor map */
 	out = new_list_intensity();
 
-	for ( i=0; i<ncom; i++ ) {
 
-		double i1, i2;
+	/* Find common reflections (taking symmetry into account) */
+	icommon = new_items();
+	ref2_transformed = new_list_intensity();
+	for ( i=0; i<num_items(i1); i++ ) {
+
 		struct refl_item *it;
 		signed int h, k, l;
+		signed int he, ke, le;
+		double val1, val2;
 
-		it = get_item(icommon, i);
+		it = get_item(i1, i);
 		h = it->h;  k = it->k;  l = it->l;
 
-		i1 = lookup_intensity(ref1, h, k, l);
-		i2 = lookup_intensity(ref2, h, k, l);
+		if ( !find_unique_equiv(i2, h, k, l, sym, &he, &ke, &le) ) {
+			continue;
+		}
 
-		set_intensity(out, h, k, l, i1/i2);
+		val1 = lookup_intensity(ref1, h, k, l);
+		val2 = lookup_intensity(ref2, he, ke, le);
+		set_intensity(ref2_transformed, h, k, l, val2);
+		set_intensity(out, h, k, l, val1/val2);
+		add_item(icommon, h, k, l);
 
 	}
+	ncom = num_items(icommon);
 
 	STATUS("%i,%i reflections: %i in common\n",
 	       num_items(i1), num_items(i2), ncom);
-	R2 = stat_r2(ref1, ref2, icommon, &scale);
+	R2 = stat_r2(ref1, ref2_transformed, icommon, &scale);
 	STATUS("R2 = %5.4f %% (scale=%5.2e)\n", R2*100.0, scale);
-	Rmerge = stat_rmerge(ref1, ref2, icommon, &scale);
+	Rmerge = stat_rmerge(ref1, ref2_transformed, icommon, &scale);
 	STATUS("Rmerge = %5.4f %% (scale=%5.2e)\n", Rmerge*100.0, scale);
-	pearson = stat_pearson(ref1, ref2, icommon);
+	pearson = stat_pearson(ref1, ref2_transformed, icommon);
 	STATUS("Pearson r = %5.4f\n", pearson);
 
 	if ( outfile != NULL ) {
diff --git a/src/get_hkl.c b/src/get_hkl.c
index 1052d786..14f79266 100644
--- a/src/get_hkl.c
+++ b/src/get_hkl.c
@@ -72,42 +72,6 @@ static void noisify_reflections(double *ref)
 }
 
 
-static void scold_user_about_symmetry(signed int h, signed int k, signed int l,
-                                      signed int he, signed int ke,
-                                      signed int le)
-{
-	ERROR("Merohedrally equivalent reflection (%i %i %i) found for "
-	      "%i %i %i.\n", he, ke, le, h, k, l);
-	ERROR("This indicates that you lied to me about the symmetry of the "
-	      "input reflections.  ");
-	ERROR("I won't be able to give you a meaningful result in this "
-	      "situation, so I'm going to give up right now.  ");
-	ERROR("Please reconsider your previous processing of the data, and "
-	      "perhaps try again with a lower symmetry for the '-y' option.\n");
-}
-
-
-static int find_unique_equiv(ReflItemList *items, signed int h, signed int k,
-                             signed int l, const char *mero, signed int *hu,
-                             signed int *ku, signed int *lu)
-{
-	int i;
-
-	for ( i=0; i<num_equivs(h, k, l, mero); i++ ) {
-
-		signed int he, ke, le;
-		get_equiv(h, k, l, &he, &ke, &le, mero, i);
-		if ( find_item(items, he, ke, le) ) {
-			*hu = he;  *ku = ke;  *lu = le;
-			return 1;
-		}
-
-	}
-
-	return 0;
-}
-
-
 static ReflItemList *twin_reflections(double *ref, ReflItemList *items,
                                       const char *holo, const char *mero)
 {
@@ -125,28 +89,14 @@ static ReflItemList *twin_reflections(double *ref, ReflItemList *items,
 		int skip;
 
 		it = get_item(items, i);
-		h = it->h;  k = it->k;  l = it->l;
 
-		/* None of the equivalent reflections should exist in the
-		 * input dataset.  That would indicate that the user lied about
-		 * the input symmetry.
-		 *
-		 * Start from j=1 to ignore the reflection itself.
-		 */
-		for ( j=1; j<num_equivs(h, k, l, mero); j++ ) {
-
-			signed int he, ke, le;
-			get_equiv(h, k, l, &he, &ke, &le, mero, j);
-			if ( !find_item(items, he, ke, le) ) continue;
-
-			scold_user_about_symmetry(h, k, l, he, ke, le);
-			abort();
-
-		}
-		/* It doesn't matter if the reflection wasn't actually the one
-		 * we define as being in the asymmetric unit cell, as long as
-		 * things aren't confused by there being more than one of it.
+		/* There is a many-to-one correspondence between reflections
+		 * in the merohedral and holohedral unit cells.  Do the
+		 * calculation only once for each reflection in the holohedral
+		 * cell, which contains fewer reflections.
 		 */
+		get_asymm(it->h, it->k, it->l, &h, &k, &l, holo);
+		if ( find_item(new, h, k, l) ) continue;
 
 		n = num_equivs(h, k, l, holo);
 
@@ -163,9 +113,6 @@ static ReflItemList *twin_reflections(double *ref, ReflItemList *items,
 			 * We might not have the particular (merohedral)
 			 * equivalent which belongs to our definition of the
 			 * asymmetric unit cell, so check them all.
-			 *
-			 * We checked earlier that there's only one of these
-			 * for each reflection.
 			 */
 			if ( !find_unique_equiv(items, he, ke, le, mero,
 			                        &hu, &ku, &lu) ) {
diff --git a/src/symmetry.c b/src/symmetry.c
index 4bc8d115..54b7f5cc 100644
--- a/src/symmetry.c
+++ b/src/symmetry.c
@@ -352,3 +352,50 @@ ReflItemList *get_twins(ReflItemList *items, const char *holo, const char *mero)
 
 	return ops;
 }
+
+
+static void scold_user_about_symmetry(signed int h, signed int k, signed int l,
+                                      signed int he, signed int ke,
+                                      signed int le)
+{
+	ERROR("Symmetrically equivalent reflection (%i %i %i) found for "
+	      "%i %i %i in the input.\n", he, ke, le, h, k, l);
+	ERROR("This indicates that you lied to me about the symmetry of the "
+	      "input reflections.  ");
+	ERROR("I won't be able to give you a meaningful result in this "
+	      "situation, so I'm going to give up right now.  ");
+	ERROR("Please reconsider your previous processing of the data, and "
+	      "perhaps try again with a lower symmetry for the '-y' option.\n");
+	abort();
+}
+
+
+int find_unique_equiv(ReflItemList *items, signed int h, signed int k,
+                      signed int l, const char *mero, signed int *hu,
+                      signed int *ku, signed int *lu)
+{
+	int i;
+	int found = 0;
+
+	for ( i=0; i<num_equivs(h, k, l, mero); i++ ) {
+
+		signed int he, ke, le;
+		int f;
+		get_equiv(h, k, l, &he, &ke, &le, mero, i);
+		f = find_item(items, he, ke, le);
+
+		/* There must only be one equivalent.  If there are more, it
+		 * indicates that the user lied about the input symmetry. */
+		if ( f && found ) {
+			scold_user_about_symmetry(he, ke, le, *hu, *ku, *lu);
+		}
+
+		if ( f && !found ) {
+			*hu = he;  *ku = ke;  *lu = le;
+			found = 1;
+		}
+
+	}
+
+	return found;
+}
diff --git a/src/symmetry.h b/src/symmetry.h
index 437066ad..bb4ffa60 100644
--- a/src/symmetry.h
+++ b/src/symmetry.h
@@ -40,5 +40,9 @@ extern const char *get_holohedral(const char *sym);
 extern ReflItemList *get_twins(ReflItemList *items,
                                const char *holo, const char *mero);
 
+extern int find_unique_equiv(ReflItemList *items, signed int h, signed int k,
+                             signed int l, const char *mero, signed int *hu,
+                             signed int *ku, signed int *lu);
+
 
 #endif	/* SYMMETRY_H */