From 4d8c846f2853ffa85c170c3e0aca5506e9c6a9a2 Mon Sep 17 00:00:00 2001 From: Helen Ginn Date: Wed, 18 Jan 2017 12:59:39 +0000 Subject: Fiddled with parameters, now nice for CPV17. Fixed a ton of bugs too --- libcrystfel/src/taketwo.c | 253 ++++++++++++++++++++++++++++++---------------- 1 file changed, 168 insertions(+), 85 deletions(-) (limited to 'libcrystfel/src/taketwo.c') diff --git a/libcrystfel/src/taketwo.c b/libcrystfel/src/taketwo.c index 028847ac..7e545fd2 100644 --- a/libcrystfel/src/taketwo.c +++ b/libcrystfel/src/taketwo.c @@ -30,6 +30,7 @@ #include #include +#include #include #include @@ -37,6 +38,7 @@ #include "index.h" #include "taketwo.h" #include "peaks.h" +#include /** * spotvec @@ -74,25 +76,25 @@ int global_nrlps; /* Maximum distance between two rlp sizes to consider info for indexing */ -#define MAX_RECIP_DISTANCE (0.15*1e10) +#define MAX_RECIP_DISTANCE (0.12*1e10) /* Tolerance for two lengths in reciprocal space to be considered the same */ -#define RECIP_TOLERANCE (0.0002*1e10) +#define RECIP_TOLERANCE (0.0005*1e10) /* Threshold for network members to consider a potential solution */ -#define NETWORK_MEMBER_THRESHOLD (50) +#define NETWORK_MEMBER_THRESHOLD (20) /* Maximum network members (obviously a solution so should stop) */ -#define MAX_NETWORK_MEMBERS (500) +#define MAX_NETWORK_MEMBERS (NETWORK_MEMBER_THRESHOLD + 5) /* Maximum dead ends for a single branch extension during indexing */ -#define MAX_DEAD_ENDS (10) +#define MAX_DEAD_ENDS (5) /* Tolerance for two angles to be considered the same */ #define ANGLE_TOLERANCE (deg2rad(1.0)) /* Tolerance for rot_mats_are_similar */ -#define TRACE_TOLERANCE (deg2rad(8.0)) +#define TRACE_TOLERANCE (deg2rad(3.0)) /** TODO: * @@ -192,8 +194,10 @@ static struct rvec rvec_cross(struct rvec a, struct rvec b) } -static void show_rvec(struct rvec r) +static void show_rvec(struct rvec r2) { + struct rvec r = r2; + normalise_rvec(&r); STATUS("[ %.3f %.3f %.3f ]\n", r.u, r.v, r.w); } @@ -295,7 +299,8 @@ static double matrix_trace(gsl_matrix *a) } -static int rot_mats_are_similar(gsl_matrix *rot1, gsl_matrix *rot2) +static int rot_mats_are_similar(gsl_matrix *rot1, gsl_matrix *rot2, + double *score) { double tr; gsl_matrix *sub; @@ -309,6 +314,8 @@ static int rot_mats_are_similar(gsl_matrix *rot1, gsl_matrix *rot2) 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); double max = sqrt(4.0*(1.0-cos(TRACE_TOLERANCE))); @@ -464,24 +471,24 @@ static int obs_vecs_match_angles(struct SpotVec *her_obs, temp_hers = realloc(*her_match_idxs, new_size); int *temp_his; 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)++; + (*match_count)++; } } } - return (match_count > 0); + return (*match_count > 0); } @@ -494,7 +501,7 @@ static int obs_angles_match_array(struct SpotVec *obs_vecs, int test_idx, * is identical - too much faff. **/ - int i; + int i, j; struct SpotVec *her_obs = &obs_vecs[test_idx]; for ( i=0; i obs_vec_count) { - return -1; + return 0; } int match_found = -1; @@ -625,7 +702,7 @@ static int grow_network(gsl_matrix *rot, struct SpotVec *obs_vecs, obs_vec_count, obs_members, match_members, start, member_num, - &match_found); + &match_found, 0); if ( member_num < 2 ) { return 0; @@ -636,8 +713,7 @@ static int grow_network(gsl_matrix *rot, struct SpotVec *obs_vecs, * on indexing altogether. **/ if ( dead_ends > MAX_DEAD_ENDS ) { - dead_ends = 0; - continue; + break; } /* We have not had too many dead ends. Try removing @@ -650,7 +726,7 @@ static int grow_network(gsl_matrix *rot, struct SpotVec *obs_vecs, } /* we have elongated membership - so reset dead_ends counter */ - // dead_ends = 0; + dead_ends = 0; obs_members[member_num] = next_index; match_members[member_num] = match_found; @@ -666,8 +742,8 @@ static int grow_network(gsl_matrix *rot, struct SpotVec *obs_vecs, if ( member_num > NETWORK_MEMBER_THRESHOLD ) break; } - /* Deal with this shit after coffee */ - /* (note: turns out there's no shit to deal with) */ + finalise_solution(rot, obs_vecs, obs_members, + match_members, member_num); return ( member_num > NETWORK_MEMBER_THRESHOLD ); } @@ -684,11 +760,10 @@ static signed int spot_idx(struct rvec *rlp) static int start_seed(struct SpotVec *obs_vecs, int obs_vec_count, int i, int j, int i_match, int j_match, gsl_matrix **rotation, - int *max_members, int accept_anything) + int *max_members) { gsl_matrix *rot_mat = gsl_matrix_calloc(3, 3); - // FIXME: go through ALL matches, not jsut first rot_mat = generate_rot_mat(obs_vecs[i].obsvec, obs_vecs[j].obsvec, obs_vecs[i].matches[i_match], @@ -697,17 +772,12 @@ static int start_seed(struct SpotVec *obs_vecs, int obs_vec_count, int i, /* Try to expand this rotation matrix to a larger network */ int success = grow_network(rot_mat, obs_vecs, obs_vec_count, - i, j, i_match, j_match, max_members, - accept_anything); - - /* return this matrix or free it and try again */ - if ( success ) { - *rotation = rot_mat; - return 1; - } else { - gsl_matrix_free(rot_mat); - return 0; - } + i, j, i_match, j_match, max_members); + + /* return this matrix and if it was immediately successful */ + *rotation = rot_mat; + + return success; } static int find_seed(struct SpotVec *obs_vecs, int obs_vec_count, @@ -715,7 +785,9 @@ static int find_seed(struct SpotVec *obs_vecs, int obs_vec_count, { /* META: Maximum achieved maximum network membership */ int max_max_members = 0; - int best_i, best_j, best_i_idx, best_j_idx; + gsl_matrix *best_rotation = NULL; + + unsigned long start_time = time(NULL); /* loop round pairs of vectors to try and find a suitable * seed to start building a self-consistent network of vectors @@ -746,42 +818,43 @@ static int find_seed(struct SpotVec *obs_vecs, int obs_vec_count, /* We have seeds! Pass each of them through the seed-starter */ /* If a seed has the highest achieved membership, make note...*/ int k; - for ( k=0; k<1; k++ ) { + for ( k=0; k max_max_members) { max_max_members = max_members; - best_i = i; - best_j = j; - best_i_idx = i_idx[k]; - best_j_idx = j_idx[k]; + gsl_matrix_free(best_rotation); + best_rotation = *rotation; + } else { + gsl_matrix_free(*rotation); + *rotation = NULL; } } + + unsigned long now_time = time(NULL); + unsigned int seconds = now_time - start_time; + + if (seconds > 60) { + /* Heading towards CrystFEL cutoff so + return your best guess and run */ + *rotation = best_rotation; + STATUS("After %i seconds, returning best guess\n", seconds); + return (best_rotation != NULL); + } } } } /* yes this } is meant to be here */ - int max_members = 0; - int success = start_seed(obs_vecs, obs_vec_count, best_i, best_j, - best_i_idx, best_j_idx, - rotation, &max_members, - 1); - - if (success) { - return success; - } else { - /* Well, shit... something went wrong. */ - return 0; - } + *rotation = best_rotation; + return (best_rotation != NULL); } @@ -1017,8 +1090,11 @@ global_nrlps = rlp_count; cleanup_taketwo_cell_vecs(cell_vecs); - find_seed(obs_vecs, obs_vec_count, &solution); - if ( solution == NULL ) return NULL; + int threshold_reached = find_seed(obs_vecs, obs_vec_count, &solution); + + if ( solution == NULL ) { + return NULL; + } result = transform_cell_gsl(cell, solution); @@ -1052,6 +1128,8 @@ int taketwo_index(struct image *image, IndexingPrivate *ipriv) rlps[n_rlps].v = 0.0; rlps[n_rlps++].w = 0.0; +// STATUS("n_rlps = %i", n_rlps); + cell = run_taketwo(tp->cell, rlps, n_rlps); if ( cell == NULL ) return 0; @@ -1099,6 +1177,11 @@ IndexingPrivate *taketwo_prepare(IndexingMethod *indm, UnitCell *cell, ERROR("TakeTwo indexing requires a unit cell.\n"); return NULL; } + + STATUS("Welcome to TakeTwo\n"); + STATUS("If you use these indexing results, please cite:\n"); + STATUS("Ginn et al., Acta Cryst. (2016). D72, 956-965\n"); + tp = malloc(sizeof(struct taketwo_private)); if ( tp == NULL ) return NULL; -- cgit v1.2.3