Finishing tweaks for uncentering stuff

Today, I have mostly been having my life made difficult by the PDB's invention of "H centering".

5 files changed, 280 insertions, 68 deletions

diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c
index 2f966466..f854efe0 100644
--- a/libcrystfel/src/cell-utils.c
+++ b/libcrystfel/src/cell-utils.c
@@ -171,21 +171,19 @@ void cell_print(UnitCell *cell)
 	double a, b, c, alpha, beta, gamma;
 	double ax, ay, az, bx, by, bz, cx, cy, cz;
 	LatticeType lt;
+	char cen;
 
 	lt = cell_get_lattice_type(cell);
+	cen = cell_get_centering(cell);
 
-	STATUS("%s %c", str_lattice(lt), cell_get_centering(cell));
+	STATUS("%s %c", str_lattice(lt), cen);
 
-	switch ( lt )
+	if ( (lt==L_MONOCLINIC) || (lt==L_TETRAGONAL) || ( lt==L_HEXAGONAL)
+	  || ( (lt==L_ORTHORHOMBIC) && (cen=='A') )
+	  || ( (lt==L_ORTHORHOMBIC) && (cen=='B') )
+	  || ( (lt==L_ORTHORHOMBIC) && (cen=='C') ) )
 	{
-		case L_MONOCLINIC :
-		case L_TETRAGONAL :
-		case L_HEXAGONAL :
 		STATUS(", unique axis %c", cell_get_unique_axis(cell));
-		break;
-
-		default :
-		break;
 	}
 
 	if ( right_handed(cell) ) {
@@ -263,6 +261,11 @@ int bravais_lattice(UnitCell *cell)
 		return 0;
 
 		case 'H' :
+		/* "Hexagonal H" is not a Bravais lattice, but rather something
+		 * invented by the PDB to make life difficult for programmers.
+		 * Accepting it as Bravais seems to be the least painful way to
+		 * handle it correctly. Yuk. */
+		if ( ua != 'c' ) return 0;
 		if ( lattice == L_HEXAGONAL ) return 1;
 		return 0;
 
@@ -294,24 +297,28 @@ static UnitCellTransformation *uncentering_transformation(UnitCell *in,
 	t = tfn_identity();
 	if ( t == NULL ) return NULL;
 
-	if ( (ua == 'a') || (cen == 'A') ) {
+	if ( ua == 'a' ) {
 		tfn_combine(t, tfn_vector(0,0,1),
 		               tfn_vector(0,1,0),
 		               tfn_vector(-1,0,0));
-		if ( cen == 'A' ) cen = 'C';
 	}
 
-	if ( (ua == 'b') || (cen == 'B') ) {
+	if ( ua == 'b' ) {
 		tfn_combine(t, tfn_vector(1,0,0),
 		               tfn_vector(0,0,1),
 		               tfn_vector(0,-1,0));
-		if ( cen == 'B' ) cen = 'C';
 	}
 
 	switch ( cen ) {
 
 		case 'P' :
+		*new_latt = lt;
+		*new_centering = 'P';
+		break;
+
 		case 'R' :
+		*new_latt = L_RHOMBOHEDRAL;
+		*new_centering = 'R';
 		break;
 
 		case 'I' :
@@ -335,7 +342,6 @@ static UnitCellTransformation *uncentering_transformation(UnitCell *in,
 		if ( lt == L_CUBIC ) {
 			*new_latt = L_RHOMBOHEDRAL;
 			*new_centering = 'R';
-
 		} else {
 			assert(lt == L_ORTHORHOMBIC);
 			*new_latt = L_TRICLINIC;
@@ -343,6 +349,8 @@ static UnitCellTransformation *uncentering_transformation(UnitCell *in,
 		}
 		break;
 
+		case 'A' :
+		case 'B' :
 		case 'C' :
 		tfn_combine(t, tfn_vector(H,H,0),
 		               tfn_vector(-H,H,0),
@@ -352,6 +360,7 @@ static UnitCellTransformation *uncentering_transformation(UnitCell *in,
 		break;
 
 		case 'H' :
+		/* Obverse setting */
 		tfn_combine(t, tfn_vector(TT,OT,OT),
 		               tfn_vector(-OT,OT,OT),
 		               tfn_vector(-OT,-TT,OT));
@@ -366,18 +375,20 @@ static UnitCellTransformation *uncentering_transformation(UnitCell *in,
 
 	}
 
-	if ( ua == 'a' ) {
-		tfn_combine(t,  tfn_vector(0,0,-1),
-		                tfn_vector(0,1,0),
-		                tfn_vector(1,0,0));
-		if ( cen == 'C' ) cen = 'A';
-	}
+	/* Reverse the axis permutation, but only if this was not an H->R
+	 * transformation */
+	if ( !((cen=='H') && (*new_latt == L_RHOMBOHEDRAL)) ) {
+		if ( ua == 'a' ) {
+			tfn_combine(t, tfn_vector(0,0,-1),
+				       tfn_vector(0,1,0),
+				       tfn_vector(1,0,0));
+		}
 
-	if ( ua == 'b' ) {
-		tfn_combine(t,  tfn_vector(1,0,0),
-		                tfn_vector(0,0,-1),
-		                tfn_vector(0,1,0));
-		if ( cen == 'C' ) cen = 'B';
+		if ( ua == 'b' ) {
+			tfn_combine(t, tfn_vector(1,0,0),
+				       tfn_vector(0,0,-1),
+				       tfn_vector(0,1,0));
+		}
 	}
 
 	return t;
@@ -406,6 +417,7 @@ UnitCell *uncenter_cell(UnitCell *in, UnitCellTransformation **t)
 
 	if ( !bravais_lattice(in) ) {
 		ERROR("Cannot uncenter: not a Bravais lattice.\n");
+		cell_print(in);
 		return NULL;
 	}
 
@@ -481,11 +493,11 @@ UnitCell *match_cell(UnitCell *cell_in, UnitCell *template_in, int verbose,
 	float angtol = deg2rad(tols[3]);
 	UnitCell *cell;
 	UnitCell *template;
-	UnitCellTransformation *to_given_cell;
+	UnitCellTransformation *uncentering;
 	UnitCell *new_cell_trans;
 
 	/* "Un-center" the template unit cell to make the comparison easier */
-	template = uncenter_cell(template_in, &to_given_cell);
+	template = uncenter_cell(template_in, &uncentering);
 
 	/* The candidate cell is also uncentered, because it might be centered
 	 * if it came from (e.g.) MOSFLM */
@@ -671,7 +683,7 @@ UnitCell *match_cell(UnitCell *cell_in, UnitCell *template_in, int verbose,
 	free(cand[2]);
 
 	/* Reverse the de-centering transformation */
-	new_cell_trans = cell_transform_inverse(new_cell, to_given_cell);
+	new_cell_trans = cell_transform_inverse(new_cell, uncentering);
 	cell_free(new_cell);
 	cell_set_lattice_type(new_cell, cell_get_lattice_type(template_in));
 	cell_set_centering(new_cell, cell_get_centering(template_in));
@@ -1115,10 +1127,13 @@ int cell_is_sensible(UnitCell *cell)
  * lattice is a conventional Bravais lattice.
  * Warnings are printied if any of the checks are failed.
  *
+ * Returns: true if cell is invalid.
+ *
  */
-void validate_cell(UnitCell *cell)
+int validate_cell(UnitCell *cell)
 {
 	int err = 0;
+	char cen, ua;
 
 	if ( !cell_is_sensible(cell) ) {
 		ERROR("Warning: Unit cell parameters are not sensible.\n");
@@ -1136,5 +1151,58 @@ void validate_cell(UnitCell *cell)
 		err = 1;
 	}
 
-	if ( err ) cell_print(cell);
+	cen = cell_get_centering(cell);
+	ua = cell_get_unique_axis(cell);
+	if ( (cen == 'A') && (ua != 'a') ) {
+		ERROR("Warning: centering doesn't match unique axis.\n");
+		err = 1;
+	}
+	if ( (cen == 'B') && (ua != 'b') ) {
+		ERROR("Warning: centering doesn't match unique axis.\n");
+		err = 1;
+	}
+	if ( (cen == 'C') && (ua != 'c') ) {
+		ERROR("Warning: centering doesn't match unique axis.\n");
+		err = 1;
+	}
+
+	return err;
+}
+
+
+/**
+ * forbidden_reflection:
+ * @cell: A %UnitCell
+ * @h: h index to check
+ * @k: k index to check
+ * @l: l index to check
+ *
+ * Returns: true if this reflection is forbidden.
+ *
+ */
+int forbidden_reflection(UnitCell *cell,
+                         signed int h, signed int k, signed int l)
+{
+	char cen;
+
+	cen = cell_get_centering(cell);
+
+	/* Reflection conditions here must match the transformation matrices
+	 * in uncentering_transformation().  tests/centering_check verifies
+	 * this (amongst other things). */
+
+	if ( cen == 'P' ) return 0;
+	if ( cen == 'R' ) return 0;
+
+	if ( cen == 'A' ) return (k+l) % 2;
+	if ( cen == 'B' ) return (h+l) % 2;
+	if ( cen == 'C' ) return (h+k) % 2;
+
+	if ( cen == 'I' ) return (h+k+l) % 2;
+	if ( cen == 'F' ) return ((h+k) % 2) || ((h+l) % 2) || ((k+l) % 2);
+
+	/* Obverse setting */
+	if ( cen == 'H' ) return (-h+k+l) % 3;
+
+	return 0;
 }
diff --git a/libcrystfel/src/cell-utils.h b/libcrystfel/src/cell-utils.h
index 8acf2f85..f92ab22d 100644
--- a/libcrystfel/src/cell-utils.h
+++ b/libcrystfel/src/cell-utils.h
@@ -53,7 +53,7 @@ extern UnitCell *load_cell_from_pdb(const char *filename);
 
 extern int cell_is_sensible(UnitCell *cell);
 
-extern void validate_cell(UnitCell *cell);
+extern int validate_cell(UnitCell *cell);
 
 extern UnitCell *uncenter_cell(UnitCell *in, UnitCellTransformation **tr);
 
@@ -63,4 +63,7 @@ extern int right_handed(UnitCell *cell);
 
 extern const char *str_lattice(LatticeType l);
 
+extern int forbidden_reflection(UnitCell *cell,
+                                signed int h, signed int k, signed int l);
+
 #endif	/* CELL_UTILS_H */
diff --git a/libcrystfel/src/cell.c b/libcrystfel/src/cell.c
index 7a62f51c..e7c9dace 100644
--- a/libcrystfel/src/cell.c
+++ b/libcrystfel/src/cell.c
@@ -260,6 +260,9 @@ UnitCell *cell_new_from_cell(UnitCell *orig)
 	cell_get_cartesian(orig, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
 	cell_set_cartesian(new, ax, ay, az, bx, by, bz, cx, cy, cz);
 	cell_set_pointgroup(new, orig->pointgroup);
+	cell_set_lattice_type(new, orig->lattice_type);
+	cell_set_centering(new, orig->centering);
+	cell_set_unique_axis(new, orig->unique_axis);
 
 	return new;
 }
@@ -691,12 +694,12 @@ UnitCell *cell_transform(UnitCell *cell, UnitCellTransformation *t)
 
 	if ( t == NULL ) return NULL;
 
-	out = cell_new();
+	out = cell_new_from_cell(cell);
 	if ( out == NULL ) return NULL;
 
-	cell_get_cartesian(cell, &ax, &ay, &az,
-	                         &bx, &by, &bz,
-	                         &cx, &cy, &cz);
+	cell_get_cartesian(out, &ax, &ay, &az,
+	                        &bx, &by, &bz,
+	                        &cx, &cy, &cz);
 
 	m = gsl_matrix_alloc(3,3);
 	a = gsl_matrix_calloc(3,3);
diff --git a/tests/centering_check.c b/tests/centering_check.c
index eca67c85..e17915d3 100644
--- a/tests/centering_check.c
+++ b/tests/centering_check.c
@@ -30,6 +30,7 @@
 #include <stdlib.h>
 #include <stdio.h>
 #include <stdarg.h>
+#include <fenv.h>
 
 #include <cell.h>
 #include <cell-utils.h>
@@ -39,24 +40,13 @@ static int check_cell(UnitCell *cell, const char *text)
 {
 	int err = 0;
 
-	if ( !cell_is_sensible(cell) ) {
-		ERROR(" %s unit cell parameters are not sensible.\n", text);
-		err = 1;
-	}
-
-	if ( !bravais_lattice(cell) ) {
-		ERROR(" %s unit cell is not a conventional Bravais"
-		      " lattice.\n", text);
-		err = 1;
-	}
+	err = validate_cell(cell);
 
-	if ( !right_handed(cell) ) {
-		ERROR(" %s unit cell is not right handed.\n", text);
-		err = 1;
+	if ( err ) {
+		ERROR("%s cell:\n", text);
+		cell_print(cell);
 	}
 
-	if ( err ) cell_print(cell);
-
 	return err;
 }
 
@@ -65,40 +55,147 @@ static int check_centering(double a, double b, double c,
                            double al, double be, double ga,
                            LatticeType latt, char cen, char ua)
 {
-	UnitCell *cell;
+	UnitCell *cell, *cref;
 	UnitCell *n;
 	UnitCellTransformation *t;
 	int fail = 0;
+	int i;
+	double asx, asy, asz;
+	double bsx, bsy, bsz;
+	double csx, csy, csz;
+	double ax, ay, az;
+	double bx, by, bz;
+	double cx, cy, cz;
 
 	STATUS("   --------------->  "
 	       "Checking %s %c (ua %c) %5.2e %5.2e %5.2e %5.2f %5.2f %5.2f\n",
 	       str_lattice(latt), cen, ua, a, b, c, al, be, ga);
 
-	cell = cell_new_from_parameters(a, b, c,
+	cref = cell_new_from_parameters(a, b, c,
 	                                deg2rad(al), deg2rad(be), deg2rad(ga));
-	cell_set_lattice_type(cell, latt);
-	cell_set_centering(cell, cen);
-	cell_set_unique_axis(cell, ua);
+	cell_set_lattice_type(cref, latt);
+	cell_set_centering(cref, cen);
+	cell_set_unique_axis(cref, ua);
+
+	cell = cell_rotate(cref, random_quaternion());
+	if ( cell == NULL ) return 1;
+	cell_free(cref);
 
-	if ( check_cell(cell, "Input") ) fail = 1;
-	//cell_print(cell);
+	check_cell(cell, "Input");
 	n = uncenter_cell(cell, &t);
 	if ( n != NULL ) {
+		STATUS("Transformation was:\n");
+		tfn_print(t);
 		if ( check_cell(n, "Output") ) fail = 1;
+		if ( !fail ) cell_print(n);
 	} else {
 		fail = 1;
 	}
 
-	STATUS("Transformation was:\n");
-	tfn_print(t);
+	cell_get_reciprocal(cell, &asx, &asy, &asz,
+	                          &bsx, &bsy, &bsz,
+	                          &csx, &csy, &csz);
+	cell_get_cartesian(n, &ax, &ay, &az,
+	                      &bx, &by, &bz,
+	                      &cx, &cy, &cz);
+
+	fesetround(1);  /* Round towards nearest */
+	for ( i=0; i<100; i++ ) {
+
+		signed int h, k, l;
+		double x, y, z;
+		double nh, nk, nl;
+		double dh, dk, dl;
+		int f = 0;
+
+		do {
+
+			h = flat_noise(0, 30);
+			k = flat_noise(0, 30);
+			l = flat_noise(0, 30);
+
+		} while ( forbidden_reflection(cell, h, k, l) );
 
-	if ( fail ) ERROR("\n");
+		x = h*asx + k*bsx + l*csx;
+		y = h*asy + k*bsy + l*csy;
+		z = h*asz + k*bsz + l*csz;
+
+		nh = x*ax + y*ay + z*az;
+		nk = x*bx + y*by + z*bz;
+		nl = x*cx + y*cy + z*cz;
+
+		dh = nh - lrint(nh);
+		dk = nk - lrint(nk);
+		dl = nl - lrint(nl);
+		if ( fabs(dh) > 0.1 ) f++;
+		if ( fabs(dk) > 0.1 ) f++;
+		if ( fabs(dl) > 0.1 ) f++;
+
+		if ( f ) {
+			STATUS("Centered %3i %3i %3i -> "
+			       "Primitive %7.2f %7.2f %7.2f\n",
+			       h, k, l, nh, nk, nl);
+			fail = 1;
+		}
+
+	}
+
+	cell_get_reciprocal(n, &asx, &asy, &asz,
+	                       &bsx, &bsy, &bsz,
+	                       &csx, &csy, &csz);
+	cell_get_cartesian(cell, &ax, &ay, &az,
+	                         &bx, &by, &bz,
+	                         &cx, &cy, &cz);
+
+	for ( i=0; i<100; i++ ) {
+
+		signed int h, k, l;
+		double x, y, z;
+		double nh, nk, nl;
+		double dh, dk, dl;
+		int f = 0;
+		long int ih, ik, il;
+
+		h = flat_noise(0, 30);
+		k = flat_noise(0, 30);
+		l = flat_noise(0, 30);
+
+		x = h*asx + k*bsx + l*csx;
+		y = h*asy + k*bsy + l*csy;
+		z = h*asz + k*bsz + l*csz;
+
+		nh = x*ax + y*ay + z*az;
+		nk = x*bx + y*by + z*bz;
+		nl = x*cx + y*cy + z*cz;
+
+		dh = nh - lrint(nh);  dk = nk - lrint(nk);  dl = nl - lrint(nl);
+
+		if ( fabs(dh) > 0.1 ) f++;
+		if ( fabs(dk) > 0.1 ) f++;
+		if ( fabs(dl) > 0.1 ) f++;
+
+		ih = lrint(nh);  ik = lrint(nk);  il = lrint(nl);
+		if ( forbidden_reflection(cell, ih, ik, il) ) {
+			STATUS("Primitive %3i %3i %3i -> "
+			       "Centered %3li %3li %3li, "
+			       "which is forbidden\n",
+			       h, k, l, ih, ik, il);
+			fail = 1;
+		}
+
+		if ( f ) {
+			STATUS("Primitive %3i %3i %3i -> "
+			       "Centered %7.2f %7.2f %7.2f\n",
+			       h, k, l, nh, nk, nl);
+			fail = 1;
+		}
+
+	}
 
 	return fail;
 }
 
 
-
 int main(int argc, char *argv[])
 {
 	int fail = 0;
@@ -133,15 +230,15 @@ int main(int argc, char *argv[])
 
 	/* Orthorhombic A */
 	fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 90.0, 90.0,
-	                        L_ORTHORHOMBIC, 'A', '*');
+	                        L_ORTHORHOMBIC, 'A', 'a');
 
 	/* Orthorhombic B */
 	fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 90.0, 90.0,
-	                        L_ORTHORHOMBIC, 'B', '*');
+	                        L_ORTHORHOMBIC, 'B', 'b');
 
 	/* Orthorhombic C */
 	fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 90.0, 90.0,
-	                        L_ORTHORHOMBIC, 'C', '*');
+	                        L_ORTHORHOMBIC, 'C', 'c');
 
 	/* Orthorhombic I */
 	fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 90.0, 90.0,
@@ -180,10 +277,6 @@ int main(int argc, char *argv[])
 	                        L_HEXAGONAL, 'P', 'c');
 
 	/* Hexagonal H (PDB-speak for rhombohedral) */
-	fail += check_centering(40e-10, 20e-10, 20e-10, 120.0, 90.0, 90.0,
-	                        L_HEXAGONAL, 'H', 'a');
-	fail += check_centering(20e-10, 40e-10, 20e-10, 90.0, 120.0, 90.0,
-	                        L_HEXAGONAL, 'H', 'b');
 	fail += check_centering(20e-10, 20e-10, 40e-10, 90.0, 90.0, 120.0,
 	                        L_HEXAGONAL, 'H', 'c');
 
diff --git a/tests/transformation_check.c b/tests/transformation_check.c
index 80324795..7d25aa04 100644
--- a/tests/transformation_check.c
+++ b/tests/transformation_check.c
@@ -73,6 +73,39 @@ static int check_transformation(UnitCell *cell, UnitCellTransformation *tfn)
 }
 
 
+static int check_identity(UnitCell *cell, UnitCellTransformation *tfn)
+{
+	UnitCell *cnew;
+	double a[9], b[9];
+	int i;
+	int fail = 0;
+
+	cnew = cell_transform(cell, tfn);
+
+	cell_get_cartesian(cell, &a[0], &a[1], &a[2],
+	                         &a[3], &a[4], &a[5],
+	                         &a[6], &a[7], &a[8]);
+	cell_get_cartesian(cnew, &b[0], &b[1], &b[2],
+	                         &b[3], &b[4], &b[5],
+	                         &b[6], &b[7], &b[8]);
+	for ( i=0; i<9; i++ ) {
+		if ( !within_tolerance(a[i], b[i], 0.1) ) {
+			fail = 1;
+			STATUS("%e %e\n", a[i], b[i]);
+		}
+	}
+
+	if ( fail ) {
+		ERROR("Original cell not recovered after transformation:\n");
+		cell_print(cell);
+		tfn_print(tfn);
+		cell_print(cnew);
+	}
+
+	return fail;
+}
+
+
 int main(int argc, char *argv[])
 {
 	int fail = 0;
@@ -125,6 +158,18 @@ int main(int argc, char *argv[])
 	fail += check_transformation(cell, tfn);
 	tfn_free(tfn);
 
+	/* Identity in two parts */
+	tfn = tfn_identity();
+	if ( tfn == NULL ) return 1;
+	tfn_combine(tfn, tfn_vector(0,0,1),
+	                 tfn_vector(0,1,0),
+	                 tfn_vector(-1,0,0));
+	tfn_combine(tfn,  tfn_vector(0,0,-1),
+	                  tfn_vector(0,1,0),
+	                  tfn_vector(1,0,0));
+	fail += check_identity(cell, tfn);
+	tfn_free(tfn);
+
 	cell_free(cell);
 
 	return fail;