From 0f4208d86bbf2da6ec04ca3c3867422e4ba90ace Mon Sep 17 00:00:00 2001 From: cppxfel Date: Fri, 30 Jun 2017 01:08:26 +0100 Subject: Removed some code, fixed some code, and so on and so forth --- libcrystfel/src/taketwo.c | 395 +++++++++++++++++----------------------------- 1 file changed, 146 insertions(+), 249 deletions(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/taketwo.c b/libcrystfel/src/taketwo.c index b37fc8f1..feb4ac36 100644 --- a/libcrystfel/src/taketwo.c +++ b/libcrystfel/src/taketwo.c @@ -102,19 +102,17 @@ struct TakeTwoCell #define RECIP_TOLERANCE (0.0010*1e10) /* Threshold for network members to consider a potential solution */ -#define NETWORK_MEMBER_THRESHOLD (25) +#define NETWORK_MEMBER_THRESHOLD (20) /* Maximum network members (obviously a solution so should stop) */ #define MAX_NETWORK_MEMBERS (NETWORK_MEMBER_THRESHOLD + 3) /* Maximum dead ends for a single branch extension during indexing */ -#define MAX_DEAD_ENDS (2) +#define MAX_DEAD_ENDS (10) /* Tolerance for two angles to be considered the same */ #define ANGLE_TOLERANCE (deg2rad(0.6)) -#define COSINE_TOLERANCE 0.010 - /* Tolerance for rot_mats_are_similar */ #define TRACE_TOLERANCE (deg2rad(4.0)) @@ -433,10 +431,8 @@ 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); + gsl_matrix *sub = gsl_matrix_calloc(3, 3); + gsl_matrix *mul = gsl_matrix_calloc(3, 3); for (i = 0; i < cell->numOps; i++) { gsl_matrix *testRot = gsl_matrix_alloc(3, 3); @@ -555,7 +551,6 @@ static int obs_shares_spot_w_array(struct SpotVec *obs_vecs, int test_idx, { int i; int total = 0; - int target = 1; struct SpotVec *her_obs = &obs_vecs[test_idx]; @@ -564,143 +559,106 @@ static int obs_shares_spot_w_array(struct SpotVec *obs_vecs, int test_idx, int shares = obs_vecs_share_spot(her_obs, his_obs); - if ( shares ) return 1;; - } - - if (total > target) { - return 1; + if ( shares ) return 1; } return 0; } -/** Note: this could potentially (and cautiously) converted to comparing - * cosines rather than angles, to lose an "acos" but different parts of the - * cosine graph are more sensitive than others, so may be a trade off... or not. - */ static int obs_vecs_match_angles(struct SpotVec *her_obs, - struct SpotVec *his_obs, - int **her_match_idxs, int **his_match_idxs, - int *match_count) + struct SpotVec *his_obs, + int **her_match_idxs, int **his_match_idxs, + int *match_count) { - int i, j; - *match_count = 0; - - double min_angle = 0.9990; - double max_angle = -0.99144; - - /* calculate angle between observed vectors */ - double obs_cos = rvec_cosine(her_obs->obsvec, his_obs->obsvec); - - /* calculate angle between all potential theoretical vectors */ - - for ( i=0; imatch_num; i++ ) { - for ( j=0; jmatch_num; j++ ) { - - struct rvec *her_match = &her_obs->matches[i]; - struct rvec *his_match = &his_obs->matches[j]; - - double theory_cos = rvec_cosine(*her_match, - *his_match); - - /* is this angle a match? */ - - double cos_diff = fabs(theory_cos - obs_cos); - - if ( cos_diff < COSINE_TOLERANCE ) { - // in the case of a brief check only - if (!her_match_idxs || !his_match_idxs) { - return 1; - } - - /* If the angles are too close to 0 - or 180, one axis ill-determined */ - 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_cos = rvec_cosine(*her_match, - theory_diff); - obs_cos = rvec_cosine(her_obs->obsvec, obs_diff); - - if (fabs(obs_cos - theory_cos) > COSINE_TOLERANCE) { - continue; - } - - theory_cos = rvec_cosine(*his_match, - theory_diff); - obs_cos = rvec_cosine(his_obs->obsvec, obs_diff); - - if (fabs(obs_cos - theory_cos) > COSINE_TOLERANCE) { - continue; - } - - size_t new_size = (*match_count + 1) * - sizeof(int); - if (her_match_idxs && his_match_idxs) - { - /* Reallocate the array to fit in another match */ - int *temp_hers; - int *temp_his; - temp_hers = realloc(*her_match_idxs, new_size); - temp_his = realloc(*his_match_idxs, new_size); - - if ( temp_hers == NULL || temp_his == NULL ) { - apologise(); - } - - (*her_match_idxs) = temp_hers; - (*his_match_idxs) = temp_his; - - (*her_match_idxs)[*match_count] = i; - (*his_match_idxs)[*match_count] = j; - } - - (*match_count)++; - } - } - } - - return (*match_count > 0); + int i, j; + *match_count = 0; + + double min_angle = deg2rad(2.5); + double max_angle = deg2rad(187.5); + + /* calculate angle between observed vectors */ + double obs_angle = rvec_angle(her_obs->obsvec, his_obs->obsvec); + + /* calculate angle between all potential theoretical vectors */ + + for ( i=0; imatch_num; i++ ) { + for ( j=0; jmatch_num; j++ ) { + + 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); + + /* is this angle a match? */ + + double angle_diff = fabs(theory_angle - obs_angle); + + if ( angle_diff < ANGLE_TOLERANCE ) { + // in the case of a brief check only + if (!her_match_idxs || !his_match_idxs) { + return 1; + } + + /* If the angles are too close to 0 + or 180, one axis ill-determined */ + if (theory_angle < min_angle || + theory_angle > 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_diff); + obs_angle = rvec_angle(her_obs->obsvec, obs_diff); + + if (fabs(obs_angle - theory_angle) > ANGLE_TOLERANCE) { + continue; + } + + theory_angle = rvec_angle(*his_match, + theory_diff); + obs_angle = rvec_angle(his_obs->obsvec, obs_diff); + + if (fabs(obs_angle - theory_angle) > ANGLE_TOLERANCE) { + continue; + } + + size_t new_size = (*match_count + 1) * + sizeof(int); + if (her_match_idxs && his_match_idxs) + { + /* Reallocate the array to fit in another match */ + int *temp_hers; + int *temp_his; + temp_hers = realloc(*her_match_idxs, new_size); + temp_his = realloc(*his_match_idxs, new_size); + + if ( temp_hers == NULL || temp_his == NULL ) { + apologise(); + } + + (*her_match_idxs) = temp_hers; + (*his_match_idxs) = temp_his; + + (*her_match_idxs)[*match_count] = i; + (*his_match_idxs)[*match_count] = j; + } + + (*match_count)++; + } + } + } + + return (*match_count > 0); } -// DELETE ME IF UNUSED -static int obs_angles_match_array(struct SpotVec *obs_vecs, int test_idx, - int *obs_members, int *match_members, int num) -{ - /* note: this is just a preliminary check to reduce unnecessary - * computation later down the line, but is not entirely accurate. - **/ - - int i = 0; - struct SpotVec *her_obs = &obs_vecs[test_idx]; - - for ( i=0; imatches[match_members[j]]; + struct rvec me_obs = me->obsvec; struct rvec you_obs = you->obsvec; - int *your_idxs = 0; - int *my_idxs = 0; + int one_is_okay = 0; - obs_vecs_match_angles(me, you, &my_idxs, - &your_idxs, &match_num); - - STATUS("2\n"); - - if (match_num == 0) { - all_ok = 0; - STATUS("shit\n"); - continue; - } - - weed_duplicate_matches(me, you, - &my_idxs, &your_idxs, - &match_num, cell); - - STATUS("3\n"); - - for ( k=0; kmatch_num; k++ ) { + gsl_matrix *test_rot; + + struct rvec me_cell = me->matches[k]; - if (my_idxs[k] < 0 || your_idxs[k] < 0) { - continue; - } - - STATUS("4\n"); - - struct rvec me_cell = me->matches[my_idxs[k]]; - struct rvec you_cell; - - you_cell = you->matches[your_idxs[k]]; - test_rot = generate_rot_mat(me_obs, you_obs, me_cell, you_cell, twiz1, twiz2); @@ -955,23 +858,34 @@ static signed int find_next_index(gsl_matrix *rot, struct SpotVec *obs_vecs, gsl_matrix_free(test_rot); - STATUS("5\n"); - - if (!ok) { - all_ok = 0; - STATUS("fuck\n"); - - break; + if (ok) { + one_is_okay = 1; + + if (matched >= 0 && k == matched) { + *match_found = k; + } else if (matched < 0) { + matched = k; + } else { + one_is_okay = 0; + break; + } } } - free(my_idxs); - free(your_idxs); + if (!one_is_okay) { + all_ok = 0; + break; + } } - + + if (all_ok) { - STATUS("phew\n"); - *match_found = i; + + for ( int j=0; j *max_members) { *max_members = member_num; } @@ -1059,19 +970,6 @@ static int grow_network(gsl_matrix *rot, struct SpotVec *obs_vecs, /* If member_num is high enough, we want to return a yes */ if ( member_num > NETWORK_MEMBER_THRESHOLD ) break; } - - /* - int i; - for (i = 0; i < member_num; i++) - { - STATUS("."); - } - - if ( member_num > NETWORK_MEMBER_THRESHOLD ) { - STATUS(" yes (%i)\n", member_num); - } else { - STATUS(" nope (%i)\n", member_num); - }*/ finish_solution(rot, obs_vecs, obs_members, match_members, member_num); @@ -1620,15 +1518,11 @@ static UnitCell *run_taketwo(UnitCell *cell, struct rvec *rlps, int rlp_count) return NULL; } -// STATUS("Returning result.\n"); - result = transform_cell_gsl(cell, solution); gsl_matrix_free(solution); cleanup_taketwo_obs_vecs(obs_vecs, obs_vec_count); cleanup_taketwo_cell(&ttCell); - STATUS("good.\n"); - return result; } @@ -1673,8 +1567,11 @@ int taketwo_index(struct image *image, IndexingPrivate *ipriv) if ( !peak_sanity_check(image, &cr, 1) ) { cell_free(cell); crystal_free(cr); + // STATUS("Rubbish!!\n"); return 0; + } else { + // STATUS("That's good!\n"); } } @@ -1711,7 +1608,7 @@ IndexingPrivate *taketwo_prepare(IndexingMethod *indm, UnitCell *cell, STATUS("***** Welcome to TakeTwo *****\n"); STATUS("************************************\n\n"); STATUS("If you use these indexing results, please keep a roof\n"); - STATUS("over my head by citing:\n"); + STATUS("over the author's head by citing:\n"); STATUS("Ginn et al., Acta Cryst. (2016). D72, 956-965\n"); STATUS("\n"); -- cgit v1.2.3