diff options
Diffstat (limited to 'tests')
| -rw-r--r-- | tests/CMakeLists.txt | 4 | ||||
| -rw-r--r-- | tests/centering_check.c | 22 | ||||
| -rw-r--r-- | tests/rational_check.c | 194 | ||||
| -rw-r--r-- | tests/transformation_check.c | 227 |
4 files changed, 355 insertions, 92 deletions
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 02f3ca02..beb141c3 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -83,3 +83,7 @@ if (HAVE_OPENCL) add_test(gpu_sim_check gpu_sim_check) endif (HAVE_OPENCL) +add_executable(rational_check rational_check.c) +target_include_directories(rational_check PRIVATE ${COMMON_INCLUDES}) +target_link_libraries(rational_check ${COMMON_LIBRARIES}) +add_test(rational_check rational_check) diff --git a/tests/centering_check.c b/tests/centering_check.c index 887311dd..cc553dd0 100644 --- a/tests/centering_check.c +++ b/tests/centering_check.c @@ -61,7 +61,8 @@ static int check_centering(double a, double b, double c, { UnitCell *cell, *cref; UnitCell *n; - UnitCellTransformation *t; + IntegerMatrix *C; + RationalMatrix *Ci; int fail = 0; int i; double asx, asy, asz; @@ -86,14 +87,19 @@ static int check_centering(double a, double b, double c, cell_free(cref); check_cell(cell, "Input"); - n = uncenter_cell(cell, &t); + n = uncenter_cell(cell, &C, &Ci); if ( n != NULL ) { - STATUS("Transformation was:\n"); - tfn_print(t); + STATUS("The back transformation is:\n"); + intmat_print(C); if ( check_cell(n, "Output") ) fail = 1; - if ( !fail ) cell_print(n); + if ( intmat_is_identity(C) && ((cen=='P') || (cen=='R')) ) { + STATUS("Identity, as expected.\n"); + } else { + cell_print(n); + } } else { - fail = 1; + STATUS("************* Could not uncenter.\n"); + return 1; } cell_get_reciprocal(cell, &asx, &asy, &asz, @@ -299,10 +305,6 @@ int main(int argc, char *argv[]) fail += check_centering(20e-10, 20e-10, 40e-10, 90.0, 90.0, 120.0, L_HEXAGONAL, 'H', 'c', rng); - /* Cubic P */ - fail += check_centering(30e-10, 30e-10, 30e-10, 90.0, 90.0, 90.0, - L_CUBIC, 'P', '*', rng); - /* Cubic I */ fail += check_centering(30e-10, 30e-10, 30e-10, 90.0, 90.0, 90.0, L_CUBIC, 'I', '*', rng); diff --git a/tests/rational_check.c b/tests/rational_check.c new file mode 100644 index 00000000..484264d9 --- /dev/null +++ b/tests/rational_check.c @@ -0,0 +1,194 @@ +/* + * rational_check.c + * + * Check that rational numbers work + * + * Copyright © 2012-2019 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. + * + * Authors: + * 2012-2019 Thomas White <taw@physics.org> + * + * 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 <http://www.gnu.org/licenses/>. + * + */ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + + +#include <stdlib.h> +#include <stdio.h> +#include <stdarg.h> +#include <gsl/gsl_rng.h> + +#include <rational.h> +#include <utils.h> + + +static Rational gen_rtnl(gsl_rng *rng, double *dptr, int sz) +{ + Rational r; + signed int num, den; + do { + num = gsl_rng_uniform_int(rng, 2*sz) - sz; + den = gsl_rng_uniform_int(rng, 2*sz) - sz; + } while ( den == 0 ); + r = rtnl(num, den); + *dptr = (double)num/den; + return r; +} + + +static int test_rational(gsl_rng *rng) +{ + Rational r1, r2, r3; + double rd1, rd2, rd3; + char op; + + r1 = gen_rtnl(rng, &rd1, 100); + r2 = gen_rtnl(rng, &rd2, 100); + + switch ( gsl_rng_uniform_int(rng, 3) ) { + + case 0: + r3 = rtnl_mul(r1, r2); + rd3 = rd1*rd2; + op = '*'; + break; + + case 1: + r3 = rtnl_div(r1, r2); + rd3 = rd1/rd2; + op = '/'; + break; + + case 2: + r3 = rtnl_add(r1, r2); + rd3 = rd1+rd2; + op = '+'; + break; + + case 3: + r3 = rtnl_sub(r1, r2); + rd3 = rd1-rd2; + op = '-'; + break; + + default: + abort(); + + } + + if ( isinf(rd3) && isinf(rtnl_as_double(r3)) ) return 0; + + if ( fabs(rtnl_as_double(r3) - rd3) > 0.001 ) { + ERROR("%10s %c %10s = %10s", rtnl_format(r1), op, + rtnl_format(r2), rtnl_format(r3)); + ERROR(" ---> %10f %c %10f = %10f\n", rd1, op, rd2, rd3); + return 1; + } + + return 0; +} + + +static int test_rational_matrix(gsl_rng *rng) +{ + const int size = 3; + RationalMatrix *rm; + Rational *rvec; + Rational *rans; + gsl_matrix *gm; + gsl_vector *gvec; + gsl_vector *gans; + int i, j; + int filt; + int fail = 0; + + rm = rtnl_mtx_new(size, size); + gm = gsl_matrix_alloc(size, size); + rvec = malloc(size*sizeof(Rational)); + gvec = gsl_vector_alloc(size); + rans = malloc(size*sizeof(Rational)); + if ( (rm==NULL) || (gm==NULL) || (rvec==NULL) + || (gvec==NULL) || (rans==NULL) ) return 1; + + /* Find a matrix equation which actually works */ + do { + for ( i=0; i<size; i++ ) { + double d; + Rational r; + for ( j=0; j<size; j++ ) { + r = gen_rtnl(rng, &d, 5); + rtnl_mtx_set(rm, i, j, r); + gsl_matrix_set(gm ,i, j, d); + } + r = gen_rtnl(rng, &d, 5); + rvec[i] = r; + gsl_vector_set(gvec ,i, d); + } + + gans = solve_svd(gvec, gm, &filt, 0); + } while ( (gans==NULL) || (filt!=0) || (isnan(gsl_vector_get(gans,0))) ); + + transform_fractional_coords_rtnl(rm, rvec, rans); + + for ( i=0; i<size; i++ ) { + if ( fabs(rtnl_as_double(rans[i]) - gsl_vector_get(gans, i) > 0.001) ) + { + STATUS("%25s %10f %10f\n", rtnl_format(rans[i]), + rtnl_as_double(rans[i]), + gsl_vector_get(gans, i)); + fail = 1; + } + } + + gsl_matrix_free(gm); + rtnl_mtx_free(rm); + free(rans); + free(rvec); + gsl_vector_free(gvec); + gsl_vector_free(gans); + return fail; +} + + +int main(int argc, char *argv[]) +{ + int fail = 0; + gsl_rng *rng; + int i; + + rng = gsl_rng_alloc(gsl_rng_mt19937); + gsl_set_error_handler_off(); + + STATUS("Arithmetic test...\n"); + for ( i=0; i<1000; i++ ) { + fail += test_rational(rng); + } + + STATUS("Matrix test...\n"); + for ( i=0; i<1000; i++ ) { + fail += test_rational_matrix(rng); + } + + gsl_rng_free(rng); + + if ( fail != 0 ) return 1; + return 0; +} diff --git a/tests/transformation_check.c b/tests/transformation_check.c index 2e7e7378..2790d6f1 100644 --- a/tests/transformation_check.c +++ b/tests/transformation_check.c @@ -37,6 +37,7 @@ #include <cell.h> #include <cell-utils.h> +#include <symmetry.h> #define MAX_REFLS (10*1024) @@ -134,52 +135,61 @@ static int compare_rvecs(struct rvec *a, int na, struct rvec *b, int nb) } -static int check_transformation(UnitCell *cell, UnitCellTransformation *tfn, - int pred_test, UnitCell *ct) +static int check_same_reflections(UnitCell *cell, UnitCell *cnew) { - UnitCell *cnew, *cback; - UnitCellTransformation *inv; - double a[9], b[9]; - int i; - int fail = 0; struct rvec *vecs; struct rvec *tvecs; int na, nb; - STATUS("-----------------------\n"); - if ( ct == NULL ) { - cnew = cell_transform(cell, tfn); + /* Check that the two cells predict the same reflections */ + vecs = all_refls(cell, 1e9, &na); + tvecs = all_refls(cnew, 1e9, &nb); + if ( compare_rvecs(vecs, na, tvecs, nb) + || compare_rvecs(tvecs, nb, vecs, na) ) + { + ERROR("********************************************** "); + ERROR("Transformed cell didn't predict the same reflections\n"); + //printf("---\n"); + //for ( i=0; i<na; i++ ) { + // printf("%e %e %e\n", vecs[i].u, vecs[i].v, vecs[i].w); + //} + //printf("---\n"); + //for ( i=0; i<nb; i++ ) { + // printf("%e %e %e\n", tvecs[i].u, tvecs[i].v, tvecs[i].w); + //} + return 1; } else { - cnew = ct; + STATUS("The cells predict the same reflections.\n"); } - cback = cell_transform_inverse(cnew, tfn); + free(vecs); + free(tvecs); + return 0; +} + + +static int check_transformation(UnitCell *cell, IntegerMatrix *tfn, + int pred_test) +{ + UnitCell *cnew, *cback; + double a[9], b[9]; + int i; + int fail = 0; + + STATUS("-----------------------\n"); + cnew = cell_transform_intmat(cell, tfn); + cback = cell_transform_intmat_inverse(cnew, tfn); + STATUS("----> Before transformation:\n"); cell_print(cell); - tfn_print(tfn); + STATUS("----> The transformation matrix:\n"); + intmat_print(tfn); + STATUS("----> After transformation:\n"); cell_print(cnew); + STATUS("----> After back transformation:\n"); + cell_print(cback); if ( pred_test ) { - /* Check that the two cells predict the same reflections */ - vecs = all_refls(cell, 1e9, &na); - tvecs = all_refls(cnew, 1e9, &nb); - if ( compare_rvecs(vecs, na, tvecs, nb) - || compare_rvecs(tvecs, nb, vecs, na) ) - { - ERROR("Transformed cell didn't predict the same reflections\n"); - //printf("---\n"); - //for ( i=0; i<na; i++ ) { - // printf("%e %e %e\n", vecs[i].u, vecs[i].v, vecs[i].w); - //} - //printf("---\n"); - //for ( i=0; i<nb; i++ ) { - // printf("%e %e %e\n", tvecs[i].u, tvecs[i].v, tvecs[i].w); - //} - return 1; - } else { - STATUS("The cells predict the same reflections.\n"); - } - free(vecs); - free(tvecs); + check_same_reflections(cell, cnew); } else { STATUS("Cells not expected to predict the same reflections.\n"); } @@ -193,18 +203,18 @@ static int check_transformation(UnitCell *cell, UnitCellTransformation *tfn, &b[6], &b[7], &b[8]); for ( i=0; i<9; i++ ) { if ( !tolerance(a[i], b[i]) ) { + //STATUS("%e %e\n", a[i], b[i]); fail = 1; - STATUS("%e %e\n", a[i], b[i]); } } if ( fail ) { + ERROR("********************************************** "); ERROR("Original cell not recovered after transformation:\n"); - cell_print(cell); - tfn_print(tfn); - inv = tfn_inverse(tfn); - tfn_print(inv); + STATUS("----> After transformation and transformation back:\n"); cell_print(cback); + } else { + STATUS("The original cell was recovered after inverse transform.\n"); } return fail; @@ -214,21 +224,73 @@ static int check_transformation(UnitCell *cell, UnitCellTransformation *tfn, static int check_uncentering(UnitCell *cell) { UnitCell *ct; - UnitCellTransformation *tr; + IntegerMatrix *C; + RationalMatrix *Ci; + UnitCell *cback; + double a[9], b[9]; + int i; + int fail = 0; + + STATUS("-----------------------\n"); + + STATUS("----> Before transformation:\n"); + cell_print_full(cell); + + ct = uncenter_cell(cell, &C, &Ci); + if ( ct == NULL ) return 1; + + STATUS("----> The primitive unit cell:\n"); + cell_print(ct); - ct = uncenter_cell(cell, &tr); - return check_transformation(cell, tr, 1, ct); + STATUS("----> The matrix to put the centering back:\n"); + intmat_print(C); + + STATUS("----> The recovered centered cell:\n"); + cback = cell_transform_intmat(ct, C); + cell_print(cback); + + cell_get_cartesian(cell, &a[0], &a[1], &a[2], + &a[3], &a[4], &a[5], + &a[6], &a[7], &a[8]); + cell_get_cartesian(cback, &b[0], &b[1], &b[2], + &b[3], &b[4], &b[5], + &b[6], &b[7], &b[8]); + for ( i=0; i<9; i++ ) { + if ( fabs(a[i] - b[i]) > 1e-12 ) { + fail = 1; + } + } + + if ( fail ) { + ERROR("********************************************** "); + ERROR("Original cell not recovered after back transformation\n"); + } + + fail += check_same_reflections(cell, ct); + cell_free(ct); + cell_free(cback); + + return fail; } -static int check_identity(UnitCell *cell, UnitCellTransformation *tfn) +static int check_identity(UnitCell *cell, IntegerMatrix *tfn) { UnitCell *cnew; double a[9], b[9]; int i; int fail = 0; - cnew = cell_transform(cell, tfn); + STATUS("-----------------------\n"); + + cnew = cell_transform_intmat(cell, tfn); + + STATUS("----> Before identity transformation:\n"); + cell_print(cell); + STATUS("----> The identity transformation matrix:\n"); + intmat_print(tfn); + STATUS("----> After identity transformation:\n"); + cell_print(cnew); cell_get_cartesian(cell, &a[0], &a[1], &a[2], &a[3], &a[4], &a[5], @@ -239,14 +301,15 @@ static int check_identity(UnitCell *cell, UnitCellTransformation *tfn) for ( i=0; i<9; i++ ) { if ( !within_tolerance(a[i], b[i], 0.1) ) { fail = 1; - STATUS("%e %e\n", a[i], b[i]); + //STATUS("%e %e\n", a[i], b[i]); } } if ( fail ) { - ERROR("Original cell not recovered after transformation:\n"); + ERROR("********************************************** "); + ERROR("Original cell not recovered after identity transformation:\n"); cell_print(cell); - tfn_print(tfn); + intmat_print(tfn); cell_print(cnew); } @@ -258,7 +321,8 @@ int main(int argc, char *argv[]) { int fail = 0; UnitCell *cell, *cref; - UnitCellTransformation *tfn; + IntegerMatrix *tfn; + IntegerMatrix *part1, *part2; gsl_rng *rng; rng = gsl_rng_alloc(gsl_rng_mt19937); @@ -273,53 +337,52 @@ int main(int argc, char *argv[]) if ( cell == NULL ) return 1; cell_free(cref); + tfn = intmat_identity(3); + /* Permutation of axes */ - tfn = tfn_identity(); if ( tfn == NULL ) return 1; - tfn_combine(tfn, tfn_vector(0,1,0), - tfn_vector(0,0,1), - tfn_vector(1,0,0)); - fail += check_transformation(cell, tfn, 1, NULL); - tfn_free(tfn); + intmat_set_all_3x3(tfn, 0,0,1, + 1,0,0, + 0,1,0); + fail += check_transformation(cell, tfn, 1); /* Doubling of cell in one direction */ - tfn = tfn_identity(); if ( tfn == NULL ) return 1; - tfn_combine(tfn, tfn_vector(2,0,0), - tfn_vector(0,1,0), - tfn_vector(0,0,1)); - fail += check_transformation(cell, tfn, 0, NULL); - tfn_free(tfn); - - /* Diagonal supercell */ - tfn = tfn_identity(); + intmat_set_all_3x3(tfn, 2,0,0, + 0,1,0, + 0,0,1); + fail += check_transformation(cell, tfn, 0); + + /* Shearing */ if ( tfn == NULL ) return 1; - tfn_combine(tfn, tfn_vector(1,1,0), - tfn_vector(0,1,0), - tfn_vector(0,0,1)); - fail += check_transformation(cell, tfn, 1, NULL); - tfn_free(tfn); + intmat_set_all_3x3(tfn, 1,0,0, + 1,1,0, + 0,0,1); + fail += check_transformation(cell, tfn, 1); /* Crazy */ - tfn = tfn_identity(); if ( tfn == NULL ) return 1; - tfn_combine(tfn, tfn_vector(1,1,0), - tfn_vector(0,1,0), - tfn_vector(0,1,1)); - fail += check_transformation(cell, tfn, 0, NULL); - tfn_free(tfn); + intmat_set_all_3x3(tfn, 1,0,0, + 1,1,1, + 0,0,1); + fail += check_transformation(cell, tfn, 0); /* Identity in two parts */ - tfn = tfn_identity(); + part1 = intmat_identity(3); + part2 = intmat_identity(3); if ( tfn == NULL ) return 1; - tfn_combine(tfn, tfn_vector(0,0,1), - tfn_vector(0,1,0), - tfn_vector(-1,0,0)); - tfn_combine(tfn, tfn_vector(0,0,-1), - tfn_vector(0,1,0), - tfn_vector(1,0,0)); + intmat_set_all_3x3(part1, 0,0,-1, + 0,1,0, + 1,0,0); + intmat_set_all_3x3(part2, 0,0,1, + 0,1,0, + -1,0,0); + tfn = intmat_intmat_mult(part1, part2); fail += check_identity(cell, tfn); - tfn_free(tfn); + intmat_free(part1); + intmat_free(part2); + + intmat_free(tfn); /* Check some uncentering transformations */ cref = cell_new_from_parameters(50e-10, 50e-10, 50e-10, |