From b1431a52affcbe1bba4121441fb4194146668641 Mon Sep 17 00:00:00 2001 From: cppxfel Date: Mon, 19 Jun 2017 17:20:07 +0200 Subject: Theoretical speed improvements by replacing acos with comparison of cosines for the quick check, but not necessarily good --- libcrystfel/src/taketwo.c | 84 +++++++++++++++++++++++++++++++---------------- 1 file changed, 56 insertions(+), 28 deletions(-) (limited to 'libcrystfel/src') diff --git a/libcrystfel/src/taketwo.c b/libcrystfel/src/taketwo.c index ff7a9f41..a3c9efb0 100644 --- a/libcrystfel/src/taketwo.c +++ b/libcrystfel/src/taketwo.c @@ -112,6 +112,7 @@ struct TakeTwoCell /* Tolerance for two angles to be considered the same */ #define ANGLE_TOLERANCE (deg2rad(0.5)) +#define COSINE_TOLERANCE 0.017 /* Tolerance for rot_mats_are_similar */ #define TRACE_TOLERANCE (deg2rad(4.0)) @@ -409,24 +410,19 @@ static SymOpList *sym_ops_for_cell(UnitCell *cell) } static int rot_mats_are_similar(gsl_matrix *rot1, gsl_matrix *rot2, + gsl_matrix *sub, gsl_matrix *mul, double *score) { double tr; - gsl_matrix *sub; - gsl_matrix *mul; - sub = gsl_matrix_calloc(3, 3); gsl_matrix_memcpy(sub, rot1); gsl_matrix_sub(sub, rot2); /* sub = rot1 - rot2 */ - mul = gsl_matrix_calloc(3, 3); gsl_blas_dgemm(CblasNoTrans, CblasTrans, 1.0, sub, sub, 0.0, mul); tr = matrix_trace(mul); if (score != NULL) *score = tr; - gsl_matrix_free(mul); - gsl_matrix_free(sub); double max = sqrt(4.0*(1.0-cos(TRACE_TOLERANCE))); return (tr < max); @@ -437,6 +433,11 @@ static int symm_rot_mats_are_similar(gsl_matrix *rot1, gsl_matrix *rot2, { int i; + gsl_matrix *sub; + gsl_matrix *mul; + sub = gsl_matrix_calloc(3, 3); + mul = gsl_matrix_calloc(3, 3); + for (i = 0; i < cell->numOps; i++) { gsl_matrix *testRot = gsl_matrix_alloc(3, 3); gsl_matrix *symOp = cell->rotSymOps[i]; @@ -444,14 +445,19 @@ static int symm_rot_mats_are_similar(gsl_matrix *rot1, gsl_matrix *rot2, gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, rot1, symOp, 0.0, testRot); - if (rot_mats_are_similar(testRot, rot2, NULL)) { - gsl_matrix_free(testRot); + if (rot_mats_are_similar(testRot, rot2, sub, mul, NULL)) { + gsl_matrix_free(testRot); + gsl_matrix_free(sub); + gsl_matrix_free(mul); return 1; } gsl_matrix_free(testRot); } + gsl_matrix_free(sub); + gsl_matrix_free(mul); + return 0; } @@ -583,11 +589,11 @@ static int obs_vecs_match_angles(struct SpotVec *her_obs, int i, j; *match_count = 0; - double min_angle = deg2rad(2.5); - double max_angle = deg2rad(187.5); + double min_angle = 0.9990; + double max_angle = -0.99144; /* calculate angle between observed vectors */ - double obs_angle = rvec_angle(her_obs->obsvec, his_obs->obsvec); + double obs_cos = rvec_cosine(her_obs->obsvec, his_obs->obsvec); /* calculate angle between all potential theoretical vectors */ @@ -597,14 +603,14 @@ static int obs_vecs_match_angles(struct SpotVec *her_obs, struct rvec *her_match = &her_obs->matches[i]; struct rvec *his_match = &his_obs->matches[j]; - double theory_angle = rvec_angle(*her_match, - *his_match); + double theory_cos = rvec_cosine(*her_match, + *his_match); /* is this angle a match? */ - double angle_diff = fabs(theory_angle - obs_angle); + double cos_diff = fabs(theory_cos - obs_cos); - if ( angle_diff < ANGLE_TOLERANCE ) { + if ( cos_diff < COSINE_TOLERANCE ) { // in the case of a brief check only if (!her_match_idxs || !his_match_idxs) { return 1; @@ -612,32 +618,32 @@ static int obs_vecs_match_angles(struct SpotVec *her_obs, /* If the angles are too close to 0 or 180, one axis ill-determined */ - if (theory_angle < min_angle || - theory_angle > max_angle) { + if (theory_cos > min_angle || + theory_cos < max_angle) { continue; } // check the third vector - + /* struct rvec theory_diff = diff_vec(*his_match, *her_match); struct rvec obs_diff = diff_vec(his_obs->obsvec, her_obs->obsvec); - theory_angle = rvec_angle(*her_match, + theory_cos = rvec_cosine(*her_match, theory_diff); - obs_angle = rvec_angle(her_obs->obsvec, obs_diff); + obs_cos = rvec_cosine(her_obs->obsvec, obs_diff); - if (fabs(obs_angle - theory_angle) > ANGLE_TOLERANCE) { + if (fabs(obs_cos - theory_cos) > COSINE_TOLERANCE) { continue; } - theory_angle = rvec_angle(*his_match, + theory_cos = rvec_cosine(*his_match, theory_diff); - obs_angle = rvec_angle(his_obs->obsvec, obs_diff); + obs_cos = rvec_cosine(his_obs->obsvec, obs_diff); - if (fabs(obs_angle - theory_angle) > ANGLE_TOLERANCE) { + if (fabs(obs_cos - theory_cos) > COSINE_TOLERANCE) { continue; - } + }*/ size_t new_size = (*match_count + 1) * sizeof(int); @@ -705,6 +711,11 @@ static signed int finish_solution(gsl_matrix *rot, struct SpotVec *obs_vecs, int *obs_members, int *match_members, int member_num) { + gsl_matrix *sub; + gsl_matrix *mul; + sub = gsl_matrix_calloc(3, 3); + mul = gsl_matrix_calloc(3, 3); + gsl_matrix **rotations = malloc(sizeof(*rotations)* pow(member_num, 2) - member_num); int i, j, count; @@ -736,6 +747,7 @@ static signed int finish_solution(gsl_matrix *rot, struct SpotVec *obs_vecs, for (j=0; j