From 3e850c6c2ed2a3dbc0fe0be27d36c0d1e6b57614 Mon Sep 17 00:00:00 2001
From: Thomas White <taw@physics.org>
Date: Thu, 28 Feb 2019 10:53:43 +0100
Subject: New way of doing cell comparisons, similar to match_cell()

---
 libcrystfel/src/cell-utils.c | 230 +++++++++++++++++++++++++++++++++++++++++++
 libcrystfel/src/cell-utils.h |   4 +
 libcrystfel/src/rational.c   |  33 +++++++
 libcrystfel/src/rational.h   |   4 +
 src/cell_tool.c              |  80 ++-------------
 5 files changed, 278 insertions(+), 73 deletions(-)

diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c
index a21ae570..2bfaa3e1 100644
--- a/libcrystfel/src/cell-utils.c
+++ b/libcrystfel/src/cell-utils.c
@@ -1793,3 +1793,233 @@ int compare_reindexed_cell_parameters_and_orientation(UnitCell *a, UnitCell *b,
 	intmat_free(m);
 	return 0;
 }
+
+
+struct cand
+{
+	Rational abc[3];
+	double fom;
+};
+
+
+static int cmpcand(const void *av, const void *bv)
+{
+	const struct cand *a = av;
+	const struct cand *b = bv;
+	return a->fom > b->fom;
+}
+
+
+static Rational *find_candidates(double len, double *a, double *b, double *c,
+                                 double ltl, int *pncand)
+{
+	Rational *r;
+	struct cand *cands;
+	const int max_cand = 1024;
+	int ncand = 0;
+	Rational *rat;
+	int nrat;
+	int nrej = 0;
+	int ia, ib, ic;
+	int i;
+
+	cands = malloc(max_cand * sizeof(struct cand));
+	if ( cands == NULL ) return NULL;
+
+	rat = rtnl_list(0, 4, -4, 4, &nrat);
+	if ( rat == NULL ) return NULL;
+
+	for ( ia=0; ia<nrat; ia++ ) {
+	for ( ib=0; ib<nrat; ib++ ) {
+	for ( ic=0; ic<nrat; ic++ ) {
+		double vec[3];
+		double abc[3];
+		double veclen;
+		abc[0] = rtnl_as_double(rat[ia]);
+		abc[1] = rtnl_as_double(rat[ib]);
+		abc[2] = rtnl_as_double(rat[ic]);
+		vec[0] = a[0]*abc[0] + b[0]*abc[1] + c[0]*abc[2];
+		vec[1] = a[1]*abc[0] + b[1]*abc[1] + c[1]*abc[2];
+		vec[2] = a[2]*abc[0] + b[2]*abc[1] + c[2]*abc[2];
+		veclen = modulus(vec[0], vec[1], vec[2]);
+		if ( within_tolerance(len, veclen, ltl*100.0) ) {
+			if ( ncand == max_cand ) {
+				nrej++;
+			} else {
+				cands[ncand].abc[0] = rat[ia];
+				cands[ncand].abc[1] = rat[ib];
+				cands[ncand].abc[2] = rat[ic];
+				cands[ncand].fom = fabs(veclen - len);
+				ncand++;
+			}
+		}
+	}
+	}
+	}
+
+	if ( nrej ) {
+		ERROR("WARNING: Too many vector candidates (%i rejected)\n", nrej);
+	}
+
+	/* Sort by difference from reference vector length */
+	qsort(cands, ncand, sizeof(struct cand), cmpcand);
+
+	r = malloc(ncand * 3 * sizeof(Rational));
+	if ( r == 0 ) return NULL;
+
+	for ( i=0; i<ncand; i++ ) {
+		r[3*i+0] = cands[i].abc[0];
+		r[3*i+1] = cands[i].abc[1];
+		r[3*i+2] = cands[i].abc[2];
+	}
+	free(cands);
+
+	*pncand = ncand;
+	return r;
+}
+
+
+/**
+ * compare_reindexed_cell_parameters:
+ * @cell_in: A %UnitCell
+ * @reference_in: Another %UnitCell
+ * @ltl: Maximum allowable fractional difference in direct-space axis lengths
+ * @atl: Maximum allowable difference in direct-space angles (in radians)
+ * @pmb: Place to store pointer to matrix
+ *
+ * Compare the @cell_in with @reference_in.  If @cell is a derivative lattice
+ * of @reference, within fractional axis length difference @ltl and absolute angle
+ * difference @atl (in radians), this function returns non-zero and stores the
+ * transformation which needs to be applied to @cell_in at @pmb.
+ *
+ * Only the cell parameters will be compared.  The relative orientations are
+ * irrelevant.
+ *
+ * Returns: non-zero if the cells match, zero for no match or error.
+ *
+ */
+int compare_reindexed_cell_parameters(UnitCell *cell_in, UnitCell *reference_in,
+                                      double ltl, double atl,
+                                      RationalMatrix **pmb)
+{
+	UnitCell *cell;
+	UnitCell *reference;
+	IntegerMatrix *centering_reference;
+	IntegerMatrix *centering_cell;
+	RationalMatrix *m;
+	double a, b, c, al, be, ga;
+	double av[3], bv[3], cv[3];
+	Rational *cand_a;
+	Rational *cand_b;
+	Rational *cand_c;
+	int ncand_a, ncand_b, ncand_c;
+	int i;
+	int ia, ib;
+
+	/* Actually compare against primitive version of reference */
+	reference = uncenter_cell(reference_in, &centering_reference);
+	if ( reference == NULL ) return 0;
+
+	/* Actually compare primitive version of cell */
+	cell = uncenter_cell(cell_in, &centering_cell);
+	if ( cell == NULL ) return 0;
+
+	/* Get target parameters */
+	cell_get_parameters(reference, &a, &b, &c, &al, &be, &ga);
+	cell_get_cartesian(cell, &av[0], &av[1], &av[2],
+	                         &bv[0], &bv[1], &bv[2],
+	                         &cv[0], &cv[1], &cv[2]);
+
+	/* Find vectors in 'cell' with lengths close to a, b and c */
+	cand_a = find_candidates(a, av, bv, cv, ltl, &ncand_a);
+	cand_b = find_candidates(b, av, bv, cv, ltl, &ncand_b);
+	cand_c = find_candidates(c, av, bv, cv, ltl, &ncand_c);
+
+	STATUS("Candidates for a: %i\n", ncand_a);
+	for ( i=0; i<10; i++ ) {
+		STATUS("%s %s %s\n", rtnl_format(cand_a[3*i+0]),
+		                     rtnl_format(cand_a[3*i+1]),
+		                     rtnl_format(cand_a[3*i+2]));
+	}
+	STATUS("Candidates for b: %i\n", ncand_b);
+	for ( i=0; i<10; i++ ) {
+		STATUS("%s %s %s\n", rtnl_format(cand_b[3*i+0]),
+		                     rtnl_format(cand_b[3*i+1]),
+		                     rtnl_format(cand_b[3*i+2]));
+	}
+	STATUS("Candidates for c: %i\n", ncand_c);
+	for ( i=0; i<10; i++ ) {
+		STATUS("%s %s %s\n", rtnl_format(cand_c[3*i+0]),
+		                     rtnl_format(cand_c[3*i+1]),
+		                     rtnl_format(cand_c[3*i+2]));
+	}
+
+	m = rtnl_mtx_new(3, 3);
+	for ( ia=0; ia<ncand_a; ia++ ) {
+		for ( ib=0; ib<ncand_b; ib++ ) {
+
+			UnitCell *test;
+			double at, bt, ct, alt, bet, gat;
+			int ic = 0;
+
+			/* Form the matrix using the first candidate for c */
+			rtnl_mtx_set(m, 0, 0, cand_a[3*ia+0]);
+			rtnl_mtx_set(m, 0, 1, cand_a[3*ia+1]);
+			rtnl_mtx_set(m, 0, 2, cand_a[3*ia+2]);
+			rtnl_mtx_set(m, 1, 0, cand_b[3*ib+0]);
+			rtnl_mtx_set(m, 1, 1, cand_b[3*ib+1]);
+			rtnl_mtx_set(m, 1, 2, cand_b[3*ib+2]);
+			rtnl_mtx_set(m, 2, 0, cand_c[3*ic+0]);
+			rtnl_mtx_set(m, 2, 1, cand_c[3*ic+1]);
+			rtnl_mtx_set(m, 2, 2, cand_c[3*ic+2]);
+
+			/* Check angle between a and b */
+			test = cell_transform_rational(cell, m);
+			cell_get_parameters(test, &at, &bt, &ct, &alt, &bet, &gat);
+			cell_free(test);
+			if ( fabs(gat - ga) > atl ) continue;
+
+			/* Gamma OK, now look for place for c axis */
+			for ( ic=0; ic<ncand_c; ic++ ) {
+
+				UnitCell *test2;
+
+				rtnl_mtx_set(m, 0, 0, cand_a[3*ia+0]);
+				rtnl_mtx_set(m, 0, 1, cand_a[3*ia+1]);
+				rtnl_mtx_set(m, 0, 2, cand_a[3*ia+2]);
+				rtnl_mtx_set(m, 1, 0, cand_b[3*ib+0]);
+				rtnl_mtx_set(m, 1, 1, cand_b[3*ib+1]);
+				rtnl_mtx_set(m, 1, 2, cand_b[3*ib+2]);
+				rtnl_mtx_set(m, 2, 0, cand_c[3*ic+0]);
+				rtnl_mtx_set(m, 2, 1, cand_c[3*ic+1]);
+				rtnl_mtx_set(m, 2, 2, cand_c[3*ic+2]);
+
+				if ( rtnl_cmp(rtnl_mtx_det(m),rtnl_zero()) == 0 ) continue;
+
+				test = cell_transform_rational(cell, m);
+				cell_get_parameters(test, &at, &bt, &ct, &alt, &bet, &gat);
+				if ( !right_handed(test) ) {
+					cell_free(test);
+					continue;
+				}
+				if ( fabs(alt - al) > atl ) {
+					cell_free(test);
+					continue;
+				}
+				if ( fabs(bet - be) > atl ) {
+					cell_free(test);
+					continue;
+				}
+				rtnl_mtx_print(m);
+				test2 = cell_transform_intmat(test, centering_reference);
+				cell_print(test2);
+				cell_free(test);
+				cell_free(test2);
+
+			}
+		}
+	}
+
+	*pmb = m;
+	return 1;
+}
diff --git a/libcrystfel/src/cell-utils.h b/libcrystfel/src/cell-utils.h
index a1a6fd46..0face43d 100644
--- a/libcrystfel/src/cell-utils.h
+++ b/libcrystfel/src/cell-utils.h
@@ -95,6 +95,10 @@ extern int compare_reindexed_cell_parameters_and_orientation(UnitCell *a,
                                                              double atl,
                                                              IntegerMatrix **pmb);
 
+extern int compare_reindexed_cell_parameters(UnitCell *cell, UnitCell *reference,
+                                             double ltl, double atl,
+                                             RationalMatrix **pmb);
+
 #ifdef __cplusplus
 }
 #endif
diff --git a/libcrystfel/src/rational.c b/libcrystfel/src/rational.c
index 5db83165..0f2bcaee 100644
--- a/libcrystfel/src/rational.c
+++ b/libcrystfel/src/rational.c
@@ -226,6 +226,39 @@ char *rtnl_format(Rational rt)
 }
 
 
+Rational *rtnl_list(signed int num_min, signed int num_max,
+                    signed int den_min, signed int den_max,
+                    int *pn)
+{
+	signed int num, den;
+	Rational *list;
+	int n = 0;
+
+	list = malloc((1+num_max-num_min)*(1+den_max-den_min)*sizeof(Rational));
+	if ( list == NULL ) return NULL;
+
+	for ( num=num_min; num<=num_max; num++ ) {
+		for ( den=den_min; den<=den_max; den++ ) {
+
+			Rational r = rtnl(num, den);
+
+			/* Denominator zero? */
+			if ( den == 0 ) continue;
+
+			/* Same as last entry? */
+			if ( (n>0) && (rtnl_cmp(list[n-1], r)==0) ) continue;
+
+			/* Can be reduced? */
+			if ( gcd(num, den) != 1 ) continue;
+
+			list[n++] = r;
+		}
+	}
+	*pn = n;
+	return list;
+}
+
+
 /**
  * SECTION:rational_matrix
  * @short_description: Rational matrices
diff --git a/libcrystfel/src/rational.h b/libcrystfel/src/rational.h
index 925f9be1..7983f56f 100644
--- a/libcrystfel/src/rational.h
+++ b/libcrystfel/src/rational.h
@@ -77,6 +77,10 @@ extern Rational rtnl_abs(Rational a);
 
 extern char *rtnl_format(Rational rt);
 
+extern Rational *rtnl_list(signed int num_min, signed int num_max,
+                           signed int den_min, signed int den_max,
+                           int *pn);
+
 extern RationalMatrix *rtnl_mtx_new(unsigned int rows, unsigned int cols);
 extern RationalMatrix *rtnl_mtx_copy(const RationalMatrix *m);
 extern Rational rtnl_mtx_get(const RationalMatrix *m, int i, int j);
diff --git a/src/cell_tool.c b/src/cell_tool.c
index 03c58f5e..d6c5b6cd 100644
--- a/src/cell_tool.c
+++ b/src/cell_tool.c
@@ -72,14 +72,9 @@ static void show_help(const char *s)
 static int comparecells(UnitCell *cell, const char *comparecell,
                         double ltl, double atl)
 {
-	signed int i[9];
-	int b[9];
-	const int maxorder = 2;
 	UnitCell *cell2;
 	RationalMatrix *m;
 
-	STATUS("Comparing with: %s\n", comparecell);
-
 	cell2 = load_cell_from_file(comparecell);
 	if ( cell2 == NULL ) {
 		ERROR("Failed to load unit cell from '%s'\n", comparecell);
@@ -92,75 +87,14 @@ static int comparecells(UnitCell *cell, const char *comparecell,
 	STATUS("------------------> The reference unit cell:\n");
 	cell_print(cell2);
 
-	STATUS("Comparing cells up to %ix each lattice length.\n", maxorder);
-	STATUS("Reciprocal axis length tolerance %f %%\n", ltl);
-	STATUS("Reciprocal angle tolerance %f degrees\n", rad2deg(atl));
-	STATUS("This will take about 30 seconds.  Please wait...\n");
-
-	m = rtnl_mtx_new(3, 3);
-	for ( i[0]=-maxorder; i[0]<=+maxorder; i[0]++ ) {
-	for ( i[1]=-maxorder; i[1]<=+maxorder; i[1]++ ) {
-	for ( i[2]=-maxorder; i[2]<=+maxorder; i[2]++ ) {
-	for ( i[3]=-maxorder; i[3]<=+maxorder; i[3]++ ) {
-	for ( i[4]=-maxorder; i[4]<=+maxorder; i[4]++ ) {
-	for ( i[5]=-maxorder; i[5]<=+maxorder; i[5]++ ) {
-	for ( i[6]=-maxorder; i[6]<=+maxorder; i[6]++ ) {
-	for ( i[7]=-maxorder; i[7]<=+maxorder; i[7]++ ) {
-	for ( i[8]=-maxorder; i[8]<=+maxorder; i[8]++ ) {
-	for ( b[0]=0; b[0]<=1; b[0]++ ) {
-	for ( b[1]=0; b[1]<=1; b[1]++ ) {
-	for ( b[2]=0; b[2]<=1; b[2]++ ) {
-	for ( b[3]=0; b[3]<=1; b[3]++ ) {
-	for ( b[4]=0; b[4]<=1; b[4]++ ) {
-	for ( b[5]=0; b[5]<=1; b[5]++ ) {
-	for ( b[6]=0; b[6]<=1; b[6]++ ) {
-	for ( b[7]=0; b[7]<=1; b[7]++ ) {
-	for ( b[8]=0; b[8]<=1; b[8]++ ) {
-
-		UnitCell *nc;
-		int j, k;
-		int l = 0;
-
-		for ( j=0; j<3; j++ ) {
-		for ( k=0; k<3; k++ ) {
-			if ( b[l] || i[l]==1 || i[l]==-1 ) {
-				rtnl_mtx_set(m, j, k, rtnl(i[l], 1));
-			} else {
-				rtnl_mtx_set(m, j, k, rtnl(1, i[l]));
-			}
-			l++;
-		}
-		}
-
-		nc = cell_transform_rational(cell, m);
-
-		if ( compare_cell_parameters(cell2, nc, ltl, atl) ) {
-			STATUS("-----------------------------------------------"
-			       "-------------------------------------------\n");
-			cell_print(nc);
-			rtnl_mtx_print(m);
-		}
-
-		cell_free(nc);
-
-	}
-	}
-	}
-	}
-	}
-	}
-	}
-	}
-	}
-	}
-	}
-	}
-	}
-	}
-	}
-	}
-	}
+	STATUS("------------------> The comparison results:\n");
+	if ( !compare_reindexed_cell_parameters(cell, cell2, ltl, atl, &m) ) {
+		STATUS("No relationship found between lattices.\n");
+		return 0;
+	} else {
+		STATUS("Relationship found:\n");
 	}
+
 	rtnl_mtx_free(m);
 
 	return 0;
-- 
cgit v1.2.3