Finish symmetry tests and update table one last time

5 files changed, 222 insertions, 20 deletions

diff --git a/doc/twin-calculator.odt b/doc/twin-calculator.odt
index 56282572..3feba7fe 100644
--- a/doc/twin-calculator.odt
+++ b/doc/twin-calculator.odt
diff --git a/doc/twin-calculator.pdf b/doc/twin-calculator.pdf
index d1bfc9fd..6d3dd582 100644
--- a/doc/twin-calculator.pdf
+++ b/doc/twin-calculator.pdf
diff --git a/src/symmetry.c b/src/symmetry.c
index e9c7d642..10c41ba7 100644
--- a/src/symmetry.c
+++ b/src/symmetry.c
@@ -456,7 +456,7 @@ static SymOpList *make_4mm()
 	SymOpList *new = new_symoplist();
 	add_symop(new, v(0,-1,0,0), v(1,0,0,0), v(0,0,0,1), 4); /* 4 // l */
 	add_symop(new, v(-1,0,0,0), v(0,1,0,0), v(0,0,0,1), 2); /* m -| l */
-	new->name = strdup("4mmm");
+	new->name = strdup("4mm");
 	return expand_ops(new);
 }
 
@@ -843,7 +843,7 @@ SymOpList *get_pointgroup(const char *sym)
 	if ( strcmp(sym, "23") == 0 ) return make_23();
 	if ( strcmp(sym, "m-3") == 0 ) return make_m3bar();
 	if ( strcmp(sym, "432") == 0 ) return make_432();
-	if ( strcmp(sym, "-432") == 0 ) return make_4bar3m();
+	if ( strcmp(sym, "-43m") == 0 ) return make_4bar3m();
 	if ( strcmp(sym, "m-3m") == 0 ) return make_m3barm();
 
 	ERROR("Unknown point group '%s'\n", sym);
@@ -1057,6 +1057,13 @@ static int ops_equal(const struct sym_op *op,
 }
 
 
+static int struct_ops_equal(const struct sym_op *op1, const struct sym_op *op2)
+{
+	return ops_equal(op1, op2->h, op2->k, op2->l);
+}
+
+
+
 /* Return true if a*b = ans */
 static int check_mult(const struct sym_op *ans,
                       const struct sym_op *a, const struct sym_op *b)
@@ -1079,6 +1086,39 @@ static int check_mult(const struct sym_op *ans,
 }
 
 
+/* Check that every operation in "target" is also in "source" */
+int is_subgroup(const SymOpList *source, const SymOpList *target)
+{
+	int n_src, n_tgt;
+	int i;
+
+	n_src = num_ops(source);
+	n_tgt = num_ops(target);
+
+	for ( i=0; i<n_tgt; i++ ) {
+
+		int j;
+		int found = 0;
+
+		for ( j=0; j<n_src; j++ ) {
+
+			if ( struct_ops_equal(&target->ops[i],
+			                      &source->ops[j] ) )
+			{
+				found = 1;
+				break;
+			}
+
+		}
+
+		if ( !found ) return 0;
+
+	}
+
+	return 1;
+}
+
+
 /**
  * get_ambiguities:
  * @source: The "source" symmetry, a %SymOpList
@@ -1100,10 +1140,23 @@ SymOpList *get_ambiguities(const SymOpList *source, const SymOpList *target)
 	SymOpList *src_reordered;
 	SymOpMask *used;
 	char *name;
+	int index;
 
 	n_src = num_ops(source);
 	n_tgt = num_ops(target);
 
+	if ( !is_subgroup(source, target) ) {
+		ERROR("'%s' is not a subgroup of '%s'\n",
+		      symmetry_name(target), symmetry_name(source));
+		return NULL;
+	}
+
+	if ( n_src % n_tgt != 0 ) {
+		ERROR("Subgroup index would be fractional.\n");
+		return NULL;
+	}
+	index = n_src / n_tgt;
+
 	src_reordered = new_symoplist();
 	used = new_symopmask(source);
 
@@ -1205,6 +1258,15 @@ SymOpList *get_ambiguities(const SymOpList *source, const SymOpList *target)
 	twins = new_symoplist();
 	for ( i=0; i<n_src; i++ ) {
 		if ( used->mask[i] == 0 ) continue;
+		if ( determinant(&src_reordered->ops[i]) < 0 ) {
+			/* A mirror or inversion turned up in the list.
+			 * That means that no pure rotational ambiguity can
+			 * account for this subgroup relationship. */
+			free_symoplist(twins);
+			free_symopmask(used);
+			free_symoplist(src_reordered);
+			return NULL;
+		}
 		add_copied_op(twins, &src_reordered->ops[i]);
 	}
 
@@ -1309,13 +1371,13 @@ void describe_symmetry(const SymOpList *s)
 
 	n = num_equivs(s, NULL);
 
-	STATUS("%10s :", symmetry_name(s));
+	STATUS("%15s :", symmetry_name(s));
 
 	for ( i=0; i<n; i++ ) {
 		char *name = name_equiv(&s->ops[i]);
 		STATUS(" %6s", name);
 		free(name);
-		if ( (i!=0) && (i%8==0) ) STATUS("\n%10s  ", "");
+		if ( (i!=0) && (i%8==0) ) STATUS("\n%15s  ", "");
 	}
 	STATUS("\n");
 }
diff --git a/src/symmetry.h b/src/symmetry.h
index 7ad93ece..c91864bb 100644
--- a/src/symmetry.h
+++ b/src/symmetry.h
@@ -47,7 +47,9 @@ extern void get_equiv(const SymOpList *ops, const SymOpMask *m, int idx,
                       signed int h, signed int k, signed int l,
                       signed int *he, signed int *ke, signed int *le);
 
-extern SymOpList *get_ambiguities(const SymOpList *source, const SymOpList *target);
+extern SymOpList *get_ambiguities(const SymOpList *source,
+                                  const SymOpList *target);
+extern int is_subgroup(const SymOpList *source, const SymOpList *target);
 
 extern int is_centrosymmetric(const SymOpList *s);
 extern const char *symmetry_name(const SymOpList *ops);
diff --git a/tests/symmetry_check.c b/tests/symmetry_check.c
index a83b0d53..2e47bf83 100644
--- a/tests/symmetry_check.c
+++ b/tests/symmetry_check.c
@@ -17,11 +17,72 @@
 
 #include <stdlib.h>
 #include <stdio.h>
+#include <stdarg.h>
 
 #include "../src/symmetry.h"
 #include "../src/utils.h"
 
 
+static void find_all_ambiguities(const char *first, ...)
+{
+	va_list vp;
+	int i;
+	const char *arg;
+	SymOpList *test[32];
+	int n = 0;
+
+	test[n++] = get_pointgroup(first);
+
+	va_start(vp, first);
+
+
+	do {
+
+		arg = va_arg(vp, const char *);
+		if ( arg != NULL ) {
+			test[n++] = get_pointgroup(arg);
+		}
+
+	} while ( arg != NULL );
+
+	for ( i=0; i<n; i++ ) {
+
+		SymOpList *holo;
+		int j;
+
+		holo = test[i];
+		STATUS("%7s :", symmetry_name(holo));
+		for ( j=0; j<n; j++ ) {
+
+			SymOpList *twins;
+
+			if ( i == j ) continue;
+
+			if ( !is_subgroup(holo, test[j]) ) continue;
+
+			twins = get_ambiguities(holo, test[j]);
+			if ( twins == NULL ) continue;
+
+			if ( num_equivs(twins, NULL) == 1 ) {
+				free_symoplist(twins);
+				continue;
+			}
+
+			if ( num_equivs(twins, NULL) == 2 ) {
+				STATUS(" %5s ", symmetry_name(test[j]));
+			} else {
+				STATUS(" %5s*", symmetry_name(test[j]));
+			}
+
+		}
+		STATUS("\n");
+
+	}
+
+	STATUS("\n");
+}
+
+
 static const char *maybenot(int v)
 {
 	if ( v ) {
@@ -66,11 +127,12 @@ static void check_pg_props(const char *pg, int answer, int centro, int *fail)
 
 
 static void check_subgroup(const char *ssource, const char *starget,
-                           int *fail)
+                           int should_be_subgroup, int should_be_ambiguity,
+                           int n_exp, int *fail)
 {
 	SymOpList *source;
 	SymOpList *target;
-	SymOpList *twins;
+	int sub;
 
 	source = get_pointgroup(ssource);
 	target = get_pointgroup(starget);
@@ -79,8 +141,56 @@ static void check_subgroup(const char *ssource, const char *starget,
 		return;
 	}
 
-	twins = get_ambiguities(source, target);
-	describe_symmetry(twins);
+	sub = is_subgroup(source, target);
+	if ( sub != should_be_subgroup ) {
+		ERROR("'%s' should%s be a subgroup of '%s', but is%s.\n",
+		      starget, maybenot(should_be_subgroup),
+		      ssource, maybenot(sub));
+		*fail = 1;
+		return;
+	}
+
+	if ( should_be_subgroup ) {
+
+		SymOpList *twins;
+		int nf;
+		int amb = 1;
+
+		twins = get_ambiguities(source, target);
+		if ( twins == NULL ) amb = 0;
+		if ( amb != should_be_ambiguity ) {
+			ERROR("'%s' should%s be a rotational subgroup of '%s'"
+			      " but is%s.\n",
+			      starget, maybenot(should_be_ambiguity),
+			      ssource, maybenot(amb));
+			*fail = 1;
+			return;
+		}
+
+		if ( amb ) {
+
+			describe_symmetry(twins);
+			nf = num_equivs(twins, NULL);
+			if ( nf != n_exp ) {
+				ERROR("Expected %i operations, found %i\n",
+				      n_exp, nf);
+				*fail = 1;
+				return;
+			}
+
+		} else {
+
+			STATUS("%15s : subgroup of %s, but no ambiguity without"
+			       " inversion or mirroring\n", starget, ssource);
+
+
+		}
+
+	} else {
+
+		STATUS("%15s : not a subgroup of %s\n", starget, ssource);
+
+	}
 
 	free_symoplist(target);
 	free_symoplist(source);
@@ -145,20 +255,48 @@ int main(int argc, char *argv[])
 	check_pg_props(   "23",  12,  0, &fail);
 	check_pg_props(  "m-3",  24,  1, &fail);
 	check_pg_props(  "432",  24,  0, &fail);
-	check_pg_props( "-432",  24,  0, &fail);
+	check_pg_props( "-43m",  24,  0, &fail);
 	check_pg_props( "m-3m",  48,  1, &fail);
 	STATUS("\n");
 
-	check_subgroup("2/m", "m", &fail);
-	check_subgroup("mmm", "mm2", &fail);
-	check_subgroup("-4m2", "-4", &fail);
-	check_subgroup("-42m", "-4", &fail);
-	check_subgroup("-3m1_H", "-3_H", &fail);
-	check_subgroup("-31m_H", "-3_H", &fail);
-	check_subgroup("6/mmm", "-3_H", &fail);
-	check_subgroup("-432", "m-3m", &fail);
-	check_subgroup("m-3", "m-3m", &fail);
-	check_subgroup("23", "432", &fail);
+	/* Check some easy subgroups */
+	check_subgroup("2/m",    "m",    1, 1,  2, &fail);
+	check_subgroup("mmm",    "mm2",  1, 1,  2, &fail);
+	check_subgroup("-4m2",   "-4",   1, 1,  2, &fail);
+	check_subgroup("-42m",   "-4",   1, 1,  2, &fail);
+	check_subgroup("-3m1_H", "-3_H", 1, 1,  2, &fail);
+	check_subgroup("-31m_H", "-3_H", 1, 1,  2, &fail);
+	check_subgroup("m-3m",   "-43m", 1, 1,  2, &fail);
+	check_subgroup("m-3m",   "m-3",  1, 1,  2, &fail);
+	check_subgroup("432",    "23",   1, 1,  2, &fail);
+	check_subgroup("6/m",    "-3_H", 1, 1,  2, &fail);
+	check_subgroup("4/m",    "-4",   1, 1,  2, &fail);
+
+	/* Tetartohedral */
+	check_subgroup("6/mmm",  "-3_H", 1, 1,  4, &fail);
+
+	/* Check some things that are valid subgroups, but no valid ambiguities
+	 * exist because inversions and mirrors are not allowed */
+	check_subgroup("-1",     "1",    1, 0, -1, &fail);
+	check_subgroup("4/mmm",  "4",    1, 0, -1, &fail);
+	check_subgroup("m-3m",   "23",   1, 0, -1, &fail);
+
+	/* Check some invalid combinations */
+	check_subgroup("432",    "-43m", 0, 0, -1, &fail);
+	check_subgroup("432",    "m-3",  0, 0, -1, &fail);
+
+	/* Derive all merohedral ambiguities */
+	STATUS("\nMerohedral ambiguities:\n\n");
+	find_all_ambiguities("1", "-1", NULL);
+	find_all_ambiguities("2", "m", "2/m", NULL);
+	find_all_ambiguities("mm2", "mmm", "222", NULL);
+	find_all_ambiguities("4", "-4", "-42m", "-4m2", "4mm",
+	                     "4/m", "422", "4/mmm", NULL);
+	find_all_ambiguities("3_R", "32_R", "-3_R", "3m_R", "-3m_R", NULL);
+	find_all_ambiguities("6", "3_H", "312_H", "321_H", "622", "-3_H",
+	                     "3m1_H", "-6", "31m_H", "-3m1_H", "-6m2",
+	                     "-62m", "-31m_H", "6/mmm", "6/m", "6mm", NULL);
+	find_all_ambiguities("23", "432", "-43m", "m-3", "m-3m", NULL);
 
 	return fail;
 }