From b4e52c888c13287eeca1a6b4222e9f3038db8b51 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Mon, 5 Aug 2019 14:06:13 +0200 Subject: indexamajig: Accept six values for unit cell tolerance --- libcrystfel/src/index.c | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) (limited to 'libcrystfel/src') diff --git a/libcrystfel/src/index.c b/libcrystfel/src/index.c index 6f0d1b8c..0f613e45 100644 --- a/libcrystfel/src/index.c +++ b/libcrystfel/src/index.c @@ -67,7 +67,7 @@ struct _indexingprivate { IndexingFlags flags; UnitCell *target_cell; - float tolerance[4]; + double tolerance[6]; struct taketwo_options *ttopts; struct xgandalf_options *xgandalf_opts; @@ -315,7 +315,7 @@ static void *prepare_method(IndexingMethod *m, UnitCell *cell, IndexingPrivate *setup_indexing(const char *method_list, UnitCell *cell, - struct detector *det, float *ltl, + struct detector *det, float *tols, IndexingFlags flags, struct taketwo_options *ttopts, struct xgandalf_options *xgandalf_opts, @@ -439,7 +439,7 @@ IndexingPrivate *setup_indexing(const char *method_list, UnitCell *cell, } else { ipriv->target_cell = NULL; } - for ( i=0; i<4; i++ ) ipriv->tolerance[i] = ltl[i]; + for ( i=0; i<6; i++ ) ipriv->tolerance[i] = tols[i]; ipriv->ttopts = ttopts; ipriv->xgandalf_opts = xgandalf_opts; -- cgit v1.2.3 From ff97bef6953ca6de47d6fe8076e2d3444e48dd1c Mon Sep 17 00:00:00 2001 From: Thomas White Date: Mon, 5 Aug 2019 14:24:11 +0200 Subject: compare_cell_parameters and compare_reindexed_cell_parameters: Accept 6 tolerances --- libcrystfel/src/cell-utils.c | 46 +++++++++++++++++++++----------------------- libcrystfel/src/cell-utils.h | 4 ++-- 2 files changed, 24 insertions(+), 26 deletions(-) (limited to 'libcrystfel/src') diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c index f6e189ad..3f3ce252 100644 --- a/libcrystfel/src/cell-utils.c +++ b/libcrystfel/src/cell-utils.c @@ -1671,12 +1671,11 @@ double cell_get_volume(UnitCell *cell) /** * \param cell1: A UnitCell * \param cell2: Another UnitCell - * \param ltl: Maximum allowable fractional difference in axis lengths - * \param atl: Maximum allowable difference in reciprocal angles (in radians) + * \param tolerance: Pointer to tolerances for a,b,c (fractional), al,be,ga (radians) * * Compare the two unit cells. If the real space parameters match to within - * fractional difference \p ltl, and the inter-axial angles match within \p atl, - * and the centering matches, this function returns 1. Otherwise 0. + * the specified tolerances, and the centering matches, this function returns 1. + * Otherwise 0. * * This function considers the cell parameters and centering, but ignores the * orientation of the cell. If you want to compare the orientation as well, @@ -1685,8 +1684,7 @@ double cell_get_volume(UnitCell *cell) * \returns non-zero if the cells match. * */ -int compare_cell_parameters(UnitCell *cell1, UnitCell *cell2, - float ltl, float atl) +int compare_cell_parameters(UnitCell *cell1, UnitCell *cell2, double *tolerance) { double a1, b1, c1, al1, be1, ga1; double a2, b2, c2, al2, be2, ga2; @@ -1698,12 +1696,12 @@ int compare_cell_parameters(UnitCell *cell1, UnitCell *cell2, cell_get_parameters(cell1, &a1, &b1, &c1, &al1, &be1, &ga1); cell_get_parameters(cell2, &a2, &b2, &c2, &al2, &be2, &ga2); - if ( !within_tolerance(a1, a2, ltl*100.0) ) return 0; - if ( !within_tolerance(b1, b2, ltl*100.0) ) return 0; - if ( !within_tolerance(c1, c2, ltl*100.0) ) return 0; - if ( fabs(al1-al2) > atl ) return 0; - if ( fabs(be1-be2) > atl ) return 0; - if ( fabs(ga1-ga2) > atl ) return 0; + if ( !within_tolerance(a1, a2, tolerance[0]*100.0) ) return 0; + if ( !within_tolerance(b1, b2, tolerance[1]*100.0) ) return 0; + if ( !within_tolerance(c1, c2, tolerance[2]*100.0) ) return 0; + if ( fabs(al1-al2) > tolerance[3] ) return 0; + if ( fabs(be1-be2) > tolerance[4] ) return 0; + if ( fabs(ga1-ga2) > tolerance[5] ) return 0; return 1; } @@ -1960,15 +1958,15 @@ static double g6_distance(double a1, double b1, double c1, /** * \param cell_in: A UnitCell * \param reference_in: Another UnitCell - * \param ltl: Maximum allowable fractional difference in direct-space axis lengths - * \param atl: Maximum allowable difference in direct-space angles (in radians) + * \param tolerance: Pointer to tolerances for a,b,c (fractional), al,be,ga (radians) * \param csl: Non-zero to look for coincidence site lattice relationships * \param pmb: Place to store pointer to matrix * * Compare the \p cell_in with \p reference_in. If \p cell is a derivative lattice - * of \p reference, within fractional axis length difference \p ltl and absolute angle - * difference \p atl (in radians), this function returns non-zero and stores the - * transformation which needs to be applied to \p cell_in at \p pmb. + * of \p reference, within fractional axis length differences \p tolerance[0..2] + * and absolute angle difference \p atl[3..5] (in radians), this function returns + * non-zero and stores the transformation which needs to be applied to + * \p cell_in at \p pmb. * * Note that the tolerances will be applied to the primitive unit cell. If * the reference cell is centered, a primitive unit cell will first be calculated. @@ -1989,7 +1987,7 @@ static double g6_distance(double a1, double b1, double c1, * */ int compare_reindexed_cell_parameters(UnitCell *cell_in, UnitCell *reference_in, - double ltl, double atl, int csl, + double *tolerance, int csl, RationalMatrix **pmb) { UnitCell *cell; @@ -2026,9 +2024,9 @@ int compare_reindexed_cell_parameters(UnitCell *cell_in, UnitCell *reference_in, &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, csl, &ncand_a); - cand_b = find_candidates(b, av, bv, cv, ltl, csl, &ncand_b); - cand_c = find_candidates(c, av, bv, cv, ltl, csl, &ncand_c); + cand_a = find_candidates(a, av, bv, cv, tolerance[0], csl, &ncand_a); + cand_b = find_candidates(b, av, bv, cv, tolerance[1], csl, &ncand_b); + cand_c = find_candidates(c, av, bv, cv, tolerance[2], csl, &ncand_c); if ( (ncand_a==0) || (ncand_b==0) || (ncand_c==0) ) { *pmb = NULL; @@ -2065,7 +2063,7 @@ int compare_reindexed_cell_parameters(UnitCell *cell_in, UnitCell *reference_in, 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; + if ( fabs(gat - ga) > tolerance[5] ) continue; /* Gamma OK, now look for place for c axis */ for ( ic=0; ic atl ) { + if ( fabs(alt - al) > tolerance[3] ) { cell_free(test); continue; } - if ( fabs(bet - be) > atl ) { + if ( fabs(bet - be) > tolerance[4] ) { cell_free(test); continue; } diff --git a/libcrystfel/src/cell-utils.h b/libcrystfel/src/cell-utils.h index be878c10..42e7214b 100644 --- a/libcrystfel/src/cell-utils.h +++ b/libcrystfel/src/cell-utils.h @@ -88,7 +88,7 @@ extern int forbidden_reflection(UnitCell *cell, extern double cell_get_volume(UnitCell *cell); extern int compare_cell_parameters(UnitCell *cell1, UnitCell *cell2, - float ltl, float atl); + double *tolerance); extern int compare_cell_parameters_and_orientation(UnitCell *cell1, @@ -103,7 +103,7 @@ extern int compare_reindexed_cell_parameters_and_orientation(UnitCell *a, IntegerMatrix **pmb); extern int compare_reindexed_cell_parameters(UnitCell *cell, UnitCell *reference, - double ltl, double atl, int csl, + double *tolerance, int csl, RationalMatrix **pmb); #ifdef __cplusplus -- cgit v1.2.3 From b9bfccba1564f365581639a1a6856ddff682a13c Mon Sep 17 00:00:00 2001 From: Thomas White Date: Mon, 5 Aug 2019 14:24:44 +0200 Subject: Add compare_permuted_cell_parameters --- libcrystfel/src/cell-utils.c | 94 ++++++++++++++++++++++++++++++++++++++++++++ libcrystfel/src/cell-utils.h | 2 + 2 files changed, 96 insertions(+) (limited to 'libcrystfel/src') diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c index 3f3ce252..b78c4831 100644 --- a/libcrystfel/src/cell-utils.c +++ b/libcrystfel/src/cell-utils.c @@ -2128,3 +2128,97 @@ int compare_reindexed_cell_parameters(UnitCell *cell_in, UnitCell *reference_in, *pmb = CiAMCB; return 1; } + + +/** + * \param cell_in: A UnitCell + * \param reference_in: Another UnitCell + * \param tolerance: Pointer to tolerances for a,b,c (fractional), al,be,ga (radians) + * \param pmb: Place to store pointer to matrix + * + * Compare \p cell_in with \p reference_in. If they are the same, + * within \p tolerance, taking into account possible permutations of the axes, + * this function returns non-zero and stores the transformation which needs to + * be applied to \p cell_in at \p pmb. + * + * Subject to the tolerances, this function will find the transformation which + * gives the best match to the reference cell, using the Euclidian norm in + * G6 [see e.g. Andrews and Bernstein, Acta Cryst. A44 (1988) p1009]. + * + * 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_permuted_cell_parameters(UnitCell *cell, UnitCell *reference, + double *tolerance, IntegerMatrix **pmb) +{ + IntegerMatrix *m; + signed int i[9]; + double a, b, c, al, be, ga; + double min_dist = +INFINITY; + + m = intmat_new(3, 3); + cell_get_parameters(reference, &a, &b, &c, &al, &be, &ga); + *pmb = NULL; + + 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 *test; + int j, k; + int l = 0; + + for ( j=0; j<3; j++ ) + for ( k=0; k<3; k++ ) + intmat_set(m, j, k, i[l++]); + + if ( intmat_det(m) != +1 ) continue; + + test = cell_transform_intmat(cell, m); + + if ( compare_cell_parameters(reference, test, tolerance) ) { + + double at, bt, ct, alt, bet, gat; + double dist; + + cell_get_parameters(test, &at, &bt, &ct, &alt, &bet, &gat); + dist = g6_distance(at, bt, ct, alt, bet, gat, + a, b, c, al, be, ga); + if ( dist < min_dist ) { + min_dist = dist; + intmat_free(*pmb); + *pmb = intmat_copy(m); + } + + } + + cell_free(test); + + } + } + } + } + } + } + } + } + } + + intmat_free(m); + + if ( isinf(min_dist) ) { + return 0; + } + + /* Solution found */ + return 1; +} diff --git a/libcrystfel/src/cell-utils.h b/libcrystfel/src/cell-utils.h index 42e7214b..a489a6d8 100644 --- a/libcrystfel/src/cell-utils.h +++ b/libcrystfel/src/cell-utils.h @@ -106,6 +106,8 @@ extern int compare_reindexed_cell_parameters(UnitCell *cell, UnitCell *reference double *tolerance, int csl, RationalMatrix **pmb); +extern int compare_permuted_cell_parameters(UnitCell *cell, UnitCell *reference, + double *tolerance, IntegerMatrix **pmb); #ifdef __cplusplus } #endif -- cgit v1.2.3 From ea74744a59728bd375dda10a9386b502eb5056c2 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Mon, 5 Aug 2019 14:25:22 +0200 Subject: Use new cell comparison functions for indexing --- libcrystfel/src/index.c | 45 ++++++++++++++++++++++++++------------------- 1 file changed, 26 insertions(+), 19 deletions(-) (limited to 'libcrystfel/src') diff --git a/libcrystfel/src/index.c b/libcrystfel/src/index.c index 0f613e45..b187e119 100644 --- a/libcrystfel/src/index.c +++ b/libcrystfel/src/index.c @@ -543,32 +543,39 @@ void map_all_peaks(struct image *image) static int check_cell(IndexingFlags flags, Crystal *cr, UnitCell *target, - float *tolerance) + double *tolerance) { - if ( (flags & INDEXING_CHECK_CELL_COMBINATIONS) - || (flags & INDEXING_CHECK_CELL_AXES) ) - { - UnitCell *out; - int reduce; - - if ( flags & INDEXING_CHECK_CELL_COMBINATIONS ) - { - reduce = 1; - } else { - reduce = 0; - } + UnitCell *out; + IntegerMatrix *im; + RationalMatrix *rm; - out = match_cell(crystal_get_cell(cr), - target, 0, tolerance, reduce); + /* Check at all? */ + if ( ! ((flags & INDEXING_CHECK_CELL_COMBINATIONS) + || (flags & INDEXING_CHECK_CELL_AXES)) ) return 0; - if ( out == NULL ) { - return 1; - } + if ( compare_permuted_cell_parameters(crystal_get_cell(cr), target, + tolerance, &im) ) + { + out = cell_transform_intmat(crystal_get_cell(cr), im); + cell_free(crystal_get_cell(cr)); + crystal_set_cell(cr, out); + intmat_free(im); + return 0; + } + if ( (flags & INDEXING_CHECK_CELL_COMBINATIONS ) + && compare_reindexed_cell_parameters(crystal_get_cell(cr), target, + tolerance, 0, &rm) ) + { + out = cell_transform_rational(crystal_get_cell(cr), rm); cell_free(crystal_get_cell(cr)); crystal_set_cell(cr, out); + rtnl_mtx_free(rm); + return 0; } - return 0; + + return 1; + } -- cgit v1.2.3 From f6c9932aa1b41e36178a2471e6728c04fb45289a Mon Sep 17 00:00:00 2001 From: Thomas White Date: Thu, 8 Aug 2019 11:05:09 +0200 Subject: Clarify cell vs reference in comparison functions. Also reindexed vs permuted --- libcrystfel/src/cell-utils.c | 61 +++++++++++++++++++++++--------------------- libcrystfel/src/cell-utils.h | 20 +++++++-------- libcrystfel/src/index.c | 7 ++--- 3 files changed, 46 insertions(+), 42 deletions(-) (limited to 'libcrystfel/src') diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c index b78c4831..2580f3fb 100644 --- a/libcrystfel/src/cell-utils.c +++ b/libcrystfel/src/cell-utils.c @@ -1669,8 +1669,8 @@ double cell_get_volume(UnitCell *cell) /** - * \param cell1: A UnitCell - * \param cell2: Another UnitCell + * \param cell: A UnitCell + * \param reference: Another UnitCell * \param tolerance: Pointer to tolerances for a,b,c (fractional), al,be,ga (radians) * * Compare the two unit cells. If the real space parameters match to within @@ -1684,17 +1684,17 @@ double cell_get_volume(UnitCell *cell) * \returns non-zero if the cells match. * */ -int compare_cell_parameters(UnitCell *cell1, UnitCell *cell2, double *tolerance) +int compare_cell_parameters(UnitCell *cell, UnitCell *reference, double *tolerance) { double a1, b1, c1, al1, be1, ga1; double a2, b2, c2, al2, be2, ga2; /* Centering must match: we don't arbitrarte primitive vs centered, * different cell choices etc */ - if ( cell_get_centering(cell1) != cell_get_centering(cell2) ) return 0; + if ( cell_get_centering(cell) != cell_get_centering(reference) ) return 0; - cell_get_parameters(cell1, &a1, &b1, &c1, &al1, &be1, &ga1); - cell_get_parameters(cell2, &a2, &b2, &c2, &al2, &be2, &ga2); + cell_get_parameters(cell, &a1, &b1, &c1, &al1, &be1, &ga1); + cell_get_parameters(reference, &a2, &b2, &c2, &al2, &be2, &ga2); if ( !within_tolerance(a1, a2, tolerance[0]*100.0) ) return 0; if ( !within_tolerance(b1, b2, tolerance[1]*100.0) ) return 0; @@ -1717,8 +1717,8 @@ static double moduli_check(double ax, double ay, double az, /** - * \param cell1: A UnitCell - * \param cell2: Another UnitCell + * \param cell: A UnitCell + * \param reference: Another UnitCell * \param ltl: Maximum allowable fractional difference in axis lengths * \param atl: Maximum allowable difference in angles (in radians) * @@ -1730,27 +1730,27 @@ static double moduli_check(double ax, double ay, double az, * If you just want to see if the parameters are the same, use * compare_cell_parameters() instead. * - * The cells \p a and \p b must have the same centering. Otherwise, this function - * always returns zero. + * \p cell and \p reference must have the same centering. Otherwise, this + * function always returns zero. * * \returns non-zero if the cells match. * */ -int compare_cell_parameters_and_orientation(UnitCell *cell1, UnitCell *cell2, +int compare_cell_parameters_and_orientation(UnitCell *cell, UnitCell *reference, const double ltl, const double atl) { double ax1, ay1, az1, bx1, by1, bz1, cx1, cy1, cz1; double ax2, ay2, az2, bx2, by2, bz2, cx2, cy2, cz2; - if ( cell_get_centering(cell1) != cell_get_centering(cell2) ) return 0; + if ( cell_get_centering(cell) != cell_get_centering(reference) ) return 0; - cell_get_cartesian(cell1, &ax1, &ay1, &az1, - &bx1, &by1, &bz1, - &cx1, &cy1, &cz1); + cell_get_cartesian(cell, &ax1, &ay1, &az1, + &bx1, &by1, &bz1, + &cx1, &cy1, &cz1); - cell_get_cartesian(cell2, &ax2, &ay2, &az2, - &bx2, &by2, &bz2, - &cx2, &cy2, &cz2); + cell_get_cartesian(reference, &ax2, &ay2, &az2, + &bx2, &by2, &bz2, + &cx2, &cy2, &cz2); if ( angle_between(ax1, ay1, az1, ax2, ay2, az2) > atl ) return 0; if ( angle_between(bx1, by1, bz1, bx2, by2, bz2) > atl ) return 0; @@ -1765,8 +1765,8 @@ int compare_cell_parameters_and_orientation(UnitCell *cell1, UnitCell *cell2, /** - * \param a: A UnitCell - * \param b: Another UnitCell + * \param cell: A UnitCell + * \param reference: Another UnitCell * \param ltl: Maximum allowable fractional difference in reciprocal axis lengths * \param atl: Maximum allowable difference in reciprocal angles (in radians) * \param pmb: Place to store pointer to matrix @@ -1774,26 +1774,27 @@ int compare_cell_parameters_and_orientation(UnitCell *cell1, UnitCell *cell2, * Compare the two unit cells. If, using any permutation of the axes, the * axes can be made to match in length (to within fractional difference \p ltl) * and the axes aligned to within \p atl radians, this function returns non-zero - * and stores the transformation to map \p b onto \p a. + * and stores the transformation which must be applied to \p cell at \p pmb. * * Note that the orientations of the cells must match, not just the parameters. * The comparison is done after reindexing using * compare_cell_parameters_and_orientation(). * - * The cells \p a and \p b must have the same centering. Otherwise, this function - * always returns zero. + * \p cell and \p reference must have the same centering. Otherwise, this + * function always returns zero. * * \returns non-zero if the cells match. * */ -int compare_reindexed_cell_parameters_and_orientation(UnitCell *a, UnitCell *b, - double ltl, double atl, - IntegerMatrix **pmb) +int compare_permuted_cell_parameters_and_orientation(UnitCell *cell, + UnitCell *reference, + double ltl, double atl, + IntegerMatrix **pmb) { IntegerMatrix *m; int i[9]; - if ( cell_get_centering(a) != cell_get_centering(b) ) return 0; + if ( cell_get_centering(cell) != cell_get_centering(reference) ) return 0; m = intmat_new(3, 3); @@ -1817,9 +1818,11 @@ int compare_reindexed_cell_parameters_and_orientation(UnitCell *a, UnitCell *b, if ( intmat_det(m) != +1 ) continue; - nc = cell_transform_intmat(b, m); + nc = cell_transform_intmat(cell, m); - if ( compare_cell_parameters_and_orientation(a, nc, ltl, atl) ) { + if ( compare_cell_parameters_and_orientation(nc, reference, + ltl, atl) ) + { if ( pmb != NULL ) *pmb = m; cell_free(nc); return 1; diff --git a/libcrystfel/src/cell-utils.h b/libcrystfel/src/cell-utils.h index a489a6d8..fa577269 100644 --- a/libcrystfel/src/cell-utils.h +++ b/libcrystfel/src/cell-utils.h @@ -87,27 +87,27 @@ extern int forbidden_reflection(UnitCell *cell, extern double cell_get_volume(UnitCell *cell); -extern int compare_cell_parameters(UnitCell *cell1, UnitCell *cell2, +extern int compare_cell_parameters(UnitCell *cell, UnitCell *reference, double *tolerance); +extern int compare_permuted_cell_parameters(UnitCell *cell, UnitCell *reference, + double *tolerance, IntegerMatrix **pmb); -extern int compare_cell_parameters_and_orientation(UnitCell *cell1, - UnitCell *cell2, +extern int compare_cell_parameters_and_orientation(UnitCell *cell, + UnitCell *reference, const double ltl, const double atl); -extern int compare_reindexed_cell_parameters_and_orientation(UnitCell *a, - UnitCell *b, - double ltl, - double atl, - IntegerMatrix **pmb); +extern int compare_permuted_cell_parameters_and_orientation(UnitCell *cell, + UnitCell *reference, + double ltl, + double atl, + IntegerMatrix **pmb); extern int compare_reindexed_cell_parameters(UnitCell *cell, UnitCell *reference, double *tolerance, int csl, RationalMatrix **pmb); -extern int compare_permuted_cell_parameters(UnitCell *cell, UnitCell *reference, - double *tolerance, IntegerMatrix **pmb); #ifdef __cplusplus } #endif diff --git a/libcrystfel/src/index.c b/libcrystfel/src/index.c index b187e119..b07dd067 100644 --- a/libcrystfel/src/index.c +++ b/libcrystfel/src/index.c @@ -542,6 +542,7 @@ void map_all_peaks(struct image *image) } +/* Return 0 for cell OK, 1 for cell incorrect */ static int check_cell(IndexingFlags flags, Crystal *cr, UnitCell *target, double *tolerance) { @@ -708,9 +709,9 @@ static int try_indexer(struct image *image, IndexingMethod indm, /* Don't do similarity check against bad crystals */ if ( crystal_get_user_flag(that_cr) ) continue; - if ( compare_reindexed_cell_parameters_and_orientation(crystal_get_cell(cr), - crystal_get_cell(that_cr), - 0.1, deg2rad(0.5), NULL) ) + if ( compare_permuted_cell_parameters_and_orientation(crystal_get_cell(cr), + crystal_get_cell(that_cr), + 0.1, deg2rad(0.5), NULL) ) { crystal_set_user_flag(cr, 1); } -- cgit v1.2.3 From fde5da00b76d0ae0d3d71b2596f9b4354ba9ff89 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Thu, 8 Aug 2019 11:05:52 +0200 Subject: compare_permuted_cell_parameters{_and_orientation}: Accept wrong-handed cell --- libcrystfel/src/cell-utils.c | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) (limited to 'libcrystfel/src') diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c index 2580f3fb..b207afe7 100644 --- a/libcrystfel/src/cell-utils.c +++ b/libcrystfel/src/cell-utils.c @@ -1811,12 +1811,14 @@ int compare_permuted_cell_parameters_and_orientation(UnitCell *cell, UnitCell *nc; int j, k; int l = 0; + signed int det; for ( j=0; j<3; j++ ) for ( k=0; k<3; k++ ) intmat_set(m, j, k, i[l++]); - if ( intmat_det(m) != +1 ) continue; + det = intmat_det(m); + if ( (det != +1) && (det != -1) ) continue; nc = cell_transform_intmat(cell, m); @@ -2179,12 +2181,14 @@ int compare_permuted_cell_parameters(UnitCell *cell, UnitCell *reference, UnitCell *test; int j, k; int l = 0; + signed int det; for ( j=0; j<3; j++ ) for ( k=0; k<3; k++ ) intmat_set(m, j, k, i[l++]); - if ( intmat_det(m) != +1 ) continue; + det = intmat_det(m); + if ( (det != +1) && (det != -1) ) continue; test = cell_transform_intmat(cell, m); -- cgit v1.2.3 From 9524362d96c79f3a396686b9400e7baed4f4967a Mon Sep 17 00:00:00 2001 From: Thomas White Date: Wed, 14 Aug 2019 14:06:22 +0200 Subject: g6_distance: Take UnitCell --- libcrystfel/src/cell-utils.c | 17 +++++++---------- 1 file changed, 7 insertions(+), 10 deletions(-) (limited to 'libcrystfel/src') diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c index b207afe7..b7662fb3 100644 --- a/libcrystfel/src/cell-utils.c +++ b/libcrystfel/src/cell-utils.c @@ -1931,9 +1931,10 @@ static Rational *find_candidates(double len, double *a, double *b, double *c, } -static void g6_components(double *g6, double a, double b, double c, - double al, double be, double ga) +void g6_components(double *g6, UnitCell *cell) { + double a, b, c, al, be, ga; + cell_get_parameters(cell, &a, &b, &c, &al, &be, &ga); g6[0] = a*a; g6[1] = b*b; g6[2] = c*c; @@ -1943,16 +1944,13 @@ static void g6_components(double *g6, double a, double b, double c, } -static double g6_distance(double a1, double b1, double c1, - double al1, double be1, double ga1, - double a2, double b2, double c2, - double al2, double be2, double ga2) +static double g6_distance(UnitCell *cell1, UnitCell *cell2) { double g1[6], g2[6]; int i; double total = 0.0; - g6_components(g1, a1, b1, c1, al1, be1, ga1); - g6_components(g2, a2, b2, c2, al2, be2, ga2); + g6_components(g1, cell1); + g6_components(g2, cell2); for ( i=0; i<6; i++ ) { total += (g1[i]-g2[i])*(g1[i]-g2[i]); } @@ -2103,8 +2101,7 @@ int compare_reindexed_cell_parameters(UnitCell *cell_in, UnitCell *reference_in, continue; } - dist = g6_distance(at, bt, ct, alt, bet, gat, - a, b, c, al, be, ga); + dist = g6_distance(test, reference); if ( dist < min_dist ) { min_dist = dist; rtnl_mtx_mtxmult(M, CB, MCB); -- cgit v1.2.3 From 0b16a4aa50bfe233d81077b1661c57b4b0ef0ef2 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Wed, 14 Aug 2019 14:06:57 +0200 Subject: G6 cell reduction Doesn't work yet --- libcrystfel/src/cell-utils.c | 203 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 203 insertions(+) (limited to 'libcrystfel/src') diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c index b7662fb3..141eaaab 100644 --- a/libcrystfel/src/cell-utils.c +++ b/libcrystfel/src/cell-utils.c @@ -2132,6 +2132,209 @@ int compare_reindexed_cell_parameters(UnitCell *cell_in, UnitCell *reference_in, } +/* Criteria from Grosse-Kunstleve et al., Acta Cryst A60 (2004) p1-6 */ +#define GT(x,y) (y < x-eps) +#define LT(x,y) (x < y-eps) +#define EQ(x,y) !(GT(x,y) || LT(x,y)) + + +static void do_R5(IntegerMatrix *R, double *g, double eps) +{ + signed int s = GT(g[3], 0.0) ? 1 : -1; + intmat_set(R, 0, 0, 1); + intmat_set(R, 1, 1, 1); + intmat_set(R, 2, 2, 1); + intmat_set(R, 1, 2, -s); + g[2] = g[1] + g[2] -s*g[3]; + g[3] = -2*s*g[1] + g[3]; + g[4] = g[4] - s*g[5]; + STATUS("R5\n"); +} + + +static void do_R6(IntegerMatrix *R, double *g, double eps) +{ + signed int s = GT(g[4], 0.0) ? 1 : -1; + intmat_set(R, 0, 0, 1); + intmat_set(R, 1, 1, 1); + intmat_set(R, 2, 2, 1); + intmat_set(R, 0, 2, -s); + g[2] = g[0] + g[2] -s*g[4]; + g[3] = g[3] - s*g[5]; + g[4] = -2*s*g[0] + g[4]; + STATUS("R6\n"); +} + + +static void do_R7(IntegerMatrix *R, double *g, double eps) +{ + signed int s = GT(g[5], 0.0) ? 1 : -1; + intmat_set(R, 0, 0, 1); + intmat_set(R, 1, 1, 1); + intmat_set(R, 2, 2, 1); + intmat_set(R, 0, 1, -s); + g[1] = g[0] + g[1] -s*g[5]; + g[3] = g[3] - s*g[4]; + g[5] = -2*s*g[0] + g[5]; + STATUS("R7\n"); +} + + +static void do_R8(IntegerMatrix *R, double *g, double eps) +{ + intmat_set(R, 0, 0, 1); + intmat_set(R, 1, 1, 1); + intmat_set(R, 2, 2, 1); + intmat_set(R, 1, 2, 1); + g[2] = g[0]+g[1]+g[2]+g[3]+g[4]+g[5]; + g[3] = 2.0*g[1] + g[3] + g[5]; + g[4] = 2.0*g[0] + g[4] + g[5]; + STATUS("R8\n"); +} + + +IntegerMatrix *reduce_g6(double *g, double eps) +{ + IntegerMatrix *T; + int finished = 0; + + T = intmat_identity(3); + STATUS("eps = %e\n", eps); + + do { + + IntegerMatrix *tmp; + IntegerMatrix *M; + IntegerMatrix *R; + + M = intmat_new(3, 3); + R = intmat_new(3, 3); + + int did_something; + do { + + did_something = 0; + + if ( GT(g[0], g[1]) || (EQ(g[0], g[1]) && GT(g[3], g[4])) ) { + + /* Swap a and b */ + double temp; + temp = g[0]; g[0] = g[1]; g[1] = temp; + temp = g[3]; g[3] = g[4]; g[4] = temp; + intmat_set(M, 1, 0, -1); + intmat_set(M, 0, 1, -1); + intmat_set(M, 2, 2, -1); + STATUS("SP1\n"); + did_something = 1; + + } else if ( GT(g[1], g[2]) + || (EQ(g[1], g[2]) && GT(g[4], g[5])) ) { + + /* Swap b and c */ + double temp; + temp = g[1]; g[1] = g[2]; g[2] = temp; + temp = g[4]; g[4] = g[5]; g[5] = temp; + intmat_set(M, 0, 0, -1); + intmat_set(M, 1, 2, -1); + intmat_set(M, 2, 1, -1); + STATUS("SP2\n"); + did_something = 1; + + } + + } while ( did_something ); + + if ( GT(g[3]*g[4]*g[5], 0.0) ) { + + intmat_set(M, 0, 0, GT(g[3], 0.0) ? 1 : -1); + intmat_set(M, 1, 1, GT(g[4], 0.0) ? 1 : -1); + intmat_set(M, 2, 2, GT(g[5], 0.0) ? 1 : -1); + g[3] = fabs(g[3]); + g[4] = fabs(g[4]); + g[5] = fabs(g[5]); + STATUS("SP3\n"); + + } else { + + int r, c; + int have_rc = 0; + intmat_set(M, 0, 0, 1); + intmat_set(M, 1, 1, 1); + intmat_set(M, 2, 2, 1); + if ( GT(g[3], 0.0) ) { + intmat_set(M, 0, 0, -1); + } else if ( !LT(g[3], 0.0) ) { + r = 0; c = 0; + have_rc = 1; + } + if ( GT(g[4], 0.0) ) { + intmat_set(M, 1, 1, -1); + } else if ( !LT(g[3], 0.0) ) { + r = 1; c = 1; + have_rc = 1; + } + if ( GT(g[5], 0.0) ) { + intmat_set(M, 2, 2, -1); + } else if ( !LT(g[3], 0.0) ) { + r = 2; c = 2; + have_rc = 1; + } + if ( LT(g[3]*g[4]*g[5], 0.0) ) { + assert(have_rc); + if ( have_rc ) { + intmat_set(M, r, c, -1); + } else { + ERROR("Don't have r,c!\n"); + abort(); + } + } + g[3] = -fabs(g[3]); + g[4] = -fabs(g[4]); + g[5] = -fabs(g[5]); + STATUS("SP4\n"); + + } + + STATUS("G6: %f %f %f %f %f %f\n", g[0], g[1], g[2], g[3], g[4], g[5]); + + if ( GT(fabs(g[3]), g[1]) ) { + do_R5(R, g, eps); + } else if ( GT(fabs(g[5]), g[0]) ) { + do_R6(R, g, eps); + } else if ( GT(fabs(g[6]), g[0]) ) { + do_R7(R, g, eps); + } else if ( LT(g[0]+g[1]+g[2]+g[3]+g[4]+g[5], 0.0) ) { + do_R8(R, g, eps); + } else if ( (EQ(g[1], g[3]) && LT(2.0*g[4], g[4])) + || (EQ(g[1], -g[3]) && LT(g[5], 0.0) ) ) { + do_R5(R, g, eps); + } else if ( (EQ(g[0], g[4]) && LT(2.0*g[3], g[5])) + || (EQ(g[0], -g[4]) && LT(g[5], 0.0) ) ) { + do_R6(R, g, eps); + } else if ( (EQ(g[0], g[5]) && LT(2.0*g[3], g[4])) + || (EQ(g[0], -g[5]) && LT(g[6], 0.0) ) ) { + do_R7(R, g, eps); + } else if ( EQ(g[0]+g[1]+g[2]+g[3]+g[4]+g[5], 0.0) + && GT(2.0*g[0] + 2.0*g[4] + g[5], 0.0) ) { + do_R8(R, g, eps); + } else { + finished = 1; + R = intmat_identity(3); + STATUS("Done.\n"); + } + + STATUS("G6: %f %f %f %f %f %f\n", g[0], g[1], g[2], g[3], g[4], g[5]); + tmp = intmat_intmat_mult(T, M); + intmat_free(T); + T = intmat_intmat_mult(tmp, R); + intmat_free(tmp); + + } while ( !finished ); + + return T; +} + + /** * \param cell_in: A UnitCell * \param reference_in: Another UnitCell -- cgit v1.2.3 From e3f4046056cf92ce5e40e5e715cbcd53df7b0621 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Wed, 14 Aug 2019 14:07:31 +0200 Subject: Framework for new unit cell comparison function --- libcrystfel/src/cell-utils.c | 105 +++++++++++++++---------------------------- libcrystfel/src/cell-utils.h | 6 +-- libcrystfel/src/index.c | 16 +------ 3 files changed, 40 insertions(+), 87 deletions(-) (limited to 'libcrystfel/src') diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c index 141eaaab..20ede189 100644 --- a/libcrystfel/src/cell-utils.c +++ b/libcrystfel/src/cell-utils.c @@ -2341,10 +2341,9 @@ IntegerMatrix *reduce_g6(double *g, double eps) * \param tolerance: Pointer to tolerances for a,b,c (fractional), al,be,ga (radians) * \param pmb: Place to store pointer to matrix * - * Compare \p cell_in with \p reference_in. If they are the same, - * within \p tolerance, taking into account possible permutations of the axes, - * this function returns non-zero and stores the transformation which needs to - * be applied to \p cell_in at \p pmb. + * Compare the \p cell_in with \p reference_in. If they represent the same + * lattice, this function returns non-zero and stores the transformation which + * needs to be applied to \p cell_in at \p pmb. * * Subject to the tolerances, this function will find the transformation which * gives the best match to the reference cell, using the Euclidian norm in @@ -2356,76 +2355,44 @@ IntegerMatrix *reduce_g6(double *g, double eps) * \returns non-zero if the cells match, zero for no match or error. * */ -int compare_permuted_cell_parameters(UnitCell *cell, UnitCell *reference, - double *tolerance, IntegerMatrix **pmb) +int compare_lattices(UnitCell *cell_in, UnitCell *reference_in, + double *tolerance, RationalMatrix **pmb) { - IntegerMatrix *m; - signed int i[9]; - double a, b, c, al, be, ga; - double min_dist = +INFINITY; - - m = intmat_new(3, 3); - cell_get_parameters(reference, &a, &b, &c, &al, &be, &ga); - *pmb = NULL; - - 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 *test; - int j, k; - int l = 0; - signed int det; - - 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) && (det != -1) ) continue; - - test = cell_transform_intmat(cell, m); - - if ( compare_cell_parameters(reference, test, tolerance) ) { - - double at, bt, ct, alt, bet, gat; - double dist; + UnitCell *cell; + UnitCell *reference; + IntegerMatrix *CBint; + RationalMatrix *CiA; + RationalMatrix *CB; + IntegerMatrix *Mcell; + IntegerMatrix *Mref; + double eps; + double G6cell[6]; + double G6ref[6]; - cell_get_parameters(test, &at, &bt, &ct, &alt, &bet, &gat); - dist = g6_distance(at, bt, ct, alt, bet, gat, - a, b, c, al, be, ga); - if ( dist < min_dist ) { - min_dist = dist; - intmat_free(*pmb); - *pmb = intmat_copy(m); - } + /* 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); - } + cell = uncenter_cell(cell_in, NULL, &CiA); + if ( cell == NULL ) return 0; - cell_free(test); + /* Calculate G6 components for both */ + g6_components(G6cell, cell); + g6_components(G6ref, reference); - } - } - } - } - } - } - } - } - } + /* Convert both to reduced basis (stably) */ + eps = pow(cell_get_volume(cell), 1.0/3.0) * 1e-5; + Mcell = reduce_g6(G6cell, eps); + eps = pow(cell_get_volume(reference), 1.0/3.0) * 1e-5; + Mref = reduce_g6(G6ref, eps); - intmat_free(m); + /* Compare cells (including nearby ones, possibly permuting + * axes if close to equality) */ - if ( isinf(min_dist) ) { - return 0; - } - - /* Solution found */ - return 1; + intmat_free(Mcell); + intmat_free(Mref); + cell_free(reference); + cell_free(cell); } diff --git a/libcrystfel/src/cell-utils.h b/libcrystfel/src/cell-utils.h index fa577269..e46d860b 100644 --- a/libcrystfel/src/cell-utils.h +++ b/libcrystfel/src/cell-utils.h @@ -90,9 +90,6 @@ extern double cell_get_volume(UnitCell *cell); extern int compare_cell_parameters(UnitCell *cell, UnitCell *reference, double *tolerance); -extern int compare_permuted_cell_parameters(UnitCell *cell, UnitCell *reference, - double *tolerance, IntegerMatrix **pmb); - extern int compare_cell_parameters_and_orientation(UnitCell *cell, UnitCell *reference, const double ltl, @@ -108,6 +105,9 @@ extern int compare_reindexed_cell_parameters(UnitCell *cell, UnitCell *reference double *tolerance, int csl, RationalMatrix **pmb); +extern int compare_lattices(UnitCell *cell_in, UnitCell *reference_in, + double *tolerance, RationalMatrix **pmb); + #ifdef __cplusplus } #endif diff --git a/libcrystfel/src/index.c b/libcrystfel/src/index.c index b07dd067..aa0cfb6a 100644 --- a/libcrystfel/src/index.c +++ b/libcrystfel/src/index.c @@ -547,26 +547,13 @@ static int check_cell(IndexingFlags flags, Crystal *cr, UnitCell *target, double *tolerance) { UnitCell *out; - IntegerMatrix *im; RationalMatrix *rm; /* Check at all? */ if ( ! ((flags & INDEXING_CHECK_CELL_COMBINATIONS) || (flags & INDEXING_CHECK_CELL_AXES)) ) return 0; - if ( compare_permuted_cell_parameters(crystal_get_cell(cr), target, - tolerance, &im) ) - { - out = cell_transform_intmat(crystal_get_cell(cr), im); - cell_free(crystal_get_cell(cr)); - crystal_set_cell(cr, out); - intmat_free(im); - return 0; - } - - if ( (flags & INDEXING_CHECK_CELL_COMBINATIONS ) - && compare_reindexed_cell_parameters(crystal_get_cell(cr), target, - tolerance, 0, &rm) ) + if ( compare_lattices(crystal_get_cell(cr), target, tolerance, &rm) ) { out = cell_transform_rational(crystal_get_cell(cr), rm); cell_free(crystal_get_cell(cr)); @@ -576,7 +563,6 @@ static int check_cell(IndexingFlags flags, Crystal *cr, UnitCell *target, } return 1; - } -- cgit v1.2.3 From da13e5c05e0762860df003812991c43225d7d379 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Thu, 15 Aug 2019 16:35:50 +0200 Subject: Add cell_get_G6() and corresponding structure --- libcrystfel/src/cell.c | 15 +++++++++++++++ libcrystfel/src/cell.h | 12 ++++++++++++ 2 files changed, 27 insertions(+) (limited to 'libcrystfel/src') diff --git a/libcrystfel/src/cell.c b/libcrystfel/src/cell.c index 6b1d58c7..e785c1c8 100644 --- a/libcrystfel/src/cell.c +++ b/libcrystfel/src/cell.c @@ -579,6 +579,21 @@ LatticeType cell_get_lattice_type(UnitCell *cell) } +struct g6 cell_get_G6(UnitCell *cell) +{ + double a, b, c, al, be, ga; + struct g6 g; + cell_get_parameters(cell, &a, &b, &c, &al, &be, &ga); + g.A = a*a; + g.B = b*b; + g.C = c*c; + g.D = 2.0*b*c*cos(al); + g.E = 2.0*a*c*cos(be); + g.F = 2.0*a*b*cos(ga); + return g; +} + + char cell_get_unique_axis(UnitCell *cell) { return cell->unique_axis; diff --git a/libcrystfel/src/cell.h b/libcrystfel/src/cell.h index 4bbc1be2..e12fc209 100644 --- a/libcrystfel/src/cell.h +++ b/libcrystfel/src/cell.h @@ -132,6 +132,18 @@ extern void cell_set_reciprocal(UnitCell *cell, extern LatticeType cell_get_lattice_type(UnitCell *cell); extern void cell_set_lattice_type(UnitCell *cell, LatticeType lattice_type); +struct g6 +{ + double A; + double B; + double C; + double D; + double E; + double F; +}; + +extern struct g6 cell_get_G6(UnitCell *cell); + extern char cell_get_centering(UnitCell *cell); extern void cell_set_centering(UnitCell *cell, char centering); -- cgit v1.2.3 From 0b8430c5401803690c8ca659b533d0b1d3b022e0 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Thu, 15 Aug 2019 16:39:10 +0200 Subject: Working Niggli reduction --- libcrystfel/src/cell-utils.c | 367 ++++++++++++++++++++++++++----------------- 1 file changed, 224 insertions(+), 143 deletions(-) (limited to 'libcrystfel/src') diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c index 20ede189..1d351910 100644 --- a/libcrystfel/src/cell-utils.c +++ b/libcrystfel/src/cell-utils.c @@ -1931,29 +1931,18 @@ static Rational *find_candidates(double len, double *a, double *b, double *c, } -void g6_components(double *g6, UnitCell *cell) -{ - double a, b, c, al, be, ga; - cell_get_parameters(cell, &a, &b, &c, &al, &be, &ga); - g6[0] = a*a; - g6[1] = b*b; - g6[2] = c*c; - g6[3] = 2.0*b*c*cos(al); - g6[4] = 2.0*a*c*cos(be); - g6[5] = 2.0*a*b*cos(ga); -} - - static double g6_distance(UnitCell *cell1, UnitCell *cell2) { - double g1[6], g2[6]; - int i; + struct g6 g1, g2; double total = 0.0; - g6_components(g1, cell1); - g6_components(g2, cell2); - for ( i=0; i<6; i++ ) { - total += (g1[i]-g2[i])*(g1[i]-g2[i]); - } + g1 = cell_get_G6(cell1); + g2 = cell_get_G6(cell2); + total += (g1.A-g2.A)*(g1.A-g2.A); + total += (g1.B-g2.B)*(g1.B-g2.B); + total += (g1.C-g2.C)*(g1.C-g2.C); + total += (g1.D-g2.D)*(g1.D-g2.D); + total += (g1.E-g2.E)*(g1.E-g2.E); + total += (g1.F-g2.F)*(g1.F-g2.F); return sqrt(total); } @@ -2133,204 +2122,296 @@ int compare_reindexed_cell_parameters(UnitCell *cell_in, UnitCell *reference_in, /* Criteria from Grosse-Kunstleve et al., Acta Cryst A60 (2004) p1-6 */ -#define GT(x,y) (y < x-eps) -#define LT(x,y) (x < y-eps) +#define GT(x,y) (y < x - eps) +#define LT(x,y) GT(y,x) #define EQ(x,y) !(GT(x,y) || LT(x,y)) +#define LTE(x,y) !(GT(x,y)) -static void do_R5(IntegerMatrix *R, double *g, double eps) +static int in_standard_presentation(struct g6 g, double eps) { - signed int s = GT(g[3], 0.0) ? 1 : -1; - intmat_set(R, 0, 0, 1); - intmat_set(R, 1, 1, 1); - intmat_set(R, 2, 2, 1); - intmat_set(R, 1, 2, -s); - g[2] = g[1] + g[2] -s*g[3]; - g[3] = -2*s*g[1] + g[3]; - g[4] = g[4] - s*g[5]; - STATUS("R5\n"); + if ( GT(g.A, g.B) ) return 0; + if ( GT(g.B, g.C) ) return 0; + + if ( EQ(g.A, g.B) && GT(fabs(g.D), fabs(g.E)) ) return 0; + if ( EQ(g.B, g.C) && GT(fabs(g.E), fabs(g.F)) ) return 0; + + if ( ( GT(g.D, 0.0) && GT(g.E, 0.0) && GT(g.F, 0.0) ) + || ( !GT(g.D, 0.0) && !GT(g.E, 0.0) && !GT(g.F, 0.0) ) ) return 1; + + return 0; } -static void do_R6(IntegerMatrix *R, double *g, double eps) +static int is_burger(struct g6 g, double eps) { - signed int s = GT(g[4], 0.0) ? 1 : -1; - intmat_set(R, 0, 0, 1); - intmat_set(R, 1, 1, 1); - intmat_set(R, 2, 2, 1); - intmat_set(R, 0, 2, -s); - g[2] = g[0] + g[2] -s*g[4]; - g[3] = g[3] - s*g[5]; - g[4] = -2*s*g[0] + g[4]; - STATUS("R6\n"); + if ( !in_standard_presentation(g, eps) ) return 0; + if ( GT(fabs(g.D), g.B) ) return 0; + if ( GT(fabs(g.E), g.A) ) return 0; + if ( GT(fabs(g.F), g.A) ) return 0; + if ( LT(g.D + g.E + g.F + g.A + g.B, 0.0) ) return 0; + return 1; } -static void do_R7(IntegerMatrix *R, double *g, double eps) +static int is_niggli(struct g6 g, double eps) { - signed int s = GT(g[5], 0.0) ? 1 : -1; - intmat_set(R, 0, 0, 1); - intmat_set(R, 1, 1, 1); - intmat_set(R, 2, 2, 1); - intmat_set(R, 0, 1, -s); - g[1] = g[0] + g[1] -s*g[5]; - g[3] = g[3] - s*g[4]; - g[5] = -2*s*g[0] + g[5]; - STATUS("R7\n"); + if ( !is_burger(g, eps) ) return 0; + + if ( EQ(g.D, g.B) && GT(g.F, 2.0*g.E) ) return 0; + if ( EQ(g.E, g.A) && GT(g.F, 2.0*g.D) ) return 0; + if ( EQ(g.F, g.A) && GT(g.E, 2.0*g.D) ) return 0; + + if ( EQ(g.D, -g.B) && !EQ(g.F, 0.0) ) return 0; + if ( EQ(g.E, -g.A) && !EQ(g.F, 0.0) ) return 0; + if ( EQ(g.F, -g.A) && !EQ(g.E, 0.0) ) return 0; + + if ( EQ(g.D + g.E + g.F + g.A + g.B, 0.0) + && GT(2.0*(g.A + g.E)+g.F, 0.0) ) return 0; + + return 1; } -static void do_R8(IntegerMatrix *R, double *g, double eps) +static void debug_lattice(struct g6 g, double eps) { - intmat_set(R, 0, 0, 1); - intmat_set(R, 1, 1, 1); - intmat_set(R, 2, 2, 1); - intmat_set(R, 1, 2, 1); - g[2] = g[0]+g[1]+g[2]+g[3]+g[4]+g[5]; - g[3] = 2.0*g[1] + g[3] + g[5]; - g[4] = 2.0*g[0] + g[4] + g[5]; - STATUS("R8\n"); + STATUS("%e %e %e %e %e %e", + 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"); + STATUS("\n"); } -IntegerMatrix *reduce_g6(double *g, double eps) +static signed int eps_sign(double v, double eps) +{ + return GT(v, 0.0) ? +1 : -1; +} + + +static void mult_in_place(IntegerMatrix *T, IntegerMatrix *M) +{ + int i, j; + IntegerMatrix *tmp = intmat_intmat_mult(T, M); + for ( i=0; i<3; i++ ) { + for ( j=0; j<3; j++ ) { + intmat_set(T, i, j, intmat_get(tmp, i, j)); + } + } + intmat_free(tmp); +} + + +/* NB The G6 components are passed by value, not reference. + * It's the caller's reponsibility to apply the matrix to the cell and + * re-calculate the G6 vector, if required. */ +IntegerMatrix *reduce_g6(struct g6 g, double eps) { IntegerMatrix *T; - int finished = 0; + IntegerMatrix *M; + int finished; - T = intmat_identity(3); STATUS("eps = %e\n", eps); + T = intmat_identity(3); + M = intmat_new(3, 3); do { - IntegerMatrix *tmp; - IntegerMatrix *M; - IntegerMatrix *R; - - M = intmat_new(3, 3); - R = intmat_new(3, 3); + int done_s1s2; - int did_something; do { + done_s1s2 = 1; - did_something = 0; - - if ( GT(g[0], g[1]) || (EQ(g[0], g[1]) && GT(g[3], g[4])) ) { + debug_lattice(g, eps); + if ( GT(g.A, g.B) + || (EQ(g.A, g.B) && GT(fabs(g.D), fabs(g.E))) ) + { /* Swap a and b */ double temp; - temp = g[0]; g[0] = g[1]; g[1] = temp; - temp = g[3]; g[3] = g[4]; g[4] = temp; + temp = g.A; g.A = g.B; g.B = temp; + temp = g.D; g.D = g.E; g.E = temp; + + intmat_zero(M); intmat_set(M, 1, 0, -1); intmat_set(M, 0, 1, -1); intmat_set(M, 2, 2, -1); - STATUS("SP1\n"); - did_something = 1; + mult_in_place(T, M); - } else if ( GT(g[1], g[2]) - || (EQ(g[1], g[2]) && GT(g[4], g[5])) ) { + STATUS("step 1\n"); + } else if ( GT(g.B, g.C) + || (EQ(g.B, g.C) && GT(fabs(g.E), fabs(g.F))) ) + { /* Swap b and c */ double temp; - temp = g[1]; g[1] = g[2]; g[2] = temp; - temp = g[4]; g[4] = g[5]; g[5] = temp; + temp = g.B; g.B = g.C; g.C = temp; + temp = g.E; g.E = g.F; g.F = temp; + + intmat_zero(M); intmat_set(M, 0, 0, -1); intmat_set(M, 1, 2, -1); intmat_set(M, 2, 1, -1); - STATUS("SP2\n"); - did_something = 1; + mult_in_place(T, M); + STATUS("step 2\n"); + done_s1s2 = 0; } - } while ( did_something ); + } while ( !done_s1s2 ); - if ( GT(g[3]*g[4]*g[5], 0.0) ) { + finished = 0; - intmat_set(M, 0, 0, GT(g[3], 0.0) ? 1 : -1); - intmat_set(M, 1, 1, GT(g[4], 0.0) ? 1 : -1); - intmat_set(M, 2, 2, GT(g[5], 0.0) ? 1 : -1); - g[3] = fabs(g[3]); - g[4] = fabs(g[4]); - g[5] = fabs(g[5]); - STATUS("SP3\n"); + debug_lattice(g, eps); - } else { + /* 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); + mult_in_place(T, M); + + g.D = fabs(g.D); + g.E = fabs(g.E); + g.F = fabs(g.F); + + STATUS("step 3\n"); + + } 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); - if ( GT(g[3], 0.0) ) { + + if ( GT(g.D, 0.0) ) { intmat_set(M, 0, 0, -1); - } else if ( !LT(g[3], 0.0) ) { + } else if ( !LT(g.D, 0.0) ) { r = 0; c = 0; have_rc = 1; } - if ( GT(g[4], 0.0) ) { + if ( GT(g.E, 0.0) ) { intmat_set(M, 1, 1, -1); - } else if ( !LT(g[3], 0.0) ) { + } else if ( !LT(g.D, 0.0) ) { r = 1; c = 1; have_rc = 1; - } - if ( GT(g[5], 0.0) ) { + } + if ( GT(g.F, 0.0) ) { intmat_set(M, 2, 2, -1); - } else if ( !LT(g[3], 0.0) ) { + } else if ( !LT(g.D, 0.0) ) { r = 2; c = 2; have_rc = 1; } - if ( LT(g[3]*g[4]*g[5], 0.0) ) { + 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(); + //abort(); } } - g[3] = -fabs(g[3]); - g[4] = -fabs(g[4]); - g[5] = -fabs(g[5]); - STATUS("SP4\n"); + mult_in_place(T, M); + + g.D = -fabs(g.D); + g.E = -fabs(g.E); + g.F = -fabs(g.F); + + STATUS("step 4\n"); } - STATUS("G6: %f %f %f %f %f %f\n", g[0], g[1], g[2], g[3], g[4], g[5]); - - if ( GT(fabs(g[3]), g[1]) ) { - do_R5(R, g, eps); - } else if ( GT(fabs(g[5]), g[0]) ) { - do_R6(R, g, eps); - } else if ( GT(fabs(g[6]), g[0]) ) { - do_R7(R, g, eps); - } else if ( LT(g[0]+g[1]+g[2]+g[3]+g[4]+g[5], 0.0) ) { - do_R8(R, g, eps); - } else if ( (EQ(g[1], g[3]) && LT(2.0*g[4], g[4])) - || (EQ(g[1], -g[3]) && LT(g[5], 0.0) ) ) { - do_R5(R, g, eps); - } else if ( (EQ(g[0], g[4]) && LT(2.0*g[3], g[5])) - || (EQ(g[0], -g[4]) && LT(g[5], 0.0) ) ) { - do_R6(R, g, eps); - } else if ( (EQ(g[0], g[5]) && LT(2.0*g[3], g[4])) - || (EQ(g[0], -g[5]) && LT(g[6], 0.0) ) ) { - do_R7(R, g, eps); - } else if ( EQ(g[0]+g[1]+g[2]+g[3]+g[4]+g[5], 0.0) - && GT(2.0*g[0] + 2.0*g[4] + g[5], 0.0) ) { - do_R8(R, g, eps); + 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)) ) + { + signed int s = g.D > 0.0 ? 1 : -1; + + intmat_zero(M); + intmat_set(M, 0, 0, 1); + intmat_set(M, 1, 1, 1); + intmat_set(M, 2, 2, 1); + intmat_set(M, 1, 2, -s); + mult_in_place(T, M); + + g.C = g.B + g.C -s*g.D; + g.D = -2*s*g.B + g.D; + g.E = g.E - s*g.F; + + STATUS("R5\n"); + + } else if ( GT(fabs(g.E), g.A) + || (EQ(g.A, g.E) && LT(2.0*g.D, g.F)) + || (EQ(g.A, -g.E) && LT(g.F, 0.0)) ) + { + signed int s = g.E > 0.0 ? 1 : -1; + + intmat_zero(M); + intmat_set(M, 0, 0, 1); + intmat_set(M, 1, 1, 1); + intmat_set(M, 2, 2, 1); + intmat_set(M, 0, 2, -s); + mult_in_place(T, M); + + g.C = g.A + g.C -s*g.E; + g.D = g.D - s*g.F; + g.E = -2*s*g.A + g.E; + + STATUS("R6\n"); + + } else if ( GT(fabs(g.F), g.A) + || (EQ(g.A, g.F) && LT(2.0*g.D, g.E)) + || (EQ(g.A, -g.F) && LT(g.E, 0.0)) ) + { + signed int s = g.F > 0.0 ? 1 : -1; + + intmat_zero(M); + intmat_set(M, 0, 0, 1); + intmat_set(M, 1, 1, 1); + intmat_set(M, 2, 2, 1); + intmat_set(M, 0, 1, -s); + mult_in_place(T, M); + + g.B = g.A + g.B -s*g.F; + g.D = g.D - s*g.E; + g.F = -2*s*g.A + g.F; + + STATUS("R7\n"); + + } 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) + && GT(2.0*g.A + 2.0*g.E + g.F, 0.0)) ) ) + { + intmat_zero(M); + intmat_set(M, 0, 0, 1); + intmat_set(M, 1, 1, 1); + intmat_set(M, 2, 2, 1); + intmat_set(M, 1, 2, 1); + mult_in_place(T, M); + + g.C = g.A+g.B+g.C+g.D+g.E+g.F; + g.D = 2.0*g.B + g.D + g.F; + g.E = 2.0*g.A + g.E + g.F; + + STATUS("R8\n"); + } else { finished = 1; - R = intmat_identity(3); - STATUS("Done.\n"); } - STATUS("G6: %f %f %f %f %f %f\n", g[0], g[1], g[2], g[3], g[4], g[5]); - tmp = intmat_intmat_mult(T, M); - intmat_free(T); - T = intmat_intmat_mult(tmp, R); - intmat_free(tmp); - } while ( !finished ); + debug_lattice(g, eps); + + intmat_free(M); return T; } @@ -2366,8 +2447,8 @@ int compare_lattices(UnitCell *cell_in, UnitCell *reference_in, IntegerMatrix *Mcell; IntegerMatrix *Mref; double eps; - double G6cell[6]; - double G6ref[6]; + struct g6 g6cell; + struct g6 g6ref; /* Un-center both cells */ reference = uncenter_cell(reference_in, &CBint, NULL); @@ -2379,14 +2460,14 @@ int compare_lattices(UnitCell *cell_in, UnitCell *reference_in, if ( cell == NULL ) return 0; /* Calculate G6 components for both */ - g6_components(G6cell, cell); - g6_components(G6ref, reference); + g6cell = cell_get_G6(cell); + g6ref = cell_get_G6(reference); /* Convert both to reduced basis (stably) */ eps = pow(cell_get_volume(cell), 1.0/3.0) * 1e-5; - Mcell = reduce_g6(G6cell, eps); + Mcell = reduce_g6(g6cell, eps); eps = pow(cell_get_volume(reference), 1.0/3.0) * 1e-5; - Mref = reduce_g6(G6ref, eps); + Mref = reduce_g6(g6ref, eps); /* Compare cells (including nearby ones, possibly permuting * axes if close to equality) */ -- cgit v1.2.3 From 29cce73cf8b97dc7e58ebc94e5ddc42a70c98233 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Fri, 16 Aug 2019 10:21:30 +0200 Subject: Determine cell reduction tolerance automatically --- libcrystfel/src/cell-utils.c | 20 +++++++++++++++----- 1 file changed, 15 insertions(+), 5 deletions(-) (limited to 'libcrystfel/src') diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c index 1d351910..5937a9c6 100644 --- a/libcrystfel/src/cell-utils.c +++ b/libcrystfel/src/cell-utils.c @@ -2203,16 +2203,28 @@ static void mult_in_place(IntegerMatrix *T, IntegerMatrix *M) } +/* Cell volume from G6 components */ +static double g6_volume(struct g6 g) +{ + return sqrt(g.A*g.B*g.C + - 0.25*(g.A*g.D*g.D + g.B*g.E*g.E + g.C*g.F*g.F - g.D*g.E*g.F)); +} + + /* NB The G6 components are passed by value, not reference. * It's the caller's reponsibility to apply the matrix to the cell and * re-calculate the G6 vector, if required. */ -IntegerMatrix *reduce_g6(struct g6 g, double eps) +IntegerMatrix *reduce_g6(struct g6 g, double epsrel) { IntegerMatrix *T; IntegerMatrix *M; int finished; + double eps; + 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); @@ -2464,10 +2476,8 @@ int compare_lattices(UnitCell *cell_in, UnitCell *reference_in, g6ref = cell_get_G6(reference); /* Convert both to reduced basis (stably) */ - eps = pow(cell_get_volume(cell), 1.0/3.0) * 1e-5; - Mcell = reduce_g6(g6cell, eps); - eps = pow(cell_get_volume(reference), 1.0/3.0) * 1e-5; - Mref = reduce_g6(g6ref, eps); + Mcell = reduce_g6(g6cell, 1e-5); + Mref = reduce_g6(g6ref, 1e-5); /* Compare cells (including nearby ones, possibly permuting * axes if close to equality) */ -- cgit v1.2.3 From 5ab96056a6d1972a248abb154ea5df0ffda44d06 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Fri, 16 Aug 2019 10:22:59 +0200 Subject: Clean up / debug --- libcrystfel/src/cell-utils.c | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) (limited to 'libcrystfel/src') diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c index 5937a9c6..097ee463 100644 --- a/libcrystfel/src/cell-utils.c +++ b/libcrystfel/src/cell-utils.c @@ -2358,7 +2358,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("R5\n"); + STATUS("step 5\n"); } else if ( GT(fabs(g.E), g.A) || (EQ(g.A, g.E) && LT(2.0*g.D, g.F)) @@ -2377,7 +2377,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("R6\n"); + STATUS("step 6\n"); } else if ( GT(fabs(g.F), g.A) || (EQ(g.A, g.F) && LT(2.0*g.D, g.E)) @@ -2396,7 +2396,7 @@ IntegerMatrix *reduce_g6(struct g6 g, double epsrel) g.D = g.D - s*g.E; g.F = -2*s*g.A + g.F; - STATUS("R7\n"); + STATUS("step 7\n"); } 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) @@ -2413,7 +2413,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("R8\n"); + STATUS("step 8\n"); } else { finished = 1; @@ -2423,6 +2423,9 @@ IntegerMatrix *reduce_g6(struct g6 g, double epsrel) debug_lattice(g, eps); + assert(is_burger(g, eps)); + assert(is_niggli(g, eps)); + intmat_free(M); return T; } @@ -2458,7 +2461,6 @@ int compare_lattices(UnitCell *cell_in, UnitCell *reference_in, RationalMatrix *CB; IntegerMatrix *Mcell; IntegerMatrix *Mref; - double eps; struct g6 g6cell; struct g6 g6ref; @@ -2481,6 +2483,7 @@ int compare_lattices(UnitCell *cell_in, UnitCell *reference_in, /* Compare cells (including nearby ones, possibly permuting * axes if close to equality) */ + STATUS("Reduced cell:\n"); intmat_free(Mcell); intmat_free(Mref); -- cgit v1.2.3 From 0b1db945a12c0e68491412baf322b761511eb016 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Wed, 21 Aug 2019 12:23:39 +0200 Subject: Tidy up comparison function definitions Especially, remove the last ltl/atl tolerance values. --- libcrystfel/src/cell-utils.c | 120 +++++++++++++++++++++++++------------------ libcrystfel/src/cell-utils.h | 20 ++++---- libcrystfel/src/index.c | 9 +++- 3 files changed, 87 insertions(+), 62 deletions(-) (limited to 'libcrystfel/src') diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c index 097ee463..961a6554 100644 --- a/libcrystfel/src/cell-utils.c +++ b/libcrystfel/src/cell-utils.c @@ -1671,7 +1671,7 @@ double cell_get_volume(UnitCell *cell) /** * \param cell: A UnitCell * \param reference: Another UnitCell - * \param tolerance: Pointer to tolerances for a,b,c (fractional), al,be,ga (radians) + * \param tols: Pointer to tolerances for a,b,c (fractional), al,be,ga (radians) * * Compare the two unit cells. If the real space parameters match to within * the specified tolerances, and the centering matches, this function returns 1. @@ -1684,7 +1684,8 @@ double cell_get_volume(UnitCell *cell) * \returns non-zero if the cells match. * */ -int compare_cell_parameters(UnitCell *cell, UnitCell *reference, double *tolerance) +int compare_cell_parameters(UnitCell *cell, UnitCell *reference, + const double *tols) { double a1, b1, c1, al1, be1, ga1; double a2, b2, c2, al2, be2, ga2; @@ -1696,12 +1697,12 @@ int compare_cell_parameters(UnitCell *cell, UnitCell *reference, double *toleran cell_get_parameters(cell, &a1, &b1, &c1, &al1, &be1, &ga1); cell_get_parameters(reference, &a2, &b2, &c2, &al2, &be2, &ga2); - if ( !within_tolerance(a1, a2, tolerance[0]*100.0) ) return 0; - if ( !within_tolerance(b1, b2, tolerance[1]*100.0) ) return 0; - if ( !within_tolerance(c1, c2, tolerance[2]*100.0) ) return 0; - if ( fabs(al1-al2) > tolerance[3] ) return 0; - if ( fabs(be1-be2) > tolerance[4] ) return 0; - if ( fabs(ga1-ga2) > tolerance[5] ) return 0; + if ( !within_tolerance(a1, a2, tols[0]*100.0) ) return 0; + if ( !within_tolerance(b1, b2, tols[1]*100.0) ) return 0; + if ( !within_tolerance(c1, c2, tols[2]*100.0) ) return 0; + if ( fabs(al1-al2) > tols[3] ) return 0; + if ( fabs(be1-be2) > tols[4] ) return 0; + if ( fabs(ga1-ga2) > tols[5] ) return 0; return 1; } @@ -1719,17 +1720,22 @@ static double moduli_check(double ax, double ay, double az, /** * \param cell: A UnitCell * \param reference: Another UnitCell - * \param ltl: Maximum allowable fractional difference in axis lengths - * \param atl: Maximum allowable difference in angles (in radians) + * \param tols: Pointer to six tolerance values (see below) * * Compare the two unit cells. If the axes match in length (to within - * fractional difference \p ltl) and the axes are aligned to within \p atl radians, - * this function returns non-zero. + * the specified tolerances), this function returns non-zero. * * This function compares the orientation of the cell as well as the parameters. * If you just want to see if the parameters are the same, use * compare_cell_parameters() instead. * + * The comparison is done by checking that the lengths of the unit cell axes are + * the same between the two cells, and that the axes have similar directions in + * 3D space. The first three tolerance values are the maximum allowable fractional + * differences between the a,b,c axis lengths (respectively) of the two cells. + * The last three tolerance values are the maximum allowable angles, in radians, + * between the directions of the a,b,c axes of the two cells. + * * \p cell and \p reference must have the same centering. Otherwise, this * function always returns zero. * @@ -1737,7 +1743,7 @@ static double moduli_check(double ax, double ay, double az, * */ int compare_cell_parameters_and_orientation(UnitCell *cell, UnitCell *reference, - const double ltl, const double atl) + const double *tols) { double ax1, ay1, az1, bx1, by1, bz1, cx1, cy1, cz1; double ax2, ay2, az2, bx2, by2, bz2, cx2, cy2, cz2; @@ -1752,13 +1758,13 @@ int compare_cell_parameters_and_orientation(UnitCell *cell, UnitCell *reference, &bx2, &by2, &bz2, &cx2, &cy2, &cz2); - if ( angle_between(ax1, ay1, az1, ax2, ay2, az2) > atl ) return 0; - if ( angle_between(bx1, by1, bz1, bx2, by2, bz2) > atl ) return 0; - if ( angle_between(cx1, cy1, cz1, cx2, cy2, cz2) > atl ) return 0; + if ( angle_between(ax1, ay1, az1, ax2, ay2, az2) > tols[3] ) return 0; + if ( angle_between(bx1, by1, bz1, bx2, by2, bz2) > tols[4] ) return 0; + if ( angle_between(cx1, cy1, cz1, cx2, cy2, cz2) > tols[5] ) return 0; - if ( moduli_check(ax1, ay1, az1, ax2, ay2, az2) > ltl ) return 0; - if ( moduli_check(bx1, by1, bz1, bx2, by2, bz2) > ltl ) return 0; - if ( moduli_check(cx1, cy1, cz1, cx2, cy2, cz2) > ltl ) return 0; + if ( moduli_check(ax1, ay1, az1, ax2, ay2, az2) > tols[0] ) return 0; + if ( moduli_check(bx1, by1, bz1, bx2, by2, bz2) > tols[1] ) return 0; + if ( moduli_check(cx1, cy1, cz1, cx2, cy2, cz2) > tols[2] ) return 0; return 1; } @@ -1767,18 +1773,23 @@ int compare_cell_parameters_and_orientation(UnitCell *cell, UnitCell *reference, /** * \param cell: A UnitCell * \param reference: Another UnitCell - * \param ltl: Maximum allowable fractional difference in reciprocal axis lengths - * \param atl: Maximum allowable difference in reciprocal angles (in radians) + * \param tols: Pointer to six tolerance values (see below) * \param pmb: Place to store pointer to matrix * * Compare the two unit cells. If, using any permutation of the axes, the - * axes can be made to match in length (to within fractional difference \p ltl) - * and the axes aligned to within \p atl radians, this function returns non-zero - * and stores the transformation which must be applied to \p cell at \p pmb. + * axes can be made to match in length and the axes aligned in space, this + * function returns non-zero and stores the transformation which must be applied + * to \p cell at \p pmb. + * + * A "permutation" means a transformation represented by a matrix with all + * elements equal to +1, 0 or -1, having determinant +1 or -1. That means that + * this function will find the relationship between a left-handed and a right- + * handed basis. * * Note that the orientations of the cells must match, not just the parameters. * The comparison is done after reindexing using - * compare_cell_parameters_and_orientation(). + * compare_cell_parameters_and_orientation(). See that function for more + * details. * * \p cell and \p reference must have the same centering. Otherwise, this * function always returns zero. @@ -1788,7 +1799,7 @@ int compare_cell_parameters_and_orientation(UnitCell *cell, UnitCell *reference, */ int compare_permuted_cell_parameters_and_orientation(UnitCell *cell, UnitCell *reference, - double ltl, double atl, + const double *tols, IntegerMatrix **pmb) { IntegerMatrix *m; @@ -1823,7 +1834,7 @@ int compare_permuted_cell_parameters_and_orientation(UnitCell *cell, nc = cell_transform_intmat(cell, m); if ( compare_cell_parameters_and_orientation(nc, reference, - ltl, atl) ) + tols) ) { if ( pmb != NULL ) *pmb = m; cell_free(nc); @@ -1950,13 +1961,13 @@ static double g6_distance(UnitCell *cell1, UnitCell *cell2) /** * \param cell_in: A UnitCell * \param reference_in: Another UnitCell - * \param tolerance: Pointer to tolerances for a,b,c (fractional), al,be,ga (radians) + * \param tols: Pointer to tolerances for a,b,c (fractional), al,be,ga (radians) * \param csl: Non-zero to look for coincidence site lattice relationships * \param pmb: Place to store pointer to matrix * * Compare the \p cell_in with \p reference_in. If \p cell is a derivative lattice - * of \p reference, within fractional axis length differences \p tolerance[0..2] - * and absolute angle difference \p atl[3..5] (in radians), this function returns + * of \p reference, within fractional axis length differences \p tols[0..2] + * and absolute angle difference \p tols[3..5] (in radians), this function returns * non-zero and stores the transformation which needs to be applied to * \p cell_in at \p pmb. * @@ -1975,12 +1986,16 @@ static double g6_distance(UnitCell *cell1, UnitCell *cell2) * meaning that the lattice points of the transformed version of \p cell_in * might not coincide with lattice points of \p reference_in. * + * This function is used by CrystFEL's cell_tool program to find non-obvious + * relationships between crystal lattices. For most routine comparisons, this + * function is probably not the one you need! + * * \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 *tolerance, int csl, - RationalMatrix **pmb) +int compare_derivative_cell_parameters(UnitCell *cell_in, UnitCell *reference_in, + const double *tols, int csl, + RationalMatrix **pmb) { UnitCell *cell; UnitCell *reference; @@ -2016,9 +2031,9 @@ int compare_reindexed_cell_parameters(UnitCell *cell_in, UnitCell *reference_in, &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, tolerance[0], csl, &ncand_a); - cand_b = find_candidates(b, av, bv, cv, tolerance[1], csl, &ncand_b); - cand_c = find_candidates(c, av, bv, cv, tolerance[2], csl, &ncand_c); + cand_a = find_candidates(a, av, bv, cv, tols[0], csl, &ncand_a); + cand_b = find_candidates(b, av, bv, cv, tols[1], csl, &ncand_b); + cand_c = find_candidates(c, av, bv, cv, tols[2], csl, &ncand_c); if ( (ncand_a==0) || (ncand_b==0) || (ncand_c==0) ) { *pmb = NULL; @@ -2055,7 +2070,7 @@ int compare_reindexed_cell_parameters(UnitCell *cell_in, UnitCell *reference_in, test = cell_transform_rational(cell, M); cell_get_parameters(test, &at, &bt, &ct, &alt, &bet, &gat); cell_free(test); - if ( fabs(gat - ga) > tolerance[5] ) continue; + if ( fabs(gat - ga) > tols[5] ) continue; /* Gamma OK, now look for place for c axis */ for ( ic=0; ic tolerance[3] ) { + if ( fabs(alt - al) > tols[3] ) { cell_free(test); continue; } - if ( fabs(bet - be) > tolerance[4] ) { + if ( fabs(bet - be) > tols[4] ) { cell_free(test); continue; } @@ -2434,25 +2449,30 @@ IntegerMatrix *reduce_g6(struct g6 g, double epsrel) /** * \param cell_in: A UnitCell * \param reference_in: Another UnitCell - * \param tolerance: Pointer to tolerances for a,b,c (fractional), al,be,ga (radians) - * \param pmb: Place to store pointer to matrix + * \param tols: Pointer to tolerances for a,b,c (fractional), al,be,ga (radians) + * \param pmb: Place to store pointer to matrix, or NULL if not needed * * Compare the \p cell_in with \p reference_in. If they represent the same - * lattice, this function returns non-zero and stores the transformation which - * needs to be applied to \p cell_in at \p pmb. + * lattice, this function returns a copy of \p cell_in transformed to look + * similar to \p reference_in. Otherwise, it returns NULL. * - * Subject to the tolerances, this function will find the transformation which - * gives the best match to the reference cell, using the Euclidian norm in - * G6 [see e.g. Andrews and Bernstein, Acta Cryst. A44 (1988) p1009]. + * If \pmb is non-NULL, the transformation which needs to be applied to + * \p cell_in will be stored there. * * Only the cell parameters will be compared. The relative orientations are - * irrelevant. + * irrelevant. The tolerances will be applied to the transformed copy of + * \p cell_in. * - * \returns non-zero if the cells match, zero for no match or error. + * This is the right function to use for deciding if an indexing solution + * matches a reference cell or not. + * + * \returns A newly allocated UnitCell, or NULL. * */ -int compare_lattices(UnitCell *cell_in, UnitCell *reference_in, - double *tolerance, RationalMatrix **pmb) +UnitCell *compare_reindexed_cell_parameters(UnitCell *cell_in, + UnitCell *reference_in, + const double *tols, + RationalMatrix **pmb) { UnitCell *cell; UnitCell *reference; diff --git a/libcrystfel/src/cell-utils.h b/libcrystfel/src/cell-utils.h index e46d860b..bb1e5cd0 100644 --- a/libcrystfel/src/cell-utils.h +++ b/libcrystfel/src/cell-utils.h @@ -88,25 +88,25 @@ extern int forbidden_reflection(UnitCell *cell, extern double cell_get_volume(UnitCell *cell); extern int compare_cell_parameters(UnitCell *cell, UnitCell *reference, - double *tolerance); + const double *tols); extern int compare_cell_parameters_and_orientation(UnitCell *cell, UnitCell *reference, - const double ltl, - const double atl); + const double *tols); extern int compare_permuted_cell_parameters_and_orientation(UnitCell *cell, UnitCell *reference, - double ltl, - double atl, + const double *tols, IntegerMatrix **pmb); -extern int compare_reindexed_cell_parameters(UnitCell *cell, UnitCell *reference, - double *tolerance, int csl, - RationalMatrix **pmb); +extern int compare_derivative_cell_parameters(UnitCell *cell, UnitCell *reference, + const double *tols, int csl, + RationalMatrix **pmb); -extern int compare_lattices(UnitCell *cell_in, UnitCell *reference_in, - double *tolerance, RationalMatrix **pmb); +extern UnitCell *compare_reindexed_cell_parameters(UnitCell *cell_in, + UnitCell *reference_in, + const double *tols, + RationalMatrix **pmb); #ifdef __cplusplus } diff --git a/libcrystfel/src/index.c b/libcrystfel/src/index.c index aa0cfb6a..45b24afb 100644 --- a/libcrystfel/src/index.c +++ b/libcrystfel/src/index.c @@ -553,7 +553,8 @@ static int check_cell(IndexingFlags flags, Crystal *cr, UnitCell *target, if ( ! ((flags & INDEXING_CHECK_CELL_COMBINATIONS) || (flags & INDEXING_CHECK_CELL_AXES)) ) return 0; - if ( compare_lattices(crystal_get_cell(cr), target, tolerance, &rm) ) + if ( compare_reindexed_cell_parameters(crystal_get_cell(cr), target, + tolerance, &rm) ) { out = cell_transform_rational(crystal_get_cell(cr), rm); cell_free(crystal_get_cell(cr)); @@ -691,13 +692,17 @@ static int try_indexer(struct image *image, IndexingMethod indm, for ( j=0; jcrystals[j]; + const double tols[] = {0.1, 0.1, 0.1, + deg2rad(5.0), + deg2rad(5.0), + deg2rad(5.0)}; /* Don't do similarity check against bad crystals */ if ( crystal_get_user_flag(that_cr) ) continue; if ( compare_permuted_cell_parameters_and_orientation(crystal_get_cell(cr), crystal_get_cell(that_cr), - 0.1, deg2rad(0.5), NULL) ) + tols, NULL) ) { crystal_set_user_flag(cr, 1); } -- cgit v1.2.3 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/src') 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 From 877bea233ed0a9150fdb900ba7fbc480aacc313c Mon Sep 17 00:00:00 2001 From: Thomas White Date: Wed, 21 Aug 2019 17:24:32 +0200 Subject: Fix matrix multiplication --- libcrystfel/src/cell-utils.c | 72 ++++++++++++++++++++++---------------------- 1 file changed, 36 insertions(+), 36 deletions(-) (limited to 'libcrystfel/src') diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c index 72225a8c..e59f1fa1 100644 --- a/libcrystfel/src/cell-utils.c +++ b/libcrystfel/src/cell-utils.c @@ -2482,18 +2482,21 @@ static IntegerMatrix *check_permutations(UnitCell *cell, UnitCell *reference, double a, b, c, al, be, ga; double min_dist = +INFINITY; IntegerMatrix *best_m = NULL; - RationalMatrix *RBCB; - RationalMatrix *rRB; + RationalMatrix *RiBCB; + RationalMatrix *rRiB; RationalMatrix *rCB; + IntegerMatrix *RiB; 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); + RiBCB = rtnl_mtx_new(3, 3); + RiB = intmat_inverse(RB); + rRiB = rtnl_mtx_from_intmat(RiB); rCB = rtnl_mtx_from_intmat(CB); - rtnl_mtx_mtxmult(rRB, rCB, RBCB); - rtnl_mtx_free(rRB); + rtnl_mtx_mtxmult(rRiB, rCB, RiBCB); + intmat_free(RiB); + rtnl_mtx_free(rRiB); rtnl_mtx_free(rCB); for ( i[0]=-1; i[0]<=+1; i[0]++ ) { @@ -2520,7 +2523,7 @@ static IntegerMatrix *check_permutations(UnitCell *cell, UnitCell *reference, if ( det != +1 ) continue; tmp = cell_transform_intmat(cell, m); - nc = cell_transform_rational(tmp, RBCB); + nc = cell_transform_rational(tmp, RiBCB); cell_free(tmp); if ( compare_cell_parameters(nc, reference, tols) ) { @@ -2544,7 +2547,7 @@ static IntegerMatrix *check_permutations(UnitCell *cell, UnitCell *reference, } } intmat_free(m); - rtnl_mtx_free(RBCB); + rtnl_mtx_free(RiBCB); return best_m; } @@ -2631,37 +2634,34 @@ UnitCell *compare_reindexed_cell_parameters(UnitCell *cell_in, P = check_permutations(cell_reduced, reference_in, RB, CB, tols); if ( P != NULL ) { - match = cell_transform_intmat(cell_reduced, P); - + //STATUS("The best permutation matrix:\n"); + //intmat_print(P); + + /* 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); + + match = cell_transform_rational(cell_in, comb); //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; - - } + if ( pmb != NULL ) *pmb = comb; } else { match = NULL; -- cgit v1.2.3 From 95cce19a7aeb6be7568ab73a5263913fe476591e Mon Sep 17 00:00:00 2001 From: Thomas White Date: Mon, 5 Aug 2019 14:25:04 +0200 Subject: Remove match_cell and match_cell_ab --- libcrystfel/src/cell-utils.c | 398 ------------------------------------------- libcrystfel/src/cell-utils.h | 5 - 2 files changed, 403 deletions(-) (limited to 'libcrystfel/src') diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c index e59f1fa1..2af3f293 100644 --- a/libcrystfel/src/cell-utils.c +++ b/libcrystfel/src/cell-utils.c @@ -574,404 +574,6 @@ UnitCell *uncenter_cell(UnitCell *in, IntegerMatrix **pC, RationalMatrix **pCi) } -#define MAX_CAND (1024) - -static int right_handed_vec(struct rvec a, struct rvec b, struct rvec c) -{ - struct rvec aCb; - double aCb_dot_c; - - /* "a" cross "b" */ - aCb.u = a.v*b.w - a.w*b.v; - aCb.v = - (a.u*b.w - a.w*b.u); - aCb.w = a.u*b.v - a.v*b.u; - - /* "a cross b" dot "c" */ - aCb_dot_c = aCb.u*c.u + aCb.v*c.v + aCb.w*c.w; - - if ( aCb_dot_c > 0.0 ) return 1; - return 0; -} - - -struct cvec { - struct rvec vec; - float na; - float nb; - float nc; - float fom; -}; - - -static int same_vector(struct cvec a, struct cvec b) -{ - if ( a.na != b.na ) return 0; - if ( a.nb != b.nb ) return 0; - if ( a.nc != b.nc ) return 0; - return 1; -} - - -/* Attempt to make 'cell' fit into 'template' somehow */ -UnitCell *match_cell(UnitCell *cell_in, UnitCell *template_in, int verbose, - const float *tols, int reduce) -{ - signed int n1l, n2l, n3l; - double asx, asy, asz; - double bsx, bsy, bsz; - double csx, csy, csz; - int i, j; - double lengths[3]; - double angles[3]; - struct cvec *cand[3]; - UnitCell *new_cell = NULL; - float best_fom = +999999999.9; /* Large number.. */ - int ncand[3] = {0,0,0}; - signed int ilow, ihigh; - float angtol = deg2rad(tols[3]); - UnitCell *cell; - UnitCell *template; - IntegerMatrix *centering; - UnitCell *new_cell_trans; - - /* "Un-center" the template unit cell to make the comparison easier */ - template = uncenter_cell(template_in, ¢ering, NULL); - if ( template == NULL ) return NULL; - - /* The candidate cell is also uncentered, because it might be centered - * if it came from (e.g.) MOSFLM */ - cell = uncenter_cell(cell_in, NULL, NULL); - if ( cell == NULL ) return NULL; - - if ( cell_get_reciprocal(template, &asx, &asy, &asz, - &bsx, &bsy, &bsz, - &csx, &csy, &csz) ) { - ERROR("Couldn't get reciprocal cell for template.\n"); - cell_free(template); - cell_free(cell); - intmat_free(centering); - return NULL; - } - - lengths[0] = modulus(asx, asy, asz); - lengths[1] = modulus(bsx, bsy, bsz); - lengths[2] = modulus(csx, csy, csz); - - angles[0] = angle_between(bsx, bsy, bsz, csx, csy, csz); - angles[1] = angle_between(asx, asy, asz, csx, csy, csz); - angles[2] = angle_between(asx, asy, asz, bsx, bsy, bsz); - - cand[0] = malloc(MAX_CAND*sizeof(struct cvec)); - cand[1] = malloc(MAX_CAND*sizeof(struct cvec)); - cand[2] = malloc(MAX_CAND*sizeof(struct cvec)); - - if ( cell_get_reciprocal(cell, &asx, &asy, &asz, - &bsx, &bsy, &bsz, - &csx, &csy, &csz) ) { - ERROR("Couldn't get reciprocal cell.\n"); - cell_free(template); - cell_free(cell); - intmat_free(centering); - return NULL; - } - - if ( reduce ) { - ilow = -2; ihigh = 4; - } else { - ilow = 0; ihigh = 1; - } - - /* Negative values mean 1/n, positive means n, zero means zero */ - for ( n1l=ilow; n1l<=ihigh; n1l++ ) { - for ( n2l=ilow; n2l<=ihigh; n2l++ ) { - for ( n3l=ilow; n3l<=ihigh; n3l++ ) { - - float n1, n2, n3; - signed int b1, b2, b3; - - n1 = (n1l>=0) ? (n1l) : (1.0/n1l); - n2 = (n2l>=0) ? (n2l) : (1.0/n2l); - n3 = (n3l>=0) ? (n3l) : (1.0/n3l); - - if ( !reduce ) { - if ( n1l + n2l + n3l > 1 ) continue; - } - - /* 'bit' values can be +1 or -1 */ - for ( b1=-1; b1<=1; b1+=2 ) { - for ( b2=-1; b2<=1; b2+=2 ) { - for ( b3=-1; b3<=1; b3+=2 ) { - - double tx, ty, tz; - double tlen; - int i; - - n1 *= b1; n2 *= b2; n3 *= b3; - - tx = n1*asx + n2*bsx + n3*csx; - ty = n1*asy + n2*bsy + n3*csy; - tz = n1*asz + n2*bsz + n3*csz; - tlen = modulus(tx, ty, tz); - - /* Test modulus for agreement with moduli of template */ - for ( i=0; i<3; i++ ) { - - if ( !within_tolerance(lengths[i], tlen, - tols[i]) ) - { - continue; - } - - if ( ncand[i] == MAX_CAND ) { - ERROR("Too many cell candidates - "); - ERROR("consider tightening the unit "); - ERROR("cell tolerances.\n"); - } else { - - double fom; - - fom = fabs(lengths[i] - tlen); - - cand[i][ncand[i]].vec.u = tx; - cand[i][ncand[i]].vec.v = ty; - cand[i][ncand[i]].vec.w = tz; - cand[i][ncand[i]].na = n1; - cand[i][ncand[i]].nb = n2; - cand[i][ncand[i]].nc = n3; - cand[i][ncand[i]].fom = fom; - - ncand[i]++; - - } - } - - } - } - } - - } - } - } - - if ( verbose ) { - STATUS("Candidates: %i %i %i\n", ncand[0], ncand[1], ncand[2]); - } - - for ( i=0; i angtol ) continue; - - fom1 = fabs(ang - angles[2]); - - for ( k=0; k angtol ) continue; - - fom2 = fom1 + fabs(ang - angles[1]); - - /* Finally, the angle between the current candidate for - * axis 1 and the kth candidate for axis 2 */ - ang = angle_between(cand[1][j].vec.u, cand[1][j].vec.v, - cand[1][j].vec.w, cand[2][k].vec.u, - cand[2][k].vec.v, cand[2][k].vec.w); - - /* ... it should be angle 0 ... */ - if ( fabs(ang - angles[0]) > angtol ) continue; - - /* Unit cell must be right-handed */ - if ( !right_handed_vec(cand[0][i].vec, cand[1][j].vec, - cand[2][k].vec) ) continue; - - fom3 = fom2 + fabs(ang - angles[0]); - fom3 += LWEIGHT * (cand[0][i].fom + cand[1][j].fom - + cand[2][k].fom); - - if ( fom3 < best_fom ) { - if ( new_cell != NULL ) free(new_cell); - new_cell = cell_new_from_reciprocal_axes( - cand[0][i].vec, cand[1][j].vec, - cand[2][k].vec); - best_fom = fom3; - } - - } - - } - } - - free(cand[0]); - free(cand[1]); - free(cand[2]); - - cell_free(cell); - - /* Reverse the de-centering transformation */ - if ( new_cell != NULL ) { - - new_cell_trans = cell_transform_intmat(new_cell, centering); - cell_free(new_cell); - cell_set_lattice_type(new_cell_trans, - cell_get_lattice_type(template_in)); - cell_set_centering(new_cell_trans, - cell_get_centering(template_in)); - cell_set_unique_axis(new_cell_trans, - cell_get_unique_axis(template_in)); - - cell_free(template); - intmat_free(centering); - - return new_cell_trans; - - } else { - cell_free(template); - intmat_free(centering); - return NULL; - } -} - - -UnitCell *match_cell_ab(UnitCell *cell_in, UnitCell *template_in) -{ - double ax, ay, az; - double bx, by, bz; - double cx, cy, cz; - int i; - double lengths[3]; - int used[3]; - struct rvec real_a, real_b, real_c; - struct rvec params[3]; - double alen, blen; - float ltl = 5.0; /* percent */ - int have_real_a; - int have_real_b; - int have_real_c; - UnitCell *cell; - UnitCell *template; - IntegerMatrix *to_given_cell; - UnitCell *new_cell; - UnitCell *new_cell_trans; - - /* "Un-center" the template unit cell to make the comparison easier */ - template = uncenter_cell(template_in, &to_given_cell, NULL); - - /* The candidate cell is also uncentered, because it might be centered - * if it came from (e.g.) MOSFLM */ - cell = uncenter_cell(cell_in, NULL, NULL); - - /* Get the lengths to match */ - if ( cell_get_cartesian(template, &ax, &ay, &az, - &bx, &by, &bz, - &cx, &cy, &cz) ) - { - ERROR("Couldn't get cell for template.\n"); - return NULL; - } - alen = modulus(ax, ay, az); - blen = modulus(bx, by, bz); - - /* Get the lengths from the cell and turn them into anonymous vectors */ - if ( cell_get_cartesian(cell, &ax, &ay, &az, - &bx, &by, &bz, - &cx, &cy, &cz) ) - { - ERROR("Couldn't get cell.\n"); - return NULL; - } - lengths[0] = modulus(ax, ay, az); - lengths[1] = modulus(bx, by, bz); - lengths[2] = modulus(cx, cy, cz); - used[0] = 0; used[1] = 0; used[2] = 0; - params[0].u = ax; params[0].v = ay; params[0].w = az; - params[1].u = bx; params[1].v = by; params[1].w = bz; - params[2].u = cx; params[2].v = cy; params[2].w = cz; - - real_a.u = 0.0; real_a.v = 0.0; real_a.w = 0.0; - real_b.u = 0.0; real_b.v = 0.0; real_b.w = 0.0; - real_c.u = 0.0; real_c.v = 0.0; real_c.w = 0.0; - - /* Check each vector against a and b */ - have_real_a = 0; - have_real_b = 0; - for ( i=0; i<3; i++ ) { - if ( within_tolerance(lengths[i], alen, ltl) - && !used[i] && !have_real_a ) - { - used[i] = 1; - memcpy(&real_a, ¶ms[i], sizeof(struct rvec)); - have_real_a = 1; - } - if ( within_tolerance(lengths[i], blen, ltl) - && !used[i] && !have_real_b ) - { - used[i] = 1; - memcpy(&real_b, ¶ms[i], sizeof(struct rvec)); - have_real_b = 1; - } - } - - /* Have we matched both a and b? */ - if ( !(have_real_a && have_real_b) ) return NULL; - - /* "c" is "the other one" */ - have_real_c = 0; - for ( i=0; i<3; i++ ) { - if ( !used[i] ) { - memcpy(&real_c, ¶ms[i], sizeof(struct rvec)); - have_real_c = 1; - } - } - - if ( !have_real_c ) { - ERROR("Huh? Couldn't find the third vector.\n"); - ERROR("Matches: %i %i %i\n", used[0], used[1], used[2]); - return NULL; - } - - /* Flip c if not right-handed */ - if ( !right_handed_vec(real_a, real_b, real_c) ) { - real_c.u = -real_c.u; - real_c.v = -real_c.v; - real_c.w = -real_c.w; - } - - new_cell = cell_new_from_direct_axes(real_a, real_b, real_c); - - /* Reverse the de-centering transformation */ - new_cell_trans = cell_transform_intmat_inverse(new_cell, to_given_cell); - 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)); - cell_set_unique_axis(new_cell, cell_get_unique_axis(template_in)); - - return new_cell_trans; -} - - /* Return sin(theta)/lambda = 1/2d. Multiply by two if you want 1/d */ double resolution(UnitCell *cell, signed int h, signed int k, signed int l) { diff --git a/libcrystfel/src/cell-utils.h b/libcrystfel/src/cell-utils.h index bb1e5cd0..5b705247 100644 --- a/libcrystfel/src/cell-utils.h +++ b/libcrystfel/src/cell-utils.h @@ -59,11 +59,6 @@ extern UnitCell *rotate_cell(UnitCell *in, double omega, double phi, extern void cell_print(UnitCell *cell); extern void cell_print_full(UnitCell *cell); -extern UnitCell *match_cell(UnitCell *cell, UnitCell *tempcell, int verbose, - const float *ltl, int reduce); - -extern UnitCell *match_cell_ab(UnitCell *cell, UnitCell *tempcell); - extern UnitCell *load_cell_from_pdb(const char *filename); extern UnitCell *load_cell_from_file(const char *filename); extern void write_cell(UnitCell *cell, FILE *fh); -- cgit v1.2.3 From 184a21098ac1325c2c02a24ac315584a464532a7 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Fri, 23 Aug 2019 12:21:00 +0200 Subject: Fix Niggli reduction logic --- libcrystfel/src/cell-utils.c | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) (limited to 'libcrystfel/src') diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c index 2af3f293..6e7cebb5 100644 --- a/libcrystfel/src/cell-utils.c +++ b/libcrystfel/src/cell-utils.c @@ -1893,9 +1893,10 @@ IntegerMatrix *reduce_g6(struct g6 g, double epsrel) debug_lattice(g, eps, 1); + } - } else if ( GT(g.B, g.C) - || (EQ(g.B, g.C) && GT(fabs(g.E), fabs(g.F))) ) + if ( GT(g.B, g.C) + || (EQ(g.B, g.C) && GT(fabs(g.E), fabs(g.F))) ) { /* Swap b and c */ double temp; @@ -1909,6 +1910,8 @@ IntegerMatrix *reduce_g6(struct g6 g, double epsrel) mult_in_place(T, M); debug_lattice(g, eps, 2); + + /* ..."and go to 1." */ done_s1s2 = 0; } -- cgit v1.2.3 From da1de864442781c2bb993fc15a5c9dca18940e38 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Wed, 28 Aug 2019 13:29:48 +0200 Subject: Avoid converting IntegerMatrix to RationalMatrix --- libcrystfel/src/cell-utils.c | 33 +++++++++++++++++++++------------ libcrystfel/src/rational.c | 23 +++++++++++++++++++++++ libcrystfel/src/rational.h | 2 ++ 3 files changed, 46 insertions(+), 12 deletions(-) (limited to 'libcrystfel/src') diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c index 6e7cebb5..dbc845bc 100644 --- a/libcrystfel/src/cell-utils.c +++ b/libcrystfel/src/cell-utils.c @@ -1844,6 +1844,23 @@ static void rtnl_mult_in_place(RationalMatrix *T, RationalMatrix *M) } +/* Replace the elements of T with those of T*M */ +static void rtnl_int_mult_in_place(RationalMatrix *T, IntegerMatrix *M) +{ + int i, j; + RationalMatrix *tmp; + + tmp = rtnl_mtx_new(3, 3); + rtnl_mtx_intmatmult(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) { @@ -2243,24 +2260,16 @@ UnitCell *compare_reindexed_cell_parameters(UnitCell *cell_in, //intmat_print(P); /* 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_int_mult_in_place(comb, RA); + rtnl_int_mult_in_place(comb, P); + rtnl_int_mult_in_place(comb, RiB); + rtnl_int_mult_in_place(comb, CB); - rtnl_mtx_free(rRA); - rtnl_mtx_free(rP); intmat_free(RiB); - rtnl_mtx_free(rRiB); - rtnl_mtx_free(rCB); match = cell_transform_rational(cell_in, comb); //STATUS("Original cell transformed to look like reference:\n"); diff --git a/libcrystfel/src/rational.c b/libcrystfel/src/rational.c index ce286d7a..f8d3cb68 100644 --- a/libcrystfel/src/rational.c +++ b/libcrystfel/src/rational.c @@ -553,6 +553,29 @@ void rtnl_mtx_mtxmult(const RationalMatrix *A, const RationalMatrix *B, } +void rtnl_mtx_intmatmult(const RationalMatrix *A, const IntegerMatrix *B, + RationalMatrix *ans) +{ + int i, j; + + assert(ans->rows == A->rows); + + for ( i=0; irows; i++ ) { + for ( j=0; jcols; j++ ) { + int k; + Rational sum = rtnl_zero(); + for ( k=0; krows; k++ ) { + Rational add; + add = rtnl_mul(rtnl_mtx_get(A, i, k), + rtnl(intmat_get(B, k, j), 1)); + sum = rtnl_add(sum, add); + } + rtnl_mtx_set(ans, i, j, sum); + } + } +} + + /* Given a "P-matrix" (see ITA chapter 5.1), calculate the fractional * coordinates of point "vec" in the original axes, given its fractional * coordinates in the transformed axes. diff --git a/libcrystfel/src/rational.h b/libcrystfel/src/rational.h index 8681bb06..012b929e 100644 --- a/libcrystfel/src/rational.h +++ b/libcrystfel/src/rational.h @@ -92,6 +92,8 @@ extern IntegerMatrix *intmat_from_rtnl_mtx(const RationalMatrix *m); extern void rtnl_mtx_free(RationalMatrix *mtx); extern void rtnl_mtx_mtxmult(const RationalMatrix *A, const RationalMatrix *B, RationalMatrix *ans); +extern void rtnl_mtx_intmatmult(const RationalMatrix *A, const IntegerMatrix *B, + RationalMatrix *ans); extern int transform_fractional_coords_rtnl(const RationalMatrix *P, const Rational *ivec, Rational *ans); -- cgit v1.2.3 From 24e2744b8b53dde715f59d98f969d44a10530675 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Wed, 28 Aug 2019 15:33:46 +0200 Subject: Rationalise matrix mutliplication compare_reindexed_cell_parameters still needs to be updated --- libcrystfel/src/cell-utils.c | 46 ++++-------------------------- libcrystfel/src/integer_matrix.c | 4 +-- libcrystfel/src/integer_matrix.h | 4 +-- libcrystfel/src/rational.c | 61 +++++++++++++++++++++++++++++++++------- libcrystfel/src/rational.h | 11 +++++--- libcrystfel/src/symmetry.c | 10 +++---- 6 files changed, 72 insertions(+), 64 deletions(-) (limited to 'libcrystfel/src') diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c index dbc845bc..de8fe7ae 100644 --- a/libcrystfel/src/cell-utils.c +++ b/libcrystfel/src/cell-utils.c @@ -1612,7 +1612,6 @@ int compare_derivative_cell_parameters(UnitCell *cell_in, UnitCell *reference_in Rational *cand_c; int ncand_a, ncand_b, ncand_c; int ia, ib; - RationalMatrix *MCB; RationalMatrix *CiAMCB; double min_dist = +INFINITY; @@ -1647,8 +1646,6 @@ int compare_derivative_cell_parameters(UnitCell *cell_in, UnitCell *reference_in } M = rtnl_mtx_new(3, 3); - MCB = rtnl_mtx_new(3, 3); - CiAMCB = rtnl_mtx_new(3, 3); for ( ia=0; iacols == B->cols); - assert(ans->rows == A->rows); - assert(A->cols == B->rows); + intmat_size(B, &B_rows, &B_cols); + assert(A->cols == B_rows); + + ans = rtnl_mtx_new(A->rows, B_cols); + if ( ans == NULL ) return NULL; for ( i=0; irows; i++ ) { for ( j=0; jcols; j++ ) { @@ -544,21 +548,27 @@ void rtnl_mtx_mtxmult(const RationalMatrix *A, const RationalMatrix *B, for ( k=0; krows; k++ ) { Rational add; add = rtnl_mul(rtnl_mtx_get(A, i, k), - rtnl_mtx_get(B, k, j)); + rtnl(intmat_get(B, k, j), 1)); sum = rtnl_add(sum, add); } rtnl_mtx_set(ans, i, j, sum); } } + + return ans; } -void rtnl_mtx_intmatmult(const RationalMatrix *A, const IntegerMatrix *B, - RationalMatrix *ans) +RationalMatrix *rtnlmtx_times_rtnlmtx(const RationalMatrix *A, + const RationalMatrix *B) { int i, j; + RationalMatrix *ans; + + assert(A->cols == B->rows); - assert(ans->rows == A->rows); + ans = rtnl_mtx_new(A->rows, B->cols); + if ( ans == NULL ) return NULL; for ( i=0; irows; i++ ) { for ( j=0; jcols; j++ ) { @@ -567,12 +577,43 @@ void rtnl_mtx_intmatmult(const RationalMatrix *A, const IntegerMatrix *B, for ( k=0; krows; k++ ) { Rational add; add = rtnl_mul(rtnl_mtx_get(A, i, k), - rtnl(intmat_get(B, k, j), 1)); + rtnl_mtx_get(B, k, j)); + sum = rtnl_add(sum, add); + } + rtnl_mtx_set(ans, i, j, sum); + } + } + + return ans; +} + + +RationalMatrix *intmat_times_rtnlmtx(const IntegerMatrix *A, const RationalMatrix *B) +{ + int i, j; + RationalMatrix *ans; + unsigned int A_rows, A_cols; + + intmat_size(A, &A_rows, &A_cols); + + ans = rtnl_mtx_new(A_rows, B->cols); + if ( ans == NULL ) return NULL; + + for ( i=0; irows; i++ ) { + for ( j=0; jcols; j++ ) { + int k; + Rational sum = rtnl_zero(); + for ( k=0; kops[i]); + r = intmat_times_intmat(P, s->ops[i]); if ( r == NULL ) { ERROR("Matrix multiplication failed.\n"); return; } - f = intmat_intmat_mult(r, Pi); + f = intmat_times_intmat(r, Pi); if ( f == NULL ) { ERROR("Matrix multiplication failed.\n"); return; @@ -1340,7 +1340,7 @@ static int check_mult(const IntegerMatrix *ans, int val; IntegerMatrix *m; - m = intmat_intmat_mult(a, b); + m = intmat_times_intmat(a, b); assert(m != NULL); val = intmat_equals(ans, m); @@ -1403,7 +1403,7 @@ static int order(const IntegerMatrix *m) IntegerMatrix *anew; - anew = intmat_intmat_mult(m, a); + anew = intmat_times_intmat(m, a); assert(anew != NULL); intmat_free(a); a = anew; -- cgit v1.2.3 From 865809d07a25b70b4ae0697d509a8f5e7a17c625 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Wed, 28 Aug 2019 16:49:55 +0200 Subject: compare_reindexed_cell_parameters: (Sort of) intelligent choice of solution --- libcrystfel/src/cell-utils.c | 134 ++++++++++++++++++++++++++++++++----------- 1 file changed, 99 insertions(+), 35 deletions(-) (limited to 'libcrystfel/src') diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c index de8fe7ae..d92bd6a9 100644 --- a/libcrystfel/src/cell-utils.c +++ b/libcrystfel/src/cell-utils.c @@ -2059,32 +2059,33 @@ static double cell_diff(UnitCell *cell, double a, double b, double c, } -static IntegerMatrix *check_permutations(UnitCell *cell, UnitCell *reference, - IntegerMatrix *RB, IntegerMatrix *CB, +static int random_int(int max) +{ + int r; + int limit = RAND_MAX; + while ( limit % max ) limit--; + do { + r = rand(); + } while ( r > limit ); + return rand() % max; +} + + +static IntegerMatrix *check_permutations(UnitCell *cell_reduced, UnitCell *reference, + RationalMatrix *CiARA, IntegerMatrix *RiBCB, const double *tols) { IntegerMatrix *m; int i[9]; double a, b, c, al, be, ga; double min_dist = +INFINITY; - IntegerMatrix *best_m = NULL; - RationalMatrix *RiBCB; - RationalMatrix *rRiB; - RationalMatrix *rCB; - IntegerMatrix *RiB; + int s, sel; + IntegerMatrix *best_m[5] = {NULL, NULL, NULL, NULL, NULL}; + int n_best = 0; m = intmat_new(3, 3); cell_get_parameters(reference, &a, &b, &c, &al, &be, &ga); - RiBCB = rtnl_mtx_new(3, 3); - RiB = intmat_inverse(RB); - rRiB = rtnl_mtx_from_intmat(RiB); - rCB = rtnl_mtx_from_intmat(CB); - rtnl_mtx_mtxmult(rRiB, rCB, RiBCB); - intmat_free(RiB); - rtnl_mtx_free(rRiB); - 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]++ ) { @@ -2108,16 +2109,34 @@ static IntegerMatrix *check_permutations(UnitCell *cell, UnitCell *reference, det = intmat_det(m); if ( det != +1 ) continue; - tmp = cell_transform_intmat(cell, m); - nc = cell_transform_rational(tmp, RiBCB); + tmp = cell_transform_intmat(cell_reduced, m); + nc = cell_transform_intmat(tmp, RiBCB); 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); + + /* If the new solution is significantly better, + * dump all the previous ones */ + for ( s=0; s