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 +++--- src/indexamajig.c | 37 +++++++++++++++++++++++++++---------- src/process_image.h | 2 +- 3 files changed, 31 insertions(+), 14 deletions(-) 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; diff --git a/src/indexamajig.c b/src/indexamajig.c index 739b4969..78f26709 100644 --- a/src/indexamajig.c +++ b/src/indexamajig.c @@ -135,7 +135,7 @@ static void show_help(const char *s) " -p, --pdb= Unit cell file (PDB or CrystFEL unit cell format)\n" " Default: 'molecule.pdb'\n" " --tolerance= Tolerances for cell comparison\n" -" Default: 5,5,5,1.5\n" +" Default: 5,5,5,1.5,1.5,1.5\n" " --no-check-cell Don't check lattice parameters against input cell\n" " --no-cell-combinations\n" " Don't use axis combinations when checking cell\n" @@ -297,10 +297,12 @@ int main(int argc, char *argv[]) iargs.noisefilter = 0; iargs.median_filter = 0; iargs.satcorr = 1; - iargs.tols[0] = 5.0; - iargs.tols[1] = 5.0; - iargs.tols[2] = 5.0; + iargs.tols[0] = 0.05; + iargs.tols[1] = 0.05; + iargs.tols[2] = 0.05; iargs.tols[3] = 1.5; + iargs.tols[4] = 1.5; + iargs.tols[5] = 1.5; iargs.threshold = 800.0; iargs.min_sq_gradient = 100000.0; iargs.min_snr = 5.0; @@ -1056,14 +1058,29 @@ int main(int argc, char *argv[]) /* Parse unit cell tolerance */ if ( toler != NULL ) { int ttt; - ttt = sscanf(toler, "%f,%f,%f,%f", - &iargs.tols[0], &iargs.tols[1], - &iargs.tols[2], &iargs.tols[3]); - if ( ttt != 4 ) { - ERROR("Invalid parameters for '--tolerance'\n"); - return 1; + ttt = sscanf(toler, "%f,%f,%f,%f,%f,%f", + &iargs.tols[0], &iargs.tols[1], &iargs.tols[2], + &iargs.tols[3], &iargs.tols[4], &iargs.tols[5]); + if ( ttt != 6 ) { + ttt = sscanf(toler, "%f,%f,%f,%f", + &iargs.tols[0], &iargs.tols[1], + &iargs.tols[2], &iargs.tols[3]); + if ( ttt != 4 ) { + ERROR("Invalid parameters for '--tolerance'\n"); + return 1; + } + iargs.tols[4] = iargs.tols[3]; + iargs.tols[5] = iargs.tols[3]; } free(toler); + + /* Percent to fraction */ + iargs.tols[0] /= 100.0; + iargs.tols[1] /= 100.0; + iargs.tols[2] /= 100.0; + iargs.tols[3] = deg2rad(iargs.tols[3]); + iargs.tols[4] = deg2rad(iargs.tols[4]); + iargs.tols[5] = deg2rad(iargs.tols[5]); } /* Parse integration radii */ diff --git a/src/process_image.h b/src/process_image.h index 8d6682a6..f8b0f4b1 100644 --- a/src/process_image.h +++ b/src/process_image.h @@ -74,7 +74,7 @@ struct index_args struct detector *det; IndexingPrivate *ipriv; int peaks; /* Peak detection method */ - float tols[4]; + float tols[6]; struct beam_params *beam; char *hdf5_peak_path; int half_pixel_shift; -- 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 ++-- src/cell_tool.c | 20 +++++++++++++++++-- 3 files changed, 42 insertions(+), 28 deletions(-) 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 diff --git a/src/cell_tool.c b/src/cell_tool.c index cd457f25..a1755f42 100644 --- a/src/cell_tool.c +++ b/src/cell_tool.c @@ -77,6 +77,7 @@ static int comparecells(UnitCell *cell, const char *comparecell, { UnitCell *cell2; RationalMatrix *m; + double tolerance[6]; cell2 = load_cell_from_file(comparecell); if ( cell2 == NULL ) { @@ -90,8 +91,15 @@ static int comparecells(UnitCell *cell, const char *comparecell, STATUS("------------------> The reference unit cell:\n"); cell_print(cell2); + tolerance[0] = ltl; + tolerance[1] = ltl; + tolerance[2] = ltl; + tolerance[3] = atl; + tolerance[4] = atl; + tolerance[5] = atl; + STATUS("------------------> The comparison results:\n"); - if ( !compare_reindexed_cell_parameters(cell, cell2, ltl, atl, csl, &m) ) { + if ( !compare_reindexed_cell_parameters(cell, cell2, tolerance, csl, &m) ) { STATUS("No relationship found between lattices.\n"); return 0; } else { @@ -227,6 +235,14 @@ static int find_ambi(UnitCell *cell, SymOpList *sym, double ltl, double atl) SymOpList *ops; signed int i[9]; const int maxorder = 3; + double tolerance[6]; + + tolerance[0] = ltl; + tolerance[1] = ltl; + tolerance[2] = ltl; + tolerance[3] = atl; + tolerance[4] = atl; + tolerance[5] = atl; ops = get_pointgroup("1"); if ( ops == NULL ) return 1; @@ -265,7 +281,7 @@ static int find_ambi(UnitCell *cell, SymOpList *sym, double ltl, double atl) nc = cell_transform_intmat(cell, m); - if ( compare_cell_parameters(cell, nc, ltl, atl) ) { + if ( compare_cell_parameters(cell, nc, tolerance) ) { if ( !intmat_is_identity(m) ) add_symop(ops, m); STATUS("-----------------------------------------------" "-------------------------------------------\n"); -- 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(+) 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(-) 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 ++--- src/whirligig.c | 12 ++++----- 4 files changed, 52 insertions(+), 48 deletions(-) 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); } diff --git a/src/whirligig.c b/src/whirligig.c index cc2b531e..a144e769 100644 --- a/src/whirligig.c +++ b/src/whirligig.c @@ -309,9 +309,9 @@ static IntegerMatrix *try_all(struct window *win, int n1, int n2, for ( i=0; in_crystals; i++ ) { for ( j=0; jn_crystals; j++ ) { - if ( compare_reindexed_cell_parameters_and_orientation(crystal_get_cell(i1->crystals[i]), - crystal_get_cell(i2->crystals[j]), - 0.1, deg2rad(5.0), &m) ) + if ( compare_permuted_cell_parameters_and_orientation(crystal_get_cell(i1->crystals[i]), + crystal_get_cell(i2->crystals[j]), + 0.1, deg2rad(5.0), &m) ) { if ( !crystal_used(win, n1, i) && !crystal_used(win, n2, j) ) @@ -374,9 +374,9 @@ static int try_join(struct window *win, int sn) for ( j=0; jimg[win->join_ptr].n_crystals; j++ ) { Crystal *cr2; cr2 = win->img[win->join_ptr].crystals[j]; - if ( compare_reindexed_cell_parameters_and_orientation(ref, crystal_get_cell(cr2), - 0.1, deg2rad(5.0), - &win->mat[sn][win->join_ptr]) ) + if ( compare_permuted_cell_parameters_and_orientation(ref, crystal_get_cell(cr2), + 0.1, deg2rad(5.0), + &win->mat[sn][win->join_ptr]) ) { win->ser[sn][win->join_ptr] = j; cell_free(ref); -- 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(-) 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 9507e87653ea18562b95cf02d66f7b285893c2af Mon Sep 17 00:00:00 2001 From: Thomas White Date: Thu, 8 Aug 2019 11:06:34 +0200 Subject: Add cellcompare_check --- tests/CMakeLists.txt | 5 + tests/cellcompare_check.c | 274 ++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 279 insertions(+) create mode 100644 tests/cellcompare_check.c diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 192913ce..bd9959ae 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -92,3 +92,8 @@ add_executable(spectrum_check spectrum_check.c) target_include_directories(spectrum_check PRIVATE ${COMMON_INCLUDES}) target_link_libraries(spectrum_check ${COMMON_LIBRARIES}) add_test(spectrum_check spectrum_check) + +add_executable(cellcompare_check cellcompare_check.c) +target_include_directories(cellcompare_check PRIVATE ${COMMON_INCLUDES}) +target_link_libraries(cellcompare_check ${COMMON_LIBRARIES}) +add_test(cellcompare_check cellcompare_check) diff --git a/tests/cellcompare_check.c b/tests/cellcompare_check.c new file mode 100644 index 00000000..68122c03 --- /dev/null +++ b/tests/cellcompare_check.c @@ -0,0 +1,274 @@ +/* + * cellcompare_check.c + * + * Check that unit cell comparison works + * + * Copyright © 2019 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. + * + * Authors: + * 2019 Thomas White + * + * This file is part of CrystFEL. + * + * CrystFEL is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * CrystFEL is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with CrystFEL. If not, see . + * + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + + +#include +#include +#include +#include + +#include +#include +#include + + +static double moduli_check(double ax, double ay, double az, + double bx, double by, double bz) +{ + double ma = modulus(ax, ay, az); + double mb = modulus(bx, by, bz); + return fabs(ma-mb)/ma; +} +static void complain(UnitCell *cell, UnitCell *cref, const char *t) +{ + STATUS("These cells should %s be the same:\n", t); + STATUS("Transformed: ----------------------------\n"); + cell_print_full(cell); + STATUS("Original: ---------------------------\n"); + cell_print_full(cref); +} + + +static RationalMatrix *random_reindexing(gsl_rng *rng) +{ + int i; + RationalMatrix *tr; + int v[] = {0, 1, 2}; + + tr = rtnl_mtx_new(3, 3); + if ( tr == NULL ) return NULL; + + gsl_ran_shuffle(rng, v, 3, sizeof(int)); + for ( i=0; i<3; i++ ) { + rtnl_mtx_set(tr, i, v[i], rtnl(1, 1)); + } + + return tr; +} + + +static IntegerMatrix *random_permutation(gsl_rng *rng) +{ + int i; + IntegerMatrix *tr; + int v[] = {0, 1, 2}; + + tr = intmat_new(3, 3); + if ( tr == NULL ) return NULL; + + gsl_ran_shuffle(rng, v, 3, sizeof(int)); + for ( i=0; i<3; i++ ) { + intmat_set(tr, i, v[i], 1); + } + + return tr; +} + + +int compare_cell_parameters_and_orientation2(UnitCell *cell1, UnitCell *cell2, + 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; + + cell_get_cartesian(cell1, &ax1, &ay1, &az1, + &bx1, &by1, &bz1, + &cx1, &cy1, &cz1); + + cell_get_cartesian(cell2, &ax2, &ay2, &az2, + &bx2, &by2, &bz2, + &cx2, &cy2, &cz2); + + cell_print_full(cell1); + cell_print_full(cell2); + + STATUS("%f\n", rad2deg(atl)); + STATUS("%f\n", rad2deg(angle_between(ax1, ay1, az1, ax2, ay2, az2))); + STATUS("%f\n", rad2deg(angle_between(bx1, by1, bz1, bx2, by2, bz2))); + STATUS("%f\n", rad2deg(angle_between(cx1, cy1, cz1, 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 ( 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; + + return 1; +} + + +int main(int argc, char *argv[]) +{ + UnitCell *cell, *cref; + gsl_rng *rng; + int i; + double tols[] = { 0.01, 0.01, 0.01, + deg2rad(1.0), deg2rad(1.0), deg2rad(1.0) }; + + rng = gsl_rng_alloc(gsl_rng_mt19937); + + cref = cell_new_from_parameters(50e-10, 55e-10, 70e-10, + deg2rad(67.0), + deg2rad(70.0), + deg2rad(77.0)); + if ( cref == NULL ) return 1; + + /* Just rotate cell */ + STATUS("Testing plain rotation...\n"); + for ( i=0; i<100; i++ ) { + cell = cell_rotate(cref, random_quaternion(rng)); + if ( cell == NULL ) return 1; + if ( !compare_cell_parameters(cell, cref, tols) ) { + complain(cell, cref, ""); + return 1; + } + if ( compare_cell_parameters_and_orientation(cell, cref, + tols[0], tols[3]) ) + { + complain(cell, cref, "NOT"); + return 1; + } + cell_free(cell); + } + + /* Permute axes but don't rotate */ + STATUS("Testing axis permutation...\n"); + for ( i=0; i<100; i++ ) { + + IntegerMatrix *m; + IntegerMatrix *tr; + + tr = random_permutation(rng); + cell = cell_transform_intmat(cref, tr); + + if ( !compare_permuted_cell_parameters_and_orientation(cell, cref, + tols[0], tols[3], &m) ) + { + complain(cell, cref, ""); + return 1; + } + + if ( compare_cell_parameters(cell, cref, tols) + && !intmat_is_identity(tr) ) + { + complain(cell, cref, "NOT"); + return 1; + } + + cell_free(cell); + intmat_free(tr); + intmat_free(m); + } + + /* Rotate cell and permute axes */ + STATUS("Testing rotation with axis permutation...\n"); + for ( i=0; i<100; i++ ) { + + IntegerMatrix *m; + IntegerMatrix *tr; + UnitCell *cell2; + + cell = cell_rotate(cref, random_quaternion(rng)); + if ( cell == NULL ) return 1; + + tr = random_permutation(rng); + cell2 = cell_transform_intmat(cell, tr); + + if ( !compare_permuted_cell_parameters(cell2, cref, tols, &m) ) { + complain(cell2, cref, ""); + return 1; + } + + if ( compare_permuted_cell_parameters_and_orientation(cell2, cref, + tols[0], tols[3], &m) ) + { + UnitCell *cc; + complain(cell2, cref, "NOT, with just permutation,"); + STATUS("Matrix was (det=%i):\n", intmat_det(m)); + intmat_print(m); + STATUS("Transformed version of cref:\n"); + cc = cell_transform_intmat(cref, m); + cell_print_full(cc); + cell_free(cc); + return 1; + } + + if ( compare_cell_parameters_and_orientation(cell2, cref, + tols[0], tols[3]) ) + { + complain(cell2, cref, "NOT, without any change,"); + return 1; + } + + cell_free(cell); + cell_free(cell2); + intmat_free(tr); + intmat_free(m); + } + + /* Reindex and rotate */ + STATUS("Testing rotation with reindexing...\n"); + for ( i=0; i<100; i++ ) { + + RationalMatrix *m; + RationalMatrix *tr; + UnitCell *cell2; + + cell = cell_rotate(cref, random_quaternion(rng)); + if ( cell == NULL ) return 1; + + tr = random_reindexing(rng); + cell2 = cell_transform_rational(cell, tr); + cell_free(cell); + + if ( !compare_reindexed_cell_parameters(cell2, cref, tols, 0, &m) ) { + complain(cell2, cref, ""); + return 1; + } + cell_free(cell2); + rtnl_mtx_free(tr); + rtnl_mtx_free(m); + } + + /* NB There's no compare_reindexed_cell_parameters_and_orientation */ + + /* Test using vectors */ + + cell_free(cref); + gsl_rng_free(rng); + + return 0; +} -- cgit v1.2.3 From 78790e418b1f43c176f3829cef1e44feafd79044 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Thu, 8 Aug 2019 14:22:02 +0200 Subject: "Final" cellcompare_check --- tests/cellcompare_check.c | 249 +++++++++++++++++++++++++++------------------- 1 file changed, 147 insertions(+), 102 deletions(-) diff --git a/tests/cellcompare_check.c b/tests/cellcompare_check.c index 68122c03..234d93a8 100644 --- a/tests/cellcompare_check.c +++ b/tests/cellcompare_check.c @@ -41,16 +41,9 @@ #include -static double moduli_check(double ax, double ay, double az, - double bx, double by, double bz) +static void complain(UnitCell *cell, UnitCell *cref, const char *t, const char *c) { - double ma = modulus(ax, ay, az); - double mb = modulus(bx, by, bz); - return fabs(ma-mb)/ma; -} -static void complain(UnitCell *cell, UnitCell *cref, const char *t) -{ - STATUS("These cells should %s be the same:\n", t); + STATUS("These cells should %sbe the same%s:\n", t, c); STATUS("Transformed: ----------------------------\n"); cell_print_full(cell); STATUS("Original: ---------------------------\n"); @@ -94,39 +87,106 @@ static IntegerMatrix *random_permutation(gsl_rng *rng) } -int compare_cell_parameters_and_orientation2(UnitCell *cell1, UnitCell *cell2, - const double ltl, const double atl) +static int check_ccp(UnitCell *cell, UnitCell *cref, double *tols, + int should_match) { - double ax1, ay1, az1, bx1, by1, bz1, cx1, cy1, cz1; - double ax2, ay2, az2, bx2, by2, bz2, cx2, cy2, cz2; + const char *a; + const char *b; + + a = should_match ? "" : "NOT "; + b = " with compare_cell_parameters"; + + if ( compare_cell_parameters(cell, cref, tols) != should_match ) + { + complain(cell, cref, a, b); + return 1; + } + return 0; +} - if ( cell_get_centering(cell1) != cell_get_centering(cell2) ) return 0; - cell_get_cartesian(cell1, &ax1, &ay1, &az1, - &bx1, &by1, &bz1, - &cx1, &cy1, &cz1); +static int check_cpcp(UnitCell *cell, UnitCell *cref, double *tols, + int should_match) +{ + IntegerMatrix *m = NULL; + const char *a; + const char *b; + + a = should_match ? "" : "NOT "; + b = " with compare_permuted_cell_parameters"; + + if ( compare_permuted_cell_parameters(cell, cref, tols, &m) != should_match ) + { + complain(cell, cref, a, b); + STATUS("Matrix was:\n"); + intmat_print(m); + intmat_free(m); + return 1; + } + intmat_free(m); + return 0; +} - cell_get_cartesian(cell2, &ax2, &ay2, &az2, - &bx2, &by2, &bz2, - &cx2, &cy2, &cz2); - cell_print_full(cell1); - cell_print_full(cell2); +static int check_ccpao(UnitCell *cell, UnitCell *cref, double *tols, + int should_match) +{ + const char *a; + const char *b; - STATUS("%f\n", rad2deg(atl)); - STATUS("%f\n", rad2deg(angle_between(ax1, ay1, az1, ax2, ay2, az2))); - STATUS("%f\n", rad2deg(angle_between(bx1, by1, bz1, bx2, by2, bz2))); - STATUS("%f\n", rad2deg(angle_between(cx1, cy1, cz1, cx2, cy2, cz2))); + a = should_match ? "" : "NOT "; + b = " with compare_cell_parameters_and_orientation"; - 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 ( compare_cell_parameters_and_orientation(cell, cref, + tols[0], tols[3]) != should_match ) + { + complain(cell, cref, a, b); + return 1; + } + 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; - return 1; +static int check_cpcpao(UnitCell *cell, UnitCell *cref, double *tols, + int should_match) +{ + IntegerMatrix *m = NULL; + const char *a; + const char *b; + + a = should_match ? "" : "NOT "; + b = " with compare_permuted_cell_parameters_and_orientation"; + + if ( compare_permuted_cell_parameters_and_orientation(cell, cref, + tols[0], tols[3], &m) != should_match ) + { + complain(cell, cref, a, b); + intmat_free(m); + return 1; + } + intmat_free(m); + return 0; +} + + +static int check_crcp(UnitCell *cell, UnitCell *cref, double *tols, + int should_match) +{ + RationalMatrix *m = NULL; + const char *a; + const char *b; + + a = should_match ? "" : "NOT "; + b = " with compare_reindexed_cell_parameters"; + + if ( compare_reindexed_cell_parameters(cell, cref, tols, 0, &m) != should_match ) + { + complain(cell, cref, a, b); + rtnl_mtx_free(m); + return 1; + } + rtnl_mtx_free(m); + return 0; } @@ -149,18 +209,16 @@ int main(int argc, char *argv[]) /* Just rotate cell */ STATUS("Testing plain rotation...\n"); for ( i=0; i<100; i++ ) { + cell = cell_rotate(cref, random_quaternion(rng)); if ( cell == NULL ) return 1; - if ( !compare_cell_parameters(cell, cref, tols) ) { - complain(cell, cref, ""); - return 1; - } - if ( compare_cell_parameters_and_orientation(cell, cref, - tols[0], tols[3]) ) - { - complain(cell, cref, "NOT"); - return 1; - } + + if ( check_ccp(cell, cref, tols, 1) ) return 1; + if ( check_cpcp(cell, cref, tols, 1) ) return 1; + if ( check_ccpao(cell, cref, tols, 0) ) return 1; + if ( check_cpcpao(cell, cref, tols, 0) ) return 1; + if ( check_crcp(cell, cref, tols, 1) ) return 1; + cell_free(cell); } @@ -168,99 +226,86 @@ int main(int argc, char *argv[]) STATUS("Testing axis permutation...\n"); for ( i=0; i<100; i++ ) { - IntegerMatrix *m; IntegerMatrix *tr; tr = random_permutation(rng); cell = cell_transform_intmat(cref, tr); - if ( !compare_permuted_cell_parameters_and_orientation(cell, cref, - tols[0], tols[3], &m) ) - { - complain(cell, cref, ""); - return 1; - } - - if ( compare_cell_parameters(cell, cref, tols) - && !intmat_is_identity(tr) ) - { - complain(cell, cref, "NOT"); - return 1; - } + if ( check_ccp(cell, cref, tols, intmat_is_identity(tr)) ) return 1; + if ( check_cpcp(cell, cref, tols, 1) ) return 1; + if ( check_ccpao(cell, cref, tols, intmat_is_identity(tr)) ) return 1; + if ( check_cpcpao(cell, cref, tols, 1) ) return 1; + if ( check_crcp(cell, cref, tols, 1) ) return 1; cell_free(cell); intmat_free(tr); - intmat_free(m); } /* Rotate cell and permute axes */ STATUS("Testing rotation with axis permutation...\n"); for ( i=0; i<100; i++ ) { - IntegerMatrix *m; IntegerMatrix *tr; UnitCell *cell2; - cell = cell_rotate(cref, random_quaternion(rng)); - if ( cell == NULL ) return 1; + cell2 = cell_rotate(cref, random_quaternion(rng)); + if ( cell2 == NULL ) return 1; tr = random_permutation(rng); - cell2 = cell_transform_intmat(cell, tr); - - if ( !compare_permuted_cell_parameters(cell2, cref, tols, &m) ) { - complain(cell2, cref, ""); - return 1; - } - - if ( compare_permuted_cell_parameters_and_orientation(cell2, cref, - tols[0], tols[3], &m) ) - { - UnitCell *cc; - complain(cell2, cref, "NOT, with just permutation,"); - STATUS("Matrix was (det=%i):\n", intmat_det(m)); - intmat_print(m); - STATUS("Transformed version of cref:\n"); - cc = cell_transform_intmat(cref, m); - cell_print_full(cc); - cell_free(cc); - return 1; - } - - if ( compare_cell_parameters_and_orientation(cell2, cref, - tols[0], tols[3]) ) - { - complain(cell2, cref, "NOT, without any change,"); - return 1; - } + cell = cell_transform_intmat(cell2, tr); + cell_free(cell2); + + if ( check_ccp(cell, cref, tols, intmat_is_identity(tr)) ) return 1; + if ( check_cpcp(cell, cref, tols, 1) ) return 1; + if ( check_ccpao(cell, cref, tols, 0) ) return 1; + if ( check_cpcpao(cell, cref, tols, 0) ) return 1; + if ( check_crcp(cell, cref, tols, 1) ) return 1; cell_free(cell); - cell_free(cell2); intmat_free(tr); - intmat_free(m); + } + + /* Reindex */ + STATUS("Testing reindexing...\n"); + for ( i=0; i<100; i++ ) { + + RationalMatrix *tr; + + tr = random_reindexing(rng); + cell = cell_transform_rational(cref, tr); + + if ( check_ccp(cell, cref, tols, rtnl_mtx_is_identity(tr)) ) return 1; + if ( check_cpcp(cell, cref, tols, rtnl_mtx_is_perm(tr)) ) return 1; + if ( check_ccpao(cell, cref, tols, rtnl_mtx_is_identity(tr)) ) return 1; + if ( check_cpcpao(cell, cref, tols, rtnl_mtx_is_perm(tr)) ) return 1; + if ( check_crcp(cell, cref, tols, 1) ) return 1; + + cell_free(cell); + rtnl_mtx_free(tr); } /* Reindex and rotate */ - STATUS("Testing rotation with reindexing...\n"); + STATUS("Testing reindexing with rotation...\n"); for ( i=0; i<100; i++ ) { - RationalMatrix *m; RationalMatrix *tr; UnitCell *cell2; - cell = cell_rotate(cref, random_quaternion(rng)); - if ( cell == NULL ) return 1; + cell2 = cell_rotate(cref, random_quaternion(rng)); + if ( cell2 == NULL ) return 1; tr = random_reindexing(rng); - cell2 = cell_transform_rational(cell, tr); - cell_free(cell); - - if ( !compare_reindexed_cell_parameters(cell2, cref, tols, 0, &m) ) { - complain(cell2, cref, ""); - return 1; - } + cell = cell_transform_rational(cell2, tr); cell_free(cell2); + + if ( check_ccp(cell, cref, tols, rtnl_mtx_is_identity(tr)) ) return 1; + if ( check_cpcp(cell, cref, tols, rtnl_mtx_is_perm(tr)) ) return 1; + if ( check_ccpao(cell, cref, tols, 0) ) return 1; + if ( check_cpcpao(cell, cref, tols, 0) ) return 1; + if ( check_crcp(cell, cref, tols, 1) ) return 1; + + cell_free(cell); rtnl_mtx_free(tr); - rtnl_mtx_free(m); } /* NB There's no compare_reindexed_cell_parameters_and_orientation */ -- cgit v1.2.3 From 41d5529aa1a6ef6c29c59d4c39767de61c57df12 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Thu, 8 Aug 2019 15:22:55 +0200 Subject: Implement random reindexing matrix --- tests/cellcompare_check.c | 42 ++++++++++++++++++++++++++++++++---------- 1 file changed, 32 insertions(+), 10 deletions(-) diff --git a/tests/cellcompare_check.c b/tests/cellcompare_check.c index 234d93a8..b3810378 100644 --- a/tests/cellcompare_check.c +++ b/tests/cellcompare_check.c @@ -53,17 +53,27 @@ static void complain(UnitCell *cell, UnitCell *cref, const char *t, const char * static RationalMatrix *random_reindexing(gsl_rng *rng) { - int i; + int i, j; RationalMatrix *tr; - int v[] = {0, 1, 2}; tr = rtnl_mtx_new(3, 3); if ( tr == NULL ) return NULL; - gsl_ran_shuffle(rng, v, 3, sizeof(int)); - for ( i=0; i<3; i++ ) { - rtnl_mtx_set(tr, i, v[i], rtnl(1, 1)); - } + do { + for ( i=0; i<3; i++ ) { + for ( j=0; j<3; j++ ) { + /* 0..8 inclusive -> -4..4 inclusive */ + signed int a = gsl_rng_uniform_int(rng, 9) - 4; + signed int b = gsl_rng_uniform_int(rng, 9) - 4; + if ( b == 0 ) { + a = 0; + b = 1; + } + rtnl_mtx_set(tr, i, j, rtnl(a, b)); + } + } + + } while ( rtnl_cmp(rtnl_mtx_det(tr), rtnl_zero()) == 0 ); return tr; } @@ -271,8 +281,14 @@ int main(int argc, char *argv[]) RationalMatrix *tr; - tr = random_reindexing(rng); - cell = cell_transform_rational(cref, tr); + cell = NULL; + tr = NULL; + do { + cell_free(cell); + rtnl_mtx_free(tr); + tr = random_reindexing(rng); + cell = cell_transform_rational(cref, tr); + } while ( cell_get_centering(cell) == '?' ); if ( check_ccp(cell, cref, tols, rtnl_mtx_is_identity(tr)) ) return 1; if ( check_cpcp(cell, cref, tols, rtnl_mtx_is_perm(tr)) ) return 1; @@ -294,8 +310,14 @@ int main(int argc, char *argv[]) cell2 = cell_rotate(cref, random_quaternion(rng)); if ( cell2 == NULL ) return 1; - tr = random_reindexing(rng); - cell = cell_transform_rational(cell2, tr); + cell = NULL; + tr = NULL; + do { + cell_free(cell); + rtnl_mtx_free(tr); + tr = random_reindexing(rng); + cell = cell_transform_rational(cell2, tr); + } while ( cell_get_centering(cell) == '?' ); cell_free(cell2); if ( check_ccp(cell, cref, tols, rtnl_mtx_is_identity(tr)) ) return 1; -- cgit v1.2.3 From 88b2fdfb4a43c6cb9104b9d9dd9b63d422846d0a Mon Sep 17 00:00:00 2001 From: Thomas White Date: Fri, 9 Aug 2019 13:59:47 +0200 Subject: Pass transformation matrix for debugging --- tests/cellcompare_check.c | 19 ++++++++++++------- 1 file changed, 12 insertions(+), 7 deletions(-) diff --git a/tests/cellcompare_check.c b/tests/cellcompare_check.c index b3810378..44ba9307 100644 --- a/tests/cellcompare_check.c +++ b/tests/cellcompare_check.c @@ -180,7 +180,7 @@ static int check_cpcpao(UnitCell *cell, UnitCell *cref, double *tols, static int check_crcp(UnitCell *cell, UnitCell *cref, double *tols, - int should_match) + RationalMatrix *tr, int should_match) { RationalMatrix *m = NULL; const char *a; @@ -189,9 +189,14 @@ static int check_crcp(UnitCell *cell, UnitCell *cref, double *tols, a = should_match ? "" : "NOT "; b = " with compare_reindexed_cell_parameters"; - if ( compare_reindexed_cell_parameters(cell, cref, tols, 0, &m) != should_match ) + if ( compare_reindexed_cell_parameters(cref, cell, tols, 1, &m) != should_match ) { complain(cell, cref, a, b); + STATUS("The transformation matrix to create the cell was:\n"); + rtnl_mtx_print(tr); + STATUS("Cell comparison says do this to the reference to " + "create the cell:\n"); + rtnl_mtx_print(m); rtnl_mtx_free(m); return 1; } @@ -227,7 +232,7 @@ int main(int argc, char *argv[]) if ( check_cpcp(cell, cref, tols, 1) ) return 1; if ( check_ccpao(cell, cref, tols, 0) ) return 1; if ( check_cpcpao(cell, cref, tols, 0) ) return 1; - if ( check_crcp(cell, cref, tols, 1) ) return 1; + if ( check_crcp(cell, cref, tols, NULL, 1) ) return 1; cell_free(cell); } @@ -245,7 +250,7 @@ int main(int argc, char *argv[]) if ( check_cpcp(cell, cref, tols, 1) ) return 1; if ( check_ccpao(cell, cref, tols, intmat_is_identity(tr)) ) return 1; if ( check_cpcpao(cell, cref, tols, 1) ) return 1; - if ( check_crcp(cell, cref, tols, 1) ) return 1; + if ( check_crcp(cell, cref, tols, NULL, 1) ) return 1; cell_free(cell); intmat_free(tr); @@ -269,7 +274,7 @@ int main(int argc, char *argv[]) if ( check_cpcp(cell, cref, tols, 1) ) return 1; if ( check_ccpao(cell, cref, tols, 0) ) return 1; if ( check_cpcpao(cell, cref, tols, 0) ) return 1; - if ( check_crcp(cell, cref, tols, 1) ) return 1; + if ( check_crcp(cell, cref, tols, NULL, 1) ) return 1; cell_free(cell); intmat_free(tr); @@ -294,7 +299,7 @@ int main(int argc, char *argv[]) if ( check_cpcp(cell, cref, tols, rtnl_mtx_is_perm(tr)) ) return 1; if ( check_ccpao(cell, cref, tols, rtnl_mtx_is_identity(tr)) ) return 1; if ( check_cpcpao(cell, cref, tols, rtnl_mtx_is_perm(tr)) ) return 1; - if ( check_crcp(cell, cref, tols, 1) ) return 1; + if ( check_crcp(cell, cref, tols, tr, 1) ) return 1; cell_free(cell); rtnl_mtx_free(tr); @@ -324,7 +329,7 @@ int main(int argc, char *argv[]) if ( check_cpcp(cell, cref, tols, rtnl_mtx_is_perm(tr)) ) return 1; if ( check_ccpao(cell, cref, tols, 0) ) return 1; if ( check_cpcpao(cell, cref, tols, 0) ) return 1; - if ( check_crcp(cell, cref, tols, 1) ) return 1; + if ( check_crcp(cell, cref, tols, tr, 1) ) return 1; cell_free(cell); rtnl_mtx_free(tr); -- cgit v1.2.3 From 81d63f562ae5fad3245b565743b2de8e6cc88048 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Fri, 9 Aug 2019 14:00:06 +0200 Subject: Set centering on input cell --- tests/cellcompare_check.c | 1 + 1 file changed, 1 insertion(+) diff --git a/tests/cellcompare_check.c b/tests/cellcompare_check.c index 44ba9307..472f3a25 100644 --- a/tests/cellcompare_check.c +++ b/tests/cellcompare_check.c @@ -219,6 +219,7 @@ int main(int argc, char *argv[]) deg2rad(67.0), deg2rad(70.0), deg2rad(77.0)); + cell_set_centering(cref, 'P'); if ( cref == NULL ) return 1; /* Just rotate cell */ -- cgit v1.2.3 From 93c6569621084680dc94e5395b39f6308ae9d11a Mon Sep 17 00:00:00 2001 From: Thomas White Date: Fri, 9 Aug 2019 14:00:27 +0200 Subject: Reduce range of reindexing matrices --- tests/cellcompare_check.c | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/tests/cellcompare_check.c b/tests/cellcompare_check.c index 472f3a25..962566a9 100644 --- a/tests/cellcompare_check.c +++ b/tests/cellcompare_check.c @@ -62,9 +62,10 @@ static RationalMatrix *random_reindexing(gsl_rng *rng) do { for ( i=0; i<3; i++ ) { for ( j=0; j<3; j++ ) { - /* 0..8 inclusive -> -4..4 inclusive */ - signed int a = gsl_rng_uniform_int(rng, 9) - 4; - signed int b = gsl_rng_uniform_int(rng, 9) - 4; + /* 0..6 inclusive -> -3..3 inclusive */ + signed int a = gsl_rng_uniform_int(rng, 7) - 3; + /* 0..2 inclusive */ + signed int b = gsl_rng_uniform_int(rng, 3); if ( b == 0 ) { a = 0; b = 1; -- cgit v1.2.3 From 7a81b3c0557f97f81fec83afd95802d1e62c1a43 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Mon, 12 Aug 2019 11:21:51 +0200 Subject: Don't use H centering for cell comparison check --- tests/cellcompare_check.c | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/tests/cellcompare_check.c b/tests/cellcompare_check.c index 962566a9..4b9ba869 100644 --- a/tests/cellcompare_check.c +++ b/tests/cellcompare_check.c @@ -295,7 +295,11 @@ int main(int argc, char *argv[]) rtnl_mtx_free(tr); tr = random_reindexing(rng); cell = cell_transform_rational(cref, tr); - } while ( cell_get_centering(cell) == '?' ); + } while ( (cell_get_centering(cell) == '?') + || (cell_get_centering(cell) == 'H' ) ); + /* H centering is no good because it needs a unique axis to + * be specified in order for uncentering in c_r_c_p to work. + * cell_transform_rational doesn't set the unique axis (yet?) */ if ( check_ccp(cell, cref, tols, rtnl_mtx_is_identity(tr)) ) return 1; if ( check_cpcp(cell, cref, tols, rtnl_mtx_is_perm(tr)) ) return 1; @@ -324,7 +328,8 @@ int main(int argc, char *argv[]) rtnl_mtx_free(tr); tr = random_reindexing(rng); cell = cell_transform_rational(cell2, tr); - } while ( cell_get_centering(cell) == '?' ); + } while ( (cell_get_centering(cell) == '?') + || (cell_get_centering(cell) == 'H' ) ); /* See above */ cell_free(cell2); if ( check_ccp(cell, cref, tols, rtnl_mtx_is_identity(tr)) ) return 1; -- cgit v1.2.3 From 9a510740e673a9b9ce82e08c2e2ccb4c3a36ebfd Mon Sep 17 00:00:00 2001 From: Thomas White Date: Mon, 12 Aug 2019 11:34:36 +0200 Subject: cellcompare_check: Progress bars --- tests/cellcompare_check.c | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/tests/cellcompare_check.c b/tests/cellcompare_check.c index 4b9ba869..48b7a9c5 100644 --- a/tests/cellcompare_check.c +++ b/tests/cellcompare_check.c @@ -224,7 +224,6 @@ int main(int argc, char *argv[]) if ( cref == NULL ) return 1; /* Just rotate cell */ - STATUS("Testing plain rotation...\n"); for ( i=0; i<100; i++ ) { cell = cell_rotate(cref, random_quaternion(rng)); @@ -237,10 +236,10 @@ int main(int argc, char *argv[]) if ( check_crcp(cell, cref, tols, NULL, 1) ) return 1; cell_free(cell); + progress_bar(i+1, 100, "Plain rotation"); } /* Permute axes but don't rotate */ - STATUS("Testing axis permutation...\n"); for ( i=0; i<100; i++ ) { IntegerMatrix *tr; @@ -256,10 +255,10 @@ int main(int argc, char *argv[]) cell_free(cell); intmat_free(tr); + progress_bar(i+1, 100, "Axis permutation"); } /* Rotate cell and permute axes */ - STATUS("Testing rotation with axis permutation...\n"); for ( i=0; i<100; i++ ) { IntegerMatrix *tr; @@ -280,10 +279,10 @@ int main(int argc, char *argv[]) cell_free(cell); intmat_free(tr); + progress_bar(i+1, 100, "Rotation with axis permutation"); } /* Reindex */ - STATUS("Testing reindexing...\n"); for ( i=0; i<100; i++ ) { RationalMatrix *tr; @@ -309,10 +308,10 @@ int main(int argc, char *argv[]) cell_free(cell); rtnl_mtx_free(tr); + progress_bar(i+1, 100, "Reindexing"); } /* Reindex and rotate */ - STATUS("Testing reindexing with rotation...\n"); for ( i=0; i<100; i++ ) { RationalMatrix *tr; @@ -340,6 +339,7 @@ int main(int argc, char *argv[]) cell_free(cell); rtnl_mtx_free(tr); + progress_bar(i+1, 100, "Reindexing with rotation"); } /* NB There's no compare_reindexed_cell_parameters_and_orientation */ -- 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(-) 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(+) 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 +------ tests/cellcompare_check.c | 79 +++++++++++++++++++------------- 4 files changed, 87 insertions(+), 119 deletions(-) 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; - } diff --git a/tests/cellcompare_check.c b/tests/cellcompare_check.c index 48b7a9c5..0e98276c 100644 --- a/tests/cellcompare_check.c +++ b/tests/cellcompare_check.c @@ -116,29 +116,6 @@ static int check_ccp(UnitCell *cell, UnitCell *cref, double *tols, } -static int check_cpcp(UnitCell *cell, UnitCell *cref, double *tols, - int should_match) -{ - IntegerMatrix *m = NULL; - const char *a; - const char *b; - - a = should_match ? "" : "NOT "; - b = " with compare_permuted_cell_parameters"; - - if ( compare_permuted_cell_parameters(cell, cref, tols, &m) != should_match ) - { - complain(cell, cref, a, b); - STATUS("Matrix was:\n"); - intmat_print(m); - intmat_free(m); - return 1; - } - intmat_free(m); - return 0; -} - - static int check_ccpao(UnitCell *cell, UnitCell *cref, double *tols, int should_match) { @@ -206,6 +183,33 @@ static int check_crcp(UnitCell *cell, UnitCell *cref, double *tols, } +static void yaro_test() +{ + UnitCell *cell; + UnitCell *reference; + UnitCell *cmatch; + float tols[] = {5, 5, 5, 1.5}; + + cell = cell_new_from_parameters(64.24e-10, 63.94e-10, 64.61e-10, + deg2rad(89.98), deg2rad(89.82), deg2rad(119.87)); + cell_set_unique_axis(cell, 'c'); + cell_set_lattice_type(cell, L_HEXAGONAL); + + reference = cell_new_from_parameters(64.7e-10, 64.7e-10, 65.2e-10, + deg2rad(90.0), deg2rad(90.0), deg2rad(120.0)); + cell_set_unique_axis(reference, 'c'); + cell_set_lattice_type(reference, L_HEXAGONAL); + + cell_print(cell); + cell_print(reference); + cmatch = match_cell(cell, reference, 0, tols, 1); + cell_print(cmatch); +} + + +extern IntegerMatrix *reduce_g6(double *g, double eps); +extern void g6_components(double *g6, UnitCell *cell); + int main(int argc, char *argv[]) { UnitCell *cell, *cref; @@ -216,13 +220,29 @@ int main(int argc, char *argv[]) rng = gsl_rng_alloc(gsl_rng_mt19937); - cref = cell_new_from_parameters(50e-10, 55e-10, 70e-10, - deg2rad(67.0), - deg2rad(70.0), - deg2rad(77.0)); + yaro_test(); + + cref = cell_new_from_parameters(3e-0, 5.196e-0, 2e-0, + deg2rad(103.9166666), + deg2rad(109.4666666), + deg2rad(134.8833333)); cell_set_centering(cref, 'P'); if ( cref == NULL ) return 1; + STATUS("The test cell:\n"); + cell_print(cref); + double g[6]; + g6_components(g, cref); + STATUS("G6: %e %e %e %e %e %e\n", g[0], g[1], g[2], g[3], g[4], g[5]); + double eps = pow(cell_get_volume(cref), 1.0/3.0) * 1e-5; + //eps = 1e-27; + IntegerMatrix *M = reduce_g6(g, eps); + STATUS("The transformation to reduce:\n"); + intmat_print(M); + STATUS("The reduced cell:\n"); + cell = cell_transform_intmat(cref, M); + cell_print(cell); + /* Just rotate cell */ for ( i=0; i<100; i++ ) { @@ -230,7 +250,6 @@ int main(int argc, char *argv[]) if ( cell == NULL ) return 1; if ( check_ccp(cell, cref, tols, 1) ) return 1; - if ( check_cpcp(cell, cref, tols, 1) ) return 1; if ( check_ccpao(cell, cref, tols, 0) ) return 1; if ( check_cpcpao(cell, cref, tols, 0) ) return 1; if ( check_crcp(cell, cref, tols, NULL, 1) ) return 1; @@ -248,7 +267,6 @@ int main(int argc, char *argv[]) cell = cell_transform_intmat(cref, tr); if ( check_ccp(cell, cref, tols, intmat_is_identity(tr)) ) return 1; - if ( check_cpcp(cell, cref, tols, 1) ) return 1; if ( check_ccpao(cell, cref, tols, intmat_is_identity(tr)) ) return 1; if ( check_cpcpao(cell, cref, tols, 1) ) return 1; if ( check_crcp(cell, cref, tols, NULL, 1) ) return 1; @@ -272,7 +290,6 @@ int main(int argc, char *argv[]) cell_free(cell2); if ( check_ccp(cell, cref, tols, intmat_is_identity(tr)) ) return 1; - if ( check_cpcp(cell, cref, tols, 1) ) return 1; if ( check_ccpao(cell, cref, tols, 0) ) return 1; if ( check_cpcpao(cell, cref, tols, 0) ) return 1; if ( check_crcp(cell, cref, tols, NULL, 1) ) return 1; @@ -301,7 +318,6 @@ int main(int argc, char *argv[]) * cell_transform_rational doesn't set the unique axis (yet?) */ if ( check_ccp(cell, cref, tols, rtnl_mtx_is_identity(tr)) ) return 1; - if ( check_cpcp(cell, cref, tols, rtnl_mtx_is_perm(tr)) ) return 1; if ( check_ccpao(cell, cref, tols, rtnl_mtx_is_identity(tr)) ) return 1; if ( check_cpcpao(cell, cref, tols, rtnl_mtx_is_perm(tr)) ) return 1; if ( check_crcp(cell, cref, tols, tr, 1) ) return 1; @@ -332,7 +348,6 @@ int main(int argc, char *argv[]) cell_free(cell2); if ( check_ccp(cell, cref, tols, rtnl_mtx_is_identity(tr)) ) return 1; - if ( check_cpcp(cell, cref, tols, rtnl_mtx_is_perm(tr)) ) return 1; if ( check_ccpao(cell, cref, tols, 0) ) return 1; if ( check_cpcpao(cell, cref, tols, 0) ) return 1; if ( check_crcp(cell, cref, tols, tr, 1) ) 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(+) 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 ++++++++++++++++++++++++++----------------- tests/cellcompare_check.c | 15 +- 2 files changed, 232 insertions(+), 150 deletions(-) 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) */ diff --git a/tests/cellcompare_check.c b/tests/cellcompare_check.c index 0e98276c..d0f51da5 100644 --- a/tests/cellcompare_check.c +++ b/tests/cellcompare_check.c @@ -207,8 +207,7 @@ static void yaro_test() } -extern IntegerMatrix *reduce_g6(double *g, double eps); -extern void g6_components(double *g6, UnitCell *cell); +extern IntegerMatrix *reduce_g6(struct g6 g, double eps); int main(int argc, char *argv[]) { @@ -222,7 +221,7 @@ int main(int argc, char *argv[]) yaro_test(); - cref = cell_new_from_parameters(3e-0, 5.196e-0, 2e-0, + cref = cell_new_from_parameters(3e-10, 5.196e-10, 2e-10, deg2rad(103.9166666), deg2rad(109.4666666), deg2rad(134.8833333)); @@ -231,11 +230,13 @@ int main(int argc, char *argv[]) STATUS("The test cell:\n"); cell_print(cref); - double g[6]; - g6_components(g, cref); - STATUS("G6: %e %e %e %e %e %e\n", g[0], g[1], g[2], g[3], g[4], g[5]); + struct g6 g; + g = cell_get_G6(cref); double eps = pow(cell_get_volume(cref), 1.0/3.0) * 1e-5; - //eps = 1e-27; + eps = eps*eps; + //eps *= 100; + //g.A = 9.0e-20; g.B = 27.0e-20; g.C = 4.0e-20; + //g.D = -5.0e-20; g.E = -4.0e-20; g.F = -22.0e-20; IntegerMatrix *M = reduce_g6(g, eps); STATUS("The transformation to reduce:\n"); intmat_print(M); -- 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 +++++++++++++++----- tests/cellcompare_check.c | 11 ++++------- 2 files changed, 19 insertions(+), 12 deletions(-) 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) */ diff --git a/tests/cellcompare_check.c b/tests/cellcompare_check.c index d0f51da5..f5787753 100644 --- a/tests/cellcompare_check.c +++ b/tests/cellcompare_check.c @@ -207,7 +207,7 @@ static void yaro_test() } -extern IntegerMatrix *reduce_g6(struct g6 g, double eps); +extern IntegerMatrix *reduce_g6(struct g6 g, double epsrel); int main(int argc, char *argv[]) { @@ -232,12 +232,9 @@ int main(int argc, char *argv[]) cell_print(cref); struct g6 g; g = cell_get_G6(cref); - double eps = pow(cell_get_volume(cref), 1.0/3.0) * 1e-5; - eps = eps*eps; - //eps *= 100; - //g.A = 9.0e-20; g.B = 27.0e-20; g.C = 4.0e-20; - //g.D = -5.0e-20; g.E = -4.0e-20; g.F = -22.0e-20; - IntegerMatrix *M = reduce_g6(g, eps); + g.A = 9.0e-20; g.B = 27.0e-20; g.C = 4.0e-20; + g.D = -5.0e-20; g.E = -4.0e-20; g.F = -22.0e-20; + IntegerMatrix *M = reduce_g6(g, 1e-5); STATUS("The transformation to reduce:\n"); intmat_print(M); STATUS("The reduced cell:\n"); -- 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(-) 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 +++- src/cell_tool.c | 2 +- src/whirligig.c | 8 ++- tests/cellcompare_check.c | 46 ++++++++++++++--- 6 files changed, 134 insertions(+), 71 deletions(-) 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); } diff --git a/src/cell_tool.c b/src/cell_tool.c index a1755f42..0942ca15 100644 --- a/src/cell_tool.c +++ b/src/cell_tool.c @@ -99,7 +99,7 @@ static int comparecells(UnitCell *cell, const char *comparecell, tolerance[5] = atl; STATUS("------------------> The comparison results:\n"); - if ( !compare_reindexed_cell_parameters(cell, cell2, tolerance, csl, &m) ) { + if ( !compare_derivative_cell_parameters(cell, cell2, tolerance, csl, &m) ) { STATUS("No relationship found between lattices.\n"); return 0; } else { diff --git a/src/whirligig.c b/src/whirligig.c index a144e769..fd94a0f2 100644 --- a/src/whirligig.c +++ b/src/whirligig.c @@ -297,6 +297,8 @@ static IntegerMatrix *try_all(struct window *win, int n1, int n2, IntegerMatrix *m; struct image *i1; struct image *i2; + const double tols[] = {0.1, 0.1, 0.1, + deg2rad(5.0), deg2rad(5.0), deg2rad(5.0)}; assert(n1 >= 0); assert(n2 >= 0); @@ -311,7 +313,7 @@ static IntegerMatrix *try_all(struct window *win, int n1, int n2, if ( compare_permuted_cell_parameters_and_orientation(crystal_get_cell(i1->crystals[i]), crystal_get_cell(i2->crystals[j]), - 0.1, deg2rad(5.0), &m) ) + tols, &m) ) { if ( !crystal_used(win, n1, i) && !crystal_used(win, n2, j) ) @@ -365,6 +367,8 @@ static int try_join(struct window *win, int sn) Crystal *cr; UnitCell *ref; const int sp = win->join_ptr - 1; + const double tols[] = {0.1, 0.1, 0.1, + deg2rad(5.0), deg2rad(5.0), deg2rad(5.0)}; /* Get the appropriately transformed cell from the last crystal in this * series */ @@ -375,7 +379,7 @@ static int try_join(struct window *win, int sn) Crystal *cr2; cr2 = win->img[win->join_ptr].crystals[j]; if ( compare_permuted_cell_parameters_and_orientation(ref, crystal_get_cell(cr2), - 0.1, deg2rad(5.0), + tols, &win->mat[sn][win->join_ptr]) ) { win->ser[sn][win->join_ptr] = j; diff --git a/tests/cellcompare_check.c b/tests/cellcompare_check.c index f5787753..0c3709ba 100644 --- a/tests/cellcompare_check.c +++ b/tests/cellcompare_check.c @@ -126,7 +126,7 @@ static int check_ccpao(UnitCell *cell, UnitCell *cref, double *tols, b = " with compare_cell_parameters_and_orientation"; if ( compare_cell_parameters_and_orientation(cell, cref, - tols[0], tols[3]) != should_match ) + tols) != should_match ) { complain(cell, cref, a, b); return 1; @@ -145,8 +145,8 @@ static int check_cpcpao(UnitCell *cell, UnitCell *cref, double *tols, a = should_match ? "" : "NOT "; b = " with compare_permuted_cell_parameters_and_orientation"; - if ( compare_permuted_cell_parameters_and_orientation(cell, cref, - tols[0], tols[3], &m) != should_match ) + if ( compare_permuted_cell_parameters_and_orientation(cell, cref, tols, &m) + != should_match ) { complain(cell, cref, a, b); intmat_free(m); @@ -157,17 +157,46 @@ static int check_cpcpao(UnitCell *cell, UnitCell *cref, double *tols, } +static int check_cdcp(UnitCell *cell, UnitCell *cref, double *tols, + RationalMatrix *tr, int should_match) +{ + RationalMatrix *m = NULL; + const char *a; + const char *b; + + a = should_match ? "" : "NOT "; + b = " with compare_derivative_cell_parameters"; + + if ( compare_derivative_cell_parameters(cref, cell, tols, 1, &m) != should_match ) + { + complain(cell, cref, a, b); + STATUS("The transformation matrix to create the cell was:\n"); + rtnl_mtx_print(tr); + STATUS("Cell comparison says do this to the reference to " + "create the cell:\n"); + rtnl_mtx_print(m); + rtnl_mtx_free(m); + return 1; + } + rtnl_mtx_free(m); + return 0; +} + + static int check_crcp(UnitCell *cell, UnitCell *cref, double *tols, RationalMatrix *tr, int should_match) { RationalMatrix *m = NULL; const char *a; const char *b; + UnitCell *match; a = should_match ? "" : "NOT "; b = " with compare_reindexed_cell_parameters"; - if ( compare_reindexed_cell_parameters(cref, cell, tols, 1, &m) != should_match ) + match = compare_reindexed_cell_parameters(cell, cref, tols, &m); + + if ( (match != NULL) != should_match ) { complain(cell, cref, a, b); STATUS("The transformation matrix to create the cell was:\n"); @@ -250,6 +279,7 @@ int main(int argc, char *argv[]) if ( check_ccp(cell, cref, tols, 1) ) return 1; if ( check_ccpao(cell, cref, tols, 0) ) return 1; if ( check_cpcpao(cell, cref, tols, 0) ) return 1; + if ( check_cdcp(cell, cref, tols, NULL, 1) ) return 1; if ( check_crcp(cell, cref, tols, NULL, 1) ) return 1; cell_free(cell); @@ -267,6 +297,7 @@ int main(int argc, char *argv[]) if ( check_ccp(cell, cref, tols, intmat_is_identity(tr)) ) return 1; if ( check_ccpao(cell, cref, tols, intmat_is_identity(tr)) ) return 1; if ( check_cpcpao(cell, cref, tols, 1) ) return 1; + if ( check_cdcp(cell, cref, tols, NULL, 1) ) return 1; if ( check_crcp(cell, cref, tols, NULL, 1) ) return 1; cell_free(cell); @@ -290,6 +321,7 @@ int main(int argc, char *argv[]) if ( check_ccp(cell, cref, tols, intmat_is_identity(tr)) ) return 1; if ( check_ccpao(cell, cref, tols, 0) ) return 1; if ( check_cpcpao(cell, cref, tols, 0) ) return 1; + if ( check_cdcp(cell, cref, tols, NULL, 1) ) return 1; if ( check_crcp(cell, cref, tols, NULL, 1) ) return 1; cell_free(cell); @@ -318,7 +350,8 @@ int main(int argc, char *argv[]) if ( check_ccp(cell, cref, tols, rtnl_mtx_is_identity(tr)) ) return 1; if ( check_ccpao(cell, cref, tols, rtnl_mtx_is_identity(tr)) ) return 1; if ( check_cpcpao(cell, cref, tols, rtnl_mtx_is_perm(tr)) ) return 1; - if ( check_crcp(cell, cref, tols, tr, 1) ) return 1; + if ( check_cdcp(cell, cref, tols, tr, 1) ) return 1; + /* check_crcp: Sometimes yes, hard to tell */ cell_free(cell); rtnl_mtx_free(tr); @@ -348,7 +381,8 @@ int main(int argc, char *argv[]) if ( check_ccp(cell, cref, tols, rtnl_mtx_is_identity(tr)) ) return 1; if ( check_ccpao(cell, cref, tols, 0) ) return 1; if ( check_cpcpao(cell, cref, tols, 0) ) return 1; - if ( check_crcp(cell, cref, tols, tr, 1) ) return 1; + if ( check_cdcp(cell, cref, tols, tr, 1) ) return 1; + /* check_crcp: Sometimes yes, hard to tell */ cell_free(cell); rtnl_mtx_free(tr); -- 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 +++++++++++++++++++++++++++++++++---------- tests/cellcompare_check.c | 36 ++--- 2 files changed, 258 insertions(+), 87 deletions(-) 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; } diff --git a/tests/cellcompare_check.c b/tests/cellcompare_check.c index 0c3709ba..4a7d4913 100644 --- a/tests/cellcompare_check.c +++ b/tests/cellcompare_check.c @@ -218,21 +218,40 @@ static void yaro_test() UnitCell *reference; UnitCell *cmatch; float tols[] = {5, 5, 5, 1.5}; + double dtols[] = {0.05, 0.05, 0.05, deg2rad(5.0), deg2rad(5.0), deg2rad(5.0)}; - cell = cell_new_from_parameters(64.24e-10, 63.94e-10, 64.61e-10, + cell = cell_new_from_parameters(63.24e-10, 63.94e-10, 64.61e-10, deg2rad(89.98), deg2rad(89.82), deg2rad(119.87)); cell_set_unique_axis(cell, 'c'); cell_set_lattice_type(cell, L_HEXAGONAL); + cell_set_centering(cell, 'P'); reference = cell_new_from_parameters(64.7e-10, 64.7e-10, 65.2e-10, deg2rad(90.0), deg2rad(90.0), deg2rad(120.0)); cell_set_unique_axis(reference, 'c'); cell_set_lattice_type(reference, L_HEXAGONAL); + cell_set_centering(reference, 'P'); + STATUS("The cell:\n"); cell_print(cell); + STATUS("The reference:\n"); cell_print(reference); cmatch = match_cell(cell, reference, 0, tols, 1); + STATUS("The match:\n"); cell_print(cmatch); + cell_free(cmatch); + + RationalMatrix *m = NULL; + cmatch = compare_reindexed_cell_parameters(cell, reference, dtols, &m); + STATUS("The new match:\n") + cell_print(cmatch); + STATUS("The matrix:\n"); + rtnl_mtx_print(m); + cell_free(cmatch); + rtnl_mtx_free(m); + + cell_free(cell); + cell_free(reference); } @@ -257,19 +276,6 @@ int main(int argc, char *argv[]) cell_set_centering(cref, 'P'); if ( cref == NULL ) return 1; - STATUS("The test cell:\n"); - cell_print(cref); - struct g6 g; - g = cell_get_G6(cref); - g.A = 9.0e-20; g.B = 27.0e-20; g.C = 4.0e-20; - g.D = -5.0e-20; g.E = -4.0e-20; g.F = -22.0e-20; - IntegerMatrix *M = reduce_g6(g, 1e-5); - STATUS("The transformation to reduce:\n"); - intmat_print(M); - STATUS("The reduced cell:\n"); - cell = cell_transform_intmat(cref, M); - cell_print(cell); - /* Just rotate cell */ for ( i=0; i<100; i++ ) { @@ -391,8 +397,6 @@ int main(int argc, char *argv[]) /* NB There's no compare_reindexed_cell_parameters_and_orientation */ - /* Test using vectors */ - cell_free(cref); gsl_rng_free(rng); -- 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(-) 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 08a8430379c52cedb4e026285aa73dd9547be17e Mon Sep 17 00:00:00 2001 From: Thomas White Date: Thu, 22 Aug 2019 13:13:08 +0200 Subject: cellcompare_check: Variable for number of trials --- tests/cellcompare_check.c | 21 +++++++++++---------- 1 file changed, 11 insertions(+), 10 deletions(-) diff --git a/tests/cellcompare_check.c b/tests/cellcompare_check.c index 4a7d4913..107702c4 100644 --- a/tests/cellcompare_check.c +++ b/tests/cellcompare_check.c @@ -262,6 +262,7 @@ int main(int argc, char *argv[]) UnitCell *cell, *cref; gsl_rng *rng; int i; + const int ntrial = 100; double tols[] = { 0.01, 0.01, 0.01, deg2rad(1.0), deg2rad(1.0), deg2rad(1.0) }; @@ -277,7 +278,7 @@ int main(int argc, char *argv[]) if ( cref == NULL ) return 1; /* Just rotate cell */ - for ( i=0; i<100; i++ ) { + for ( i=0; i 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 - tests/cellcompare_check.c | 10 +- 3 files changed, 5 insertions(+), 408 deletions(-) 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); diff --git a/tests/cellcompare_check.c b/tests/cellcompare_check.c index 107702c4..7aa59fe4 100644 --- a/tests/cellcompare_check.c +++ b/tests/cellcompare_check.c @@ -217,7 +217,7 @@ static void yaro_test() UnitCell *cell; UnitCell *reference; UnitCell *cmatch; - float tols[] = {5, 5, 5, 1.5}; + //float tols[] = {5, 5, 5, 1.5}; double dtols[] = {0.05, 0.05, 0.05, deg2rad(5.0), deg2rad(5.0), deg2rad(5.0)}; cell = cell_new_from_parameters(63.24e-10, 63.94e-10, 64.61e-10, @@ -236,10 +236,10 @@ static void yaro_test() cell_print(cell); STATUS("The reference:\n"); cell_print(reference); - cmatch = match_cell(cell, reference, 0, tols, 1); - STATUS("The match:\n"); - cell_print(cmatch); - cell_free(cmatch); + //cmatch = match_cell(cell, reference, 0, tols, 1); + //STATUS("The match:\n"); + //cell_print(cmatch); + //cell_free(cmatch); RationalMatrix *m = NULL; cmatch = compare_reindexed_cell_parameters(cell, reference, dtols, &m); -- cgit v1.2.3 From 4efec1fa1dc08f4472531ef5b24de4d4e1b74aa9 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Thu, 22 Aug 2019 16:26:14 +0200 Subject: cellcompare_check: Fix terminology, and reduce number of trials --- tests/cellcompare_check.c | 18 ++++++++---------- 1 file changed, 8 insertions(+), 10 deletions(-) diff --git a/tests/cellcompare_check.c b/tests/cellcompare_check.c index 7aa59fe4..f674a7f2 100644 --- a/tests/cellcompare_check.c +++ b/tests/cellcompare_check.c @@ -51,7 +51,7 @@ static void complain(UnitCell *cell, UnitCell *cref, const char *t, const char * } -static RationalMatrix *random_reindexing(gsl_rng *rng) +static RationalMatrix *random_derivative(gsl_rng *rng) { int i, j; RationalMatrix *tr; @@ -262,7 +262,7 @@ int main(int argc, char *argv[]) UnitCell *cell, *cref; gsl_rng *rng; int i; - const int ntrial = 100; + const int ntrial = 10; double tols[] = { 0.01, 0.01, 0.01, deg2rad(1.0), deg2rad(1.0), deg2rad(1.0) }; @@ -336,7 +336,7 @@ int main(int argc, char *argv[]) progress_bar(i+1, ntrial, "Rotation with axis permutation"); } - /* Reindex */ + /* Derivative lattice */ for ( i=0; i 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(-) 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 f289a626ccc0e0fb74f28acd9bd714f44343eacb Mon Sep 17 00:00:00 2001 From: Thomas White Date: Mon, 26 Aug 2019 16:57:04 +0200 Subject: cell_tool: Just do all the comparisons, get rid of --csl --- src/cell_tool.c | 52 ++++++++++++++++++++++++++++++++++++++++++---------- 1 file changed, 42 insertions(+), 10 deletions(-) diff --git a/src/cell_tool.c b/src/cell_tool.c index 0942ca15..5fef5d9f 100644 --- a/src/cell_tool.c +++ b/src/cell_tool.c @@ -67,13 +67,12 @@ static void show_help(const char *s) " --tolerance= Set the tolerances for cell comparison.\n" " Default: 5,1.5 (axis percentage, angle deg).\n" " --highres=n Resolution limit (Angstroms) for --rings\n" -" --csl Allow --compare to find coincidence site lattice relationships.\n" ); } static int comparecells(UnitCell *cell, const char *comparecell, - double ltl, double atl, int csl) + double ltl, double atl) { UnitCell *cell2; RationalMatrix *m; @@ -98,10 +97,9 @@ static int comparecells(UnitCell *cell, const char *comparecell, tolerance[4] = atl; tolerance[5] = atl; - STATUS("------------------> The comparison results:\n"); - if ( !compare_derivative_cell_parameters(cell, cell2, tolerance, csl, &m) ) { + STATUS("------------------> Reindexed (strictly the same lattice):\n"); + if ( !compare_reindexed_cell_parameters(cell, cell2, tolerance, &m) ) { STATUS("No relationship found between lattices.\n"); - return 0; } else { UnitCell *trans; STATUS("Relationship found. To become similar to the reference" @@ -117,7 +115,44 @@ static int comparecells(UnitCell *cell, const char *comparecell, } - rtnl_mtx_free(m); + STATUS("------------------> Derivative lattice " + "(strictly the same lattice):\n"); + if ( !compare_derivative_cell_parameters(cell, cell2, tolerance, 0, &m) ) { + STATUS("No relationship found between lattices.\n"); + } else { + UnitCell *trans; + STATUS("Relationship found. To become similar to the reference" + " cell, the input cell should be transformed by:\n"); + rtnl_mtx_print(m); + STATUS("Transformed version of input unit cell:\n"); + trans = cell_transform_rational(cell, m); + cell_print(trans); + cell_free(trans); + STATUS("NB transformed cell might not really be triclinic, " + "it's just that I don't (yet) know how to work out what " + "it is.\n"); + rtnl_mtx_free(m); + } + + STATUS("------------------> Coincidence site lattice " + "(not strictly the same lattice):\n"); + if ( !compare_derivative_cell_parameters(cell, cell2, tolerance, 1, &m) ) { + STATUS("No relationship found between lattices.\n"); + return 0; + } else { + UnitCell *trans; + STATUS("Relationship found. To become similar to the reference" + " cell, the input cell should be transformed by:\n"); + rtnl_mtx_print(m); + STATUS("Transformed version of input unit cell:\n"); + trans = cell_transform_rational(cell, m); + cell_print(trans); + cell_free(trans); + STATUS("NB transformed cell might not really be triclinic, " + "it's just that I don't (yet) know how to work out what " + "it is.\n"); + rtnl_mtx_free(m); + } return 0; } @@ -442,7 +477,6 @@ int main(int argc, char *argv[]) float highres; double rmax = 1/(2.0e-10); char *trans_str = NULL; - int csl = 0; /* Long options */ const struct option longopts[] = { @@ -461,7 +495,6 @@ int main(int argc, char *argv[]) {"transform", 1, NULL, 4}, {"highres", 1, NULL, 5}, - {"csl", 0, &csl, 1}, {0, 0, NULL, 0} }; @@ -588,8 +621,7 @@ int main(int argc, char *argv[]) if ( mode == CT_FINDAMBI ) return find_ambi(cell, sym, ltl, atl); if ( mode == CT_UNCENTER ) return uncenter(cell, out_file); if ( mode == CT_RINGS ) return all_rings(cell, sym, rmax); - if ( mode == CT_COMPARE ) return comparecells(cell, comparecell, - ltl, atl, csl); + if ( mode == CT_COMPARE ) return comparecells(cell, comparecell, ltl, atl); if ( mode == CT_TRANSFORM ) return transform(cell, trans_str, out_file); if ( mode == CT_CHOICES ) return cell_choices(cell); -- cgit v1.2.3 From 71e7acb7668073b8f37565aa59bfd67d3b4cd166 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Tue, 27 Aug 2019 10:26:36 +0200 Subject: cell_tool: Clarify interpretation of tolerances --- doc/man/cell_tool.1 | 2 +- src/cell_tool.c | 3 +++ 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/doc/man/cell_tool.1 b/doc/man/cell_tool.1 index 576943d0..86dd9ae0 100644 --- a/doc/man/cell_tool.1 +++ b/doc/man/cell_tool.1 @@ -51,7 +51,7 @@ There are an infinite number of primitive unit cell for any lattice. This progr .PP The program will compare the two cells, and report if \fImy_structure.cell\fR can be made to look similar to \fIreference.cell\fR applying any transformation matrix. .PP -The tolerance \fItols\fR is given as lengthtol,angtol, in percent and degrees respectively, which will be applied to the real-space unit cell axis lengths and angles. +The tolerance \fItols\fR is given as lengthtol,angtol, in percent and degrees respectively, which will be applied to the real-space unit cell axis lengths and angles. If either of the unit cells are centered, a primitive version of will be used \fIfor some of the comparison results\fR. Read the output carefully. .SH TRANSFORMING A UNIT CELL .PP diff --git a/src/cell_tool.c b/src/cell_tool.c index 5fef5d9f..83056a6f 100644 --- a/src/cell_tool.c +++ b/src/cell_tool.c @@ -98,6 +98,7 @@ static int comparecells(UnitCell *cell, const char *comparecell, tolerance[5] = atl; STATUS("------------------> Reindexed (strictly the same lattice):\n"); + STATUS("Tolerances applied directly to the unit cells\n"); if ( !compare_reindexed_cell_parameters(cell, cell2, tolerance, &m) ) { STATUS("No relationship found between lattices.\n"); } else { @@ -117,6 +118,7 @@ static int comparecells(UnitCell *cell, const char *comparecell, STATUS("------------------> Derivative lattice " "(strictly the same lattice):\n"); + STATUS("Tolerances applied to primitive versions of the unit cells\n"); if ( !compare_derivative_cell_parameters(cell, cell2, tolerance, 0, &m) ) { STATUS("No relationship found between lattices.\n"); } else { @@ -136,6 +138,7 @@ static int comparecells(UnitCell *cell, const char *comparecell, STATUS("------------------> Coincidence site lattice " "(not strictly the same lattice):\n"); + STATUS("Tolerances applied to primitive versions of the unit cells\n"); if ( !compare_derivative_cell_parameters(cell, cell2, tolerance, 1, &m) ) { STATUS("No relationship found between lattices.\n"); return 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(-) 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 941889a4d9fc9663d70f80eb63c79c73eafb06c9 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Wed, 28 Aug 2019 13:30:11 +0200 Subject: cell_explorer: Add a warning --- src/cell_explorer.c | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/cell_explorer.c b/src/cell_explorer.c index 446aec5d..812e16d5 100644 --- a/src/cell_explorer.c +++ b/src/cell_explorer.c @@ -2040,6 +2040,9 @@ int main(int argc, char *argv[]) } w.cells[w.n_cells] = crystal_get_cell(cr); + if ( !right_handed(w.cells[w.n_cells]) ) { + ERROR("WARNING: Left-handed cell encountered\n"); + } w.indms[w.n_cells] = cur.indexed_by; w.n_cells++; -- 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 +++---- tests/transformation_check.c | 2 +- 7 files changed, 73 insertions(+), 65 deletions(-) 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; diff --git a/tests/transformation_check.c b/tests/transformation_check.c index 33f44340..8c8f7fa0 100644 --- a/tests/transformation_check.c +++ b/tests/transformation_check.c @@ -396,7 +396,7 @@ int main(int argc, char *argv[]) intmat_set_all_3x3(part2, 0,0,1, 0,1,0, -1,0,0); - tfn = intmat_intmat_mult(part1, part2); + tfn = intmat_times_intmat(part1, part2); fail += check_identity(cell, tfn); intmat_free(part1); intmat_free(part2); -- 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(-) 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