From beedf1a19a31ecd83f887ceba6e7f08eba6e930b Mon Sep 17 00:00:00 2001 From: Thomas White Date: Wed, 21 Aug 2019 16:27:34 +0200 Subject: Final implementation of Niggli-based cell comparison --- libcrystfel/src/cell-utils.c | 309 +++++++++++++++++++++++++++++++++---------- 1 file changed, 238 insertions(+), 71 deletions(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c index 961a6554..72225a8c 100644 --- a/libcrystfel/src/cell-utils.c +++ b/libcrystfel/src/cell-utils.c @@ -2169,7 +2169,7 @@ static int is_burger(struct g6 g, double eps) } -static int is_niggli(struct g6 g, double eps) +static int UNUSED is_niggli(struct g6 g, double eps) { if ( !is_burger(g, eps) ) return 0; @@ -2188,20 +2188,27 @@ static int is_niggli(struct g6 g, double eps) } -static void debug_lattice(struct g6 g, double eps) +static signed int eps_sign(double v, double eps) +{ + return GT(v, 0.0) ? +1 : -1; +} + + +static void debug_lattice(struct g6 g, double eps, int step) { - STATUS("%e %e %e %e %e %e", +#if 0 + STATUS("After step %i: %e %e %e %e %e %e --", step, g.A/1e-0, g.B/1e-0, g.C/1e-0, g.D/1e-0, g.E/1e-0, g.F/1e-0); if ( is_burger(g, eps) ) STATUS(" B"); if ( is_niggli(g, eps) ) STATUS(" N"); + if ( eps_sign(g.D, eps)*eps_sign(g.E, eps)*eps_sign(g.F, eps) > 0 ) { + STATUS(" I"); + } else { + STATUS(" II"); + } STATUS("\n"); -} - - -static signed int eps_sign(double v, double eps) -{ - return GT(v, 0.0) ? +1 : -1; +#endif } @@ -2218,6 +2225,23 @@ static void mult_in_place(IntegerMatrix *T, IntegerMatrix *M) } +/* Replace the elements of T with those of T*M */ +static void rtnl_mult_in_place(RationalMatrix *T, RationalMatrix *M) +{ + int i, j; + RationalMatrix *tmp; + + tmp = rtnl_mtx_new(3, 3); + rtnl_mtx_mtxmult(T, M, tmp); + for ( i=0; i<3; i++ ) { + for ( j=0; j<3; j++ ) { + rtnl_mtx_set(T, i, j, rtnl_mtx_get(tmp, i, j)); + } + } + rtnl_mtx_free(tmp); +} + + /* Cell volume from G6 components */ static double g6_volume(struct g6 g) { @@ -2238,11 +2262,12 @@ IntegerMatrix *reduce_g6(struct g6 g, double epsrel) eps = pow(g6_volume(g), 1.0/3.0) * epsrel; eps = eps*eps; - STATUS("eps = %e\n", eps); T = intmat_identity(3); M = intmat_new(3, 3); + debug_lattice(g, eps, 0); + do { int done_s1s2; @@ -2250,8 +2275,6 @@ IntegerMatrix *reduce_g6(struct g6 g, double epsrel) do { done_s1s2 = 1; - debug_lattice(g, eps); - if ( GT(g.A, g.B) || (EQ(g.A, g.B) && GT(fabs(g.D), fabs(g.E))) ) { @@ -2266,7 +2289,8 @@ IntegerMatrix *reduce_g6(struct g6 g, double epsrel) intmat_set(M, 2, 2, -1); mult_in_place(T, M); - STATUS("step 1\n"); + debug_lattice(g, eps, 1); + } else if ( GT(g.B, g.C) || (EQ(g.B, g.C) && GT(fabs(g.E), fabs(g.F))) ) @@ -2282,7 +2306,7 @@ IntegerMatrix *reduce_g6(struct g6 g, double epsrel) intmat_set(M, 2, 1, -1); mult_in_place(T, M); - STATUS("step 2\n"); + debug_lattice(g, eps, 2); done_s1s2 = 0; } @@ -2290,72 +2314,62 @@ IntegerMatrix *reduce_g6(struct g6 g, double epsrel) finished = 0; - debug_lattice(g, eps); - /* K-G paper says g3*g4*g3, which I assume is a misprint */ if ( eps_sign(g.D, eps) * eps_sign(g.E, eps) * eps_sign(g.F, eps) > 0 ) { intmat_zero(M); - intmat_set(M, 0, 0, GT(g.D, 0.0) ? 1 : -1); - intmat_set(M, 1, 1, GT(g.E, 0.0) ? 1 : -1); - intmat_set(M, 2, 2, GT(g.F, 0.0) ? 1 : -1); + intmat_set(M, 0, 0, LT(g.D, 0.0) ? -1 : 1); + intmat_set(M, 1, 1, LT(g.E, 0.0) ? -1 : 1); + intmat_set(M, 2, 2, LT(g.F, 0.0) ? -1 : 1); mult_in_place(T, M); g.D = fabs(g.D); g.E = fabs(g.E); g.F = fabs(g.F); - STATUS("step 3\n"); + debug_lattice(g, eps, 3); } else { - int r, c; - int have_rc = 0; - - intmat_zero(M); - intmat_set(M, 0, 0, 1); - intmat_set(M, 1, 1, 1); - intmat_set(M, 2, 2, 1); + int z = 4; + signed int ijk[3] = { 1, 1, 1 }; if ( GT(g.D, 0.0) ) { - intmat_set(M, 0, 0, -1); + ijk[0] = -1; } else if ( !LT(g.D, 0.0) ) { - r = 0; c = 0; - have_rc = 1; + z = 0; } if ( GT(g.E, 0.0) ) { - intmat_set(M, 1, 1, -1); + ijk[1] = -1; } else if ( !LT(g.D, 0.0) ) { - r = 1; c = 1; - have_rc = 1; - } + z = 1; + } if ( GT(g.F, 0.0) ) { - intmat_set(M, 2, 2, -1); + ijk[2] = -1; } else if ( !LT(g.D, 0.0) ) { - r = 2; c = 2; - have_rc = 1; + z = 2; } - if ( LT(g.D*g.E*g.F, 0.0) ) { - assert(have_rc); - if ( have_rc ) { - intmat_set(M, r, c, -1); - } else { - ERROR("Don't have r,c!\n"); - //abort(); + if ( ijk[0]*ijk[1]*ijk[2] < 0 ) { + if ( z == 4 ) { + ERROR("No element in reduction step 4"); + abort(); } + ijk[z] = -1; } + intmat_zero(M); + intmat_set(M, 0, 0, ijk[0]); + intmat_set(M, 1, 1, ijk[1]); + intmat_set(M, 2, 2, ijk[2]); mult_in_place(T, M); g.D = -fabs(g.D); g.E = -fabs(g.E); g.F = -fabs(g.F); - STATUS("step 4\n"); + debug_lattice(g, eps, 4); } - debug_lattice(g, eps); - if ( GT(fabs(g.D), g.B) || (EQ(g.B, g.D) && LT(2.0*g.E, g.F)) || (EQ(g.B, -g.D) && LT(g.F, 0.0)) ) @@ -2373,7 +2387,7 @@ IntegerMatrix *reduce_g6(struct g6 g, double epsrel) g.D = -2*s*g.B + g.D; g.E = g.E - s*g.F; - STATUS("step 5\n"); + debug_lattice(g, eps, 5); } else if ( GT(fabs(g.E), g.A) || (EQ(g.A, g.E) && LT(2.0*g.D, g.F)) @@ -2392,7 +2406,7 @@ IntegerMatrix *reduce_g6(struct g6 g, double epsrel) g.D = g.D - s*g.F; g.E = -2*s*g.A + g.E; - STATUS("step 6\n"); + debug_lattice(g, eps, 6); } else if ( GT(fabs(g.F), g.A) || (EQ(g.A, g.F) && LT(2.0*g.D, g.E)) @@ -2411,10 +2425,10 @@ IntegerMatrix *reduce_g6(struct g6 g, double epsrel) g.D = g.D - s*g.E; g.F = -2*s*g.A + g.F; - STATUS("step 7\n"); + debug_lattice(g, eps, 7); - } else if ( LT(g.A+g.B+g.C+g.D+g.E+g.F, 0.0) - || ( (EQ(g.A+g.B+g.C+g.D+g.E+g.F, 0.0) + } else if ( LT(g.A+g.B+g.D+g.E+g.F, 0.0) /* not g.C */ + || ( (EQ(g.A+g.B+g.D+g.E+g.F, 0.0) && GT(2.0*g.A + 2.0*g.E + g.F, 0.0)) ) ) { intmat_zero(M); @@ -2428,7 +2442,7 @@ IntegerMatrix *reduce_g6(struct g6 g, double epsrel) g.D = 2.0*g.B + g.D + g.F; g.E = 2.0*g.A + g.E + g.F; - STATUS("step 8\n"); + debug_lattice(g, eps, 8); } else { finished = 1; @@ -2436,7 +2450,7 @@ IntegerMatrix *reduce_g6(struct g6 g, double epsrel) } while ( !finished ); - debug_lattice(g, eps); + debug_lattice(g, eps, 99); assert(is_burger(g, eps)); assert(is_niggli(g, eps)); @@ -2446,6 +2460,96 @@ IntegerMatrix *reduce_g6(struct g6 g, double epsrel) } +static double cell_diff(UnitCell *cell, double a, double b, double c, + double al, double be, double ga) +{ + double ta, tb, tc, tal, tbe, tga; + double diff = 0.0; + cell_get_parameters(cell, &ta, &tb, &tc, &tal, &tbe, &tga); + diff += fabs(a - ta); + diff += fabs(b - tb); + diff += fabs(c - tc); + return diff; +} + + +static IntegerMatrix *check_permutations(UnitCell *cell, UnitCell *reference, + IntegerMatrix *RB, IntegerMatrix *CB, + const double *tols) +{ + IntegerMatrix *m; + int i[9]; + double a, b, c, al, be, ga; + double min_dist = +INFINITY; + IntegerMatrix *best_m = NULL; + RationalMatrix *RBCB; + RationalMatrix *rRB; + RationalMatrix *rCB; + + m = intmat_new(3, 3); + cell_get_parameters(reference, &a, &b, &c, &al, &be, &ga); + + RBCB = rtnl_mtx_new(3, 3); + rRB = rtnl_mtx_from_intmat(RB); + rCB = rtnl_mtx_from_intmat(CB); + rtnl_mtx_mtxmult(rRB, rCB, RBCB); + rtnl_mtx_free(rRB); + rtnl_mtx_free(rCB); + + for ( i[0]=-1; i[0]<=+1; i[0]++ ) { + for ( i[1]=-1; i[1]<=+1; i[1]++ ) { + for ( i[2]=-1; i[2]<=+1; i[2]++ ) { + for ( i[3]=-1; i[3]<=+1; i[3]++ ) { + for ( i[4]=-1; i[4]<=+1; i[4]++ ) { + for ( i[5]=-1; i[5]<=+1; i[5]++ ) { + for ( i[6]=-1; i[6]<=+1; i[6]++ ) { + for ( i[7]=-1; i[7]<=+1; i[7]++ ) { + for ( i[8]=-1; i[8]<=+1; i[8]++ ) { + + UnitCell *nc; + int j, k; + int l = 0; + signed int det; + UnitCell *tmp; + + for ( j=0; j<3; j++ ) + for ( k=0; k<3; k++ ) + intmat_set(m, j, k, i[l++]); + + det = intmat_det(m); + if ( det != +1 ) continue; + + tmp = cell_transform_intmat(cell, m); + nc = cell_transform_rational(tmp, RBCB); + cell_free(tmp); + + if ( compare_cell_parameters(nc, reference, tols) ) { + double dist = cell_diff(nc, a, b, c, al, be, ga); + if ( dist < min_dist ) { + intmat_free(best_m); + min_dist = dist; + best_m = intmat_copy(m); + } + } + + cell_free(nc); + + } + } + } + } + } + } + } + } + } + intmat_free(m); + rtnl_mtx_free(RBCB); + + return best_m; +} + + /** * \param cell_in: A UnitCell * \param reference_in: Another UnitCell @@ -2461,7 +2565,9 @@ IntegerMatrix *reduce_g6(struct g6 g, double epsrel) * * Only the cell parameters will be compared. The relative orientations are * irrelevant. The tolerances will be applied to the transformed copy of - * \p cell_in. + * \p cell_in, i.e. the version of the input cell which looks similar to + * \p reference_in. Subject to the tolerances, the cell will be chosen which + * has the lowest total absolute error in unit cell axis lengths. * * This is the right function to use for deciding if an indexing solution * matches a reference cell or not. @@ -2476,37 +2582,98 @@ UnitCell *compare_reindexed_cell_parameters(UnitCell *cell_in, { UnitCell *cell; UnitCell *reference; - IntegerMatrix *CBint; + IntegerMatrix *CB; RationalMatrix *CiA; - RationalMatrix *CB; - IntegerMatrix *Mcell; - IntegerMatrix *Mref; + IntegerMatrix *RA; + IntegerMatrix *RB; + IntegerMatrix *P; + UnitCell *cell_reduced; + UnitCell *reference_reduced; + UnitCell *match; struct g6 g6cell; struct g6 g6ref; + //STATUS("The input cell:\n"); + //cell_print(cell_in); + + //STATUS("The reference cell:\n"); + //cell_print(reference_in); + /* Un-center both cells */ - reference = uncenter_cell(reference_in, &CBint, NULL); - if ( reference == NULL ) return 0; - CB = rtnl_mtx_from_intmat(CBint); - intmat_free(CBint); + reference = uncenter_cell(reference_in, &CB, NULL); + if ( reference == NULL ) return NULL; cell = uncenter_cell(cell_in, NULL, &CiA); - if ( cell == NULL ) return 0; + if ( cell == NULL ) return NULL; /* Calculate G6 components for both */ g6cell = cell_get_G6(cell); g6ref = cell_get_G6(reference); /* Convert both to reduced basis (stably) */ - Mcell = reduce_g6(g6cell, 1e-5); - Mref = reduce_g6(g6ref, 1e-5); + RA = reduce_g6(g6cell, 1e-5); + //STATUS("------------------------------------------\n"); + RB = reduce_g6(g6ref, 1e-5); + + cell_reduced = cell_transform_intmat(cell, RA); - /* Compare cells (including nearby ones, possibly permuting - * axes if close to equality) */ - STATUS("Reduced cell:\n"); + //STATUS("Reduced cell:\n"); + //cell_print(cell_reduced); + //STATUS("Reduced reference:\n"); + reference_reduced = cell_transform_intmat(reference, RB); + //cell_print(reference_reduced); - intmat_free(Mcell); - intmat_free(Mref); + /* The primitive (non-reduced) cells are no longer needed */ cell_free(reference); cell_free(cell); + + /* Within tolerance? */ + P = check_permutations(cell_reduced, reference_in, RB, CB, tols); + if ( P != NULL ) { + + match = cell_transform_intmat(cell_reduced, P); + + //STATUS("Original cell transformed to look like reference:\n"); + //cell_print(match); + + /* Match found. Combined matrix requested? */ + if ( pmb != NULL ) { + + /* Calculate combined matrix: CB.RiB.RA.CiA */ + RationalMatrix *rRA = rtnl_mtx_from_intmat(RA); + RationalMatrix *rP = rtnl_mtx_from_intmat(P); + IntegerMatrix *RiB = intmat_inverse(RB); + RationalMatrix *rRiB = rtnl_mtx_from_intmat(RiB); + RationalMatrix *rCB = rtnl_mtx_from_intmat(CB); + RationalMatrix *comb = rtnl_mtx_identity(3); + + rtnl_mult_in_place(comb, CiA); + rtnl_mult_in_place(comb, rRA); + rtnl_mult_in_place(comb, rP); + rtnl_mult_in_place(comb, rRiB); + rtnl_mult_in_place(comb, rCB); + + rtnl_mtx_free(rRA); + rtnl_mtx_free(rP); + intmat_free(RiB); + rtnl_mtx_free(rRiB); + rtnl_mtx_free(rCB); + + *pmb = comb; + + } + + } else { + match = NULL; + } + + rtnl_mtx_free(CiA); + intmat_free(RA); + intmat_free(RB); + intmat_free(CB); + intmat_free(P); + cell_free(reference_reduced); + cell_free(cell_reduced); + + return match; } -- cgit v1.2.3