From d0adbfaf3b7bbadc372f162451d53fa4d2bb43ea Mon Sep 17 00:00:00 2001 From: Thomas White Date: Mon, 1 Oct 2018 17:00:17 +0200 Subject: cell_tool: Use libcrystfel cells_are_similar() function --- libcrystfel/src/cell-utils.c | 4 ++-- libcrystfel/src/cell-utils.h | 9 ++++++--- 2 files changed, 8 insertions(+), 5 deletions(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c index 7b1984bb..852887a6 100644 --- a/libcrystfel/src/cell-utils.c +++ b/libcrystfel/src/cell-utils.c @@ -1703,8 +1703,8 @@ static double moduli_check(double ax, double ay, double az, } -static int cells_are_similar(UnitCell *cell1, UnitCell *cell2, - const double ltl, const double atl) +int cells_are_similar(UnitCell *cell1, UnitCell *cell2, + const double ltl, const double atl) { double asx1, asy1, asz1, bsx1, bsy1, bsz1, csx1, csy1, csz1; double asx2, asy2, asz2, bsx2, bsy2, bsz2, csx2, csy2, csz2; diff --git a/libcrystfel/src/cell-utils.h b/libcrystfel/src/cell-utils.h index 5e2b2825..cc721634 100644 --- a/libcrystfel/src/cell-utils.h +++ b/libcrystfel/src/cell-utils.h @@ -3,13 +3,13 @@ * * Unit Cell utility functions * - * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2012-2018 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * Copyright © 2012 Lorenzo Galli * * Authors: - * 2009-2013,2014,2017 Thomas White - * 2012 Lorenzo Galli + * 2009-2018 Thomas White + * 2012 Lorenzo Galli * * This file is part of CrystFEL. * @@ -53,6 +53,9 @@ extern UnitCell *transform_cell_gsl(UnitCell *in, gsl_matrix *m); extern void cell_print(UnitCell *cell); +extern int cells_are_similar(UnitCell *cell1, UnitCell *cell2, + const double ltl, const double atl); + extern UnitCell *match_cell(UnitCell *cell, UnitCell *tempcell, int verbose, const float *ltl, int reduce); -- cgit v1.2.3 From e06952e9eefb0affb65d5361d1919d28bee5974d Mon Sep 17 00:00:00 2001 From: Thomas White Date: Mon, 15 Oct 2018 16:01:47 +0200 Subject: Remove cells_are_similar() from API There are two functions with this name. Avoid confusion. --- libcrystfel/src/cell-utils.c | 4 ++-- libcrystfel/src/cell-utils.h | 3 --- 2 files changed, 2 insertions(+), 5 deletions(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c index 852887a6..7b1984bb 100644 --- a/libcrystfel/src/cell-utils.c +++ b/libcrystfel/src/cell-utils.c @@ -1703,8 +1703,8 @@ static double moduli_check(double ax, double ay, double az, } -int cells_are_similar(UnitCell *cell1, UnitCell *cell2, - const double ltl, const double atl) +static int cells_are_similar(UnitCell *cell1, UnitCell *cell2, + const double ltl, const double atl) { double asx1, asy1, asz1, bsx1, bsy1, bsz1, csx1, csy1, csz1; double asx2, asy2, asz2, bsx2, bsy2, bsz2, csx2, csy2, csz2; diff --git a/libcrystfel/src/cell-utils.h b/libcrystfel/src/cell-utils.h index cc721634..d0af717a 100644 --- a/libcrystfel/src/cell-utils.h +++ b/libcrystfel/src/cell-utils.h @@ -53,9 +53,6 @@ extern UnitCell *transform_cell_gsl(UnitCell *in, gsl_matrix *m); extern void cell_print(UnitCell *cell); -extern int cells_are_similar(UnitCell *cell1, UnitCell *cell2, - const double ltl, const double atl); - extern UnitCell *match_cell(UnitCell *cell, UnitCell *tempcell, int verbose, const float *ltl, int reduce); -- cgit v1.2.3 From 481c13eac53cce18272ff3ef3e4994170c7d2e62 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Tue, 23 Oct 2018 14:18:40 +0200 Subject: Change horribly confusing names of unit cell comparison functions cells_are_similar -> compare_cell_parameters_and_orientation compare_cells -> compare_reindexed_cell_parameters_and_orientation cell_tool.c:cells_the_same -> cellutils.c:compare_cell_parameters All comparisons now done in real space, checking that centering is the same, and without uncentering anything. --- libcrystfel/src/cell-utils.c | 114 ++++++++++++++++++++++++++++++++++--------- libcrystfel/src/cell-utils.h | 16 +++++- libcrystfel/src/index.c | 6 +-- 3 files changed, 107 insertions(+), 29 deletions(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c index 7b1984bb..3d63249d 100644 --- a/libcrystfel/src/cell-utils.c +++ b/libcrystfel/src/cell-utils.c @@ -1694,6 +1694,48 @@ double cell_get_volume(UnitCell *cell) } +/** + * compare_cell_parameters: + * @a: A %UnitCell + * @b: Another %UnitCell + * @ltl: Maximum allowable fractional difference in axis lengths + * @atl: Maximum allowable difference in reciprocal angles (in radians) + * + * Compare the two unit cells. If the real space parameters match to within + * fractional difference @ltl, and the inter-axial angles match within @atl, + * 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, + * use compare_cell_parameters_and_orientation() instead. + * + * Returns: non-zero if the cells match. + * + */ +int compare_cell_parameters(UnitCell *cell1, UnitCell *cell2, + float ltl, float atl) +{ + 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; + + 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) ) return 0; + if ( !within_tolerance(b1, b2, ltl) ) return 0; + if ( !within_tolerance(c1, c2, ltl) ) return 0; + if ( fabs(al1-al2) > atl ) return 0; + if ( fabs(be1-be2) > atl ) return 0; + if ( fabs(ga1-ga2) > atl ) return 0; + + return 1; +} + + static double moduli_check(double ax, double ay, double az, double bx, double by, double bz) { @@ -1703,28 +1745,42 @@ static double moduli_check(double ax, double ay, double az, } -static int cells_are_similar(UnitCell *cell1, UnitCell *cell2, - const double ltl, const double atl) +/** + * compare_cell_parameters_and_orientation: + * @a: A %UnitCell + * @b: Another %UnitCell + * @ltl: Maximum allowable fractional difference in reciprocal axis lengths + * @atl: Maximum allowable difference in reciprocal angles (in radians) + * + * Compare the two unit cells. If the axes match in length (to within + * fractional difference @ltl) and the axes are aligned to within @atl radians, + * 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 cells @a and @b 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, + const double ltl, const double atl) { double asx1, asy1, asz1, bsx1, bsy1, bsz1, csx1, csy1, csz1; double asx2, asy2, asz2, bsx2, bsy2, bsz2, csx2, csy2, csz2; - UnitCell *pcell1, *pcell2; - - /* Compare primitive cells, not centered */ - pcell1 = uncenter_cell(cell1, NULL); - pcell2 = uncenter_cell(cell2, NULL); - - cell_get_reciprocal(pcell1, &asx1, &asy1, &asz1, - &bsx1, &bsy1, &bsz1, - &csx1, &csy1, &csz1); - cell_get_reciprocal(pcell2, &asx2, &asy2, &asz2, - &bsx2, &bsy2, &bsz2, - &csx2, &csy2, &csz2); + if ( cell_get_centering(cell1) != cell_get_centering(cell2) ) return 0; + cell_get_cartesian(cell1, &asx1, &asy1, &asz1, + &bsx1, &bsy1, &bsz1, + &csx1, &csy1, &csz1); - cell_free(pcell1); - cell_free(pcell2); + cell_get_cartesian(cell2, &asx2, &asy2, &asz2, + &bsx2, &bsy2, &bsz2, + &csx2, &csy2, &csz2); if ( angle_between(asx1, asy1, asz1, asx2, asy2, asz2) > atl ) return 0; if ( angle_between(bsx1, bsy1, bsz1, bsx2, bsy2, bsz2) > atl ) return 0; @@ -1738,28 +1794,38 @@ static int cells_are_similar(UnitCell *cell1, UnitCell *cell2, } - /** - * compare_cells: + * compare_reindexed_cell_parameters_and_orientation: * @a: A %UnitCell * @b: Another %UnitCell * @ltl: Maximum allowable fractional difference in reciprocal axis lengths * @atl: Maximum allowable difference in reciprocal angles (in radians) * @pmb: Place to store pointer to matrix * - * Compare the two units cells. If they agree to within @ltl and @atl, using - * any change of axes, returns non-zero and stores the transformation to map @b - * onto @a. + * 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 @ltl) + * and the axes aligned to within @atl radians, this function returns non-zero + * and stores the transformation to map @b onto @a. + * + * 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 @a and @b must have the same centering. Otherwise, this function + * always returns zero. * * Returns: non-zero if the cells match. * */ -int compare_cells(UnitCell *a, UnitCell *b, double ltl, double atl, - IntegerMatrix **pmb) +int compare_reindexed_cell_parameters_and_orientation(UnitCell *a, UnitCell *b, + double ltl, double atl, + IntegerMatrix **pmb) { IntegerMatrix *m; int i[9]; + if ( cell_get_centering(a) != cell_get_centering(b) ) return 0; + m = intmat_new(3, 3); for ( i[0]=-1; i[0]<=+1; i[0]++ ) { @@ -1786,7 +1852,7 @@ int compare_cells(UnitCell *a, UnitCell *b, double ltl, double atl, tfn = tfn_from_intmat(m); nc = cell_transform(b, tfn); - if ( cells_are_similar(a, nc, ltl, atl) ) { + if ( compare_cell_parameters_and_orientation(a, nc, ltl, atl) ) { if ( pmb != NULL ) *pmb = m; tfn_free(tfn); cell_free(nc); diff --git a/libcrystfel/src/cell-utils.h b/libcrystfel/src/cell-utils.h index d0af717a..cf736e4a 100644 --- a/libcrystfel/src/cell-utils.h +++ b/libcrystfel/src/cell-utils.h @@ -80,8 +80,20 @@ extern int forbidden_reflection(UnitCell *cell, extern double cell_get_volume(UnitCell *cell); -extern int compare_cells(UnitCell *a, UnitCell *b, double ltl, double atl, - IntegerMatrix **pmb); +extern int compare_cell_parameters(UnitCell *cell1, UnitCell *cell2, + float ltl, float atl); + + +extern int compare_cell_parameters_and_orientation(UnitCell *cell1, + UnitCell *cell2, + const double ltl, + const double atl); + +extern int compare_reindexed_cell_parameters_and_orientation(UnitCell *a, + UnitCell *b, + double ltl, + double atl, + IntegerMatrix **pmb); #ifdef __cplusplus } diff --git a/libcrystfel/src/index.c b/libcrystfel/src/index.c index f72334c4..ad9c0de3 100644 --- a/libcrystfel/src/index.c +++ b/libcrystfel/src/index.c @@ -698,9 +698,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_cells(crystal_get_cell(cr), - crystal_get_cell(that_cr), - 0.1, deg2rad(0.5), NULL) ) + if ( compare_reindexed_cell_parameters_and_orientation(crystal_get_cell(cr), + crystal_get_cell(that_cr), + 0.1, deg2rad(0.5), NULL) ) { crystal_set_user_flag(cr, 1); } -- cgit v1.2.3 From eaff24431149de7708b524e857b3c0807cb17c96 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Wed, 20 Feb 2019 17:06:53 +0100 Subject: Simplify the output of cell_print() --- libcrystfel/src/cell-utils.c | 26 +++++++++++++++++++------- libcrystfel/src/cell-utils.h | 1 + 2 files changed, 20 insertions(+), 7 deletions(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c index 3d63249d..39e893cf 100644 --- a/libcrystfel/src/cell-utils.c +++ b/libcrystfel/src/cell-utils.c @@ -219,11 +219,6 @@ int right_handed(UnitCell *cell) void cell_print(UnitCell *cell) { - double asx, asy, asz; - double bsx, bsy, bsz; - double csx, csy, csz; - double a, b, c, alpha, beta, gamma; - double ax, ay, az, bx, by, bz, cx, cy, cz; LatticeType lt; char cen; @@ -251,12 +246,31 @@ void cell_print(UnitCell *cell) } if ( cell_has_parameters(cell) ) { + + double a, b, c, alpha, beta, gamma; cell_get_parameters(cell, &a, &b, &c, &alpha, &beta, &gamma); STATUS("a b c alpha beta gamma\n"); STATUS("%6.2f %6.2f %6.2f A %6.2f %6.2f %6.2f deg\n", a*1e10, b*1e10, c*1e10, rad2deg(alpha), rad2deg(beta), rad2deg(gamma)); + } else { + STATUS("Unit cell parameters are not specified.\n"); + } +} + + +void cell_print_full(UnitCell *cell) +{ + + cell_print(cell); + + if ( cell_has_parameters(cell) ) { + + double asx, asy, asz; + double bsx, bsy, bsz; + double csx, csy, csz; + double ax, ay, az, bx, by, bz, cx, cy, cz; cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz); @@ -283,8 +297,6 @@ void cell_print(UnitCell *cell) STATUS("Cell representation is %s.\n", cell_rep(cell)); - } else { - STATUS("Unit cell parameters are not specified.\n"); } } diff --git a/libcrystfel/src/cell-utils.h b/libcrystfel/src/cell-utils.h index cf736e4a..47145481 100644 --- a/libcrystfel/src/cell-utils.h +++ b/libcrystfel/src/cell-utils.h @@ -52,6 +52,7 @@ extern UnitCell *rotate_cell(UnitCell *in, double omega, double phi, extern UnitCell *transform_cell_gsl(UnitCell *in, gsl_matrix *m); 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); -- cgit v1.2.3 From 4337cafe052c4ad238c969dfa4cb7c7ac52f5e07 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Thu, 25 Oct 2018 16:48:36 +0200 Subject: Use IntegerMatrix for all unit cell transformations Get rid of UnitCellTransformation, a thin wrapper which didn't do anything. --- libcrystfel/src/cell-utils.c | 260 ++++++++---------------- libcrystfel/src/cell-utils.h | 3 +- libcrystfel/src/cell.c | 414 +++++++++++++-------------------------- libcrystfel/src/cell.h | 24 +-- libcrystfel/src/integer_matrix.c | 39 ++++ libcrystfel/src/integer_matrix.h | 9 + libcrystfel/src/taketwo.c | 3 +- libcrystfel/src/xgandalf.c | 12 +- 8 files changed, 284 insertions(+), 480 deletions(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c index 39e893cf..507d4cc2 100644 --- a/libcrystfel/src/cell-utils.c +++ b/libcrystfel/src/cell-utils.c @@ -361,157 +361,114 @@ int bravais_lattice(UnitCell *cell) } -static UnitCellTransformation *uncentering_transformation(UnitCell *in, - char *new_centering, - LatticeType *new_latt) +/* Given a centered cell @in, return the integer transformation matrix which + * turns a primitive cell into @in. Set new_centering and new_latt to the + * centering and lattice type of the primitive cell (usually aP, sometimes rR, + * rarely mP) */ +static IntegerMatrix *centering_transformation(UnitCell *in, + char *new_centering, + LatticeType *new_latt, + char *new_ua) { - UnitCellTransformation *t; - const double OT = 1.0/3.0; - const double TT = 2.0/3.0; - const double H = 0.5; LatticeType lt; char ua, cen; + IntegerMatrix *t = NULL; lt = cell_get_lattice_type(in); ua = cell_get_unique_axis(in); cen = cell_get_centering(in); - t = tfn_identity(); - if ( t == NULL ) return NULL; - - if ( ua == 'a' ) { - tfn_combine(t, tfn_vector(0,1,0), - tfn_vector(0,0,1), - tfn_vector(1,0,0)); - if ( lt == L_MONOCLINIC ) { - assert(cen != 'A'); - switch ( cen ) { - case 'B' : cen = 'A'; break; - case 'C' : cen = 'B'; break; - case 'I' : cen = 'I'; break; - } - } - } - - if ( ua == 'b' ) { - tfn_combine(t, tfn_vector(0,0,1), - tfn_vector(1,0,0), - tfn_vector(0,1,0)); - if ( lt == L_MONOCLINIC ) { - assert(cen != 'B'); - switch ( cen ) { - case 'C' : cen = 'A'; break; - case 'A' : cen = 'B'; break; - case 'I' : cen = 'I'; break; - } - } - } - - switch ( cen ) { - - case 'P' : + if ( (cen=='P') || (cen=='R') ) { + *new_centering = cen; *new_latt = lt; - *new_centering = 'P'; - break; - - case 'R' : - *new_latt = L_RHOMBOHEDRAL; - *new_centering = 'R'; - break; + *new_ua = ua; + t = intmat_identity(3); + } - case 'I' : - tfn_combine(t, tfn_vector(-H,H,H), - tfn_vector(H,-H,H), - tfn_vector(H,H,-H)); + if ( cen == 'I' ) { + t = intmat_create_3x3(0, 1, 1, + 1, 0, 1, + 1, 1, 0); if ( lt == L_CUBIC ) { *new_latt = L_RHOMBOHEDRAL; *new_centering = 'R'; + *new_ua = '*'; } else { - /* Tetragonal or orthorhombic */ *new_latt = L_TRICLINIC; *new_centering = 'P'; + *new_ua = '*'; } - break; + } - case 'F' : - tfn_combine(t, tfn_vector(0,H,H), - tfn_vector(H,0,H), - tfn_vector(H,H,0)); + if ( cen == 'F' ) { + t = intmat_create_3x3(-1, 1, 1, + 1, -1, 1, + 1, 1, -1); if ( lt == L_CUBIC ) { *new_latt = L_RHOMBOHEDRAL; *new_centering = 'R'; + *new_ua = '*'; } else { - assert(lt == L_ORTHORHOMBIC); *new_latt = L_TRICLINIC; *new_centering = 'P'; + *new_ua = '*'; } - break; + } - case 'A' : - tfn_combine(t, tfn_vector( 1, 0, 0), - tfn_vector( 0, H, H), - tfn_vector( 0,-H, H)); + if ( (lt == L_HEXAGONAL) && (cen == 'H') && (ua == 'c') ) { + /* Obverse setting */ + t = intmat_create_3x3(1, -1, 0, + 0, 1, -1, + 1, 1, 1); + assert(lt == L_HEXAGONAL); + assert(ua == 'c'); + *new_latt = L_RHOMBOHEDRAL; + *new_centering = 'R'; + *new_ua = '*'; + } + + if ( cen == 'A' ) { + t = intmat_create_3x3(1, 0, 0, + 0, 1, -1, + 0, 1, 1); if ( lt == L_ORTHORHOMBIC ) { *new_latt = L_MONOCLINIC; + *new_centering = 'P'; + *new_ua = 'a'; } else { *new_latt = L_TRICLINIC; + *new_centering = 'P'; + *new_ua = '*'; } - *new_centering = 'P'; - break; + } - case 'B' : - tfn_combine(t, tfn_vector( H, 0, H), - tfn_vector( 0, 1, 0), - tfn_vector(-H, 0, H)); + if ( cen == 'B' ) { + t = intmat_create_3x3(1, 0, -1, + 0, 1, 0, + 1, 0, 1); if ( lt == L_ORTHORHOMBIC ) { *new_latt = L_MONOCLINIC; + *new_centering = 'P'; + *new_ua = 'b'; } else { *new_latt = L_TRICLINIC; + *new_centering = 'P'; + *new_ua = '*'; } - *new_centering = 'P'; - break; + } - case 'C' : - tfn_combine(t, tfn_vector( H, H, 0), - tfn_vector(-H, H, 0), - tfn_vector( 0, 0, 1)); + if ( cen == 'C' ) { + t = intmat_create_3x3(1, -1, 0, + 1, 1, 0, + 0, 0, 1); if ( lt == L_ORTHORHOMBIC ) { *new_latt = L_MONOCLINIC; + *new_centering = 'P'; + *new_ua = 'c'; } else { *new_latt = L_TRICLINIC; - } - *new_centering = 'P'; - break; - - case 'H' : - /* Obverse setting */ - tfn_combine(t, tfn_vector(TT,OT,OT), - tfn_vector(-OT,OT,OT), - tfn_vector(-OT,-TT,OT)); - assert(lt == L_HEXAGONAL); - *new_latt = L_RHOMBOHEDRAL; - *new_centering = 'R'; - break; - - default : - ERROR("Invalid centering '%c'\n", cell_get_centering(in)); - return NULL; - - } - - /* Reverse the axis permutation, but only if this was not an H->R - * transformation */ - if ( !((cen=='H') && (*new_latt == L_RHOMBOHEDRAL)) ) { - if ( ua == 'a' ) { - tfn_combine(t, tfn_vector(0,0,1), - tfn_vector(1,0,0), - tfn_vector(0,1,0)); - } - - if ( ua == 'b' ) { - tfn_combine(t, tfn_vector(0,1,0), - tfn_vector(0,0,1), - tfn_vector(1,0,0)); + *new_centering = 'P'; + *new_ua = '*'; } } @@ -532,32 +489,28 @@ static UnitCellTransformation *uncentering_transformation(UnitCell *in, * setting. * */ -UnitCell *uncenter_cell(UnitCell *in, UnitCellTransformation **t) +UnitCell *uncenter_cell(UnitCell *in, IntegerMatrix **t) { - UnitCellTransformation *tt; + IntegerMatrix *tt; char new_centering; LatticeType new_latt; + char new_ua; UnitCell *out; - if ( !bravais_lattice(in) ) { - ERROR("Cannot uncenter: not a Bravais lattice.\n"); - cell_print(in); - return NULL; - } - - tt = uncentering_transformation(in, &new_centering, &new_latt); + tt = centering_transformation(in, &new_centering, &new_latt, &new_ua); if ( tt == NULL ) return NULL; - out = cell_transform(in, tt); + out = cell_transform_inverse(in, tt); if ( out == NULL ) return NULL; cell_set_lattice_type(out, new_latt); cell_set_centering(out, new_centering); + cell_set_unique_axis(out, new_ua); if ( t != NULL ) { *t = tt; } else { - tfn_free(tt); + intmat_free(tt); } return out; @@ -621,11 +574,11 @@ UnitCell *match_cell(UnitCell *cell_in, UnitCell *template_in, int verbose, float angtol = deg2rad(tols[3]); UnitCell *cell; UnitCell *template; - UnitCellTransformation *uncentering; + IntegerMatrix *centering; UnitCell *new_cell_trans; /* "Un-center" the template unit cell to make the comparison easier */ - template = uncenter_cell(template_in, &uncentering); + template = uncenter_cell(template_in, ¢ering); if ( template == NULL ) return NULL; /* The candidate cell is also uncentered, because it might be centered @@ -639,7 +592,7 @@ UnitCell *match_cell(UnitCell *cell_in, UnitCell *template_in, int verbose, ERROR("Couldn't get reciprocal cell for template.\n"); cell_free(template); cell_free(cell); - tfn_free(uncentering); + intmat_free(centering); return NULL; } @@ -661,7 +614,7 @@ UnitCell *match_cell(UnitCell *cell_in, UnitCell *template_in, int verbose, ERROR("Couldn't get reciprocal cell.\n"); cell_free(template); cell_free(cell); - tfn_free(uncentering); + intmat_free(centering); return NULL; } @@ -823,7 +776,7 @@ UnitCell *match_cell(UnitCell *cell_in, UnitCell *template_in, int verbose, /* Reverse the de-centering transformation */ if ( new_cell != NULL ) { - new_cell_trans = cell_transform_inverse(new_cell, uncentering); + new_cell_trans = cell_transform(new_cell, centering); cell_free(new_cell); cell_set_lattice_type(new_cell_trans, cell_get_lattice_type(template_in)); @@ -833,13 +786,13 @@ UnitCell *match_cell(UnitCell *cell_in, UnitCell *template_in, int verbose, cell_get_unique_axis(template_in)); cell_free(template); - tfn_free(uncentering); + intmat_free(centering); return new_cell_trans; } else { cell_free(template); - tfn_free(uncentering); + intmat_free(centering); return NULL; } } @@ -862,7 +815,7 @@ UnitCell *match_cell_ab(UnitCell *cell_in, UnitCell *template_in) int have_real_c; UnitCell *cell; UnitCell *template; - UnitCellTransformation *to_given_cell; + IntegerMatrix *to_given_cell; UnitCell *new_cell; UnitCell *new_cell_trans; @@ -1455,49 +1408,6 @@ void cell_fudge_gslcblas() } -UnitCell *transform_cell_gsl(UnitCell *in, gsl_matrix *m) -{ - gsl_matrix *c; - double asx, asy, asz; - double bsx, bsy, bsz; - double csx, csy, csz; - gsl_matrix *res; - UnitCell *out; - - cell_get_reciprocal(in, &asx, &asy, &asz, &bsx, &bsy, - &bsz, &csx, &csy, &csz); - - c = gsl_matrix_alloc(3, 3); - gsl_matrix_set(c, 0, 0, asx); - gsl_matrix_set(c, 1, 0, asy); - gsl_matrix_set(c, 2, 0, asz); - gsl_matrix_set(c, 0, 1, bsx); - gsl_matrix_set(c, 1, 1, bsy); - gsl_matrix_set(c, 2, 1, bsz); - gsl_matrix_set(c, 0, 2, csx); - gsl_matrix_set(c, 1, 2, csy); - gsl_matrix_set(c, 2, 2, csz); - - res = gsl_matrix_calloc(3, 3); - gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, m, c, 0.0, res); - - out = cell_new_from_cell(in); - cell_set_reciprocal(out, gsl_matrix_get(res, 0, 0), - gsl_matrix_get(res, 1, 0), - gsl_matrix_get(res, 2, 0), - gsl_matrix_get(res, 0, 1), - gsl_matrix_get(res, 1, 1), - gsl_matrix_get(res, 2, 1), - gsl_matrix_get(res, 0, 2), - gsl_matrix_get(res, 1, 2), - gsl_matrix_get(res, 2, 2)); - - gsl_matrix_free(res); - gsl_matrix_free(c); - return out; -} - - /** * rotate_cell: * @in: A %UnitCell to rotate @@ -1658,7 +1568,7 @@ int forbidden_reflection(UnitCell *cell, cen = cell_get_centering(cell); /* Reflection conditions here must match the transformation matrices - * in uncentering_transformation(). tests/centering_check verifies + * in centering_transformation(). tests/centering_check verifies * this (amongst other things). */ if ( cen == 'P' ) return 0; @@ -1850,7 +1760,6 @@ int compare_reindexed_cell_parameters_and_orientation(UnitCell *a, UnitCell *b, for ( i[7]=-1; i[7]<=+1; i[7]++ ) { for ( i[8]=-1; i[8]<=+1; i[8]++ ) { - UnitCellTransformation *tfn; UnitCell *nc; int j, k; int l = 0; @@ -1861,17 +1770,14 @@ int compare_reindexed_cell_parameters_and_orientation(UnitCell *a, UnitCell *b, if ( intmat_det(m) != +1 ) continue; - tfn = tfn_from_intmat(m); - nc = cell_transform(b, tfn); + nc = cell_transform(b, m); if ( compare_cell_parameters_and_orientation(a, nc, ltl, atl) ) { if ( pmb != NULL ) *pmb = m; - tfn_free(tfn); cell_free(nc); return 1; } - tfn_free(tfn); cell_free(nc); } diff --git a/libcrystfel/src/cell-utils.h b/libcrystfel/src/cell-utils.h index 47145481..1a4f3bcc 100644 --- a/libcrystfel/src/cell-utils.h +++ b/libcrystfel/src/cell-utils.h @@ -49,7 +49,6 @@ extern double resolution(UnitCell *cell, extern UnitCell *cell_rotate(UnitCell *in, struct quaternion quat); extern UnitCell *rotate_cell(UnitCell *in, double omega, double phi, double rot); -extern UnitCell *transform_cell_gsl(UnitCell *in, gsl_matrix *m); extern void cell_print(UnitCell *cell); extern void cell_print_full(UnitCell *cell); @@ -67,7 +66,7 @@ extern int cell_is_sensible(UnitCell *cell); extern int validate_cell(UnitCell *cell); -extern UnitCell *uncenter_cell(UnitCell *in, UnitCellTransformation **t); +extern UnitCell *uncenter_cell(UnitCell *in, IntegerMatrix **m); extern int bravais_lattice(UnitCell *cell); diff --git a/libcrystfel/src/cell.c b/libcrystfel/src/cell.c index cc18b49d..242094f6 100644 --- a/libcrystfel/src/cell.c +++ b/libcrystfel/src/cell.c @@ -45,6 +45,7 @@ #include "cell.h" #include "utils.h" #include "image.h" +#include "integer_matrix.h" /** @@ -600,70 +601,88 @@ const char *cell_rep(UnitCell *cell) } -struct _unitcelltransformation +UnitCell *cell_transform_gsl_direct(UnitCell *in, gsl_matrix *m) { - gsl_matrix *m; -}; - - -/** - * tfn_inverse: - * @t: A %UnitCellTransformation. - * - * Calculates the inverse of @t. That is, if you apply cell_transform() to a - * %UnitCell using @t, and then apply cell_transform() to the result using - * tfn_inverse(@t) instead of @t, you will recover the same lattice vectors - * (but note that the lattice type, centering and unique axis information will - * be lost). - * - * Returns: The inverse of @t. - * - */ -UnitCellTransformation *tfn_inverse(UnitCellTransformation *t) -{ - int s; - gsl_matrix *m; - gsl_matrix *inv; - gsl_permutation *perm; - UnitCellTransformation *out; - - m = gsl_matrix_alloc(3, 3); - if ( m == NULL ) return NULL; + gsl_matrix *c; + double asx, asy, asz; + double bsx, bsy, bsz; + double csx, csy, csz; + gsl_matrix *res; + UnitCell *out; - out = tfn_identity(); - if ( out == NULL ) { - gsl_matrix_free(m); - return NULL; - } + cell_get_cartesian(in, &asx, &asy, &asz, &bsx, &bsy, + &bsz, &csx, &csy, &csz); + + c = gsl_matrix_alloc(3, 3); + gsl_matrix_set(c, 0, 0, asx); + gsl_matrix_set(c, 0, 1, asy); + gsl_matrix_set(c, 0, 2, asz); + gsl_matrix_set(c, 1, 0, bsx); + gsl_matrix_set(c, 1, 1, bsy); + gsl_matrix_set(c, 1, 2, bsz); + gsl_matrix_set(c, 2, 0, csx); + gsl_matrix_set(c, 2, 1, csy); + gsl_matrix_set(c, 2, 2, csz); + + res = gsl_matrix_calloc(3, 3); + gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, m, c, 0.0, res); + + out = cell_new_from_cell(in); + cell_set_cartesian(out, gsl_matrix_get(res, 0, 0), + gsl_matrix_get(res, 0, 1), + gsl_matrix_get(res, 0, 2), + gsl_matrix_get(res, 1, 0), + gsl_matrix_get(res, 1, 1), + gsl_matrix_get(res, 1, 2), + gsl_matrix_get(res, 2, 0), + gsl_matrix_get(res, 2, 1), + gsl_matrix_get(res, 2, 2)); + + gsl_matrix_free(res); + gsl_matrix_free(c); + return out; +} - gsl_matrix_memcpy(m, t->m); - perm = gsl_permutation_alloc(m->size1); - if ( perm == NULL ) { - ERROR("Couldn't allocate permutation\n"); - return NULL; - } - inv = gsl_matrix_alloc(m->size1, m->size2); - if ( inv == NULL ) { - ERROR("Couldn't allocate inverse\n"); - gsl_permutation_free(perm); - return NULL; - } - if ( gsl_linalg_LU_decomp(m, perm, &s) ) { - ERROR("Couldn't decompose matrix\n"); - gsl_permutation_free(perm); - return NULL; - } - if ( gsl_linalg_LU_invert(m, perm, inv) ) { - ERROR("Couldn't invert transformation matrix\n"); - gsl_permutation_free(perm); - return NULL; - } - gsl_permutation_free(perm); +UnitCell *cell_transform_gsl_reciprocal(UnitCell *in, gsl_matrix *m) +{ + gsl_matrix *c; + double asx, asy, asz; + double bsx, bsy, bsz; + double csx, csy, csz; + gsl_matrix *res; + UnitCell *out; - gsl_matrix_free(out->m); - gsl_matrix_free(m); - out->m = inv; + cell_get_reciprocal(in, &asx, &asy, &asz, &bsx, &bsy, + &bsz, &csx, &csy, &csz); + + c = gsl_matrix_alloc(3, 3); + gsl_matrix_set(c, 0, 0, asx); + gsl_matrix_set(c, 0, 1, asy); + gsl_matrix_set(c, 0, 2, asz); + gsl_matrix_set(c, 1, 0, bsx); + gsl_matrix_set(c, 1, 1, bsy); + gsl_matrix_set(c, 1, 2, bsz); + gsl_matrix_set(c, 2, 0, csx); + gsl_matrix_set(c, 2, 1, csy); + gsl_matrix_set(c, 2, 2, csz); + + res = gsl_matrix_calloc(3, 3); + gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, m, c, 0.0, res); + + out = cell_new_from_cell(in); + cell_set_reciprocal(out, gsl_matrix_get(res, 0, 0), + gsl_matrix_get(res, 0, 1), + gsl_matrix_get(res, 0, 2), + gsl_matrix_get(res, 1, 0), + gsl_matrix_get(res, 1, 1), + gsl_matrix_get(res, 1, 2), + gsl_matrix_get(res, 2, 0), + gsl_matrix_get(res, 2, 1), + gsl_matrix_get(res, 2, 2)); + + gsl_matrix_free(res); + gsl_matrix_free(c); return out; } @@ -671,63 +690,40 @@ UnitCellTransformation *tfn_inverse(UnitCellTransformation *t) /** * cell_transform: * @cell: A %UnitCell. - * @t: A %UnitCellTransformation. + * @t: An %IntegerMatrix. * - * Applies @t to @cell. Note that the lattice type, centering and unique axis - * information will not be preserved. + * Applies @t to @cell. * * Returns: Transformed copy of @cell. * */ -UnitCell *cell_transform(UnitCell *cell, UnitCellTransformation *t) +UnitCell *cell_transform(UnitCell *cell, IntegerMatrix *m) { UnitCell *out; - double ax, ay, az; - double bx, by, bz; - double cx, cy, cz; - gsl_matrix *m; - gsl_matrix *a; + gsl_matrix *tm; - if ( t == NULL ) return NULL; - - out = cell_new_from_cell(cell); - if ( out == NULL ) return NULL; - - if ( cell_get_cartesian(out, &ax, &ay, &az, - &bx, &by, &bz, - &cx, &cy, &cz) ) return NULL; + if ( m == NULL ) return NULL; - m = gsl_matrix_alloc(3,3); - a = gsl_matrix_calloc(3,3); - if ( (m == NULL) || (a == NULL) ) { - cell_free(out); + tm = gsl_matrix_alloc(3,3); + if ( tm == NULL ) { return NULL; } - gsl_matrix_set(m, 0, 0, ax); - gsl_matrix_set(m, 0, 1, ay); - gsl_matrix_set(m, 0, 2, az); - gsl_matrix_set(m, 1, 0, bx); - gsl_matrix_set(m, 1, 1, by); - gsl_matrix_set(m, 1, 2, bz); - gsl_matrix_set(m, 2, 0, cx); - gsl_matrix_set(m, 2, 1, cy); - gsl_matrix_set(m, 2, 2, cz); + gsl_matrix_set(tm, 0, 0, intmat_get(m, 0, 0)); + gsl_matrix_set(tm, 0, 1, intmat_get(m, 0, 1)); + gsl_matrix_set(tm, 0, 2, intmat_get(m, 0, 2)); + gsl_matrix_set(tm, 1, 0, intmat_get(m, 1, 0)); + gsl_matrix_set(tm, 1, 1, intmat_get(m, 1, 1)); + gsl_matrix_set(tm, 1, 2, intmat_get(m, 1, 2)); + gsl_matrix_set(tm, 2, 0, intmat_get(m, 2, 0)); + gsl_matrix_set(tm, 2, 1, intmat_get(m, 2, 1)); + gsl_matrix_set(tm, 2, 2, intmat_get(m, 2, 2)); - gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, t->m, m, 0.0, a); + out = cell_transform_gsl_direct(cell, tm); - cell_set_cartesian(out, gsl_matrix_get(a, 0, 0), - gsl_matrix_get(a, 0, 1), - gsl_matrix_get(a, 0, 2), - gsl_matrix_get(a, 1, 0), - gsl_matrix_get(a, 1, 1), - gsl_matrix_get(a, 1, 2), - gsl_matrix_get(a, 2, 0), - gsl_matrix_get(a, 2, 1), - gsl_matrix_get(a, 2, 2)); + /* FIXME: Update centering, unique axis, lattice type */ - gsl_matrix_free(a); - gsl_matrix_free(m); + gsl_matrix_free(tm); return out; } @@ -736,197 +732,67 @@ UnitCell *cell_transform(UnitCell *cell, UnitCellTransformation *t) /** * cell_transform_inverse: * @cell: A %UnitCell. - * @t: A %UnitCellTransformation. + * @m: An %IntegerMatrix * - * Applies the inverse of @t to @cell. + * Applies the inverse of @m to @cell. * * Returns: Transformed copy of @cell. * */ -UnitCell *cell_transform_inverse(UnitCell *cell, UnitCellTransformation *t) +UnitCell *cell_transform_inverse(UnitCell *cell, IntegerMatrix *m) { - UnitCellTransformation *inv; UnitCell *out; + gsl_matrix *tm; + gsl_matrix *inv; + gsl_permutation *perm; + int s; - inv = tfn_inverse(t); - out = cell_transform(cell, inv); - tfn_free(inv); - return out; -} - - -/** - * tfn_identity: - * - * Returns: A %UnitCellTransformation corresponding to an identity operation. - * - */ -UnitCellTransformation *tfn_identity() -{ - UnitCellTransformation *tfn; - - tfn = calloc(1, sizeof(UnitCellTransformation)); - if ( tfn == NULL ) return NULL; + if ( m == NULL ) return NULL; - tfn->m = gsl_matrix_alloc(3, 3); - if ( tfn->m == NULL ) { - free(tfn); + tm = gsl_matrix_alloc(3,3); + if ( tm == NULL ) { return NULL; } - gsl_matrix_set_identity(tfn->m); - - return tfn; -} - - - -/** - * tfn_from_intmat: - * @m: An %IntegerMatrix - * - * Returns: A %UnitCellTransformation corresponding to @m. - * - */ -UnitCellTransformation *tfn_from_intmat(IntegerMatrix *m) -{ - UnitCellTransformation *tfn; - int i, j; - - tfn = tfn_identity(); - if ( tfn == NULL ) return NULL; - - for ( i=0; i<3; i++ ) { - for ( j=0; j<3; j++ ) { - gsl_matrix_set(tfn->m, i, j, intmat_get(m, i, j)); - } - } - - return tfn; -} - - -/** - * tfn_combine: - * @t: A %UnitCellTransformation - * @na: Pointer to three doubles representing naa, nab, nac - * @nb: Pointer to three doubles representing nba, nbb, nbc - * @nc: Pointer to three doubles representing nca, ncb, ncc - * - * Updates @t such that it represents its previous transformation followed by - * a new transformation, corresponding to letting a = naa*a + nab*b + nac*c. - * Likewise, a = nba*a + nbb*b + nbc*c and c = nca*a + ncb*b + ncc*c. - * - */ -void tfn_combine(UnitCellTransformation *t, double *na, double *nb, double *nc) -{ - gsl_matrix *a; - gsl_matrix *n; - - n = gsl_matrix_alloc(3, 3); - a = gsl_matrix_calloc(3, 3); - if ( (n == NULL) || (a == NULL) ) { - return; + gsl_matrix_set(tm, 0, 0, intmat_get(m, 0, 0)); + gsl_matrix_set(tm, 0, 1, intmat_get(m, 0, 1)); + gsl_matrix_set(tm, 0, 2, intmat_get(m, 0, 2)); + gsl_matrix_set(tm, 1, 0, intmat_get(m, 1, 0)); + gsl_matrix_set(tm, 1, 1, intmat_get(m, 1, 1)); + gsl_matrix_set(tm, 1, 2, intmat_get(m, 1, 2)); + gsl_matrix_set(tm, 2, 0, intmat_get(m, 2, 0)); + gsl_matrix_set(tm, 2, 1, intmat_get(m, 2, 1)); + gsl_matrix_set(tm, 2, 2, intmat_get(m, 2, 2)); + + perm = gsl_permutation_alloc(3); + if ( perm == NULL ) { + ERROR("Couldn't allocate permutation\n"); + return NULL; } - gsl_matrix_set(n, 0, 0, na[0]); - gsl_matrix_set(n, 0, 1, na[1]); - gsl_matrix_set(n, 0, 2, na[2]); - gsl_matrix_set(n, 1, 0, nb[0]); - gsl_matrix_set(n, 1, 1, nb[1]); - gsl_matrix_set(n, 1, 2, nb[2]); - gsl_matrix_set(n, 2, 0, nc[0]); - gsl_matrix_set(n, 2, 1, nc[1]); - gsl_matrix_set(n, 2, 2, nc[2]); - - free(na); - free(nb); - free(nc); - - gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, n, t->m, 0.0, a); - - gsl_matrix_free(t->m); - t->m = a; - gsl_matrix_free(n); -} - - -/** - * tfn_vector: - * @a: Amount of "a" to include in new vector - * @b: Amount of "b" to include in new vector - * @c: Amount of "c" to include in new vector - * - * This is a convenience function to use when sending vectors to tfn_combine(): - * tfn_combine(tfn, tfn_vector(1,0,0), - * tfn_vector(0,2,0), - * tfn_vector(0,0,1)); - * - */ -double *tfn_vector(double a, double b, double c) -{ - double *vec = malloc(3*sizeof(double)); - if ( vec == NULL ) return NULL; - vec[0] = a; vec[1] = b; vec[2] = c; - return vec; -} - - -/** - * tfn_print: - * @t: A %UnitCellTransformation - * - * Prints information about @t to stderr. - * - */ -void tfn_print(UnitCellTransformation *t) -{ - gsl_permutation *perm; - gsl_matrix *lu; - int s; - - STATUS("New a = %+.2fa %+.2fb %+.2fc\n", gsl_matrix_get(t->m, 0, 0), - gsl_matrix_get(t->m, 0, 1), - gsl_matrix_get(t->m, 0, 2)); - STATUS("New b = %+.2fa %+.2fb %+.2fc\n", gsl_matrix_get(t->m, 1, 0), - gsl_matrix_get(t->m, 1, 1), - gsl_matrix_get(t->m, 1, 2)); - STATUS("New c = %+.2fa %+.2fb %+.2fc\n", gsl_matrix_get(t->m, 2, 0), - gsl_matrix_get(t->m, 2, 1), - gsl_matrix_get(t->m, 2, 2)); - lu = gsl_matrix_alloc(3, 3); - if ( lu == NULL ) { - ERROR("Couldn't allocate LU decomposition.\n"); - return; + inv = gsl_matrix_alloc(3, 3); + if ( inv == NULL ) { + ERROR("Couldn't allocate inverse\n"); + gsl_permutation_free(perm); + return NULL; } - - gsl_matrix_memcpy(lu, t->m); - - perm = gsl_permutation_alloc(t->m->size1); - if ( perm == NULL ) { - ERROR("Couldn't allocate permutation.\n"); - gsl_matrix_free(lu); - return; + if ( gsl_linalg_LU_decomp(tm, perm, &s) ) { + ERROR("Couldn't decompose matrix\n"); + gsl_permutation_free(perm); + return NULL; } - if ( gsl_linalg_LU_decomp(lu, perm, &s) ) { - ERROR("LU decomposition failed.\n"); + if ( gsl_linalg_LU_invert(tm, perm, inv) ) { + ERROR("Couldn't invert transformation matrix\n"); gsl_permutation_free(perm); - gsl_matrix_free(lu); - return; + return NULL; } + gsl_permutation_free(perm); - STATUS("Transformation determinant = %.2f\n", gsl_linalg_LU_det(lu, s)); -} + out = cell_transform_gsl_direct(cell, inv); + /* FIXME: Update centering, unique axis, lattice type */ -/** - * tfn_free: - * @t: A %UnitCellTransformation - * - * Frees all resources associated with @t. - * - */ -void tfn_free(UnitCellTransformation *t) -{ - gsl_matrix_free(t->m); - free(t); + gsl_matrix_free(tm); + gsl_matrix_free(inv); + + return out; } diff --git a/libcrystfel/src/cell.h b/libcrystfel/src/cell.h index 39e6a1ed..e6ccbdf7 100644 --- a/libcrystfel/src/cell.h +++ b/libcrystfel/src/cell.h @@ -91,14 +91,6 @@ typedef enum typedef struct _unitcell UnitCell; -/** - * UnitCellTransformation: - * - * This opaque data structure represents a tranformation of a unit cell, such - * as a rotation or a centering operation. - **/ -typedef struct _unitcelltransformation UnitCellTransformation; - #ifdef __cplusplus extern "C" { #endif @@ -156,18 +148,10 @@ extern void cell_set_unique_axis(UnitCell *cell, char unique_axis); extern const char *cell_rep(UnitCell *cell); -extern UnitCell *cell_transform(UnitCell *cell, UnitCellTransformation *t); -extern UnitCell *cell_transform_inverse(UnitCell *cell, - UnitCellTransformation *t); - -extern UnitCellTransformation *tfn_identity(void); -extern UnitCellTransformation *tfn_from_intmat(IntegerMatrix *m); -extern void tfn_combine(UnitCellTransformation *t, - double *na, double *nb, double *nc); -extern void tfn_print(UnitCellTransformation *t); -extern UnitCellTransformation *tfn_inverse(UnitCellTransformation *t); -extern double *tfn_vector(double a, double b, double c); -extern void tfn_free(UnitCellTransformation *t); +extern UnitCell *cell_transform_gsl_direct(UnitCell *in, gsl_matrix *m); +extern UnitCell *cell_transform_gsl_reciprocal(UnitCell *in, gsl_matrix *m); +extern UnitCell *cell_transform(UnitCell *cell, IntegerMatrix *m); +extern UnitCell *cell_transform_inverse(UnitCell *cell, IntegerMatrix *m); #ifdef __cplusplus } diff --git a/libcrystfel/src/integer_matrix.c b/libcrystfel/src/integer_matrix.c index c31072b4..1e616e43 100644 --- a/libcrystfel/src/integer_matrix.c +++ b/libcrystfel/src/integer_matrix.c @@ -342,6 +342,9 @@ static IntegerMatrix *intmat_cofactors(const IntegerMatrix *m) * * Calculates the inverse of @m. Inefficiently. * + * Works only if the inverse of the matrix is also an integer matrix, + * i.e. if the determinant of @m is +/- 1. + * * Returns: the inverse of @m, or NULL on error. **/ IntegerMatrix *intmat_inverse(const IntegerMatrix *m) @@ -532,3 +535,39 @@ IntegerMatrix *intmat_identity(int size) return m; } + + +/** + * intmat_set_all_3x3 + * @m: An %IntegerMatrix + * + * Returns: an identity %IntegerMatrix with side length @size, or NULL on error. + * + */ +void intmat_set_all_3x3(IntegerMatrix *m, + signed int m11, signed int m12, signed int m13, + signed int m21, signed int m22, signed int m23, + signed int m31, signed int m32, signed int m33) +{ + intmat_set(m, 0, 0, m11); + intmat_set(m, 0, 1, m12); + intmat_set(m, 0, 2, m13); + + intmat_set(m, 1, 0, m21); + intmat_set(m, 1, 1, m22); + intmat_set(m, 1, 2, m23); + + intmat_set(m, 2, 0, m31); + intmat_set(m, 2, 1, m32); + intmat_set(m, 2, 2, m33); +} + + +IntegerMatrix *intmat_create_3x3(signed int m11, signed int m12, signed int m13, + signed int m21, signed int m22, signed int m23, + signed int m31, signed int m32, signed int m33) +{ + IntegerMatrix *t = intmat_new(3, 3); + intmat_set_all_3x3(t, m11, m12, m13, m21, m22, m23, m31, m32, m33); + return t; +} diff --git a/libcrystfel/src/integer_matrix.h b/libcrystfel/src/integer_matrix.h index 6616b5e8..c4983aef 100644 --- a/libcrystfel/src/integer_matrix.h +++ b/libcrystfel/src/integer_matrix.h @@ -56,6 +56,15 @@ extern void intmat_set(IntegerMatrix *m, unsigned int i, unsigned int j, extern signed int intmat_get(const IntegerMatrix *m, unsigned int i, unsigned int j); +extern void intmat_set_all_3x3(IntegerMatrix *m, + signed int m11, signed int m12, signed int m13, + signed int m21, signed int m22, signed int m23, + signed int m31, signed int m32, signed int m33); + +extern IntegerMatrix *intmat_create_3x3(signed int m11, signed int m12, signed int m13, + signed int m21, signed int m22, signed int m23, + signed int m31, signed int m32, signed int m33); + /* Matrix-(int)vector multiplication */ extern signed int *intmat_intvec_mult(const IntegerMatrix *m, const signed int *vec); diff --git a/libcrystfel/src/taketwo.c b/libcrystfel/src/taketwo.c index 203e5a05..37a29722 100644 --- a/libcrystfel/src/taketwo.c +++ b/libcrystfel/src/taketwo.c @@ -98,6 +98,7 @@ #include #include +#include "cell.h" #include "cell-utils.h" #include "index.h" #include "taketwo.h" @@ -2058,7 +2059,7 @@ static UnitCell *run_taketwo(UnitCell *cell, const struct taketwo_options *opts, tp->numPrevs++; /* Prepare the solution for CrystFEL friendliness */ - result = transform_cell_gsl(cell, solution); + result = cell_transform_gsl_reciprocal(cell, solution); cleanup_taketwo_cell(&ttCell); return result; diff --git a/libcrystfel/src/xgandalf.c b/libcrystfel/src/xgandalf.c index 5f31df33..e6618560 100644 --- a/libcrystfel/src/xgandalf.c +++ b/libcrystfel/src/xgandalf.c @@ -46,7 +46,7 @@ struct xgandalf_private_data { IndexingMethod indm; UnitCell *cellTemplate; Lattice_t sampleRealLattice_A; //same as cellTemplate - UnitCellTransformation *uncenteringTransformation; + IntegerMatrix *centeringTransformation; LatticeTransform_t latticeReductionTransform; }; @@ -114,7 +114,7 @@ int run_xgandalf(struct image *image, void *ipriv) if(xgandalf_private_data->cellTemplate != NULL){ restoreCell(uc, &xgandalf_private_data->latticeReductionTransform); - UnitCell *new_cell_trans = cell_transform_inverse(uc, xgandalf_private_data->uncenteringTransformation); + UnitCell *new_cell_trans = cell_transform(uc, xgandalf_private_data->centeringTransformation); cell_free(uc); uc = new_cell_trans; @@ -147,7 +147,7 @@ void *xgandalf_prepare(IndexingMethod *indm, UnitCell *cell, allocReciprocalPeaks(&(xgandalf_private_data->reciprocalPeaks_1_per_A)); xgandalf_private_data->indm = *indm; xgandalf_private_data->cellTemplate = NULL; - xgandalf_private_data->uncenteringTransformation = NULL; + xgandalf_private_data->centeringTransformation = NULL; float tolerance = xgandalf_opts->tolerance; samplingPitch_t samplingPitch = xgandalf_opts->sampling_pitch; @@ -157,7 +157,7 @@ void *xgandalf_prepare(IndexingMethod *indm, UnitCell *cell, xgandalf_private_data->cellTemplate = cell; - UnitCell* primitiveCell = uncenter_cell(cell, &xgandalf_private_data->uncenteringTransformation); + UnitCell* primitiveCell = uncenter_cell(cell, &xgandalf_private_data->centeringTransformation); UnitCell *uc = cell_new_from_cell(primitiveCell); reduceCell(primitiveCell, &xgandalf_private_data->latticeReductionTransform); @@ -250,8 +250,8 @@ void xgandalf_cleanup(void *pp) freeReciprocalPeaks(xgandalf_private_data->reciprocalPeaks_1_per_A); IndexerPlain_delete(xgandalf_private_data->indexer); - if(xgandalf_private_data->uncenteringTransformation != NULL){ - tfn_free(xgandalf_private_data->uncenteringTransformation); + if(xgandalf_private_data->centeringTransformation != NULL){ + intmat_free(xgandalf_private_data->centeringTransformation); } free(xgandalf_private_data); } -- cgit v1.2.3 From 9a3ef1de0b661085a26d9be60bcff9e261783acd Mon Sep 17 00:00:00 2001 From: Thomas White Date: Mon, 29 Oct 2018 10:05:24 +0100 Subject: Add new rational number library --- libcrystfel/CMakeLists.txt | 2 + libcrystfel/src/integer_matrix.c | 107 ++++++++- libcrystfel/src/integer_matrix.h | 16 +- libcrystfel/src/rational.c | 472 +++++++++++++++++++++++++++++++++++++++ libcrystfel/src/rational.h | 94 ++++++++ 5 files changed, 688 insertions(+), 3 deletions(-) create mode 100644 libcrystfel/src/rational.c create mode 100644 libcrystfel/src/rational.h (limited to 'libcrystfel') diff --git a/libcrystfel/CMakeLists.txt b/libcrystfel/CMakeLists.txt index 247b925e..bd4f44ca 100644 --- a/libcrystfel/CMakeLists.txt +++ b/libcrystfel/CMakeLists.txt @@ -62,6 +62,7 @@ set(LIBCRYSTFEL_SOURCES src/peakfinder8.c src/taketwo.c src/xgandalf.c + src/rational.c ) if (HAVE_FFTW) @@ -101,6 +102,7 @@ set(LIBCRYSTFEL_HEADERS src/peakfinder8.h src/taketwo.h src/xgandalf.h + src/rational.h ) add_library(${PROJECT_NAME} SHARED diff --git a/libcrystfel/src/integer_matrix.c b/libcrystfel/src/integer_matrix.c index 1e616e43..6e2ade6a 100644 --- a/libcrystfel/src/integer_matrix.c +++ b/libcrystfel/src/integer_matrix.c @@ -35,7 +35,9 @@ #include #include +#include "rational.h" #include "integer_matrix.h" +#include "utils.h" /** @@ -95,7 +97,7 @@ IntegerMatrix *intmat_new(unsigned int rows, unsigned int cols) * * Returns: a newly allocated copy of @m, or NULL on error/ **/ -IntegerMatrix *intmat_copy(IntegerMatrix *m) +IntegerMatrix *intmat_copy(const IntegerMatrix *m) { IntegerMatrix *p; int i, j; @@ -127,6 +129,19 @@ void intmat_free(IntegerMatrix *m) } +void intmat_size(const IntegerMatrix *m, unsigned int *rows, unsigned int *cols) +{ + if ( m == NULL ) { + *rows = 0; + *cols = 0; + return; + } + + *rows = m->rows; + *cols = m->cols; +} + + /** * intmat_set: * @m: An %IntegerMatrix @@ -162,6 +177,96 @@ signed int intmat_get(const IntegerMatrix *m, unsigned int i, unsigned int j) } +/* intmat_intvec_mult: + * @m: An %IntegerMatrix + * @vec: An array of floats + * @ans: An array of floats in which to store the result + * + * Multiplies the matrix @m by the vector @vec. The size of @vec must equal the + * number of columns in @m, and the size of the result equals the number of rows + * in @m. + * + * Returns: non-zero on error + **/ +int intmat_floatvec_mult(const IntegerMatrix *m, const float *vec, float *ans) +{ + unsigned int i; + + for ( i=0; irows; i++ ) { + + unsigned int j; + + ans[i] = 0.0; + for ( j=0; jcols; j++ ) { + ans[i] += intmat_get(m, i, j) * vec[j]; + } + + } + + return 0; +} + + +/* intmat_rationalvec_mult: + * @m: An %IntegerMatrix + * @vec: An array of %Rational + * @ans: An array of %Rational in which to store the result + * + * Multiplies the matrix @m by the vector @vec. The size of @vec must equal the + * number of columns in @m, and the size of the result equals the number of rows + * in @m. + * + * Returns: non-zero on error + **/ +int intmat_rationalvec_mult(const IntegerMatrix *m, const Rational *vec, + Rational *ans) +{ + unsigned int i; + + + for ( i=0; irows; i++ ) { + + unsigned int j; + + ans[i] = rtnl_zero(); + for ( j=0; jcols; j++ ) { + Rational t; + t = rtnl_mul(vec[j], rtnl(intmat_get(m, i, j), 1)); + ans[i] = rtnl_add(ans[i], t); + } + + } + + return 0; +} + + +/* intmat_solve_rational: + * @m: An %IntegerMatrix + * @vec: An array of %Rational + * @ans: An array of %Rational in which to store the result + * + * Solves the matrix equation m*ans = vec, where @ans and @vec are + * vectors of %Rational. + * + * This is just a convenience function which creates a %RationalMatrix out of + * @m and then calls rtnl_mtx_solve(). + * + * Returns: non-zero on error + **/ +int intmat_solve_rational(const IntegerMatrix *m, const Rational *vec, + Rational *ans) +{ + RationalMatrix *rm; + int r; + + rm = rtnl_mtx_from_intmat(m); + r = rtnl_mtx_solve(rm, vec, ans); + rtnl_mtx_free(rm); + return r; +} + + /** * intmat_intvec_mult: * @m: An %IntegerMatrix diff --git a/libcrystfel/src/integer_matrix.h b/libcrystfel/src/integer_matrix.h index c4983aef..aeea38dd 100644 --- a/libcrystfel/src/integer_matrix.h +++ b/libcrystfel/src/integer_matrix.h @@ -40,17 +40,22 @@ **/ typedef struct _integermatrix IntegerMatrix; +#include "rational.h" + #ifdef __cplusplus extern "C" { #endif /* Alloc/dealloc */ extern IntegerMatrix *intmat_new(unsigned int rows, unsigned int cols); -extern IntegerMatrix *intmat_copy(IntegerMatrix *m); +extern IntegerMatrix *intmat_copy(const IntegerMatrix *m); extern IntegerMatrix *intmat_identity(int size); extern void intmat_free(IntegerMatrix *m); /* Get/set */ +extern void intmat_size(const IntegerMatrix *m, unsigned int *rows, + unsigned int *cols); + extern void intmat_set(IntegerMatrix *m, unsigned int i, unsigned int j, signed int v); extern signed int intmat_get(const IntegerMatrix *m, @@ -65,14 +70,21 @@ extern IntegerMatrix *intmat_create_3x3(signed int m11, signed int m12, signed i signed int m21, signed int m22, signed int m23, signed int m31, signed int m32, signed int m33); -/* Matrix-(int)vector multiplication */ +/* Matrix-vector multiplication */ +extern int intmat_floatvec_mult(const IntegerMatrix *m, const float *vec, + float *ans); extern signed int *intmat_intvec_mult(const IntegerMatrix *m, const signed int *vec); +extern int intmat_rationalvec_mult(const IntegerMatrix *m, const Rational *vec, + Rational *ans); /* Matrix-matrix multiplication */ extern IntegerMatrix *intmat_intmat_mult(const IntegerMatrix *a, const IntegerMatrix *b); +extern int intmat_solve_rational(const IntegerMatrix *m, const Rational *vec, + Rational *ans); + /* Inverse */ extern IntegerMatrix *intmat_inverse(const IntegerMatrix *m); diff --git a/libcrystfel/src/rational.c b/libcrystfel/src/rational.c new file mode 100644 index 00000000..bc783711 --- /dev/null +++ b/libcrystfel/src/rational.c @@ -0,0 +1,472 @@ +/* + * rational.c + * + * A small rational number library + * + * 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 "rational.h" +#include "integer_matrix.h" +#include "utils.h" + + +/** + * SECTION:rational + * @short_description: Rational numbers + * @title: Rational numbers + * @section_id: + * @see_also: + * @include: "rational.h" + * @Image: + * + * A rational number library + */ + + +/* Eucliden algorithm for finding greatest common divisor */ +static signed int gcd(signed int a, signed int b) +{ + while ( b != 0 ) { + signed int t = b; + b = a % b; + a = t; + } + return a; +} + + +static void squish(Rational *rt) +{ + signed int g; + + if ( rt->num == 0 ) { + rt->den = 1; + return; + } + + g = gcd(rt->num, rt->den); + assert(g != 0); + rt->num /= g; + rt->den /= g; + + if ( rt->den < 0 ) { + rt->num = -rt->num; + rt->den = -rt->den; + } +} + + +Rational rtnl_zero() +{ + Rational r; + r.num = 0; + r.den = 1; + return r; +} + + +Rational rtnl(signed long long int num, signed long long int den) +{ + Rational r; + r.num = num; + r.den = den; + squish(&r); + return r; +} + + +double rtnl_as_double(Rational r) +{ + return (double)r.num/r.den; +} + + +static void overflow(long long int c, long long int a, long long int b) +{ + setlocale(LC_ALL, ""); + ERROR("Overflow detected in rational number library.\n"); + ERROR("%'lli < %'lli * %'lli\n", c, a, b); + abort(); +} + + +static void check_overflow(long long int c, long long int a, long long int b) +{ + if ( (a==0) || (b==0) ) { + if ( c != 0 ) overflow(c,a,b); + } else if ( (llabs(c) < llabs(a)) || (llabs(c) < llabs(b)) ) { + overflow(c,a,b); + } +} + + +Rational rtnl_mul(Rational a, Rational b) +{ + Rational r; + r.num = a.num * b.num; + r.den = a.den * b.den; + check_overflow(r.num, a.num, b.num); + check_overflow(r.den, a.den, b.den); + squish(&r); + return r; +} + + +Rational rtnl_div(Rational a, Rational b) +{ + signed int t = b.num; + b.num = b.den; + b.den = t; + return rtnl_mul(a, b); +} + + +Rational rtnl_add(Rational a, Rational b) +{ + Rational r, trt1, trt2; + + trt1.num = a.num * b.den; + trt2.num = b.num * a.den; + check_overflow(trt1.num, a.num, b.den); + check_overflow(trt2.num, b.num, a.den); + + trt1.den = a.den * b.den; + trt2.den = trt1.den; + check_overflow(trt1.den, a.den, b.den); + + r.num = trt1.num + trt2.num; + r.den = trt1.den; + squish(&r); + return r; +} + + +Rational rtnl_sub(Rational a, Rational b) +{ + b.num = -b.num; + return rtnl_add(a, b); +} + + +/* -1, 0 +1 respectively for ab */ +signed int rtnl_cmp(Rational a, Rational b) +{ + Rational trt1, trt2; + + trt1.num = a.num * b.den; + trt2.num = b.num * a.den; + + trt1.den = a.den * b.den; + trt2.den = a.den * b.den; + + if ( trt1.num > trt2.num ) return +1; + if ( trt1.num < trt2.num ) return -1; + return 0; +} + + +Rational rtnl_abs(Rational a) +{ + Rational r = a; + squish(&r); + if ( r.num < 0 ) r.num = -r.num; + return r; +} + + +/** + * rtnl_format + * @rt: A %Rational + * + * Formats @rt as a string + * + * Returns: a string which should be freed by the caller + */ +char *rtnl_format(Rational rt) +{ + char *v = malloc(32); + if ( v == NULL ) return NULL; + if ( rt.den == 1 ) { + snprintf(v, 31, "%lli", rt.num); + } else { + snprintf(v, 31, "%lli/%lli", rt.num, rt.den); + } + return v; +} + + +/** + * SECTION:rational_matrix + * @short_description: Rational matrices + * @title: Rational matrices + * @section_id: + * @see_also: + * @include: "rational.h" + * @Image: + * + * A rational matrix library + */ + + +struct _rationalmatrix +{ + unsigned int rows; + unsigned int cols; + Rational *v; +}; + + +/** + * rtnl_mtx_new: + * @rows: Number of rows that the new matrix is to have + * @cols: Number of columns that the new matrix is to have + * + * Allocates a new %RationalMatrix with all elements set to zero. + * + * Returns: a new %RationalMatrix, or NULL on error. + **/ +RationalMatrix *rtnl_mtx_new(unsigned int rows, unsigned int cols) +{ + RationalMatrix *m; + + m = malloc(sizeof(RationalMatrix)); + if ( m == NULL ) return NULL; + + m->v = calloc(rows*cols, sizeof(Rational)); + if ( m->v == NULL ) { + free(m); + return NULL; + } + + m->rows = rows; + m->cols = cols; + + return m; +} + + +RationalMatrix *rtnl_mtx_copy(const RationalMatrix *m) +{ + RationalMatrix *n; + int i; + + n = rtnl_mtx_new(m->rows, m->cols); + if ( n == NULL ) return NULL; + + for ( i=0; irows*m->cols; i++ ) { + n->v[i] = m->v[i]; + } + + return n; +} + + +Rational rtnl_mtx_get(const RationalMatrix *m, int i, int j) +{ + assert(m != NULL); + return m->v[j+m->cols*i]; +} + + +void rtnl_mtx_set(const RationalMatrix *m, int i, int j, Rational v) +{ + assert(m != NULL); + m->v[j+m->cols*i] = v; +} + + +RationalMatrix *rtnl_mtx_from_intmat(const IntegerMatrix *m) +{ + RationalMatrix *n; + unsigned int rows, cols; + int i, j; + + intmat_size(m, &rows, &cols); + n = rtnl_mtx_new(rows, cols); + if ( n == NULL ) return NULL; + + for ( i=0; iv[j+cols*i] = rtnl(intmat_get(m, i, j), 1); + } + } + + return n; +} + + +void rtnl_mtx_free(RationalMatrix *mtx) +{ + if ( mtx == NULL ) return; + free(mtx->v); + free(mtx); +} + + +/* rtnl_mtx_solve: + * @m: A %RationalMatrix + * @vec: An array of %Rational + * @ans: An array of %Rational in which to store the result + * + * Solves the matrix equation m*ans = vec, where @ans and @vec are + * vectors of %Rational. + * + * In this version, @m must be square. + * + * The number of columns in @m must equal the length of @ans, and the number + * of rows in @m must equal the length of @vec, but note that there is no way + * for this function to check that this is the case. + * + * Returns: non-zero on error + **/ +int rtnl_mtx_solve(const RationalMatrix *m, const Rational *ivec, Rational *ans) +{ + RationalMatrix *cm; + Rational *vec; + int h, k; + int i; + + if ( m->rows != m->cols ) return 1; + + /* Copy the matrix and vector because the calculation will + * be done in-place */ + cm = rtnl_mtx_copy(m); + if ( cm == NULL ) return 1; + + vec = malloc(m->rows*sizeof(Rational)); + if ( vec == NULL ) return 1; + for ( h=0; hrows; h++ ) vec[h] = ivec[h]; + + /* Gaussian elimination with partial pivoting */ + h = 0; + k = 0; + while ( h<=m->rows && k<=m->cols ) { + + int prow = 0; + Rational pval = rtnl_zero(); + Rational t; + + /* Find the row with the largest value in column k */ + for ( i=h; irows; i++ ) { + if ( rtnl_cmp(rtnl_abs(rtnl_mtx_get(cm, i, k)), pval) > 0 ) { + pval = rtnl_abs(rtnl_mtx_get(cm, i, k)); + prow = i; + } + } + + if ( rtnl_cmp(pval, rtnl_zero()) == 0 ) { + k++; + continue; + } + + /* Swap 'prow' with row h */ + for ( i=0; icols; i++ ) { + t = rtnl_mtx_get(cm, h, i); + rtnl_mtx_set(cm, h, i, rtnl_mtx_get(cm, prow, i)); + rtnl_mtx_set(cm, prow, i, t); + } + t = vec[prow]; + vec[prow] = vec[h]; + vec[h] = t; + + /* Divide and subtract rows below */ + for ( i=h+1; irows; i++ ) { + + int j; + Rational dval; + + dval = rtnl_div(rtnl_mtx_get(cm, i, k), + rtnl_mtx_get(cm, h, k)); + + for ( j=0; jcols; j++ ) { + Rational t = rtnl_mtx_get(cm, i, j); + Rational p = rtnl_mul(dval, rtnl_mtx_get(cm, h, j)); + t = rtnl_sub(t, p); + rtnl_mtx_set(cm, i, j, t); + } + + /* Divide the right hand side as well */ + Rational t = vec[i]; + Rational p = rtnl_mul(dval, vec[h]); + vec[i] = rtnl_sub(t, p); + } + + h++; + k++; + + + } + + /* Back-substitution */ + for ( i=m->rows-1; i>=0; i-- ) { + int j; + Rational sum = rtnl_zero(); + for ( j=i+1; jcols; j++ ) { + Rational av; + av = rtnl_mul(rtnl_mtx_get(cm, i, j), ans[j]); + sum = rtnl_add(sum, av); + } + sum = rtnl_sub(vec[i], sum); + ans[i] = rtnl_div(sum, rtnl_mtx_get(cm, i, i)); + } + + free(vec); + rtnl_mtx_free(cm); + + return 0; +} + + +/** + * rtnl_mtx_print + * @m: A %RationalMatrix + * + * Prints @m to stderr. + * + */ +void rtnl_mtx_print(const RationalMatrix *m) +{ + unsigned int i, j; + + for ( i=0; irows; i++ ) { + + fprintf(stderr, "[ "); + for ( j=0; jcols; j++ ) { + char *v = rtnl_format(m->v[j+m->cols*i]); + fprintf(stderr, "%4s ", v); + free(v); + } + fprintf(stderr, "]\n"); + } +} diff --git a/libcrystfel/src/rational.h b/libcrystfel/src/rational.h new file mode 100644 index 00000000..8b182163 --- /dev/null +++ b/libcrystfel/src/rational.h @@ -0,0 +1,94 @@ +/* + * rational.h + * + * A small rational number library + * + * 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 . + * + */ + +#ifndef RATIONAL_H +#define RATIONAL_H + +#ifdef HAVE_CONFIG_H +#include +#endif + +/** + * Rational + * + * The Rational is an opaque-ish data structure representing a rational number. + * + * "Opaque-ish" means that the structure isn't technically opaque, allowing you + * to assign and allocate them easily. But you shouldn't look at or set its + * contents, except by using the accessor functions. + **/ +typedef struct { + /* Private, don't modify */ + signed long long int num; + signed long long int den; +} Rational; + + +/** + * RationalMatrix + * + * The RationalMatrix is an opaque data structure representing a matrix of + * rational numbers. + **/ +typedef struct _rationalmatrix RationalMatrix; + +#include "integer_matrix.h" + +#ifdef __cplusplus +extern "C" { +#endif + +extern Rational rtnl_zero(void); +extern Rational rtnl(signed long long int num, signed long long int den); +extern double rtnl_as_double(Rational r); + +extern Rational rtnl_mul(Rational a, Rational b); +extern Rational rtnl_div(Rational a, Rational b); +extern Rational rtnl_add(Rational a, Rational b); +extern Rational rtnl_sub(Rational a, Rational b); + +extern signed int rtnl_cmp(Rational a, Rational b); +extern Rational rtnl_abs(Rational a); + +extern char *rtnl_format(Rational rt); + +extern RationalMatrix *rtnl_mtx_new(unsigned int rows, unsigned int cols); +extern RationalMatrix *rtnl_mtx_copy(const RationalMatrix *m); +extern Rational rtnl_mtx_get(const RationalMatrix *m, int i, int j); +extern void rtnl_mtx_set(const RationalMatrix *m, int i, int j, Rational v); +extern RationalMatrix *rtnl_mtx_from_intmat(const IntegerMatrix *m); +extern void rtnl_mtx_free(RationalMatrix *mtx); +extern int rtnl_mtx_solve(const RationalMatrix *m, const Rational *vec, + Rational *ans); +extern void rtnl_mtx_print(const RationalMatrix *m); + +#ifdef __cplusplus +} +#endif + +#endif /* RATIONAL_H */ -- cgit v1.2.3 From 24ce9e4ac098d7744fb23f535eb97f3375425fd8 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Mon, 11 Feb 2019 18:15:32 +0100 Subject: Symmetry operation parser using Flex/Bison --- libcrystfel/CMakeLists.txt | 10 +++++ libcrystfel/src/symmetry.c | 35 ++++++++------- libcrystfel/src/symop.l | 51 ++++++++++++++++++++++ libcrystfel/src/symop.y | 106 +++++++++++++++++++++++++++++++++++++++++++++ 4 files changed, 184 insertions(+), 18 deletions(-) create mode 100644 libcrystfel/src/symop.l create mode 100644 libcrystfel/src/symop.y (limited to 'libcrystfel') diff --git a/libcrystfel/CMakeLists.txt b/libcrystfel/CMakeLists.txt index bd4f44ca..014d2a5a 100644 --- a/libcrystfel/CMakeLists.txt +++ b/libcrystfel/CMakeLists.txt @@ -7,6 +7,8 @@ find_package(PINKINDEXER) find_package(NBP) find_package(FDIP) find_package(ZLIB REQUIRED) +find_package(FLEX REQUIRED) +find_package(BISON REQUIRED) pkg_search_module(FFTW fftw3) set(HAVE_CURSES ${CURSES_FOUND}) @@ -30,6 +32,12 @@ endif() configure_file(config.h.cmake.in config.h) +bison_target(symopp src/symop.y ${CMAKE_CURRENT_BINARY_DIR}/symop-parse.c COMPILE_FLAGS --report=all) +flex_target(symopl src/symop.l ${CMAKE_CURRENT_BINARY_DIR}/symop-lex.c + DEFINES_FILE ${CMAKE_CURRENT_BINARY_DIR}/symop-lex.h) +add_flex_bison_dependency(symopl symopp) +include_directories(${PROJECT_SOURCE_DIR}/src) + set(LIBCRYSTFEL_SOURCES src/reflist.c src/utils.c @@ -63,6 +71,8 @@ set(LIBCRYSTFEL_SOURCES src/taketwo.c src/xgandalf.c src/rational.c + ${BISON_symopp_OUTPUTS} + ${FLEX_symopl_OUTPUTS} ) if (HAVE_FFTW) diff --git a/libcrystfel/src/symmetry.c b/libcrystfel/src/symmetry.c index 06cc368c..19c754ea 100644 --- a/libcrystfel/src/symmetry.c +++ b/libcrystfel/src/symmetry.c @@ -41,6 +41,8 @@ #include "symmetry.h" #include "utils.h" #include "integer_matrix.h" +#include "symop-parse.h" +#include "symop-lex.h" /** @@ -1713,28 +1715,25 @@ static IntegerMatrix *parse_symmetry_operation(const char *s) SymOpList *parse_symmetry_operations(const char *s) { - SymOpList *sol; - char **ops; - int n, i; + YY_BUFFER_STATE b; + RationalMatrix *m; + int r; - sol = new_symoplist(); - if ( sol == NULL ) return NULL; + m = rtnl_mtx_new(3, 3); + b = symop_scan_string(s); + r = symopparse(m); + symop_delete_buffer(b); - n = assplode(s, ";:", &ops, ASSPLODE_NONE); - for ( i=0; i + * + * 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 . + * + */ + +%{ + #define YYDEBUG 1 + #include "rational.h" + #include "symop-parse.h" +%} + +%option prefix="symop" +%option noyywrap nounput noinput + +%% + +[,] { return COMMA; } +[0-9]+ { symoplval.n = atoi(yytext); return NUMBER; } +[/] { return DIVIDE; } +[+] { return PLUS; } +[-] { return MINUS; } +[ahx] { return H; } +[bky] { return K; } +[clz] { return L; } +[(] { return OPENB; } +[)] { return CLOSEB; } + +%% diff --git a/libcrystfel/src/symop.y b/libcrystfel/src/symop.y new file mode 100644 index 00000000..4aec791a --- /dev/null +++ b/libcrystfel/src/symop.y @@ -0,0 +1,106 @@ +/* + * symop.y + * + * Parser for symmetry operations + * + * 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 . + * + */ + +%{ + #include + + #include "rational.h" + + extern int symoplex(); + extern int symopparse(RationalMatrix *m); + void symoperror(RationalMatrix *m, const char *s); +%} + +%define api.prefix {symop} + +%parse-param {RationalMatrix *m} + +%union { + RationalMatrix *m; /* Full rational matrix */ + Rational rv[3]; /* Rational vector, e.g. '1/2h+3k' */ + Rational r; /* Rational number */ + int n; /* Just a number */ +} + +%token COMMA +%token NUMBER +%token OPENB CLOSEB +%token H K L + +%left PLUS MINUS +%left DIVIDE +%precedence MUL +%precedence NEG + +%type symop +%type axexpr +%type part +%type NUMBER +%type fraction + +%% + +symop: + axexpr COMMA axexpr COMMA axexpr { rtnl_mtx_set(m, 0, 0, $1[0]); + rtnl_mtx_set(m, 0, 1, $1[1]); + rtnl_mtx_set(m, 0, 2, $1[2]); + rtnl_mtx_set(m, 1, 0, $3[0]); + rtnl_mtx_set(m, 1, 1, $3[1]); + rtnl_mtx_set(m, 1, 2, $3[2]); + rtnl_mtx_set(m, 2, 0, $5[0]); + rtnl_mtx_set(m, 2, 1, $5[1]); + rtnl_mtx_set(m, 2, 2, $5[2]); + } +; + +axexpr: + part { int i; for ( i=0; i<3; i++ ) $$[i] = $1[i]; } +| axexpr PLUS axexpr { int i; for ( i=0; i<3; i++ ) $$[i] = rtnl_add($1[i], $3[i]); } +| axexpr MINUS axexpr { int i; for ( i=0; i<3; i++ ) $$[i] = rtnl_sub($1[i], $3[i]); } +| MINUS axexpr %prec NEG { int i; for ( i=0; i<3; i++ ) $$[i] = rtnl_sub(rtnl_zero(), $2[i]); } +| OPENB axexpr CLOSEB { int i; for ( i=0; i<3; i++ ) $$[i] = $2[i]; } +| axexpr DIVIDE NUMBER { int i; for ( i=0; i<3; i++ ) $$[i] = rtnl_div($1[i], rtnl($3, 1)); } +| NUMBER axexpr %prec MUL { int i; for ( i=0; i<3; i++ ) $$[i] = rtnl_mul($2[i], rtnl($1, 1)); } +| fraction axexpr %prec MUL { int i; for ( i=0; i<3; i++ ) $$[i] = rtnl_mul($2[i], $1); } +; + +part: + H { $$[0] = rtnl(1, 1); $$[1] = rtnl_zero(); $$[2] = rtnl_zero(); } +| K { $$[1] = rtnl(1, 1); $$[0] = rtnl_zero(); $$[2] = rtnl_zero(); } +| L { $$[2] = rtnl(1, 1); $$[0] = rtnl_zero(); $$[1] = rtnl_zero(); } +; + +fraction: + NUMBER DIVIDE NUMBER { $$ = rtnl($1, $3); } +; + +%% + +void symoperror(RationalMatrix *m, const char *s) { + printf("Error\n"); +} -- cgit v1.2.3 From a99bb041b12d6f10a00e75a4d76083767199a7a7 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Mon, 18 Feb 2019 15:35:01 +0100 Subject: Add rtnl_mtx_det() and rtnl_mtx_mult() --- libcrystfel/src/rational.c | 91 ++++++++++++++++++++++++++++++++++++++++++++++ libcrystfel/src/rational.h | 3 ++ 2 files changed, 94 insertions(+) (limited to 'libcrystfel') diff --git a/libcrystfel/src/rational.c b/libcrystfel/src/rational.c index bc783711..b6ca5dd4 100644 --- a/libcrystfel/src/rational.c +++ b/libcrystfel/src/rational.c @@ -470,3 +470,94 @@ void rtnl_mtx_print(const RationalMatrix *m) fprintf(stderr, "]\n"); } } + + +void rtnl_mtx_mult(const RationalMatrix *m, const Rational *vec, Rational *ans) +{ + int i, j; + + for ( i=0; irows; i++ ) { + ans[i] = rtnl_zero(); + for ( j=0; jcols; j++ ) { + Rational add; + add = rtnl_mul(rtnl_mtx_get(m, i, j), vec[j]); + ans[i] = rtnl_add(ans[i], add); + } + } +} + + +static RationalMatrix *delete_row_and_column(const RationalMatrix *m, + unsigned int di, unsigned int dj) +{ + RationalMatrix *n; + unsigned int i, j; + + n = rtnl_mtx_new(m->rows-1, m->cols-1); + if ( n == NULL ) return NULL; + + for ( i=0; irows; i++ ) { + for ( j=0; jcols; j++ ) { + + Rational val; + unsigned int gi, gj; + + gi = (i>=di) ? i+1 : i; + gj = (j>=dj) ? j+1 : j; + val = rtnl_mtx_get(m, gi, gj); + rtnl_mtx_set(n, i, j, val); + + } + } + + return n; +} + + +static Rational cofactor(const RationalMatrix *m, + unsigned int i, unsigned int j) +{ + RationalMatrix *n; + Rational t, C; + + n = delete_row_and_column(m, i, j); + if ( n == NULL ) { + fprintf(stderr, "Failed to allocate matrix.\n"); + return rtnl_zero(); + } + + /* -1 if odd, +1 if even */ + t = (i+j) & 0x1 ? rtnl(-1, 1) : rtnl(1, 1); + + C = rtnl_mul(t, rtnl_mtx_det(n)); + rtnl_mtx_free(n); + + return C; +} + + + +Rational rtnl_mtx_det(const RationalMatrix *m) +{ + unsigned int i, j; + Rational det; + + assert(m->rows == m->cols); /* Otherwise determinant doesn't exist */ + + if ( m->rows == 2 ) { + Rational a, b; + a = rtnl_mul(rtnl_mtx_get(m, 0, 0), rtnl_mtx_get(m, 1, 1)); + b = rtnl_mul(rtnl_mtx_get(m, 0, 1), rtnl_mtx_get(m, 1, 0)); + return rtnl_sub(a, b); + } + + i = 0; /* Fixed */ + det = rtnl_zero(); + for ( j=0; jcols; j++ ) { + Rational a; + a = rtnl_mul(rtnl_mtx_get(m, i, j), cofactor(m, i, j)); + det = rtnl_add(det, a); + } + + return det; +} diff --git a/libcrystfel/src/rational.h b/libcrystfel/src/rational.h index 8b182163..e2d3e5bf 100644 --- a/libcrystfel/src/rational.h +++ b/libcrystfel/src/rational.h @@ -83,9 +83,12 @@ extern Rational rtnl_mtx_get(const RationalMatrix *m, int i, int j); extern void rtnl_mtx_set(const RationalMatrix *m, int i, int j, Rational v); extern RationalMatrix *rtnl_mtx_from_intmat(const IntegerMatrix *m); extern void rtnl_mtx_free(RationalMatrix *mtx); +extern void rtnl_mtx_mult(const RationalMatrix *m, const Rational *vec, + Rational *ans); extern int rtnl_mtx_solve(const RationalMatrix *m, const Rational *vec, Rational *ans); extern void rtnl_mtx_print(const RationalMatrix *m); +extern Rational rtnl_mtx_det(const RationalMatrix *m); #ifdef __cplusplus } -- cgit v1.2.3 From 169e7c5677ffc9c296c0a7eeddb0b77e024a4a55 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Mon, 18 Feb 2019 16:04:28 +0100 Subject: Be clear about whether functions need rational or integer transformations --- libcrystfel/src/cell-utils.c | 8 ++-- libcrystfel/src/cell.c | 88 ++++++++++++++++++++++++++++++++------------ libcrystfel/src/cell.h | 8 +++- libcrystfel/src/xgandalf.c | 2 +- 4 files changed, 75 insertions(+), 31 deletions(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c index 507d4cc2..da1e6fd5 100644 --- a/libcrystfel/src/cell-utils.c +++ b/libcrystfel/src/cell-utils.c @@ -500,7 +500,7 @@ UnitCell *uncenter_cell(UnitCell *in, IntegerMatrix **t) tt = centering_transformation(in, &new_centering, &new_latt, &new_ua); if ( tt == NULL ) return NULL; - out = cell_transform_inverse(in, tt); + out = cell_transform_intmat_inverse(in, tt); if ( out == NULL ) return NULL; cell_set_lattice_type(out, new_latt); @@ -776,7 +776,7 @@ UnitCell *match_cell(UnitCell *cell_in, UnitCell *template_in, int verbose, /* Reverse the de-centering transformation */ if ( new_cell != NULL ) { - new_cell_trans = cell_transform(new_cell, centering); + 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)); @@ -905,7 +905,7 @@ UnitCell *match_cell_ab(UnitCell *cell_in, UnitCell *template_in) new_cell = cell_new_from_direct_axes(real_a, real_b, real_c); /* Reverse the de-centering transformation */ - new_cell_trans = cell_transform_inverse(new_cell, to_given_cell); + 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)); @@ -1770,7 +1770,7 @@ int compare_reindexed_cell_parameters_and_orientation(UnitCell *a, UnitCell *b, if ( intmat_det(m) != +1 ) continue; - nc = cell_transform(b, m); + nc = cell_transform_intmat(b, m); if ( compare_cell_parameters_and_orientation(a, nc, ltl, atl) ) { if ( pmb != NULL ) *pmb = m; diff --git a/libcrystfel/src/cell.c b/libcrystfel/src/cell.c index 242094f6..5eb13349 100644 --- a/libcrystfel/src/cell.c +++ b/libcrystfel/src/cell.c @@ -688,16 +688,16 @@ UnitCell *cell_transform_gsl_reciprocal(UnitCell *in, gsl_matrix *m) /** - * cell_transform: + * cell_transform_rational: * @cell: A %UnitCell. - * @t: An %IntegerMatrix. + * @t: A %RationalMatrix. * * Applies @t to @cell. * * Returns: Transformed copy of @cell. * */ -UnitCell *cell_transform(UnitCell *cell, IntegerMatrix *m) +UnitCell *cell_transform_rational(UnitCell *cell, RationalMatrix *m) { UnitCell *out; gsl_matrix *tm; @@ -709,15 +709,15 @@ UnitCell *cell_transform(UnitCell *cell, IntegerMatrix *m) return NULL; } - gsl_matrix_set(tm, 0, 0, intmat_get(m, 0, 0)); - gsl_matrix_set(tm, 0, 1, intmat_get(m, 0, 1)); - gsl_matrix_set(tm, 0, 2, intmat_get(m, 0, 2)); - gsl_matrix_set(tm, 1, 0, intmat_get(m, 1, 0)); - gsl_matrix_set(tm, 1, 1, intmat_get(m, 1, 1)); - gsl_matrix_set(tm, 1, 2, intmat_get(m, 1, 2)); - gsl_matrix_set(tm, 2, 0, intmat_get(m, 2, 0)); - gsl_matrix_set(tm, 2, 1, intmat_get(m, 2, 1)); - gsl_matrix_set(tm, 2, 2, intmat_get(m, 2, 2)); + gsl_matrix_set(tm, 0, 0, rtnl_as_double(rtnl_mtx_get(m, 0, 0))); + gsl_matrix_set(tm, 0, 1, rtnl_as_double(rtnl_mtx_get(m, 0, 1))); + gsl_matrix_set(tm, 0, 2, rtnl_as_double(rtnl_mtx_get(m, 0, 2))); + gsl_matrix_set(tm, 1, 0, rtnl_as_double(rtnl_mtx_get(m, 1, 0))); + gsl_matrix_set(tm, 1, 1, rtnl_as_double(rtnl_mtx_get(m, 1, 1))); + gsl_matrix_set(tm, 1, 2, rtnl_as_double(rtnl_mtx_get(m, 1, 2))); + gsl_matrix_set(tm, 2, 0, rtnl_as_double(rtnl_mtx_get(m, 2, 0))); + gsl_matrix_set(tm, 2, 1, rtnl_as_double(rtnl_mtx_get(m, 2, 1))); + gsl_matrix_set(tm, 2, 2, rtnl_as_double(rtnl_mtx_get(m, 2, 2))); out = cell_transform_gsl_direct(cell, tm); @@ -730,16 +730,36 @@ UnitCell *cell_transform(UnitCell *cell, IntegerMatrix *m) /** - * cell_transform_inverse: + * cell_transform_intmat: * @cell: A %UnitCell. - * @m: An %IntegerMatrix + * @t: An %IntegerMatrix. + * + * Applies @t to @cell. + * + * Returns: Transformed copy of @cell. + * + */ +UnitCell *cell_transform_intmat(UnitCell *cell, IntegerMatrix *m) +{ + UnitCell *ans; + RationalMatrix *mtx = rtnl_mtx_from_intmat(m); + ans = cell_transform_rational(cell, mtx); + rtnl_mtx_free(mtx); + return ans; +} + + +/** + * cell_transform_rational_inverse: + * @cell: A %UnitCell. + * @m: A %RationalMatrix * * Applies the inverse of @m to @cell. * * Returns: Transformed copy of @cell. * */ -UnitCell *cell_transform_inverse(UnitCell *cell, IntegerMatrix *m) +UnitCell *cell_transform_rational_inverse(UnitCell *cell, RationalMatrix *m) { UnitCell *out; gsl_matrix *tm; @@ -754,15 +774,15 @@ UnitCell *cell_transform_inverse(UnitCell *cell, IntegerMatrix *m) return NULL; } - gsl_matrix_set(tm, 0, 0, intmat_get(m, 0, 0)); - gsl_matrix_set(tm, 0, 1, intmat_get(m, 0, 1)); - gsl_matrix_set(tm, 0, 2, intmat_get(m, 0, 2)); - gsl_matrix_set(tm, 1, 0, intmat_get(m, 1, 0)); - gsl_matrix_set(tm, 1, 1, intmat_get(m, 1, 1)); - gsl_matrix_set(tm, 1, 2, intmat_get(m, 1, 2)); - gsl_matrix_set(tm, 2, 0, intmat_get(m, 2, 0)); - gsl_matrix_set(tm, 2, 1, intmat_get(m, 2, 1)); - gsl_matrix_set(tm, 2, 2, intmat_get(m, 2, 2)); + gsl_matrix_set(tm, 0, 0, rtnl_as_double(rtnl_mtx_get(m, 0, 0))); + gsl_matrix_set(tm, 0, 1, rtnl_as_double(rtnl_mtx_get(m, 0, 1))); + gsl_matrix_set(tm, 0, 2, rtnl_as_double(rtnl_mtx_get(m, 0, 2))); + gsl_matrix_set(tm, 1, 0, rtnl_as_double(rtnl_mtx_get(m, 1, 0))); + gsl_matrix_set(tm, 1, 1, rtnl_as_double(rtnl_mtx_get(m, 1, 1))); + gsl_matrix_set(tm, 1, 2, rtnl_as_double(rtnl_mtx_get(m, 1, 2))); + gsl_matrix_set(tm, 2, 0, rtnl_as_double(rtnl_mtx_get(m, 2, 0))); + gsl_matrix_set(tm, 2, 1, rtnl_as_double(rtnl_mtx_get(m, 2, 1))); + gsl_matrix_set(tm, 2, 2, rtnl_as_double(rtnl_mtx_get(m, 2, 2))); perm = gsl_permutation_alloc(3); if ( perm == NULL ) { @@ -796,3 +816,23 @@ UnitCell *cell_transform_inverse(UnitCell *cell, IntegerMatrix *m) return out; } + + +/** + * cell_transform_intmat_inverse: + * @cell: A %UnitCell. + * @m: An %IntegerMatrix + * + * Applies the inverse of @m to @cell. + * + * Returns: Transformed copy of @cell. + * + */ +UnitCell *cell_transform_intmat_inverse(UnitCell *cell, IntegerMatrix *m) +{ + UnitCell *ans; + RationalMatrix *mtx = rtnl_mtx_from_intmat(m); + ans = cell_transform_rational_inverse(cell, mtx); + rtnl_mtx_free(mtx); + return ans; +} diff --git a/libcrystfel/src/cell.h b/libcrystfel/src/cell.h index e6ccbdf7..8e6809b5 100644 --- a/libcrystfel/src/cell.h +++ b/libcrystfel/src/cell.h @@ -150,8 +150,12 @@ extern const char *cell_rep(UnitCell *cell); extern UnitCell *cell_transform_gsl_direct(UnitCell *in, gsl_matrix *m); extern UnitCell *cell_transform_gsl_reciprocal(UnitCell *in, gsl_matrix *m); -extern UnitCell *cell_transform(UnitCell *cell, IntegerMatrix *m); -extern UnitCell *cell_transform_inverse(UnitCell *cell, IntegerMatrix *m); + +extern UnitCell *cell_transform_rational(UnitCell *cell, RationalMatrix *m); +extern UnitCell *cell_transform_rational_inverse(UnitCell *cell, RationalMatrix *m); + +extern UnitCell *cell_transform_intmat(UnitCell *cell, IntegerMatrix *m); +extern UnitCell *cell_transform_intmat_inverse(UnitCell *cell, IntegerMatrix *m); #ifdef __cplusplus } diff --git a/libcrystfel/src/xgandalf.c b/libcrystfel/src/xgandalf.c index e6618560..b122cb4c 100644 --- a/libcrystfel/src/xgandalf.c +++ b/libcrystfel/src/xgandalf.c @@ -114,7 +114,7 @@ int run_xgandalf(struct image *image, void *ipriv) if(xgandalf_private_data->cellTemplate != NULL){ restoreCell(uc, &xgandalf_private_data->latticeReductionTransform); - UnitCell *new_cell_trans = cell_transform(uc, xgandalf_private_data->centeringTransformation); + UnitCell *new_cell_trans = cell_transform_intmat(uc, xgandalf_private_data->centeringTransformation); cell_free(uc); uc = new_cell_trans; -- cgit v1.2.3 From a2f4977e0f8bd9becd50ab5a2ef903038273133c Mon Sep 17 00:00:00 2001 From: Thomas White Date: Mon, 18 Feb 2019 16:56:52 +0100 Subject: Implement parse_symmetry_operations --- libcrystfel/src/rational.c | 25 ++++++++++ libcrystfel/src/rational.h | 1 + libcrystfel/src/symmetry.c | 122 ++++++++++++++++++--------------------------- libcrystfel/src/symmetry.h | 9 ++-- libcrystfel/src/symop.l | 1 + libcrystfel/src/symop.y | 41 +++++++++++++-- 6 files changed, 117 insertions(+), 82 deletions(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/rational.c b/libcrystfel/src/rational.c index b6ca5dd4..5db83165 100644 --- a/libcrystfel/src/rational.c +++ b/libcrystfel/src/rational.c @@ -326,6 +326,31 @@ RationalMatrix *rtnl_mtx_from_intmat(const IntegerMatrix *m) } +IntegerMatrix *intmat_from_rtnl_mtx(const RationalMatrix *m) +{ + IntegerMatrix *n; + int i, j; + + n = intmat_new(m->rows, m->cols); + if ( n == NULL ) return NULL; + + for ( i=0; irows; i++ ) { + for ( j=0; jcols; j++ ) { + Rational v = rtnl_mtx_get(m, i, j); + squish(&v); + if ( v.den != 1 ) { + ERROR("Rational matrix can't be converted to integers\n"); + intmat_free(n); + return NULL; + } + intmat_set(n, i, j, v.num); + } + } + + return n; +} + + void rtnl_mtx_free(RationalMatrix *mtx) { if ( mtx == NULL ) return; diff --git a/libcrystfel/src/rational.h b/libcrystfel/src/rational.h index e2d3e5bf..925f9be1 100644 --- a/libcrystfel/src/rational.h +++ b/libcrystfel/src/rational.h @@ -82,6 +82,7 @@ extern RationalMatrix *rtnl_mtx_copy(const RationalMatrix *m); extern Rational rtnl_mtx_get(const RationalMatrix *m, int i, int j); extern void rtnl_mtx_set(const RationalMatrix *m, int i, int j, Rational v); extern RationalMatrix *rtnl_mtx_from_intmat(const IntegerMatrix *m); +extern IntegerMatrix *intmat_from_rtnl_mtx(const RationalMatrix *m); extern void rtnl_mtx_free(RationalMatrix *mtx); extern void rtnl_mtx_mult(const RationalMatrix *m, const Rational *vec, Rational *ans); diff --git a/libcrystfel/src/symmetry.c b/libcrystfel/src/symmetry.c index 19c754ea..9ee80a16 100644 --- a/libcrystfel/src/symmetry.c +++ b/libcrystfel/src/symmetry.c @@ -3,12 +3,12 @@ * * Symmetry * - * Copyright © 2012-2016 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2012-2019 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * * Authors: - * 2010-2014,2016 Thomas White - * 2014 Kenneth Beyerlein + * 2010-2019 Thomas White + * 2014 Kenneth Beyerlein * * This file is part of CrystFEL. * @@ -1638,102 +1638,76 @@ SymOpList *get_ambiguities(const SymOpList *source, const SymOpList *target) } -static IntegerMatrix *parse_symmetry_operation(const char *s) +/* Parse a single symmetry operation, e.g. 'h,-2k,(h+l)/3' */ +RationalMatrix *parse_symmetry_operation(const char *s) { - IntegerMatrix *m; - char **els; - int n, i; + YY_BUFFER_STATE b; + RationalMatrix *m; + int r; + m = rtnl_mtx_new(3, 3); + b = symop_scan_string(s); + r = symopparse(m, NULL); + symop_delete_buffer(b); - n = assplode(s, ",", &els, ASSPLODE_NONE); - if ( n != 3 ) { - for ( i=0; i 0 ) { - mult *= 10; - mult += els[i][c] - '0'; - } else { - mult *= els[i][c] - '0'; - } - ndigit++; - } - } - - intmat_set(m, i, 0, nh); - intmat_set(m, i, 1, nk); - intmat_set(m, i, 2, nl); - - free(els[i]); + return m; +} - } - free(els); - return m; +/** + * parse_cell_transformation + * @s: Textual representation of cell transformation + * + * Parses @s, for example 'a,(b+a)/2,c', and returns the corresponding + * %RationalMatrix. + * + * Returns: A %RationalMatrix describing the transformation, or NULL on error. + * + */ +RationalMatrix *parse_cell_transformation(const char *s) +{ + return parse_symmetry_operation(s); } +/** + * parse_symmetry_operations + * @s: Textual representation of a list of symmetry operations + * + * Parses @s, for example 'h,k,l;k,h,-l', and returns the corresponding + * %SymOpList + * + * Returns: A %SymOpList, or NULL on error. + * + */ SymOpList *parse_symmetry_operations(const char *s) { YY_BUFFER_STATE b; RationalMatrix *m; + SymOpList *list; int r; - m = rtnl_mtx_new(3, 3); + m = rtnl_mtx_new(3, 3); /* Scratch space for parser */ + list = new_symoplist(); /* The result we want */ + b = symop_scan_string(s); - r = symopparse(m); + r = symopparse(m, list); symop_delete_buffer(b); + rtnl_mtx_free(m); if ( r ) { ERROR("Failed to parse '%s'\n", s); - rtnl_mtx_free(m); + free_symoplist(list); return NULL; } - STATUS("Parsed OK\n"); - rtnl_mtx_print(m); - - return NULL; + return list; } diff --git a/libcrystfel/src/symmetry.h b/libcrystfel/src/symmetry.h index a61ca96f..ab2aa934 100644 --- a/libcrystfel/src/symmetry.h +++ b/libcrystfel/src/symmetry.h @@ -3,12 +3,12 @@ * * Symmetry * - * Copyright © 2012-2016 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2012-2019 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * * Authors: - * 2010-2014,2016 Thomas White - * 2014 Kenneth Beyerlein + * 2010-2019 Thomas White + * 2014 Kenneth Beyerlein * * This file is part of CrystFEL. * @@ -36,6 +36,7 @@ #include "integer_matrix.h" +#include "rational.h" /** * SymOpList @@ -93,7 +94,9 @@ extern int is_centric(signed int h, signed int k, signed int l, extern void pointgroup_warning(const char *sym); extern void add_symop(SymOpList *ops, IntegerMatrix *m); +extern RationalMatrix *parse_symmetry_operation(const char *s); extern SymOpList *parse_symmetry_operations(const char *s); +extern RationalMatrix *parse_cell_transformation(const char *s); extern char *get_matrix_name(const IntegerMatrix *m, int row); #ifdef __cplusplus diff --git a/libcrystfel/src/symop.l b/libcrystfel/src/symop.l index b00b35b4..ed38a65e 100644 --- a/libcrystfel/src/symop.l +++ b/libcrystfel/src/symop.l @@ -47,5 +47,6 @@ [clz] { return L; } [(] { return OPENB; } [)] { return CLOSEB; } +[;] { return SEMICOLON; } %% diff --git a/libcrystfel/src/symop.y b/libcrystfel/src/symop.y index 4aec791a..a304ea37 100644 --- a/libcrystfel/src/symop.y +++ b/libcrystfel/src/symop.y @@ -30,15 +30,20 @@ #include #include "rational.h" + #include "symmetry.h" extern int symoplex(); - extern int symopparse(RationalMatrix *m); - void symoperror(RationalMatrix *m, const char *s); + extern int symopparse(RationalMatrix *m, SymOpList *list); + void symoperror(RationalMatrix *m, SymOpList *list, const char *s); %} %define api.prefix {symop} -%parse-param {RationalMatrix *m} +%parse-param {RationalMatrix *m} {SymOpList *list} + +%code requires { + #include "symmetry.h" +} %union { RationalMatrix *m; /* Full rational matrix */ @@ -47,6 +52,7 @@ int n; /* Just a number */ } +%token SEMICOLON %token COMMA %token NUMBER %token OPENB CLOSEB @@ -63,8 +69,33 @@ %type NUMBER %type fraction +%{ +static int try_add_symop(SymOpList *list, RationalMatrix *m) +{ + if ( list == NULL ) { + yyerror(m, list, "Must be a single symmetry operation"); + return 1; + } else { + IntegerMatrix *im; + im = intmat_from_rtnl_mtx(m); + if ( im == NULL ) { + yyerror(m, list, "Symmetry operations must all be integer"); + return 1; + } else { + add_symop(list, im); + } + } + return 0; +} +%} + %% +symoplist: + symop { if ( try_add_symop(list, m) ) YYERROR; } +| symoplist SEMICOLON symop { if ( try_add_symop(list, m) ) YYERROR; } +; + symop: axexpr COMMA axexpr COMMA axexpr { rtnl_mtx_set(m, 0, 0, $1[0]); rtnl_mtx_set(m, 0, 1, $1[1]); @@ -101,6 +132,6 @@ fraction: %% -void symoperror(RationalMatrix *m, const char *s) { - printf("Error\n"); +void symoperror(RationalMatrix *m, SymOpList *list, const char *s) { + printf("Error: %s\n", s); } -- cgit v1.2.3 From 830868c3edd603fc2466248ed363066f312c9ac8 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Mon, 18 Feb 2019 17:07:08 +0100 Subject: Better handling of multiple symop error --- libcrystfel/src/symop.y | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/symop.y b/libcrystfel/src/symop.y index a304ea37..bebf823f 100644 --- a/libcrystfel/src/symop.y +++ b/libcrystfel/src/symop.y @@ -70,10 +70,13 @@ %type fraction %{ -static int try_add_symop(SymOpList *list, RationalMatrix *m) +static int try_add_symop(SymOpList *list, RationalMatrix *m, int complain) { if ( list == NULL ) { - yyerror(m, list, "Must be a single symmetry operation"); + /* Only complain if this isn't the only operation provided */ + if ( complain ) { + yyerror(m, list, "Must be a single symmetry operation"); + } return 1; } else { IntegerMatrix *im; @@ -92,8 +95,8 @@ static int try_add_symop(SymOpList *list, RationalMatrix *m) %% symoplist: - symop { if ( try_add_symop(list, m) ) YYERROR; } -| symoplist SEMICOLON symop { if ( try_add_symop(list, m) ) YYERROR; } + symop { try_add_symop(list, m, 0); } +| symoplist SEMICOLON symop { if ( try_add_symop(list, m, 1) ) YYERROR; } ; symop: -- cgit v1.2.3 From 2dad9da2661389cf5c34a49e933dab20bd6d943c Mon Sep 17 00:00:00 2001 From: Thomas White Date: Thu, 28 Feb 2019 10:53:11 +0100 Subject: Fix tolerances (again) --- libcrystfel/src/cell-utils.c | 12 ++++++------ libcrystfel/src/cell-utils.h | 2 +- 2 files changed, 7 insertions(+), 7 deletions(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c index da1e6fd5..a21ae570 100644 --- a/libcrystfel/src/cell-utils.c +++ b/libcrystfel/src/cell-utils.c @@ -3,13 +3,13 @@ * * Unit Cell utility functions * - * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2012-2019 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * Copyright © 2012 Lorenzo Galli * * Authors: - * 2009-2012,2014-2017 Thomas White - * 2012 Lorenzo Galli + * 2009-2019 Thomas White + * 2012 Lorenzo Galli * * This file is part of CrystFEL. * @@ -1647,9 +1647,9 @@ 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) ) return 0; - if ( !within_tolerance(b1, b2, ltl) ) return 0; - if ( !within_tolerance(c1, c2, ltl) ) return 0; + 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; diff --git a/libcrystfel/src/cell-utils.h b/libcrystfel/src/cell-utils.h index 1a4f3bcc..a1a6fd46 100644 --- a/libcrystfel/src/cell-utils.h +++ b/libcrystfel/src/cell-utils.h @@ -3,7 +3,7 @@ * * Unit Cell utility functions * - * Copyright © 2012-2018 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2012-2019 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * Copyright © 2012 Lorenzo Galli * -- cgit v1.2.3 From 3e850c6c2ed2a3dbc0fe0be27d36c0d1e6b57614 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Thu, 28 Feb 2019 10:53:43 +0100 Subject: New way of doing cell comparisons, similar to match_cell() --- libcrystfel/src/cell-utils.c | 230 +++++++++++++++++++++++++++++++++++++++++++ libcrystfel/src/cell-utils.h | 4 + libcrystfel/src/rational.c | 33 +++++++ libcrystfel/src/rational.h | 4 + 4 files changed, 271 insertions(+) (limited to 'libcrystfel') diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c index a21ae570..2bfaa3e1 100644 --- a/libcrystfel/src/cell-utils.c +++ b/libcrystfel/src/cell-utils.c @@ -1793,3 +1793,233 @@ int compare_reindexed_cell_parameters_and_orientation(UnitCell *a, UnitCell *b, intmat_free(m); return 0; } + + +struct cand +{ + Rational abc[3]; + double fom; +}; + + +static int cmpcand(const void *av, const void *bv) +{ + const struct cand *a = av; + const struct cand *b = bv; + return a->fom > b->fom; +} + + +static Rational *find_candidates(double len, double *a, double *b, double *c, + double ltl, int *pncand) +{ + Rational *r; + struct cand *cands; + const int max_cand = 1024; + int ncand = 0; + Rational *rat; + int nrat; + int nrej = 0; + int ia, ib, ic; + int i; + + cands = malloc(max_cand * sizeof(struct cand)); + if ( cands == NULL ) return NULL; + + rat = rtnl_list(0, 4, -4, 4, &nrat); + if ( rat == NULL ) return NULL; + + for ( ia=0; ia atl ) continue; + + /* Gamma OK, now look for place for c axis */ + for ( ic=0; ic atl ) { + cell_free(test); + continue; + } + if ( fabs(bet - be) > atl ) { + cell_free(test); + continue; + } + rtnl_mtx_print(m); + test2 = cell_transform_intmat(test, centering_reference); + cell_print(test2); + cell_free(test); + cell_free(test2); + + } + } + } + + *pmb = m; + return 1; +} diff --git a/libcrystfel/src/cell-utils.h b/libcrystfel/src/cell-utils.h index a1a6fd46..0face43d 100644 --- a/libcrystfel/src/cell-utils.h +++ b/libcrystfel/src/cell-utils.h @@ -95,6 +95,10 @@ extern int compare_reindexed_cell_parameters_and_orientation(UnitCell *a, double atl, IntegerMatrix **pmb); +extern int compare_reindexed_cell_parameters(UnitCell *cell, UnitCell *reference, + double ltl, double atl, + RationalMatrix **pmb); + #ifdef __cplusplus } #endif diff --git a/libcrystfel/src/rational.c b/libcrystfel/src/rational.c index 5db83165..0f2bcaee 100644 --- a/libcrystfel/src/rational.c +++ b/libcrystfel/src/rational.c @@ -226,6 +226,39 @@ char *rtnl_format(Rational rt) } +Rational *rtnl_list(signed int num_min, signed int num_max, + signed int den_min, signed int den_max, + int *pn) +{ + signed int num, den; + Rational *list; + int n = 0; + + list = malloc((1+num_max-num_min)*(1+den_max-den_min)*sizeof(Rational)); + if ( list == NULL ) return NULL; + + for ( num=num_min; num<=num_max; num++ ) { + for ( den=den_min; den<=den_max; den++ ) { + + Rational r = rtnl(num, den); + + /* Denominator zero? */ + if ( den == 0 ) continue; + + /* Same as last entry? */ + if ( (n>0) && (rtnl_cmp(list[n-1], r)==0) ) continue; + + /* Can be reduced? */ + if ( gcd(num, den) != 1 ) continue; + + list[n++] = r; + } + } + *pn = n; + return list; +} + + /** * SECTION:rational_matrix * @short_description: Rational matrices diff --git a/libcrystfel/src/rational.h b/libcrystfel/src/rational.h index 925f9be1..7983f56f 100644 --- a/libcrystfel/src/rational.h +++ b/libcrystfel/src/rational.h @@ -77,6 +77,10 @@ extern Rational rtnl_abs(Rational a); extern char *rtnl_format(Rational rt); +extern Rational *rtnl_list(signed int num_min, signed int num_max, + signed int den_min, signed int den_max, + int *pn); + extern RationalMatrix *rtnl_mtx_new(unsigned int rows, unsigned int cols); extern RationalMatrix *rtnl_mtx_copy(const RationalMatrix *m); extern Rational rtnl_mtx_get(const RationalMatrix *m, int i, int j); -- cgit v1.2.3 From 68061d0e3c42f61fa7664e0f0996cade13057391 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Tue, 5 Mar 2019 14:16:26 +0100 Subject: Add rtnl_mtx_identity() --- libcrystfel/src/rational.c | 16 ++++++++++++++++ libcrystfel/src/rational.h | 1 + 2 files changed, 17 insertions(+) (limited to 'libcrystfel') diff --git a/libcrystfel/src/rational.c b/libcrystfel/src/rational.c index 0f2bcaee..58178559 100644 --- a/libcrystfel/src/rational.c +++ b/libcrystfel/src/rational.c @@ -292,6 +292,7 @@ struct _rationalmatrix RationalMatrix *rtnl_mtx_new(unsigned int rows, unsigned int cols) { RationalMatrix *m; + int i; m = malloc(sizeof(RationalMatrix)); if ( m == NULL ) return NULL; @@ -305,6 +306,21 @@ RationalMatrix *rtnl_mtx_new(unsigned int rows, unsigned int cols) m->rows = rows; m->cols = cols; + for ( i=0; irows*m->cols; i++ ) { + m->v[i] = rtnl_zero(); + } + + return m; +} + + +RationalMatrix *rtnl_mtx_identity(int rows) +{ + int i; + RationalMatrix *m = rtnl_mtx_new(rows, rows); + for ( i=0; i Date: Wed, 6 Mar 2019 11:40:37 +0100 Subject: Keep track of the "un-centering" matrix, as well as the "centering" This makes it easy to reverse the transformation, if required, which it is when comparing centered cells. --- libcrystfel/src/cell-utils.c | 157 ++++++++++++++++++++++++++++++------------- libcrystfel/src/cell-utils.h | 3 +- libcrystfel/src/rational.c | 25 +++++++ libcrystfel/src/rational.h | 2 + libcrystfel/src/xgandalf.c | 2 +- 5 files changed, 142 insertions(+), 47 deletions(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c index 2bfaa3e1..3a5ef782 100644 --- a/libcrystfel/src/cell-utils.c +++ b/libcrystfel/src/cell-utils.c @@ -361,18 +361,45 @@ int bravais_lattice(UnitCell *cell) } +static RationalMatrix *create_rtnl_mtx(signed int a1, signed int a2, + signed int b1, signed int b2, + signed int c1, signed int c2, + signed int d1, signed int d2, + signed int e1, signed int e2, + signed int f1, signed int f2, + signed int g1, signed int g2, + signed int h1, signed int h2, + signed int i1, signed int i2) +{ + RationalMatrix *m = rtnl_mtx_new(3, 3); + if ( m == NULL ) return NULL; + rtnl_mtx_set(m, 0, 0, rtnl(a1, a2)); + rtnl_mtx_set(m, 0, 1, rtnl(b1, b2)); + rtnl_mtx_set(m, 0, 2, rtnl(c1, c2)); + rtnl_mtx_set(m, 1, 0, rtnl(d1, d2)); + rtnl_mtx_set(m, 1, 1, rtnl(e1, e2)); + rtnl_mtx_set(m, 1, 2, rtnl(f1, f2)); + rtnl_mtx_set(m, 2, 0, rtnl(g1, g2)); + rtnl_mtx_set(m, 2, 1, rtnl(h1, h2)); + rtnl_mtx_set(m, 2, 2, rtnl(i1, h2)); + return m; +} + + /* Given a centered cell @in, return the integer transformation matrix which * turns a primitive cell into @in. Set new_centering and new_latt to the * centering and lattice type of the primitive cell (usually aP, sometimes rR, - * rarely mP) */ + * rarely mP). Store the inverse matrix at pCi */ static IntegerMatrix *centering_transformation(UnitCell *in, char *new_centering, LatticeType *new_latt, - char *new_ua) + char *new_ua, + RationalMatrix **pCi) { LatticeType lt; char ua, cen; - IntegerMatrix *t = NULL; + IntegerMatrix *C = NULL; + RationalMatrix *Ci = NULL; lt = cell_get_lattice_type(in); ua = cell_get_unique_axis(in); @@ -382,13 +409,17 @@ static IntegerMatrix *centering_transformation(UnitCell *in, *new_centering = cen; *new_latt = lt; *new_ua = ua; - t = intmat_identity(3); + C = intmat_identity(3); + Ci = rtnl_mtx_identity(3); } if ( cen == 'I' ) { - t = intmat_create_3x3(0, 1, 1, + C = intmat_create_3x3(0, 1, 1, 1, 0, 1, 1, 1, 0); + Ci = create_rtnl_mtx(-1,2, 1,2, 1,2, + 1,2, -1,2, 1,2, + 1,2, 1,2, -1,2); if ( lt == L_CUBIC ) { *new_latt = L_RHOMBOHEDRAL; *new_centering = 'R'; @@ -401,9 +432,12 @@ static IntegerMatrix *centering_transformation(UnitCell *in, } if ( cen == 'F' ) { - t = intmat_create_3x3(-1, 1, 1, + C = intmat_create_3x3(-1, 1, 1, 1, -1, 1, 1, 1, -1); + Ci = create_rtnl_mtx( 0,1, 1,2, 1,2, + 1,2, 0,1, 1,2, + 1,2, 1,2, 0,1); if ( lt == L_CUBIC ) { *new_latt = L_RHOMBOHEDRAL; *new_centering = 'R'; @@ -417,9 +451,12 @@ static IntegerMatrix *centering_transformation(UnitCell *in, if ( (lt == L_HEXAGONAL) && (cen == 'H') && (ua == 'c') ) { /* Obverse setting */ - t = intmat_create_3x3(1, -1, 0, + C = intmat_create_3x3(1, -1, 0, 0, 1, -1, 1, 1, 1); + Ci = create_rtnl_mtx( 2,3, 1,3, 1,3, + -1,3, 1,3, 1,3, + -1,3, -2,3, 1,3); assert(lt == L_HEXAGONAL); assert(ua == 'c'); *new_latt = L_RHOMBOHEDRAL; @@ -428,9 +465,12 @@ static IntegerMatrix *centering_transformation(UnitCell *in, } if ( cen == 'A' ) { - t = intmat_create_3x3(1, 0, 0, - 0, 1, -1, - 0, 1, 1); + C = intmat_create_3x3(0, 0, -1, + 1, 1, 0, + 1, -1, 0); + Ci = create_rtnl_mtx( 0,1, 1,2, 1,2, + 0,1, 1,2, -1,2, + -1,1, 0,1, 0,1); if ( lt == L_ORTHORHOMBIC ) { *new_latt = L_MONOCLINIC; *new_centering = 'P'; @@ -443,9 +483,12 @@ static IntegerMatrix *centering_transformation(UnitCell *in, } if ( cen == 'B' ) { - t = intmat_create_3x3(1, 0, -1, - 0, 1, 0, - 1, 0, 1); + C = intmat_create_3x3(1, 1, 0, + 0, 0, 1, + 1, -1, 0); + Ci = create_rtnl_mtx( 1,2, 0,1, 1,2, + 1,2, 0,1, -1,2, + 0,1, 1,1, 0,1); if ( lt == L_ORTHORHOMBIC ) { *new_latt = L_MONOCLINIC; *new_centering = 'P'; @@ -458,9 +501,12 @@ static IntegerMatrix *centering_transformation(UnitCell *in, } if ( cen == 'C' ) { - t = intmat_create_3x3(1, -1, 0, + C = intmat_create_3x3(1, -1, 0, 1, 1, 0, 0, 0, 1); + Ci = create_rtnl_mtx( 1,2, 1,2, 0,1, + -1,2, 1,2, 0,1, + 0,1, 0,1, 1,1); if ( lt == L_ORTHORHOMBIC ) { *new_latt = L_MONOCLINIC; *new_centering = 'P'; @@ -472,45 +518,57 @@ static IntegerMatrix *centering_transformation(UnitCell *in, } } - return t; + *pCi = Ci; + return C; } /** * uncenter_cell: * @in: A %UnitCell - * @t: Location at which to store the transformation which was used. + * @C: Location at which to store the centering transformation + * @Ci: Location at which to store the inverse centering transformation + * + * Turns any cell into a primitive one, e.g. for comparison purposes. * - * Turns any cell into a primitive one, e.g. for comparison purposes. The - * transformation which was used is stored at @t, which can be NULL if the - * transformation is not required. + * The transformation which was used is stored at @Ci. The centering + * transformation, which is the transformation you should apply if you want to + * get back the original cell, will be stored at @C. Either or both of these + * can be NULL if you don't need that information. * - * Returns: a primitive version of @in in a conventional (unique axis c) - * setting. + * Returns: a primitive version of @in. * */ -UnitCell *uncenter_cell(UnitCell *in, IntegerMatrix **t) +UnitCell *uncenter_cell(UnitCell *in, IntegerMatrix **pC, RationalMatrix **pCi) { - IntegerMatrix *tt; + IntegerMatrix *C; + RationalMatrix *Ci; char new_centering; LatticeType new_latt; char new_ua; UnitCell *out; - tt = centering_transformation(in, &new_centering, &new_latt, &new_ua); - if ( tt == NULL ) return NULL; + C = centering_transformation(in, &new_centering, &new_latt, + &new_ua, &Ci); + if ( C == NULL ) return NULL; - out = cell_transform_intmat_inverse(in, tt); + out = cell_transform_rational(in, Ci); if ( out == NULL ) return NULL; cell_set_lattice_type(out, new_latt); cell_set_centering(out, new_centering); cell_set_unique_axis(out, new_ua); - if ( t != NULL ) { - *t = tt; + if ( pC != NULL ) { + *pC = C; } else { - intmat_free(tt); + intmat_free(C); + } + + if ( pCi != NULL ) { + *pCi = Ci; + } else { + rtnl_mtx_free(Ci); } return out; @@ -578,12 +636,12 @@ UnitCell *match_cell(UnitCell *cell_in, UnitCell *template_in, int verbose, UnitCell *new_cell_trans; /* "Un-center" the template unit cell to make the comparison easier */ - template = uncenter_cell(template_in, ¢ering); + 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); + cell = uncenter_cell(cell_in, NULL, NULL); if ( cell == NULL ) return NULL; if ( cell_get_reciprocal(template, &asx, &asy, &asz, @@ -820,11 +878,11 @@ UnitCell *match_cell_ab(UnitCell *cell_in, UnitCell *template_in) UnitCell *new_cell_trans; /* "Un-center" the template unit cell to make the comparison easier */ - template = uncenter_cell(template_in, &to_given_cell); + template = uncenter_cell(template_in, &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); + cell = uncenter_cell(cell_in, NULL, NULL); /* Get the lengths to match */ if ( cell_get_cartesian(template, &ax, &ay, &az, @@ -1904,8 +1962,9 @@ int compare_reindexed_cell_parameters(UnitCell *cell_in, UnitCell *reference_in, { UnitCell *cell; UnitCell *reference; - IntegerMatrix *centering_reference; - IntegerMatrix *centering_cell; + IntegerMatrix *CBint; + RationalMatrix *CiA; + RationalMatrix *CB; RationalMatrix *m; double a, b, c, al, be, ga; double av[3], bv[3], cv[3]; @@ -1915,13 +1974,18 @@ int compare_reindexed_cell_parameters(UnitCell *cell_in, UnitCell *reference_in, int ncand_a, ncand_b, ncand_c; int i; int ia, ib; + RationalMatrix *MCiA; + RationalMatrix *CBMCiA; + int ret = 0; /* Actually compare against primitive version of reference */ - reference = uncenter_cell(reference_in, ¢ering_reference); + reference = uncenter_cell(reference_in, &CBint, NULL); if ( reference == NULL ) return 0; + CB = rtnl_mtx_from_intmat(CBint); + intmat_free(CBint); /* Actually compare primitive version of cell */ - cell = uncenter_cell(cell_in, ¢ering_cell); + cell = uncenter_cell(cell_in, NULL, &CiA); if ( cell == NULL ) return 0; /* Get target parameters */ @@ -1955,6 +2019,8 @@ int compare_reindexed_cell_parameters(UnitCell *cell_in, UnitCell *reference_in, } m = rtnl_mtx_new(3, 3); + MCiA = rtnl_mtx_new(3, 3); + CBMCiA = rtnl_mtx_new(3, 3); for ( ia=0; iacols == A->cols); + assert(ans->rows == B->rows); + assert(A->cols == B->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_mtx_get(B, k, j)); + sum = rtnl_add(sum, add); + } + rtnl_mtx_set(ans, i, j, sum); + } + } +} + + void rtnl_mtx_mult(const RationalMatrix *m, const Rational *vec, Rational *ans) { int i, j; diff --git a/libcrystfel/src/rational.h b/libcrystfel/src/rational.h index 8590e920..9ed20bfd 100644 --- a/libcrystfel/src/rational.h +++ b/libcrystfel/src/rational.h @@ -89,6 +89,8 @@ extern RationalMatrix *rtnl_mtx_from_intmat(const IntegerMatrix *m); extern RationalMatrix *rtnl_mtx_identity(int rows); 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_mult(const RationalMatrix *m, const Rational *vec, Rational *ans); extern int rtnl_mtx_solve(const RationalMatrix *m, const Rational *vec, diff --git a/libcrystfel/src/xgandalf.c b/libcrystfel/src/xgandalf.c index b122cb4c..51819956 100644 --- a/libcrystfel/src/xgandalf.c +++ b/libcrystfel/src/xgandalf.c @@ -157,7 +157,7 @@ void *xgandalf_prepare(IndexingMethod *indm, UnitCell *cell, xgandalf_private_data->cellTemplate = cell; - UnitCell* primitiveCell = uncenter_cell(cell, &xgandalf_private_data->centeringTransformation); + UnitCell* primitiveCell = uncenter_cell(cell, &xgandalf_private_data->centeringTransformation, NULL); UnitCell *uc = cell_new_from_cell(primitiveCell); reduceCell(primitiveCell, &xgandalf_private_data->latticeReductionTransform); -- cgit v1.2.3 From f231e9e1a8293f103acffced3c295116e9f17a52 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Wed, 6 Mar 2019 12:00:27 +0100 Subject: Tidy up cell comparison --- libcrystfel/src/cell-utils.c | 38 ++++++++++++++------------------------ 1 file changed, 14 insertions(+), 24 deletions(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c index 3a5ef782..7e622023 100644 --- a/libcrystfel/src/cell-utils.c +++ b/libcrystfel/src/cell-utils.c @@ -1972,11 +1972,9 @@ int compare_reindexed_cell_parameters(UnitCell *cell_in, UnitCell *reference_in, Rational *cand_b; Rational *cand_c; int ncand_a, ncand_b, ncand_c; - int i; int ia, ib; RationalMatrix *MCiA; RationalMatrix *CBMCiA; - int ret = 0; /* Actually compare against primitive version of reference */ reference = uncenter_cell(reference_in, &CBint, NULL); @@ -1999,25 +1997,6 @@ int compare_reindexed_cell_parameters(UnitCell *cell_in, UnitCell *reference_in, cand_b = find_candidates(b, av, bv, cv, ltl, &ncand_b); cand_c = find_candidates(c, av, bv, cv, ltl, &ncand_c); - STATUS("Candidates for a: %i\n", ncand_a); - for ( i=0; i<10; i++ ) { - STATUS("%s %s %s\n", rtnl_format(cand_a[3*i+0]), - rtnl_format(cand_a[3*i+1]), - rtnl_format(cand_a[3*i+2])); - } - STATUS("Candidates for b: %i\n", ncand_b); - for ( i=0; i<10; i++ ) { - STATUS("%s %s %s\n", rtnl_format(cand_b[3*i+0]), - rtnl_format(cand_b[3*i+1]), - rtnl_format(cand_b[3*i+2])); - } - STATUS("Candidates for c: %i\n", ncand_c); - for ( i=0; i<10; i++ ) { - STATUS("%s %s %s\n", rtnl_format(cand_c[3*i+0]), - rtnl_format(cand_c[3*i+1]), - rtnl_format(cand_c[3*i+2])); - } - m = rtnl_mtx_new(3, 3); MCiA = rtnl_mtx_new(3, 3); CBMCiA = rtnl_mtx_new(3, 3); @@ -2080,13 +2059,24 @@ int compare_reindexed_cell_parameters(UnitCell *cell_in, UnitCell *reference_in, rtnl_mtx_mtxmult(m, CiA, MCiA); rtnl_mtx_mtxmult(CB, MCiA, CBMCiA); *pmb = CBMCiA; - ret = 1; - cell_free(test); + rtnl_mtx_free(m); + rtnl_mtx_free(MCiA); + /* Not CBMCiA because we are returning it */ + free(cand_a); + free(cand_b); + free(cand_c); + return 1; } } } - return ret; + rtnl_mtx_free(m); + rtnl_mtx_free(MCiA); + rtnl_mtx_free(CBMCiA); + free(cand_a); + free(cand_b); + free(cand_c); + return 0; } -- cgit v1.2.3 From 66e7b71f0acc1f13c9a97c1c63aea9cdf26700cd Mon Sep 17 00:00:00 2001 From: Thomas White Date: Wed, 6 Mar 2019 14:03:23 +0100 Subject: Find the best cell match, not just the first one --- libcrystfel/src/cell-utils.c | 67 ++++++++++++++++++++++++++++++++++---------- 1 file changed, 52 insertions(+), 15 deletions(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c index 7e622023..b9ccb4f6 100644 --- a/libcrystfel/src/cell-utils.c +++ b/libcrystfel/src/cell-utils.c @@ -1937,6 +1937,35 @@ 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) +{ + 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(double a1, double b1, double c1, + double al1, double be1, double ga1, + double a2, double b2, double c2, + double al2, double be2, double ga2) +{ + 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); + for ( i=0; i<6; i++ ) { + total += (g1[i]-g2[i])*(g1[i]-g2[i]); + } + return sqrt(total); +} + + /** * compare_reindexed_cell_parameters: * @cell_in: A %UnitCell @@ -1975,6 +2004,7 @@ int compare_reindexed_cell_parameters(UnitCell *cell_in, UnitCell *reference_in, int ia, ib; RationalMatrix *MCiA; RationalMatrix *CBMCiA; + double min_dist = +INFINITY; /* Actually compare against primitive version of reference */ reference = uncenter_cell(reference_in, &CBint, NULL); @@ -2005,6 +2035,7 @@ int compare_reindexed_cell_parameters(UnitCell *cell_in, UnitCell *reference_in, UnitCell *test; double at, bt, ct, alt, bet, gat; + double dist; int ic = 0; /* Form the matrix using the first candidate for c */ @@ -2054,29 +2085,35 @@ int compare_reindexed_cell_parameters(UnitCell *cell_in, UnitCell *reference_in, continue; } - /* m is a rational matrix which turns Ap into Bp - * We need to return CB.M.CiA, which turns A into B */ - rtnl_mtx_mtxmult(m, CiA, MCiA); - rtnl_mtx_mtxmult(CB, MCiA, CBMCiA); - *pmb = CBMCiA; + dist = g6_distance(at, bt, ct, alt, bet, gat, + a, b, c, al, be, ga); + if ( dist < min_dist ) { + min_dist = dist; + rtnl_mtx_mtxmult(m, CiA, MCiA); + } + cell_free(test); - rtnl_mtx_free(m); - rtnl_mtx_free(MCiA); - /* Not CBMCiA because we are returning it */ - free(cand_a); - free(cand_b); - free(cand_c); - return 1; } } } rtnl_mtx_free(m); - rtnl_mtx_free(MCiA); - rtnl_mtx_free(CBMCiA); free(cand_a); free(cand_b); free(cand_c); - return 0; + + if ( isinf(min_dist) ) { + rtnl_mtx_free(CBMCiA); + rtnl_mtx_free(MCiA); + *pmb = NULL; + return 0; + } + + /* Solution found */ + STATUS("minimum norm = %e\n", min_dist); + rtnl_mtx_mtxmult(CB, MCiA, CBMCiA); + rtnl_mtx_free(MCiA); + *pmb = CBMCiA; + return 1; } -- cgit v1.2.3 From 0f4dbbc0a8f9b05bcaff1b771a1ee4f03d070681 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Wed, 6 Mar 2019 14:42:13 +0100 Subject: More tidying up of cell comparison --- libcrystfel/src/cell-utils.c | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c index b9ccb4f6..819f4926 100644 --- a/libcrystfel/src/cell-utils.c +++ b/libcrystfel/src/cell-utils.c @@ -1884,7 +1884,7 @@ static Rational *find_candidates(double len, double *a, double *b, double *c, cands = malloc(max_cand * sizeof(struct cand)); if ( cands == NULL ) return NULL; - rat = rtnl_list(0, 4, -4, 4, &nrat); + rat = rtnl_list(-5, 5, 1, 4, &nrat); if ( rat == NULL ) return NULL; for ( ia=0; ia Date: Fri, 22 Feb 2019 11:42:29 +0100 Subject: Initial centering determination --- libcrystfel/src/cell.c | 37 +++++++++++++++++++++++++++++++++++-- 1 file changed, 35 insertions(+), 2 deletions(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/cell.c b/libcrystfel/src/cell.c index 5eb13349..a26e26a6 100644 --- a/libcrystfel/src/cell.c +++ b/libcrystfel/src/cell.c @@ -46,6 +46,7 @@ #include "utils.h" #include "image.h" #include "integer_matrix.h" +#include "rational.h" /** @@ -687,6 +688,31 @@ UnitCell *cell_transform_gsl_reciprocal(UnitCell *in, gsl_matrix *m) } +static char determine_centering(IntegerMatrix *m, char cen) +{ + Rational c[3]; + Rational nc[3]; + + c[0] = rtnl(1, 2); + c[1] = rtnl(1, 2); + c[2] = rtnl_zero(); + intmat_rationalvec_mult(m, c, nc); + STATUS("%s,%s,%s -> %s,%s,%s\n", + rtnl_format(c[0]), rtnl_format(c[1]), rtnl_format(c[2]), + rtnl_format(nc[0]), rtnl_format(nc[1]), rtnl_format(nc[2])); + + c[0] = rtnl_zero(); + c[1] = rtnl_zero(); + c[2] = rtnl_zero(); + intmat_solve_rational(m, nc, c); + STATUS("%s,%s,%s <- %s,%s,%s\n", + rtnl_format(c[0]), rtnl_format(c[1]), rtnl_format(c[2]), + rtnl_format(nc[0]), rtnl_format(nc[1]), rtnl_format(nc[2])); + + return 'A'; +} + + /** * cell_transform_rational: * @cell: A %UnitCell. @@ -701,6 +727,7 @@ UnitCell *cell_transform_rational(UnitCell *cell, RationalMatrix *m) { UnitCell *out; gsl_matrix *tm; + char ncen; if ( m == NULL ) return NULL; @@ -720,10 +747,16 @@ UnitCell *cell_transform_rational(UnitCell *cell, RationalMatrix *m) gsl_matrix_set(tm, 2, 2, rtnl_as_double(rtnl_mtx_get(m, 2, 2))); out = cell_transform_gsl_direct(cell, tm); + gsl_matrix_free(tm); - /* FIXME: Update centering, unique axis, lattice type */ + ncen = determine_centering(m, cell_get_centering(cell)); + if ( ncen == '*' ) { + cell_free(out); + return NULL; + } + cell_set_centering(out, ncen); - gsl_matrix_free(tm); + /* FIXME: Update unique axis, lattice type */ return out; } -- cgit v1.2.3 From d0e9fd4cb3407916e6ceb455a20e01c09847316a Mon Sep 17 00:00:00 2001 From: Thomas White Date: Thu, 24 Jan 2019 11:58:22 +0100 Subject: Automatic centering determination after transformation --- libcrystfel/src/cell.c | 191 ++++++++++++++++++++++++++++++++++++++++++++----- 1 file changed, 172 insertions(+), 19 deletions(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/cell.c b/libcrystfel/src/cell.c index a26e26a6..131e1ac8 100644 --- a/libcrystfel/src/cell.c +++ b/libcrystfel/src/cell.c @@ -97,6 +97,20 @@ struct _unitcell { char unique_axis; }; +typedef enum { + CMASK_P = 1<<0, + CMASK_A = 1<<1, + CMASK_B = 1<<2, + CMASK_C = 1<<3, + CMASK_I = 1<<4, + CMASK_F = 1<<5, + CMASK_H = 1<<6, + CMASK_R = 1<<7 +} CenteringMask; + +#define CMASK_ALL (CMASK_P | CMASK_A | CMASK_B | CMASK_C | CMASK_I \ + | CMASK_F | CMASK_H | CMASK_R) + /************************** Setters and Constructors **************************/ @@ -688,28 +702,167 @@ UnitCell *cell_transform_gsl_reciprocal(UnitCell *in, gsl_matrix *m) } -static char determine_centering(IntegerMatrix *m, char cen) +static int centering_has_point(char cen, Rational *p) +{ + /* First, put the point into the range 0..1 */ + while ( rtnl_cmp(p[0], rtnl_zero()) < 0 ) p[0] = rtnl_add(p[0], rtnl(1, 1)); + while ( rtnl_cmp(p[1], rtnl_zero()) < 0 ) p[1] = rtnl_add(p[1], rtnl(1, 1)); + while ( rtnl_cmp(p[2], rtnl_zero()) < 0 ) p[2] = rtnl_add(p[2], rtnl(1, 1)); + while ( rtnl_cmp(p[0], rtnl(1, 1)) >= 0 ) p[0] = rtnl_sub(p[0], rtnl(1, 1)); + while ( rtnl_cmp(p[1], rtnl(1, 1)) >= 0 ) p[1] = rtnl_sub(p[1], rtnl(1, 1)); + while ( rtnl_cmp(p[2], rtnl(1, 1)) >= 0 ) p[2] = rtnl_sub(p[2], rtnl(1, 1)); + + /* 0,0,0 is present in all centerings */ + if ( (rtnl_cmp(p[0], rtnl_zero()) == 0) + && (rtnl_cmp(p[1], rtnl_zero()) == 0) + && (rtnl_cmp(p[2], rtnl_zero()) == 0) ) return 1; + + /* Only I has 1/2 , 1/2, 1/2 */ + if ( (rtnl_cmp(p[0], rtnl(1,2)) == 0) + && (rtnl_cmp(p[1], rtnl(1,2)) == 0) + && (rtnl_cmp(p[2], rtnl(1,2)) == 0) + && (cen == 'I') ) return 1; + + /* A or F has 0 , 1/2, 1/2 */ + if ( (rtnl_cmp(p[0], rtnl_zero()) == 0) + && (rtnl_cmp(p[1], rtnl(1,2)) == 0) + && (rtnl_cmp(p[2], rtnl(1,2)) == 0) + && ((cen == 'A') || (cen == 'F')) ) return 1; + + /* B or F has 1/2 , 0 , 1/2 */ + if ( (rtnl_cmp(p[0], rtnl(1,2)) == 0) + && (rtnl_cmp(p[1], rtnl_zero()) == 0) + && (rtnl_cmp(p[2], rtnl(1,2)) == 0) + && ((cen == 'B') || (cen == 'F')) ) return 1; + + /* C or F has 1/2 , 1/2 , 0 */ + if ( (rtnl_cmp(p[0], rtnl(1,2)) == 0) + && (rtnl_cmp(p[1], rtnl(1,2)) == 0) + && (rtnl_cmp(p[2], rtnl_zero()) == 0) + && ((cen == 'C') || (cen == 'F')) ) return 1; + + /* H has 2/3 , 1/3 , 1/3 */ + if ( (rtnl_cmp(p[0], rtnl(2,3)) == 0) + && (rtnl_cmp(p[1], rtnl(1,3)) == 0) + && (rtnl_cmp(p[2], rtnl(1,3)) == 0) + && (cen == 'H') ) return 1; + + /* H has 1/3 , 2/3 , 2/3 */ + if ( (rtnl_cmp(p[0], rtnl(1,3)) == 0) + && (rtnl_cmp(p[1], rtnl(2,3)) == 0) + && (rtnl_cmp(p[2], rtnl(2,3)) == 0) + && (cen == 'H') ) return 1; + + return 0; +} + + +static void maybe_eliminate(CenteringMask c, CenteringMask *cmask, Rational *nc, + char cen) +{ + /* Skip test if this centering isn't even a candidate */ + if ( !(*cmask & c) ) return; + + if ( !centering_has_point(cen, nc) ) { + *cmask |= c; + *cmask ^= c; + } +} + + +/* Check if the point x,y,z in the original cell matches any lattice point + * in the transformed cell */ +static void check_point_fwd(RationalMatrix *m, CenteringMask *cmask, + Rational x, Rational y, Rational z) +{ + Rational c[3] = {x, y, z}; + Rational nc[3]; + + /* Transform the lattice point */ + rtnl_mtx_solve(m, c, nc); + + /* Eliminate any centerings which don't include the transformed point */ + maybe_eliminate(CMASK_P, cmask, nc, 'P'); + maybe_eliminate(CMASK_R, cmask, nc, 'R'); + maybe_eliminate(CMASK_A, cmask, nc, 'A'); + maybe_eliminate(CMASK_B, cmask, nc, 'B'); + maybe_eliminate(CMASK_C, cmask, nc, 'C'); + maybe_eliminate(CMASK_I, cmask, nc, 'I'); + maybe_eliminate(CMASK_F, cmask, nc, 'F'); + maybe_eliminate(CMASK_H, cmask, nc, 'H'); +} + + +/* Check if the point x,y,z in the transformed cell matches any lattice point + * in the original cell. If not, eliminate "exclude" from "*mask". */ +static void check_point_bwd(RationalMatrix *m, CenteringMask *mask, + char cen, CenteringMask exclude, + Rational x, Rational y, Rational z) { - Rational c[3]; Rational nc[3]; + Rational c[3] = {x, y, z}; + + rtnl_mtx_mult(m, c, nc); + + if ( !centering_has_point(cen, nc) ) { + *mask |= exclude; + *mask ^= exclude; /* Unset bits */ + } +} + + +static char cmask_decode(CenteringMask mask) +{ + char res[32]; + + res[0] = '\0'; + + if ( mask & CMASK_H ) strcat(res, "H"); + if ( mask & CMASK_F ) strcat(res, "F"); + if ( mask & CMASK_I ) strcat(res, "I"); + if ( mask & CMASK_A ) strcat(res, "A"); + if ( mask & CMASK_B ) strcat(res, "B"); + if ( mask & CMASK_C ) strcat(res, "C"); + if ( mask & CMASK_P ) strcat(res, "P"); + if ( mask & CMASK_R ) strcat(res, "R"); + + STATUS("possible centerings: %s\n", res); + return res[0]; +} + + +static char determine_centering(RationalMatrix *m, char cen) +{ + CenteringMask cmask = CMASK_ALL; + + /* Check whether the current centering can provide all the lattice + * points for the transformed cell. Eliminate any centerings for which + * it can't. */ + check_point_bwd(m, &cmask, cen, CMASK_A | CMASK_F, rtnl_zero(), rtnl(1,2), rtnl(1,2)); + check_point_bwd(m, &cmask, cen, CMASK_B | CMASK_F, rtnl(1,2), rtnl_zero(), rtnl(1,2)); + check_point_bwd(m, &cmask, cen, CMASK_C | CMASK_F, rtnl(1,2), rtnl(1,2), rtnl_zero()); + check_point_bwd(m, &cmask, cen, CMASK_I, rtnl(1,2), rtnl(1,2), rtnl(1,2)); + check_point_bwd(m, &cmask, cen, CMASK_H, rtnl(2,3), rtnl(1,3), rtnl(1,3)); + check_point_bwd(m, &cmask, cen, CMASK_H, rtnl(1,3), rtnl(2,3), rtnl(2,3)); + + /* Check whether the current centering's lattice points will all + * coincide with lattice points in the new centering. Eliminate any + * centerings for which they don't (they give "excess lattice points"). */ + switch ( cen ) { + + case 'P' : + case 'R' : + break; + + case 'C' : + check_point_fwd(m, &cmask, rtnl(1,2), rtnl(1,2), rtnl_zero()); + break; + + /* FIXME: Everything else */ + + } - c[0] = rtnl(1, 2); - c[1] = rtnl(1, 2); - c[2] = rtnl_zero(); - intmat_rationalvec_mult(m, c, nc); - STATUS("%s,%s,%s -> %s,%s,%s\n", - rtnl_format(c[0]), rtnl_format(c[1]), rtnl_format(c[2]), - rtnl_format(nc[0]), rtnl_format(nc[1]), rtnl_format(nc[2])); - - c[0] = rtnl_zero(); - c[1] = rtnl_zero(); - c[2] = rtnl_zero(); - intmat_solve_rational(m, nc, c); - STATUS("%s,%s,%s <- %s,%s,%s\n", - rtnl_format(c[0]), rtnl_format(c[1]), rtnl_format(c[2]), - rtnl_format(nc[0]), rtnl_format(nc[1]), rtnl_format(nc[2])); - - return 'A'; + return cmask_decode(cmask); } -- cgit v1.2.3 From 8de701d8137fe64ee2d19e89061befd17b31484d Mon Sep 17 00:00:00 2001 From: Thomas White Date: Wed, 6 Mar 2019 17:04:47 +0100 Subject: Add remaining centering types and finalise centering determination --- libcrystfel/src/cell.c | 43 +++++++++++++++++++++++++++++++++++++------ 1 file changed, 37 insertions(+), 6 deletions(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/cell.c b/libcrystfel/src/cell.c index 131e1ac8..32e45856 100644 --- a/libcrystfel/src/cell.c +++ b/libcrystfel/src/cell.c @@ -826,7 +826,7 @@ static char cmask_decode(CenteringMask mask) if ( mask & CMASK_P ) strcat(res, "P"); if ( mask & CMASK_R ) strcat(res, "R"); - STATUS("possible centerings: %s\n", res); + if ( strlen(res) == 0 ) return '?'; return res[0]; } @@ -854,11 +854,32 @@ static char determine_centering(RationalMatrix *m, char cen) case 'R' : break; + case 'A' : + check_point_fwd(m, &cmask, rtnl_zero(), rtnl(1,2), rtnl(1,2)); + break; + + case 'B' : + check_point_fwd(m, &cmask, rtnl(1,2), rtnl_zero(), rtnl(1,2)); + break; + case 'C' : check_point_fwd(m, &cmask, rtnl(1,2), rtnl(1,2), rtnl_zero()); break; - /* FIXME: Everything else */ + case 'I' : + check_point_fwd(m, &cmask, rtnl(1,2), rtnl(1,2), rtnl(1,2)); + break; + + case 'F' : + check_point_fwd(m, &cmask, rtnl_zero(), rtnl(1,2), rtnl(1,2)); + check_point_fwd(m, &cmask, rtnl(1,2), rtnl_zero(), rtnl(1,2)); + check_point_fwd(m, &cmask, rtnl(1,2), rtnl(1,2), rtnl_zero()); + break; + + case 'H' : + check_point_fwd(m, &cmask, rtnl(2,3), rtnl(1,3), rtnl(1,3)); + check_point_fwd(m, &cmask, rtnl(1,3), rtnl(2,3), rtnl(2,3)); + break; } @@ -873,6 +894,14 @@ static char determine_centering(RationalMatrix *m, char cen) * * Applies @t to @cell. * + * This function will determine the centering of the resulting unit cell, + * producing '?' if any lattice points cannot be accounted for. Note that if + * there are 'excess' lattice points in the transformed cell, the centering + * will still be '?' even if the lattice points for another centering are + * all present. + * + * The lattice type will be set to triclinic, and the unique axis to '?'. + * * Returns: Transformed copy of @cell. * */ @@ -903,13 +932,11 @@ UnitCell *cell_transform_rational(UnitCell *cell, RationalMatrix *m) gsl_matrix_free(tm); ncen = determine_centering(m, cell_get_centering(cell)); - if ( ncen == '*' ) { - cell_free(out); - return NULL; - } cell_set_centering(out, ncen); /* FIXME: Update unique axis, lattice type */ + cell_set_lattice_type(out, L_TRICLINIC); + cell_set_unique_axis(out, '?'); return out; } @@ -922,6 +949,10 @@ UnitCell *cell_transform_rational(UnitCell *cell, RationalMatrix *m) * * Applies @t to @cell. * + * This is just a convenience function which turns @m into a %RationalMatrix + * and then calls cell_transform_rational(). See the documentation for that + * function for some important information. + * * Returns: Transformed copy of @cell. * */ -- cgit v1.2.3 From e94501d025e8d6b4f4e446004e80a1badfbba30f Mon Sep 17 00:00:00 2001 From: Thomas White Date: Sat, 9 Mar 2019 11:35:48 +0100 Subject: Fix typo --- libcrystfel/src/cell-utils.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c index 819f4926..cc17f108 100644 --- a/libcrystfel/src/cell-utils.c +++ b/libcrystfel/src/cell-utils.c @@ -381,7 +381,7 @@ static RationalMatrix *create_rtnl_mtx(signed int a1, signed int a2, rtnl_mtx_set(m, 1, 2, rtnl(f1, f2)); rtnl_mtx_set(m, 2, 0, rtnl(g1, g2)); rtnl_mtx_set(m, 2, 1, rtnl(h1, h2)); - rtnl_mtx_set(m, 2, 2, rtnl(i1, h2)); + rtnl_mtx_set(m, 2, 2, rtnl(i1, i2)); return m; } -- cgit v1.2.3 From 252f4073f0041027b5e9e6773fe892ffe465ed23 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Sat, 9 Mar 2019 11:37:34 +0100 Subject: Single point of truth for matrix memory layout --- libcrystfel/src/cell.c | 32 ++++++++++++++------------------ libcrystfel/src/rational.c | 4 ++-- 2 files changed, 16 insertions(+), 20 deletions(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/cell.c b/libcrystfel/src/cell.c index 32e45856..130a24b1 100644 --- a/libcrystfel/src/cell.c +++ b/libcrystfel/src/cell.c @@ -910,6 +910,7 @@ UnitCell *cell_transform_rational(UnitCell *cell, RationalMatrix *m) UnitCell *out; gsl_matrix *tm; char ncen; + int i, j; if ( m == NULL ) return NULL; @@ -918,15 +919,12 @@ UnitCell *cell_transform_rational(UnitCell *cell, RationalMatrix *m) return NULL; } - gsl_matrix_set(tm, 0, 0, rtnl_as_double(rtnl_mtx_get(m, 0, 0))); - gsl_matrix_set(tm, 0, 1, rtnl_as_double(rtnl_mtx_get(m, 0, 1))); - gsl_matrix_set(tm, 0, 2, rtnl_as_double(rtnl_mtx_get(m, 0, 2))); - gsl_matrix_set(tm, 1, 0, rtnl_as_double(rtnl_mtx_get(m, 1, 0))); - gsl_matrix_set(tm, 1, 1, rtnl_as_double(rtnl_mtx_get(m, 1, 1))); - gsl_matrix_set(tm, 1, 2, rtnl_as_double(rtnl_mtx_get(m, 1, 2))); - gsl_matrix_set(tm, 2, 0, rtnl_as_double(rtnl_mtx_get(m, 2, 0))); - gsl_matrix_set(tm, 2, 1, rtnl_as_double(rtnl_mtx_get(m, 2, 1))); - gsl_matrix_set(tm, 2, 2, rtnl_as_double(rtnl_mtx_get(m, 2, 2))); + for ( i=0; i<3; i++ ) { + for ( j=0; j<3; j++ ) { + gsl_matrix_set(tm, i, j, + rtnl_as_double(rtnl_mtx_get(m, i, j))); + } + } out = cell_transform_gsl_direct(cell, tm); gsl_matrix_free(tm); @@ -983,6 +981,7 @@ UnitCell *cell_transform_rational_inverse(UnitCell *cell, RationalMatrix *m) gsl_matrix *inv; gsl_permutation *perm; int s; + int i, j; if ( m == NULL ) return NULL; @@ -991,15 +990,12 @@ UnitCell *cell_transform_rational_inverse(UnitCell *cell, RationalMatrix *m) return NULL; } - gsl_matrix_set(tm, 0, 0, rtnl_as_double(rtnl_mtx_get(m, 0, 0))); - gsl_matrix_set(tm, 0, 1, rtnl_as_double(rtnl_mtx_get(m, 0, 1))); - gsl_matrix_set(tm, 0, 2, rtnl_as_double(rtnl_mtx_get(m, 0, 2))); - gsl_matrix_set(tm, 1, 0, rtnl_as_double(rtnl_mtx_get(m, 1, 0))); - gsl_matrix_set(tm, 1, 1, rtnl_as_double(rtnl_mtx_get(m, 1, 1))); - gsl_matrix_set(tm, 1, 2, rtnl_as_double(rtnl_mtx_get(m, 1, 2))); - gsl_matrix_set(tm, 2, 0, rtnl_as_double(rtnl_mtx_get(m, 2, 0))); - gsl_matrix_set(tm, 2, 1, rtnl_as_double(rtnl_mtx_get(m, 2, 1))); - gsl_matrix_set(tm, 2, 2, rtnl_as_double(rtnl_mtx_get(m, 2, 2))); + for ( i=0; i<3; i++ ) { + for ( j=0; j<3; j++ ) { + gsl_matrix_set(tm, i, j, + rtnl_as_double(rtnl_mtx_get(m, i, j))); + } + } perm = gsl_permutation_alloc(3); if ( perm == NULL ) { diff --git a/libcrystfel/src/rational.c b/libcrystfel/src/rational.c index 4f02c8b1..5f17b0ea 100644 --- a/libcrystfel/src/rational.c +++ b/libcrystfel/src/rational.c @@ -367,7 +367,7 @@ RationalMatrix *rtnl_mtx_from_intmat(const IntegerMatrix *m) for ( i=0; iv[j+cols*i] = rtnl(intmat_get(m, i, j), 1); + rtnl_mtx_set(n, i, j, rtnl(intmat_get(m, i, j), 1)); } } @@ -537,7 +537,7 @@ void rtnl_mtx_print(const RationalMatrix *m) fprintf(stderr, "[ "); for ( j=0; jcols; j++ ) { - char *v = rtnl_format(m->v[j+m->cols*i]); + char *v = rtnl_format(rtnl_mtx_get(m, i, j)); fprintf(stderr, "%4s ", v); free(v); } -- cgit v1.2.3 From 555a0319ff558123aaabb8214190c53b621ffe95 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Sat, 9 Mar 2019 11:38:49 +0100 Subject: Add a comment --- libcrystfel/src/rational.c | 1 + 1 file changed, 1 insertion(+) (limited to 'libcrystfel') diff --git a/libcrystfel/src/rational.c b/libcrystfel/src/rational.c index 5f17b0ea..20a73bc8 100644 --- a/libcrystfel/src/rational.c +++ b/libcrystfel/src/rational.c @@ -571,6 +571,7 @@ void rtnl_mtx_mtxmult(const RationalMatrix *A, const RationalMatrix *B, } +/* Calculate ans = m.vec, where 'ans' and 'vec' are column vectors */ void rtnl_mtx_mult(const RationalMatrix *m, const Rational *vec, Rational *ans) { int i, j; -- cgit v1.2.3 From cb51ba0e53425ae663db313d730c1df644e4f76b Mon Sep 17 00:00:00 2001 From: Thomas White Date: Sat, 9 Mar 2019 11:42:55 +0100 Subject: Change matrix notation to match ITA chapter 5.1 --- libcrystfel/src/cell-utils.c | 72 +++++++++++++++++++++------------------- libcrystfel/src/cell.c | 60 ++++++++++++++++----------------- libcrystfel/src/integer_matrix.c | 15 +++++---- libcrystfel/src/rational.c | 4 +-- libcrystfel/src/symmetry.c | 6 ++-- libcrystfel/src/symop.y | 12 +++---- 6 files changed, 88 insertions(+), 81 deletions(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c index cc17f108..605aa449 100644 --- a/libcrystfel/src/cell-utils.c +++ b/libcrystfel/src/cell-utils.c @@ -405,6 +405,10 @@ static IntegerMatrix *centering_transformation(UnitCell *in, ua = cell_get_unique_axis(in); cen = cell_get_centering(in); + /* Write the matrices exactly as they appear in ITA Table 5.1.3.1. + * C is "P", and Ci is "Q=P^-1". Vice-versa if the transformation + * should go the opposite way to what's written in the first column. */ + if ( (cen=='P') || (cen=='R') ) { *new_centering = cen; *new_latt = lt; @@ -451,12 +455,12 @@ static IntegerMatrix *centering_transformation(UnitCell *in, if ( (lt == L_HEXAGONAL) && (cen == 'H') && (ua == 'c') ) { /* Obverse setting */ - C = intmat_create_3x3(1, -1, 0, - 0, 1, -1, - 1, 1, 1); - Ci = create_rtnl_mtx( 2,3, 1,3, 1,3, - -1,3, 1,3, 1,3, - -1,3, -2,3, 1,3); + C = intmat_create_3x3( 1, 0, 1, + -1, 1, 1, + 0, -1, 1); + Ci = create_rtnl_mtx( 2,3, -1,3, -1,3, + 1,3, 1,3, -2,3, + 1,3, 1,3, 1,3); assert(lt == L_HEXAGONAL); assert(ua == 'c'); *new_latt = L_RHOMBOHEDRAL; @@ -465,12 +469,12 @@ static IntegerMatrix *centering_transformation(UnitCell *in, } if ( cen == 'A' ) { - C = intmat_create_3x3(0, 0, -1, - 1, 1, 0, - 1, -1, 0); - Ci = create_rtnl_mtx( 0,1, 1,2, 1,2, + C = intmat_create_3x3( 1, 0, 0, + 0, 1, 1, + 0, -1, 1); + Ci = create_rtnl_mtx( 1,1, 0,1, 0,1, 0,1, 1,2, -1,2, - -1,1, 0,1, 0,1); + 0,1, 1,2, 1,2); if ( lt == L_ORTHORHOMBIC ) { *new_latt = L_MONOCLINIC; *new_centering = 'P'; @@ -483,12 +487,12 @@ static IntegerMatrix *centering_transformation(UnitCell *in, } if ( cen == 'B' ) { - C = intmat_create_3x3(1, 1, 0, - 0, 0, 1, - 1, -1, 0); - Ci = create_rtnl_mtx( 1,2, 0,1, 1,2, - 1,2, 0,1, -1,2, - 0,1, 1,1, 0,1); + C = intmat_create_3x3( 1, 0, 1, + 0, 1, 0, + -1, 0, 1); + Ci = create_rtnl_mtx( 1,2, 0,1, -1,2, + 0,1, 1,1, 0,1, + 1,2, 0,1, 1,2); if ( lt == L_ORTHORHOMBIC ) { *new_latt = L_MONOCLINIC; *new_centering = 'P'; @@ -501,11 +505,11 @@ static IntegerMatrix *centering_transformation(UnitCell *in, } if ( cen == 'C' ) { - C = intmat_create_3x3(1, -1, 0, - 1, 1, 0, - 0, 0, 1); - Ci = create_rtnl_mtx( 1,2, 1,2, 0,1, - -1,2, 1,2, 0,1, + C = intmat_create_3x3( 1, 1, 0, + -1, 1, 0, + 0, 0, 1); + Ci = create_rtnl_mtx( 1,2, -1,2, 0,1, + 1,2, 1,2, 0,1, 0,1, 0,1, 1,1); if ( lt == L_ORTHORHOMBIC ) { *new_latt = L_MONOCLINIC; @@ -2047,13 +2051,13 @@ int compare_reindexed_cell_parameters(UnitCell *cell_in, UnitCell *reference_in, /* Form the matrix using the first candidate for c */ rtnl_mtx_set(m, 0, 0, cand_a[3*ia+0]); - rtnl_mtx_set(m, 0, 1, cand_a[3*ia+1]); - rtnl_mtx_set(m, 0, 2, cand_a[3*ia+2]); - rtnl_mtx_set(m, 1, 0, cand_b[3*ib+0]); + rtnl_mtx_set(m, 1, 0, cand_a[3*ia+1]); + rtnl_mtx_set(m, 2, 0, cand_a[3*ia+2]); + rtnl_mtx_set(m, 0, 1, cand_b[3*ib+0]); rtnl_mtx_set(m, 1, 1, cand_b[3*ib+1]); - rtnl_mtx_set(m, 1, 2, cand_b[3*ib+2]); - rtnl_mtx_set(m, 2, 0, cand_c[3*ic+0]); - rtnl_mtx_set(m, 2, 1, cand_c[3*ic+1]); + rtnl_mtx_set(m, 2, 1, cand_b[3*ib+2]); + rtnl_mtx_set(m, 0, 2, cand_c[3*ic+0]); + rtnl_mtx_set(m, 1, 2, cand_c[3*ic+1]); rtnl_mtx_set(m, 2, 2, cand_c[3*ic+2]); /* Check angle between a and b */ @@ -2066,13 +2070,13 @@ int compare_reindexed_cell_parameters(UnitCell *cell_in, UnitCell *reference_in, for ( ic=0; icrows * sizeof(signed int)); if ( ans == NULL ) return NULL; - for ( i=0; irows; i++ ) { - - unsigned int j; + for ( j=0; jcols; j++ ) { - ans[i] = 0; - for ( j=0; jcols; j++ ) { + unsigned int i; + ans[j] = 0; + for ( i=0; irows; i++ ) { ans[i] += intmat_get(m, i, j) * vec[j]; } diff --git a/libcrystfel/src/rational.c b/libcrystfel/src/rational.c index 20a73bc8..ab1eff91 100644 --- a/libcrystfel/src/rational.c +++ b/libcrystfel/src/rational.c @@ -551,8 +551,8 @@ void rtnl_mtx_mtxmult(const RationalMatrix *A, const RationalMatrix *B, { int i, j; - assert(ans->cols == A->cols); - assert(ans->rows == B->rows); + assert(ans->cols == B->cols); + assert(ans->rows == A->rows); assert(A->cols == B->rows); for ( i=0; irows; i++ ) { diff --git a/libcrystfel/src/symmetry.c b/libcrystfel/src/symmetry.c index 9ee80a16..05ee7c6b 100644 --- a/libcrystfel/src/symmetry.c +++ b/libcrystfel/src/symmetry.c @@ -197,9 +197,9 @@ static void add_symop_v(SymOpList *ops, m = intmat_new(3, 3); assert(m != NULL); - for ( i=0; i<3; i++ ) intmat_set(m, 0, i, h[i]); - for ( i=0; i<3; i++ ) intmat_set(m, 1, i, k[i]); - for ( i=0; i<3; i++ ) intmat_set(m, 2, i, l[i]); + for ( i=0; i<3; i++ ) intmat_set(m, i, 0, h[i]); + for ( i=0; i<3; i++ ) intmat_set(m, i, 1, k[i]); + for ( i=0; i<3; i++ ) intmat_set(m, i, 2, l[i]); free(h); free(k); diff --git a/libcrystfel/src/symop.y b/libcrystfel/src/symop.y index bebf823f..be31c21d 100644 --- a/libcrystfel/src/symop.y +++ b/libcrystfel/src/symop.y @@ -101,13 +101,13 @@ symoplist: symop: axexpr COMMA axexpr COMMA axexpr { rtnl_mtx_set(m, 0, 0, $1[0]); - rtnl_mtx_set(m, 0, 1, $1[1]); - rtnl_mtx_set(m, 0, 2, $1[2]); - rtnl_mtx_set(m, 1, 0, $3[0]); + rtnl_mtx_set(m, 1, 0, $1[1]); + rtnl_mtx_set(m, 2, 0, $1[2]); + rtnl_mtx_set(m, 0, 1, $3[0]); rtnl_mtx_set(m, 1, 1, $3[1]); - rtnl_mtx_set(m, 1, 2, $3[2]); - rtnl_mtx_set(m, 2, 0, $5[0]); - rtnl_mtx_set(m, 2, 1, $5[1]); + rtnl_mtx_set(m, 2, 1, $3[2]); + rtnl_mtx_set(m, 0, 2, $5[0]); + rtnl_mtx_set(m, 1, 2, $5[1]); rtnl_mtx_set(m, 2, 2, $5[2]); } ; -- cgit v1.2.3 From c04d93ef0e6faf50ad243932da158239f83cc619 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Mon, 11 Mar 2019 11:53:59 +0100 Subject: Get rid of cell_transform_gsl_reciprocal() It's confusing to have matrices acting on both real and reciprocal vectors. Let's standardise instead on transformations always applying to real-space cells. It was only used in one place. --- libcrystfel/src/cell.c | 43 ------------------------------------------- libcrystfel/src/cell.h | 1 - libcrystfel/src/taketwo.c | 2 +- 3 files changed, 1 insertion(+), 45 deletions(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/cell.c b/libcrystfel/src/cell.c index 5e00ac9c..9ae77a6d 100644 --- a/libcrystfel/src/cell.c +++ b/libcrystfel/src/cell.c @@ -659,49 +659,6 @@ UnitCell *cell_transform_gsl_direct(UnitCell *in, gsl_matrix *m) } -UnitCell *cell_transform_gsl_reciprocal(UnitCell *in, gsl_matrix *m) -{ - gsl_matrix *c; - double asx, asy, asz; - double bsx, bsy, bsz; - double csx, csy, csz; - gsl_matrix *res; - UnitCell *out; - - cell_get_reciprocal(in, &asx, &asy, &asz, &bsx, &bsy, - &bsz, &csx, &csy, &csz); - - c = gsl_matrix_alloc(3, 3); - gsl_matrix_set(c, 0, 0, asx); - gsl_matrix_set(c, 1, 0, asy); - gsl_matrix_set(c, 2, 0, asz); - gsl_matrix_set(c, 0, 1, bsx); - gsl_matrix_set(c, 1, 1, bsy); - gsl_matrix_set(c, 2, 1, bsz); - gsl_matrix_set(c, 0, 2, csx); - gsl_matrix_set(c, 1, 2, csy); - gsl_matrix_set(c, 2, 2, csz); - - res = gsl_matrix_calloc(3, 3); - gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, c, m, 0.0, res); - - out = cell_new_from_cell(in); - cell_set_reciprocal(out, gsl_matrix_get(res, 0, 0), - gsl_matrix_get(res, 1, 0), - gsl_matrix_get(res, 2, 0), - gsl_matrix_get(res, 0, 1), - gsl_matrix_get(res, 1, 1), - gsl_matrix_get(res, 2, 1), - gsl_matrix_get(res, 0, 2), - gsl_matrix_get(res, 1, 2), - gsl_matrix_get(res, 2, 2)); - - gsl_matrix_free(res); - gsl_matrix_free(c); - return out; -} - - static int centering_has_point(char cen, Rational *p) { /* First, put the point into the range 0..1 */ diff --git a/libcrystfel/src/cell.h b/libcrystfel/src/cell.h index 8e6809b5..c868cd29 100644 --- a/libcrystfel/src/cell.h +++ b/libcrystfel/src/cell.h @@ -149,7 +149,6 @@ extern void cell_set_unique_axis(UnitCell *cell, char unique_axis); extern const char *cell_rep(UnitCell *cell); extern UnitCell *cell_transform_gsl_direct(UnitCell *in, gsl_matrix *m); -extern UnitCell *cell_transform_gsl_reciprocal(UnitCell *in, gsl_matrix *m); extern UnitCell *cell_transform_rational(UnitCell *cell, RationalMatrix *m); extern UnitCell *cell_transform_rational_inverse(UnitCell *cell, RationalMatrix *m); diff --git a/libcrystfel/src/taketwo.c b/libcrystfel/src/taketwo.c index 37a29722..572ef793 100644 --- a/libcrystfel/src/taketwo.c +++ b/libcrystfel/src/taketwo.c @@ -2059,7 +2059,7 @@ static UnitCell *run_taketwo(UnitCell *cell, const struct taketwo_options *opts, tp->numPrevs++; /* Prepare the solution for CrystFEL friendliness */ - result = cell_transform_gsl_reciprocal(cell, solution); + result = cell_transform_gsl_direct(cell, solution); cleanup_taketwo_cell(&ttCell); return result; -- cgit v1.2.3 From 81b02525b7843fdb4d158d818bc7c5031c265456 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Mon, 11 Mar 2019 16:49:17 +0100 Subject: Remove commented code --- libcrystfel/src/taketwo.c | 29 ----------------------------- 1 file changed, 29 deletions(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/taketwo.c b/libcrystfel/src/taketwo.c index 572ef793..678b5c72 100644 --- a/libcrystfel/src/taketwo.c +++ b/libcrystfel/src/taketwo.c @@ -464,35 +464,6 @@ static void closest_rot_mat(struct rvec vec1, struct rvec vec2, rotation_around_axis(axis, bestAngle, twizzle); } -/* -static double matrix_angle(gsl_matrix *m) -{ - double a = gsl_matrix_get(m, 0, 0); - double b = gsl_matrix_get(m, 1, 1); - double c = gsl_matrix_get(m, 2, 2); - - double cos_t = (a + b + c - 1) / 2; - double theta = acos(cos_t); - - return theta; -} - -static struct rvec matrix_axis(gsl_matrix *a) -{ - double ang = matrix_angle(a); - double cos_t = cos(ang); - double p = gsl_matrix_get(a, 0, 0); - double q = gsl_matrix_get(a, 1, 1); - double r = gsl_matrix_get(a, 2, 2); - double x = sqrt((p - cos_t) / (1 - cos_t)); - double y = sqrt((q - cos_t) / (1 - cos_t)); - double z = sqrt((r - cos_t) / (1 - cos_t)); - - struct rvec v = new_rvec(x, y, z); - return v; -} -*/ - static double matrix_trace(gsl_matrix *a) { int i; -- cgit v1.2.3 From 2fa8ae0c7b042891d3f1dccedbd124e840dda68e Mon Sep 17 00:00:00 2001 From: Thomas White Date: Mon, 11 Mar 2019 17:49:11 +0100 Subject: TakeTwo: take change in IntegerMatrix notation into account --- libcrystfel/src/taketwo.c | 23 ++++++++++++++++------- 1 file changed, 16 insertions(+), 7 deletions(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/taketwo.c b/libcrystfel/src/taketwo.c index 678b5c72..458f90f9 100644 --- a/libcrystfel/src/taketwo.c +++ b/libcrystfel/src/taketwo.c @@ -1573,7 +1573,8 @@ static unsigned int start_seeds(gsl_matrix **rotation, struct TakeTwoCell *cell) } -static void set_gsl_matrix(gsl_matrix *mat, double asx, double asy, double asz, +static void set_gsl_matrix(gsl_matrix *mat, + double asx, double asy, double asz, double bsx, double bsy, double bsz, double csx, double csy, double csz) { @@ -1602,15 +1603,23 @@ static int generate_rotation_sym_ops(struct TakeTwoCell *ttCell) gsl_matrix *recip = gsl_matrix_alloc(3, 3); gsl_matrix *cart = gsl_matrix_alloc(3, 3); - cell_get_reciprocal(ttCell->cell, &asx, &asy, &asz, &bsx, &bsy, - &bsz, &csx, &csy, &csz); - set_gsl_matrix(recip, asx, asy, asz, bsx, bsy, bsz, csx, csy, csz); + cell_get_reciprocal(ttCell->cell, &asx, &asy, &asz, + &bsx, &bsy, &bsz, + &csx, &csy, &csz); - cell_get_cartesian(ttCell->cell, &asx, &asy, &asz, &bsx, &bsy, - &bsz, &csx, &csy, &csz); + set_gsl_matrix(recip, asx, asy, asz, + asx, bsy, bsz, + csx, csy, csz); + + cell_get_cartesian(ttCell->cell, &asx, &asy, &asz, + &bsx, &bsy, &bsz, + &csx, &csy, &csz); + + set_gsl_matrix(cart, asx, bsx, csx, + asy, bsy, csy, + asz, bsz, csz); - set_gsl_matrix(cart, asx, asy, asz, bsx, bsy, bsz, csx, csy, csz); int i, j, k; int numOps = num_equivs(rawList, NULL); -- cgit v1.2.3 From f32d46e54f5b9d645fb65accec52126099b49589 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Mon, 11 Mar 2019 18:26:16 +0100 Subject: TakeTwo: Restore correct matrix multiplication for solution --- libcrystfel/src/taketwo.c | 50 ++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 49 insertions(+), 1 deletion(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/taketwo.c b/libcrystfel/src/taketwo.c index 458f90f9..933fa76a 100644 --- a/libcrystfel/src/taketwo.c +++ b/libcrystfel/src/taketwo.c @@ -356,6 +356,54 @@ static void show_rvec(struct rvec r2) * functions called under the core functions, still specialised (Level 3) * ------------------------------------------------------------------------*/ +/* cell_transform_gsl_direct() doesn't do quite what we want here. + * The matrix m should be post-multiplied by a matrix of real or reciprocal + * basis vectors (it doesn't matter which because it's just a rotation). + * M contains the basis vectors written in columns: M' = mM */ +static UnitCell *cell_post_smiley_face(UnitCell *in, gsl_matrix *m) +{ + gsl_matrix *c; + double asx, asy, asz; + double bsx, bsy, bsz; + double csx, csy, csz; + gsl_matrix *res; + UnitCell *out; + + cell_get_cartesian(in, &asx, &asy, &asz, + &bsx, &bsy, &bsz, + &csx, &csy, &csz); + + c = gsl_matrix_alloc(3, 3); + gsl_matrix_set(c, 0, 0, asx); + gsl_matrix_set(c, 1, 0, asy); + gsl_matrix_set(c, 2, 0, asz); + gsl_matrix_set(c, 0, 1, bsx); + gsl_matrix_set(c, 1, 1, bsy); + gsl_matrix_set(c, 2, 1, bsz); + gsl_matrix_set(c, 0, 2, csx); + gsl_matrix_set(c, 1, 2, csy); + gsl_matrix_set(c, 2, 2, csz); + + res = gsl_matrix_calloc(3, 3); + gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, m, c, 0.0, res); + + out = cell_new_from_cell(in); + cell_set_cartesian(out, gsl_matrix_get(res, 0, 0), + gsl_matrix_get(res, 1, 0), + gsl_matrix_get(res, 2, 0), + gsl_matrix_get(res, 0, 1), + gsl_matrix_get(res, 1, 1), + gsl_matrix_get(res, 2, 1), + gsl_matrix_get(res, 0, 2), + gsl_matrix_get(res, 1, 2), + gsl_matrix_get(res, 2, 2)); + + gsl_matrix_free(res); + gsl_matrix_free(c); + return out; +} + + static void rotation_around_axis(struct rvec c, double th, gsl_matrix *res) { @@ -2039,7 +2087,7 @@ static UnitCell *run_taketwo(UnitCell *cell, const struct taketwo_options *opts, tp->numPrevs++; /* Prepare the solution for CrystFEL friendliness */ - result = cell_transform_gsl_direct(cell, solution); + result = cell_post_smiley_face(cell, solution); cleanup_taketwo_cell(&ttCell); return result; -- cgit v1.2.3 From ee77fdf67c9bb5430c7a8eb5005673e507964931 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Tue, 12 Mar 2019 14:31:26 +0100 Subject: Get rid of two unused functions --- libcrystfel/src/integer_matrix.c | 64 ---------------------------------------- libcrystfel/src/integer_matrix.h | 4 --- 2 files changed, 68 deletions(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/integer_matrix.c b/libcrystfel/src/integer_matrix.c index 1b448771..b6f55ac2 100644 --- a/libcrystfel/src/integer_matrix.c +++ b/libcrystfel/src/integer_matrix.c @@ -177,70 +177,6 @@ signed int intmat_get(const IntegerMatrix *m, unsigned int i, unsigned int j) } -/* intmat_intvec_mult: - * @m: An %IntegerMatrix - * @vec: An array of floats - * @ans: An array of floats in which to store the result - * - * Multiplies the matrix @m by the vector @vec. The size of @vec must equal the - * number of columns in @m, and the size of the result equals the number of rows - * in @m. - * - * Returns: non-zero on error - **/ -int intmat_floatvec_mult(const IntegerMatrix *m, const float *vec, float *ans) -{ - unsigned int i; - - for ( i=0; irows; i++ ) { - - unsigned int j; - - ans[i] = 0.0; - for ( j=0; jcols; j++ ) { - ans[i] += intmat_get(m, i, j) * vec[j]; - } - - } - - return 0; -} - - -/* intmat_rationalvec_mult: - * @m: An %IntegerMatrix - * @vec: An array of %Rational - * @ans: An array of %Rational in which to store the result - * - * Multiplies the matrix @m by the vector @vec. The size of @vec must equal the - * number of columns in @m, and the size of the result equals the number of rows - * in @m. - * - * Returns: non-zero on error - **/ -int intmat_rationalvec_mult(const IntegerMatrix *m, const Rational *vec, - Rational *ans) -{ - unsigned int i; - - - for ( i=0; irows; i++ ) { - - unsigned int j; - - ans[i] = rtnl_zero(); - for ( j=0; jcols; j++ ) { - Rational t; - t = rtnl_mul(vec[j], rtnl(intmat_get(m, i, j), 1)); - ans[i] = rtnl_add(ans[i], t); - } - - } - - return 0; -} - - /* intmat_solve_rational: * @m: An %IntegerMatrix * @vec: An array of %Rational diff --git a/libcrystfel/src/integer_matrix.h b/libcrystfel/src/integer_matrix.h index aeea38dd..2da2328a 100644 --- a/libcrystfel/src/integer_matrix.h +++ b/libcrystfel/src/integer_matrix.h @@ -71,12 +71,8 @@ extern IntegerMatrix *intmat_create_3x3(signed int m11, signed int m12, signed i signed int m31, signed int m32, signed int m33); /* Matrix-vector multiplication */ -extern int intmat_floatvec_mult(const IntegerMatrix *m, const float *vec, - float *ans); extern signed int *intmat_intvec_mult(const IntegerMatrix *m, const signed int *vec); -extern int intmat_rationalvec_mult(const IntegerMatrix *m, const Rational *vec, - Rational *ans); /* Matrix-matrix multiplication */ extern IntegerMatrix *intmat_intmat_mult(const IntegerMatrix *a, -- cgit v1.2.3 From 4fdfdfa0eef645e290a8ab98c9ff6507c2bc1f5f Mon Sep 17 00:00:00 2001 From: Thomas White Date: Tue, 12 Mar 2019 14:31:49 +0100 Subject: Rename intmat_intvec_mult to transform_indices --- libcrystfel/src/integer_matrix.c | 29 ++++++++++++++++------------- libcrystfel/src/integer_matrix.h | 3 +-- libcrystfel/src/symmetry.c | 2 +- 3 files changed, 18 insertions(+), 16 deletions(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/integer_matrix.c b/libcrystfel/src/integer_matrix.c index b6f55ac2..e74e0ea4 100644 --- a/libcrystfel/src/integer_matrix.c +++ b/libcrystfel/src/integer_matrix.c @@ -183,7 +183,7 @@ signed int intmat_get(const IntegerMatrix *m, unsigned int i, unsigned int j) * @ans: An array of %Rational in which to store the result * * Solves the matrix equation m*ans = vec, where @ans and @vec are - * vectors of %Rational. + * column vectors of %Rational numbers. * * This is just a convenience function which creates a %RationalMatrix out of * @m and then calls rtnl_mtx_solve(). @@ -204,35 +204,38 @@ int intmat_solve_rational(const IntegerMatrix *m, const Rational *vec, /** - * intmat_intvec_mult: - * @m: An %IntegerMatrix - * @vec: An array of signed integers + * transform_indices: + * @P: An %IntegerMatrix + * @hkl: An array of signed integers + * + * Apply transformation matrix P to a set of reciprocal space Miller indices. * - * Multiplies the matrix @m by the vector @vec. The size of @vec must equal the - * number of columns in @m, and the size of the result equals the number of rows - * in @m. + * In other words: + * Multiplies the matrix @P by the row vector @hkl. The size of @vec must equal + * the number of columns in @P, and the size of the result equals the number of + * rows in @P. * * The multiplication looks like this: - * (a1, a2, a3) = (vec1, vec2, vec3) m + * (a1, a2, a3) = (hkl1, hkl2, hkl3) P * Therefore matching the notation in ITA chapter 5.1. * * Returns: a newly allocated array of signed integers containing the answer, * or NULL on error. **/ -signed int *intmat_intvec_mult(const IntegerMatrix *m, const signed int *vec) +signed int *transform_indices(const IntegerMatrix *P, const signed int *hkl) { signed int *ans; unsigned int j; - ans = malloc(m->rows * sizeof(signed int)); + ans = malloc(P->rows * sizeof(signed int)); if ( ans == NULL ) return NULL; - for ( j=0; jcols; j++ ) { + for ( j=0; jcols; j++ ) { unsigned int i; ans[j] = 0; - for ( i=0; irows; i++ ) { - ans[i] += intmat_get(m, i, j) * vec[j]; + for ( i=0; irows; i++ ) { + ans[i] += intmat_get(P, i, j) * hkl[j]; } } diff --git a/libcrystfel/src/integer_matrix.h b/libcrystfel/src/integer_matrix.h index 2da2328a..9a2e4ce2 100644 --- a/libcrystfel/src/integer_matrix.h +++ b/libcrystfel/src/integer_matrix.h @@ -71,8 +71,7 @@ extern IntegerMatrix *intmat_create_3x3(signed int m11, signed int m12, signed i signed int m31, signed int m32, signed int m33); /* Matrix-vector multiplication */ -extern signed int *intmat_intvec_mult(const IntegerMatrix *m, - const signed int *vec); +extern signed int *transform_indices(const IntegerMatrix *P, const signed int *hkl); /* Matrix-matrix multiplication */ extern IntegerMatrix *intmat_intmat_mult(const IntegerMatrix *a, diff --git a/libcrystfel/src/symmetry.c b/libcrystfel/src/symmetry.c index 05ee7c6b..9b5c683d 100644 --- a/libcrystfel/src/symmetry.c +++ b/libcrystfel/src/symmetry.c @@ -1126,7 +1126,7 @@ static void do_op(const IntegerMatrix *op, v[0] = h; v[1] = k; v[2] = l; - ans = intmat_intvec_mult(op, v); + ans = transform_indices(op, v); assert(ans != NULL); *he = ans[0]; *ke = ans[1]; *le = ans[2]; -- cgit v1.2.3 From 15e772078538ea071b873526d2ed64f965731bcb Mon Sep 17 00:00:00 2001 From: Thomas White Date: Wed, 13 Mar 2019 14:22:10 +0100 Subject: Fix matrix notation in symmetry operations --- libcrystfel/src/symmetry.c | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/symmetry.c b/libcrystfel/src/symmetry.c index 9b5c683d..69dada16 100644 --- a/libcrystfel/src/symmetry.c +++ b/libcrystfel/src/symmetry.c @@ -371,13 +371,13 @@ static void expand_ops(SymOpList *s) /* Transform all the operations in a SymOpList by a given matrix. * The matrix must have a determinant of +/- 1 (otherwise its inverse would * not also be an integer matrix). */ -static void transform_ops(SymOpList *s, IntegerMatrix *t) +static void transform_ops(SymOpList *s, IntegerMatrix *P) { int n, i; - IntegerMatrix *inv; + IntegerMatrix *Pi; signed int det; - det = intmat_det(t); + det = intmat_det(P); if ( det == -1 ) { ERROR("WARNING: mirrored SymOpList.\n"); } else if ( det != 1 ) { @@ -385,8 +385,8 @@ static void transform_ops(SymOpList *s, IntegerMatrix *t) return; } - inv = intmat_inverse(t); - if ( inv == NULL ) { + Pi = intmat_inverse(P); + if ( Pi == NULL ) { ERROR("Failed to invert matrix.\n"); return; } @@ -396,13 +396,13 @@ static void transform_ops(SymOpList *s, IntegerMatrix *t) IntegerMatrix *r, *f; - r = intmat_intmat_mult(s->ops[i], t); + r = intmat_intmat_mult(P, s->ops[i]); if ( r == NULL ) { ERROR("Matrix multiplication failed.\n"); return; } - f = intmat_intmat_mult(inv, r); + f = intmat_intmat_mult(r, Pi); if ( f == NULL ) { ERROR("Matrix multiplication failed.\n"); return; @@ -415,7 +415,7 @@ static void transform_ops(SymOpList *s, IntegerMatrix *t) } - intmat_free(inv); + intmat_free(Pi); } @@ -1014,14 +1014,14 @@ static SymOpList *getpg_arbitrary_ua(const char *sym, size_t s) case 'a' : intmat_set(t, 0, 2, 1); - intmat_set(t, 1, 1, 1); - intmat_set(t, 2, 0, -1); + intmat_set(t, 1, 0, 1); + intmat_set(t, 2, 1, 1); break; case 'b' : - intmat_set(t, 0, 0, 1); + intmat_set(t, 0, 1, 1); intmat_set(t, 1, 2, 1); - intmat_set(t, 2, 1, -1); + intmat_set(t, 2, 0, 1); break; -- cgit v1.2.3 From bc84a4f8659244e11fddcb509664a2121bc40279 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Wed, 13 Mar 2019 15:41:05 +0100 Subject: Rename rtnl_mtx_solve and rtnl_mtx_mult to make their purpose clearer --- libcrystfel/src/cell.c | 4 ++-- libcrystfel/src/rational.c | 8 +++++--- libcrystfel/src/rational.h | 10 ++++++---- 3 files changed, 13 insertions(+), 9 deletions(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/cell.c b/libcrystfel/src/cell.c index 9ae77a6d..63b84fc3 100644 --- a/libcrystfel/src/cell.c +++ b/libcrystfel/src/cell.c @@ -736,7 +736,7 @@ static void check_point_fwd(RationalMatrix *m, CenteringMask *cmask, Rational nc[3]; /* Transform the lattice point */ - rtnl_mtx_solve(m, c, nc); + transform_fractional_coords_rtnl(m, c, nc); /* Eliminate any centerings which don't include the transformed point */ maybe_eliminate(CMASK_P, cmask, nc, 'P'); @@ -759,7 +759,7 @@ static void check_point_bwd(RationalMatrix *m, CenteringMask *mask, Rational nc[3]; Rational c[3] = {x, y, z}; - rtnl_mtx_mult(m, c, nc); + transform_fractional_coords_rtnl_inverse(m, c, nc); if ( !centering_has_point(cen, nc) ) { *mask |= exclude; diff --git a/libcrystfel/src/rational.c b/libcrystfel/src/rational.c index ab1eff91..c1512b8d 100644 --- a/libcrystfel/src/rational.c +++ b/libcrystfel/src/rational.c @@ -424,7 +424,8 @@ void rtnl_mtx_free(RationalMatrix *mtx) * * Returns: non-zero on error **/ -int rtnl_mtx_solve(const RationalMatrix *m, const Rational *ivec, Rational *ans) +int transform_fractional_coords_rtnl(const RationalMatrix *m, + const Rational *ivec, Rational *ans) { RationalMatrix *cm; Rational *vec; @@ -571,8 +572,9 @@ void rtnl_mtx_mtxmult(const RationalMatrix *A, const RationalMatrix *B, } -/* Calculate ans = m.vec, where 'ans' and 'vec' are column vectors */ -void rtnl_mtx_mult(const RationalMatrix *m, const Rational *vec, Rational *ans) +void transform_fractional_coords_rtnl_inverse(const RationalMatrix *m, + const Rational *vec, + Rational *ans) { int i, j; diff --git a/libcrystfel/src/rational.h b/libcrystfel/src/rational.h index 9ed20bfd..23c918cf 100644 --- a/libcrystfel/src/rational.h +++ b/libcrystfel/src/rational.h @@ -91,10 +91,12 @@ 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_mult(const RationalMatrix *m, const Rational *vec, - Rational *ans); -extern int rtnl_mtx_solve(const RationalMatrix *m, const Rational *vec, - Rational *ans); +extern int transform_fractional_coords_rtnl(const RationalMatrix *P, + const Rational *ivec, + Rational *ans); +extern void transform_fractional_coords_rtnl_inverse(const RationalMatrix *P, + const Rational *vec, + Rational *ans); extern void rtnl_mtx_print(const RationalMatrix *m); extern Rational rtnl_mtx_det(const RationalMatrix *m); -- cgit v1.2.3 From 74fe33b767d9c94b5c22e274a7518937e8a7adf4 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Wed, 13 Mar 2019 15:41:38 +0100 Subject: Rename m to P, to match documentation (including ITA) --- libcrystfel/src/cell.c | 40 ++++++++++++++++++++-------------------- libcrystfel/src/rational.c | 44 +++++++++++++++++++++++++++----------------- 2 files changed, 47 insertions(+), 37 deletions(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/cell.c b/libcrystfel/src/cell.c index 63b84fc3..1dcb6bfb 100644 --- a/libcrystfel/src/cell.c +++ b/libcrystfel/src/cell.c @@ -729,14 +729,14 @@ static void maybe_eliminate(CenteringMask c, CenteringMask *cmask, Rational *nc, /* Check if the point x,y,z in the original cell matches any lattice point * in the transformed cell */ -static void check_point_fwd(RationalMatrix *m, CenteringMask *cmask, +static void check_point_fwd(RationalMatrix *P, CenteringMask *cmask, Rational x, Rational y, Rational z) { Rational c[3] = {x, y, z}; Rational nc[3]; /* Transform the lattice point */ - transform_fractional_coords_rtnl(m, c, nc); + transform_fractional_coords_rtnl(P, c, nc); /* Eliminate any centerings which don't include the transformed point */ maybe_eliminate(CMASK_P, cmask, nc, 'P'); @@ -752,14 +752,14 @@ static void check_point_fwd(RationalMatrix *m, CenteringMask *cmask, /* Check if the point x,y,z in the transformed cell matches any lattice point * in the original cell. If not, eliminate "exclude" from "*mask". */ -static void check_point_bwd(RationalMatrix *m, CenteringMask *mask, +static void check_point_bwd(RationalMatrix *P, CenteringMask *mask, char cen, CenteringMask exclude, Rational x, Rational y, Rational z) { Rational nc[3]; Rational c[3] = {x, y, z}; - transform_fractional_coords_rtnl_inverse(m, c, nc); + transform_fractional_coords_rtnl_inverse(P, c, nc); if ( !centering_has_point(cen, nc) ) { *mask |= exclude; @@ -788,19 +788,19 @@ static char cmask_decode(CenteringMask mask) } -static char determine_centering(RationalMatrix *m, char cen) +static char determine_centering(RationalMatrix *P, char cen) { CenteringMask cmask = CMASK_ALL; /* Check whether the current centering can provide all the lattice * points for the transformed cell. Eliminate any centerings for which * it can't. */ - check_point_bwd(m, &cmask, cen, CMASK_A | CMASK_F, rtnl_zero(), rtnl(1,2), rtnl(1,2)); - check_point_bwd(m, &cmask, cen, CMASK_B | CMASK_F, rtnl(1,2), rtnl_zero(), rtnl(1,2)); - check_point_bwd(m, &cmask, cen, CMASK_C | CMASK_F, rtnl(1,2), rtnl(1,2), rtnl_zero()); - check_point_bwd(m, &cmask, cen, CMASK_I, rtnl(1,2), rtnl(1,2), rtnl(1,2)); - check_point_bwd(m, &cmask, cen, CMASK_H, rtnl(2,3), rtnl(1,3), rtnl(1,3)); - check_point_bwd(m, &cmask, cen, CMASK_H, rtnl(1,3), rtnl(2,3), rtnl(2,3)); + check_point_bwd(P, &cmask, cen, CMASK_A | CMASK_F, rtnl_zero(), rtnl(1,2), rtnl(1,2)); + check_point_bwd(P, &cmask, cen, CMASK_B | CMASK_F, rtnl(1,2), rtnl_zero(), rtnl(1,2)); + check_point_bwd(P, &cmask, cen, CMASK_C | CMASK_F, rtnl(1,2), rtnl(1,2), rtnl_zero()); + check_point_bwd(P, &cmask, cen, CMASK_I, rtnl(1,2), rtnl(1,2), rtnl(1,2)); + check_point_bwd(P, &cmask, cen, CMASK_H, rtnl(2,3), rtnl(1,3), rtnl(1,3)); + check_point_bwd(P, &cmask, cen, CMASK_H, rtnl(1,3), rtnl(2,3), rtnl(2,3)); /* Check whether the current centering's lattice points will all * coincide with lattice points in the new centering. Eliminate any @@ -812,30 +812,30 @@ static char determine_centering(RationalMatrix *m, char cen) break; case 'A' : - check_point_fwd(m, &cmask, rtnl_zero(), rtnl(1,2), rtnl(1,2)); + check_point_fwd(P, &cmask, rtnl_zero(), rtnl(1,2), rtnl(1,2)); break; case 'B' : - check_point_fwd(m, &cmask, rtnl(1,2), rtnl_zero(), rtnl(1,2)); + check_point_fwd(P, &cmask, rtnl(1,2), rtnl_zero(), rtnl(1,2)); break; case 'C' : - check_point_fwd(m, &cmask, rtnl(1,2), rtnl(1,2), rtnl_zero()); + check_point_fwd(P, &cmask, rtnl(1,2), rtnl(1,2), rtnl_zero()); break; case 'I' : - check_point_fwd(m, &cmask, rtnl(1,2), rtnl(1,2), rtnl(1,2)); + check_point_fwd(P, &cmask, rtnl(1,2), rtnl(1,2), rtnl(1,2)); break; case 'F' : - check_point_fwd(m, &cmask, rtnl_zero(), rtnl(1,2), rtnl(1,2)); - check_point_fwd(m, &cmask, rtnl(1,2), rtnl_zero(), rtnl(1,2)); - check_point_fwd(m, &cmask, rtnl(1,2), rtnl(1,2), rtnl_zero()); + check_point_fwd(P, &cmask, rtnl_zero(), rtnl(1,2), rtnl(1,2)); + check_point_fwd(P, &cmask, rtnl(1,2), rtnl_zero(), rtnl(1,2)); + check_point_fwd(P, &cmask, rtnl(1,2), rtnl(1,2), rtnl_zero()); break; case 'H' : - check_point_fwd(m, &cmask, rtnl(2,3), rtnl(1,3), rtnl(1,3)); - check_point_fwd(m, &cmask, rtnl(1,3), rtnl(2,3), rtnl(2,3)); + check_point_fwd(P, &cmask, rtnl(2,3), rtnl(1,3), rtnl(1,3)); + check_point_fwd(P, &cmask, rtnl(1,3), rtnl(2,3), rtnl(2,3)); break; } diff --git a/libcrystfel/src/rational.c b/libcrystfel/src/rational.c index c1512b8d..65b209ff 100644 --- a/libcrystfel/src/rational.c +++ b/libcrystfel/src/rational.c @@ -409,7 +409,7 @@ void rtnl_mtx_free(RationalMatrix *mtx) /* rtnl_mtx_solve: - * @m: A %RationalMatrix + * @P: A %RationalMatrix * @vec: An array of %Rational * @ans: An array of %Rational in which to store the result * @@ -422,9 +422,13 @@ void rtnl_mtx_free(RationalMatrix *mtx) * of rows in @m must equal the length of @vec, but note that there is no way * for this function to check that this is the case. * + * Given that P is transformation of the unit cell axes (see ITA chapter 5.1), + * this function calculates the fractional coordinates of a point in the + * transformed axes, given its fractional coordinates in the original axes. + * * Returns: non-zero on error **/ -int transform_fractional_coords_rtnl(const RationalMatrix *m, +int transform_fractional_coords_rtnl(const RationalMatrix *P, const Rational *ivec, Rational *ans) { RationalMatrix *cm; @@ -432,28 +436,28 @@ int transform_fractional_coords_rtnl(const RationalMatrix *m, int h, k; int i; - if ( m->rows != m->cols ) return 1; + if ( P->rows != P->cols ) return 1; /* Copy the matrix and vector because the calculation will * be done in-place */ - cm = rtnl_mtx_copy(m); + cm = rtnl_mtx_copy(P); if ( cm == NULL ) return 1; - vec = malloc(m->rows*sizeof(Rational)); + vec = malloc(cm->rows*sizeof(Rational)); if ( vec == NULL ) return 1; - for ( h=0; hrows; h++ ) vec[h] = ivec[h]; + for ( h=0; hrows; h++ ) vec[h] = ivec[h]; /* Gaussian elimination with partial pivoting */ h = 0; k = 0; - while ( h<=m->rows && k<=m->cols ) { + while ( h<=cm->rows && k<=cm->cols ) { int prow = 0; Rational pval = rtnl_zero(); Rational t; /* Find the row with the largest value in column k */ - for ( i=h; irows; i++ ) { + for ( i=h; irows; i++ ) { if ( rtnl_cmp(rtnl_abs(rtnl_mtx_get(cm, i, k)), pval) > 0 ) { pval = rtnl_abs(rtnl_mtx_get(cm, i, k)); prow = i; @@ -466,7 +470,7 @@ int transform_fractional_coords_rtnl(const RationalMatrix *m, } /* Swap 'prow' with row h */ - for ( i=0; icols; i++ ) { + for ( i=0; icols; i++ ) { t = rtnl_mtx_get(cm, h, i); rtnl_mtx_set(cm, h, i, rtnl_mtx_get(cm, prow, i)); rtnl_mtx_set(cm, prow, i, t); @@ -476,7 +480,7 @@ int transform_fractional_coords_rtnl(const RationalMatrix *m, vec[h] = t; /* Divide and subtract rows below */ - for ( i=h+1; irows; i++ ) { + for ( i=h+1; irows; i++ ) { int j; Rational dval; @@ -484,7 +488,7 @@ int transform_fractional_coords_rtnl(const RationalMatrix *m, dval = rtnl_div(rtnl_mtx_get(cm, i, k), rtnl_mtx_get(cm, h, k)); - for ( j=0; jcols; j++ ) { + for ( j=0; jcols; j++ ) { Rational t = rtnl_mtx_get(cm, i, j); Rational p = rtnl_mul(dval, rtnl_mtx_get(cm, h, j)); t = rtnl_sub(t, p); @@ -504,10 +508,10 @@ int transform_fractional_coords_rtnl(const RationalMatrix *m, } /* Back-substitution */ - for ( i=m->rows-1; i>=0; i-- ) { + for ( i=cm->rows-1; i>=0; i-- ) { int j; Rational sum = rtnl_zero(); - for ( j=i+1; jcols; j++ ) { + for ( j=i+1; jcols; j++ ) { Rational av; av = rtnl_mul(rtnl_mtx_get(cm, i, j), ans[j]); sum = rtnl_add(sum, av); @@ -572,17 +576,23 @@ void rtnl_mtx_mtxmult(const RationalMatrix *A, const RationalMatrix *B, } -void transform_fractional_coords_rtnl_inverse(const RationalMatrix *m, +/* 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. + * To transform in the opposite direction, we would multiply by Q = P^-1. + * Therefore, this direction is the "easy way" and needs just a matrix + * multiplication. */ +void transform_fractional_coords_rtnl_inverse(const RationalMatrix *P, const Rational *vec, Rational *ans) { int i, j; - for ( i=0; irows; i++ ) { + for ( i=0; irows; i++ ) { ans[i] = rtnl_zero(); - for ( j=0; jcols; j++ ) { + for ( j=0; jcols; j++ ) { Rational add; - add = rtnl_mul(rtnl_mtx_get(m, i, j), vec[j]); + add = rtnl_mul(rtnl_mtx_get(P, i, j), vec[j]); ans[i] = rtnl_add(ans[i], add); } } -- cgit v1.2.3 From e2640716fb422a9de8b1611a464ec6bfe28dfe6f Mon Sep 17 00:00:00 2001 From: Thomas White Date: Wed, 13 Mar 2019 15:41:55 +0100 Subject: Remove intmat_solve_rational Not used anywhere --- libcrystfel/src/integer_matrix.c | 26 -------------------------- libcrystfel/src/integer_matrix.h | 3 --- 2 files changed, 29 deletions(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/integer_matrix.c b/libcrystfel/src/integer_matrix.c index e74e0ea4..90d58b35 100644 --- a/libcrystfel/src/integer_matrix.c +++ b/libcrystfel/src/integer_matrix.c @@ -177,32 +177,6 @@ signed int intmat_get(const IntegerMatrix *m, unsigned int i, unsigned int j) } -/* intmat_solve_rational: - * @m: An %IntegerMatrix - * @vec: An array of %Rational - * @ans: An array of %Rational in which to store the result - * - * Solves the matrix equation m*ans = vec, where @ans and @vec are - * column vectors of %Rational numbers. - * - * This is just a convenience function which creates a %RationalMatrix out of - * @m and then calls rtnl_mtx_solve(). - * - * Returns: non-zero on error - **/ -int intmat_solve_rational(const IntegerMatrix *m, const Rational *vec, - Rational *ans) -{ - RationalMatrix *rm; - int r; - - rm = rtnl_mtx_from_intmat(m); - r = rtnl_mtx_solve(rm, vec, ans); - rtnl_mtx_free(rm); - return r; -} - - /** * transform_indices: * @P: An %IntegerMatrix diff --git a/libcrystfel/src/integer_matrix.h b/libcrystfel/src/integer_matrix.h index 9a2e4ce2..6fb3e399 100644 --- a/libcrystfel/src/integer_matrix.h +++ b/libcrystfel/src/integer_matrix.h @@ -77,9 +77,6 @@ extern signed int *transform_indices(const IntegerMatrix *P, const signed int *h extern IntegerMatrix *intmat_intmat_mult(const IntegerMatrix *a, const IntegerMatrix *b); -extern int intmat_solve_rational(const IntegerMatrix *m, const Rational *vec, - Rational *ans); - /* Inverse */ extern IntegerMatrix *intmat_inverse(const IntegerMatrix *m); -- cgit v1.2.3 From fc5e8a45f99c7927af26bccf341144ef556673da Mon Sep 17 00:00:00 2001 From: Thomas White Date: Wed, 13 Mar 2019 15:53:31 +0100 Subject: validate_cell: Allow to tell the difference between serious problems and warnings --- libcrystfel/src/cell-utils.c | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c index 605aa449..a36ebe1c 100644 --- a/libcrystfel/src/cell-utils.c +++ b/libcrystfel/src/cell-utils.c @@ -1570,7 +1570,9 @@ int cell_is_sensible(UnitCell *cell) * lattice is a conventional Bravais lattice. * Warnings are printied if any of the checks are failed. * - * Returns: true if cell is invalid. + * Returns: zero if the cell is fine, 1 if it is unconventional but otherwise + * OK (e.g. left-handed or not a Bravais lattice), and 2 if there is a serious + * problem such as the parameters being physically impossible. * */ int validate_cell(UnitCell *cell) @@ -1580,7 +1582,7 @@ int validate_cell(UnitCell *cell) if ( cell_has_parameters(cell) && !cell_is_sensible(cell) ) { ERROR("WARNING: Unit cell parameters are not sensible.\n"); - err = 1; + err = 2; } if ( !bravais_lattice(cell) ) { @@ -1604,7 +1606,7 @@ int validate_cell(UnitCell *cell) || ((cen == 'C') && (ua == 'c')) ) { ERROR("WARNING: A, B or C centering matches unique" " axis.\n"); - err = 1; + err = 2; } } -- cgit v1.2.3 From c229a577b5c5578f1e83f8a6ee569bb5681524b9 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Wed, 13 Mar 2019 16:10:08 +0100 Subject: determine_centering: Check 1,1,1 in both directions The original centering must be able to provide a lattice point for 1,1,1 in the new centering. The new centering's 1,1,1 point must also map to a lattice point in the original centering. Note that it's no use checking 0,0,0, because it'll just be a load of zeroes. --- libcrystfel/src/cell.c | 8 ++++++++ 1 file changed, 8 insertions(+) (limited to 'libcrystfel') diff --git a/libcrystfel/src/cell.c b/libcrystfel/src/cell.c index 1dcb6bfb..1e463bf7 100644 --- a/libcrystfel/src/cell.c +++ b/libcrystfel/src/cell.c @@ -801,6 +801,7 @@ static char determine_centering(RationalMatrix *P, char cen) check_point_bwd(P, &cmask, cen, CMASK_I, rtnl(1,2), rtnl(1,2), rtnl(1,2)); check_point_bwd(P, &cmask, cen, CMASK_H, rtnl(2,3), rtnl(1,3), rtnl(1,3)); check_point_bwd(P, &cmask, cen, CMASK_H, rtnl(1,3), rtnl(2,3), rtnl(2,3)); + check_point_bwd(P, &cmask, cen, CMASK_ALL, rtnl(1,1), rtnl(1,1), rtnl(1,1)); /* Check whether the current centering's lattice points will all * coincide with lattice points in the new centering. Eliminate any @@ -809,31 +810,38 @@ static char determine_centering(RationalMatrix *P, char cen) case 'P' : case 'R' : + check_point_fwd(P, &cmask, rtnl(1,1), rtnl(1,1), rtnl(1,1)); break; case 'A' : + check_point_fwd(P, &cmask, rtnl(1,1), rtnl(1,1), rtnl(1,1)); check_point_fwd(P, &cmask, rtnl_zero(), rtnl(1,2), rtnl(1,2)); break; case 'B' : + check_point_fwd(P, &cmask, rtnl(1,1), rtnl(1,1), rtnl(1,1)); check_point_fwd(P, &cmask, rtnl(1,2), rtnl_zero(), rtnl(1,2)); break; case 'C' : + check_point_fwd(P, &cmask, rtnl(1,1), rtnl(1,1), rtnl(1,1)); check_point_fwd(P, &cmask, rtnl(1,2), rtnl(1,2), rtnl_zero()); break; case 'I' : + check_point_fwd(P, &cmask, rtnl(1,1), rtnl(1,1), rtnl(1,1)); check_point_fwd(P, &cmask, rtnl(1,2), rtnl(1,2), rtnl(1,2)); break; case 'F' : + check_point_fwd(P, &cmask, rtnl(1,1), rtnl(1,1), rtnl(1,1)); check_point_fwd(P, &cmask, rtnl_zero(), rtnl(1,2), rtnl(1,2)); check_point_fwd(P, &cmask, rtnl(1,2), rtnl_zero(), rtnl(1,2)); check_point_fwd(P, &cmask, rtnl(1,2), rtnl(1,2), rtnl_zero()); break; case 'H' : + check_point_fwd(P, &cmask, rtnl(1,1), rtnl(1,1), rtnl(1,1)); check_point_fwd(P, &cmask, rtnl(2,3), rtnl(1,3), rtnl(1,3)); check_point_fwd(P, &cmask, rtnl(1,3), rtnl(2,3), rtnl(2,3)); break; -- cgit v1.2.3 From a527e38ccc907a64420a2bd73245a0512c8baa87 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Wed, 13 Mar 2019 16:56:54 +0100 Subject: cell_transform_rational: Complain if matrix is singular --- libcrystfel/src/cell.c | 4 ++++ 1 file changed, 4 insertions(+) (limited to 'libcrystfel') diff --git a/libcrystfel/src/cell.c b/libcrystfel/src/cell.c index 1e463bf7..ce9d37bb 100644 --- a/libcrystfel/src/cell.c +++ b/libcrystfel/src/cell.c @@ -876,9 +876,13 @@ UnitCell *cell_transform_rational(UnitCell *cell, RationalMatrix *m) gsl_matrix *tm; char ncen; int i, j; + Rational det; if ( m == NULL ) return NULL; + det = rtnl_mtx_det(m); + if ( rtnl_cmp(det, rtnl_zero()) == 0 ) return NULL; + tm = gsl_matrix_alloc(3,3); if ( tm == NULL ) { return NULL; -- cgit v1.2.3