From c24bf8f1f75c18beec6a90df4e7fa620b5ce8d82 Mon Sep 17 00:00:00 2001 From: Helen Ginn Date: Sat, 21 Apr 2018 12:47:56 +0100 Subject: Fiddling with comments in text --- libcrystfel/src/taketwo.c | 88 ++++++++++++++++++++++++++++++++++++++++++----- 1 file changed, 79 insertions(+), 9 deletions(-) diff --git a/libcrystfel/src/taketwo.c b/libcrystfel/src/taketwo.c index 0857f813..b2c3fdef 100644 --- a/libcrystfel/src/taketwo.c +++ b/libcrystfel/src/taketwo.c @@ -28,6 +28,67 @@ * */ +/** + * \class TakeTwo + * Code outline. + * --- Get ready for calculation --- + * Pre-calculate symmertry operations (generate_rotation_symops()) + * Pre-calculate theoretical vectors from unit cell dimensions + * (gen_theoretical_vecs()) + * Generate observed vectors from data (gen_observed_vecs()) + * Match observed vectors to theoretical vectors (match_obs_to_cell_vecs()) + * + * --- Business bit --- + * Find starting seeds (find_seed()): + * - Loop through pairs of observed vectors + * - If they share a spot, find matching pairs of theoretical vectors + * - Remove all duplicate matches due to symmetry operations + * - For the remainder, loop through the matches and extend the seeds + * (start_seed()). + * - If it returns a membership greater than the highest member threshold, + * return the matrix to CrystFEL. + * + * Extending a seed (start_seed()): + * - Generate a rotation matrix which matches the chosen start seed. + * - Loop through all observed vectors starting from 0. + * - Find another vector to add to the network, if available + * (find_next_index()). + * - If the index is not available, then give up if there were too many dead + * ends. Otherwise, remove the last member and pretend like that didn't + * happen, so it won't happen again. + * - Add the vector to increment the membership list. + * - If the membership number exceeds the maximum, tidy up the solution and + * return a success. + * - If the membership does not, then resume the loop and search for the + * next vector. + * + * Finding the next member (find_next_index()): + * - Go through the observed vectors, starting from the last index + 1 to + * explore only the "new" vectors. + * - If the vector does not share a spot with the current array of vectors, + * then skip it. + * - We must loop through all the current vectors in the network, and see if + * they match the newcomer for a given matching theoretical vector. + * - We only accept a match if the rotation matrix matches the seed matrix + * for a single matching theoretical vector. + * - If it does match, we can return a success. + * + * Tidying the solution (finish_solution()): + * - This chooses the most common rotation matrix of the bunch to choose to + * send to CrystFEL. But this should probably take the average instead, + * which is very possible. + * + * Clean up the mess (cleanup_taketwo_obs_vecs()) + */ + +/** + * Helen's to-do list + * - + * + * + * - Improve the final solution + */ + #include #include #include @@ -73,24 +134,30 @@ struct taketwo_private UnitCell *cell; }; -// These rotation symmetry operators +/** + * Internal structure which gets passed the various functions looking after + * the indexing bits and bobs. */ struct TakeTwoCell { - UnitCell *cell; + UnitCell *cell; /**< Contains unit cell dimensions */ gsl_matrix **rotSymOps; unsigned int numOps; - struct SpotVec **obs_vecs; // Pointer to an array + struct SpotVec **obs_vecs; int obs_vec_count; /* Options */ int member_thresh; - double len_tol; /* In reciprocal metres */ - double angle_tol; /* In radians */ - double trace_tol; /* Contains sqrt(4*(1-cos(angle))) */ + double len_tol; /**< In reciprocal metres */ + double angle_tol; /**< In radians */ + double trace_tol; /**< Contains sqrt(4*(1-cos(angle))) */ + /**< Temporary memory always allocated for calculations */ gsl_matrix *twiz1Tmp; + /**< Temporary memory always allocated for calculations */ gsl_matrix *twiz2Tmp; + /**< Temporary memory always allocated for calculations */ gsl_vector *vec1Tmp; + /**< Temporary memory always allocated for calculations */ gsl_vector *vec2Tmp; }; @@ -844,6 +911,11 @@ static signed int find_next_index(gsl_matrix *rot, int *obs_members, if (ok) { one_is_okay = 1; + /* We are only happy if the vector + * matches for only one kind of + * theoretical vector. We don't want to + * accept mixtures of theoretical vector + * matches. */ if (matched >= 0 && k == matched) { *match_found = k; } else if (matched < 0) { @@ -1314,9 +1386,7 @@ static int gen_observed_vecs(struct rvec *rlps, int rlp_count, } } - /* Sort such that the shortest and least error-prone distances - are searched first. - */ + /* Sort such that the shortest distances are searched first. */ qsort(*obs_vecs, count, sizeof(struct SpotVec), compare_spot_vecs); *obs_vec_count = count; -- cgit v1.2.3 From 32d4f204e33f0ad1d26545654d893928c5e0a142 Mon Sep 17 00:00:00 2001 From: Helen Ginn Date: Sat, 21 Apr 2018 12:49:55 +0100 Subject: Theoretical vectors have their own structure including whether they are in the asymmetric unit --- libcrystfel/src/taketwo.c | 111 +++++++++++++++++++--------------------------- 1 file changed, 46 insertions(+), 65 deletions(-) diff --git a/libcrystfel/src/taketwo.c b/libcrystfel/src/taketwo.c index b2c3fdef..2776a352 100644 --- a/libcrystfel/src/taketwo.c +++ b/libcrystfel/src/taketwo.c @@ -118,15 +118,22 @@ struct SpotVec { struct rvec obsvec; - struct rvec *matches; + struct TheoryVec *matches; int match_num; - struct rvec *asym_matches; - int asym_match_num; double distance; struct rvec *her_rlp; struct rvec *his_rlp; }; +/** + * theoryvec + */ + +struct TheoryVec +{ + struct rvec vec; + int asym; +}; struct taketwo_private { @@ -639,8 +646,8 @@ static int obs_vecs_match_angles(struct SpotVec *her_obs, 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]; + struct rvec *her_match = &her_obs->matches[i].vec; + struct rvec *his_match = &his_obs->matches[j].vec; double theory_angle = rvec_angle(*her_match, *his_match); @@ -739,8 +746,8 @@ static signed int finish_solution(gsl_matrix *rot, struct SpotVec *obs_vecs, struct rvec i_obsvec = i_vec.obsvec; struct rvec j_obsvec = j_vec.obsvec; - struct rvec i_cellvec = i_vec.matches[match_members[i]]; - struct rvec j_cellvec = j_vec.matches[match_members[j]]; + struct rvec i_cellvec = i_vec.matches[match_members[i]].vec; + struct rvec j_cellvec = j_vec.matches[match_members[j]].vec; rotations[count] = generate_rot_mat(i_obsvec, j_obsvec, i_cellvec, j_cellvec, @@ -806,8 +813,8 @@ static int weed_duplicate_matches(struct SpotVec *her_obs, struct rvec i_obsvec = her_obs->obsvec; struct rvec j_obsvec = his_obs->obsvec; - struct rvec i_cellvec = her_obs->matches[her_match]; - struct rvec j_cellvec = his_obs->matches[his_match]; + struct rvec i_cellvec = her_obs->matches[her_match].vec; + struct rvec j_cellvec = his_obs->matches[his_match].vec; gsl_matrix *mat = generate_rot_mat(i_obsvec, j_obsvec, i_cellvec, j_cellvec, cell); @@ -885,18 +892,21 @@ static signed int find_next_index(gsl_matrix *rot, int *obs_members, struct SpotVec *me = &obs_vecs[i]; struct SpotVec *you = &obs_vecs[obs_members[j]]; - struct rvec you_cell = you->matches[match_members[j]]; + struct rvec you_cell = you->matches[match_members[j]].vec; struct rvec me_obs = me->obsvec; struct rvec you_obs = you->obsvec; int one_is_okay = 0; + /* Loop through all possible theoretical vector + * matches for the newcomer.. */ + for ( k=0; kmatch_num; k++ ) { gsl_matrix *test_rot; - struct rvec me_cell = me->matches[k]; + struct rvec me_cell = me->matches[k].vec; test_rot = generate_rot_mat(me_obs, you_obs, me_cell, you_cell, @@ -1058,8 +1068,8 @@ static int start_seed(int i, int j, int i_match, int j_match, rot_mat = generate_rot_mat(obs_vecs[i].obsvec, obs_vecs[j].obsvec, - obs_vecs[i].matches[i_match], - obs_vecs[j].matches[j_match], + obs_vecs[i].matches[i_match].vec, + obs_vecs[j].matches[j_match].vec, cell); /* Try to expand this rotation matrix to a larger network */ @@ -1256,10 +1266,9 @@ static int generate_rotation_sym_ops(struct TakeTwoCell *ttCell) return 1; } - struct sortme { - struct rvec v; + struct TheoryVec v; double dist; }; @@ -1271,9 +1280,9 @@ static int sort_func(const void *av, const void *bv) } -static int match_obs_to_cell_vecs(struct rvec *cell_vecs, int cell_vec_count, +static int match_obs_to_cell_vecs(struct TheoryVec *cell_vecs, int cell_vec_count, struct SpotVec *obs_vecs, int obs_vec_count, - int is_asymmetric, struct TakeTwoCell *cell) + struct TakeTwoCell *cell) { int i, j; @@ -1284,7 +1293,7 @@ static int match_obs_to_cell_vecs(struct rvec *cell_vecs, int cell_vec_count, for ( j=0; jtrace_tol))); } - success = match_obs_to_cell_vecs(asym_vecs, asym_vec_count, - obs_vecs, obs_vec_count, 1, &ttCell); - - success = match_obs_to_cell_vecs(cell_vecs, cell_vec_count, - obs_vecs, obs_vec_count, 0, &ttCell); + success = match_obs_to_cell_vecs(theory_vecs, cell_vec_count, + obs_vecs, obs_vec_count, &ttCell); - free(cell_vecs); - free(asym_vecs); + free(theory_vecs); if ( !success ) return NULL; -- cgit v1.2.3 From 2bfe47fef2e1d36575d38bef7da1fc29afecd354 Mon Sep 17 00:00:00 2001 From: Helen Ginn Date: Sat, 21 Apr 2018 20:21:39 +0100 Subject: Beginning of support for ignoring previous solutions --- libcrystfel/src/taketwo.c | 79 ++++++++++++++++++++++++++++++++++++++++++----- 1 file changed, 72 insertions(+), 7 deletions(-) diff --git a/libcrystfel/src/taketwo.c b/libcrystfel/src/taketwo.c index 2776a352..ea449b8b 100644 --- a/libcrystfel/src/taketwo.c +++ b/libcrystfel/src/taketwo.c @@ -149,6 +149,8 @@ struct TakeTwoCell UnitCell *cell; /**< Contains unit cell dimensions */ gsl_matrix **rotSymOps; unsigned int numOps; + gsl_matrix **prevSols; /**< Previous solutions to be ignored */ + unsigned int numPrevs; /**< Previous solution count */ struct SpotVec **obs_vecs; int obs_vec_count; @@ -790,6 +792,21 @@ static signed int finish_solution(gsl_matrix *rot, struct SpotVec *obs_vecs, return 1; } +gsl_matrix *rot_mat_from_indices(struct SpotVec *her_obs, struct SpotVec *his_obs, + int her_match, int his_match, + struct TakeTwoCell *cell) +{ + struct rvec i_obsvec = her_obs->obsvec; + struct rvec j_obsvec = his_obs->obsvec; + struct rvec i_cellvec = her_obs->matches[her_match].vec; + struct rvec j_cellvec = his_obs->matches[his_match].vec; + + gsl_matrix *mat = generate_rot_mat(i_obsvec, j_obsvec, + i_cellvec, j_cellvec, cell); + + return mat; +} + static int weed_duplicate_matches(struct SpotVec *her_obs, struct SpotVec *his_obs, int **her_match_idxs, int **his_match_idxs, @@ -807,17 +824,47 @@ static int weed_duplicate_matches(struct SpotVec *her_obs, signed int i, j; int duplicates = 0; + /* First we weed out duplicates with previous solutions */ + for (i = *match_count - 1; i >= 0; i--) { int her_match = (*her_match_idxs)[i]; int his_match = (*his_match_idxs)[i]; - struct rvec i_obsvec = her_obs->obsvec; - struct rvec j_obsvec = his_obs->obsvec; - struct rvec i_cellvec = her_obs->matches[her_match].vec; - struct rvec j_cellvec = his_obs->matches[his_match].vec; + gsl_matrix *mat; + mat = rot_mat_from_indices(her_obs, his_obs, + her_match, his_match, cell); + + if (mat == NULL) + { + continue; + } + + for (j = 0; j < cell->numPrevs; j++) + { + int sim = symm_rot_mats_are_similar(cell->prevSols[j], + mat, cell); - gsl_matrix *mat = generate_rot_mat(i_obsvec, j_obsvec, - i_cellvec, j_cellvec, cell); + /* Found a duplicate with a previous solution */ + if (sim) + { + (*her_match_idxs)[i] = -1; + (*his_match_idxs)[i] = -1; + break; + } + } + + gsl_matrix_free(mat); + } + + /* Now we weed out the self-duplicates from the remaining batch */ + + for (i = *match_count - 1; i >= 0; i--) { + int her_match = (*her_match_idxs)[i]; + int his_match = (*his_match_idxs)[i]; + + gsl_matrix *mat; + mat = rot_mat_from_indices(her_obs, his_obs, + her_match, his_match, cell); int found = 0; @@ -1171,6 +1218,22 @@ static int find_seed(gsl_matrix **rotation, struct TakeTwoCell *cell) } /* yes this } is meant to be here */ *rotation = best_rotation; + + if (best_rotation != NULL) + { + int num = cell->numPrevs + 1; + gsl_matrix **tmp = malloc(sizeof(gsl_matrix *) * num); + + if (tmp) { + cell->prevSols = tmp; + } else { + return 0; + } + + cell->prevSols[num - 1] = best_rotation; + cell->numPrevs = num; + } + return (best_rotation != NULL); } @@ -1530,11 +1593,13 @@ static UnitCell *run_taketwo(UnitCell *cell, const struct taketwo_options *opts, struct TakeTwoCell ttCell; ttCell.cell = cell; ttCell.rotSymOps = NULL; + ttCell.prevSols = NULL; ttCell.twiz1Tmp = gsl_matrix_calloc(3, 3); ttCell.twiz2Tmp = gsl_matrix_calloc(3, 3); ttCell.vec1Tmp = gsl_vector_calloc(3); ttCell.vec2Tmp = gsl_vector_calloc(3); ttCell.numOps = 0; + ttCell.numPrevs = 0; success = generate_rotation_sym_ops(&ttCell); @@ -1586,8 +1651,8 @@ static UnitCell *run_taketwo(UnitCell *cell, const struct taketwo_options *opts, } result = transform_cell_gsl(cell, solution); - gsl_matrix_free(solution); cleanup_taketwo_obs_vecs(obs_vecs, obs_vec_count); + gsl_matrix_free(solution); cleanup_taketwo_cell(&ttCell); return result; -- cgit v1.2.3 From 2153f3e863fa042a3dba68a5b2c612ed16b004ba Mon Sep 17 00:00:00 2001 From: Helen Ginn Date: Sat, 21 Apr 2018 20:24:28 +0100 Subject: Create a new Seed struct instead of passing around tiny little things to various functions --- libcrystfel/src/taketwo.c | 160 +++++++++++++++++++++++++--------------------- 1 file changed, 87 insertions(+), 73 deletions(-) diff --git a/libcrystfel/src/taketwo.c b/libcrystfel/src/taketwo.c index ea449b8b..f92506cc 100644 --- a/libcrystfel/src/taketwo.c +++ b/libcrystfel/src/taketwo.c @@ -135,6 +135,20 @@ struct TheoryVec int asym; }; + +/** + * seed + */ + +struct Seed +{ + struct SpotVec *obs1; + struct SpotVec *obs2; + int idx1; + int idx2; + double score; +}; + struct taketwo_private { IndexingMethod indm; @@ -631,89 +645,97 @@ static int obs_shares_spot_w_array(struct SpotVec *obs_vecs, int test_idx, 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 TakeTwoCell *cell) + struct Seed *seeds, int *match_count, + struct TakeTwoCell *cell) { - int i, j; - *match_count = 0; + int i, j; + *match_count = 0; - double min_angle = deg2rad(2.5); - double max_angle = deg2rad(187.5); + 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 observed vectors */ + double obs_angle = rvec_angle(her_obs->obsvec, his_obs->obsvec); - /* calculate angle between all potential theoretical vectors */ + /* calculate angle between all potential theoretical vectors */ - for ( i=0; imatch_num; i++ ) { - for ( j=0; jmatch_num; j++ ) { + for ( i=0; imatch_num; i++ ) { + for ( j=0; jmatch_num; j++ ) { + double score = 0; - struct rvec *her_match = &her_obs->matches[i].vec; - struct rvec *his_match = &his_obs->matches[j].vec; + struct rvec *her_match = &her_obs->matches[i].vec; + struct rvec *his_match = &his_obs->matches[j].vec; - double theory_angle = rvec_angle(*her_match, - *his_match); + double theory_angle = rvec_angle(*her_match, + *his_match); - /* is this angle a match? */ + /* is this angle a match? */ - double angle_diff = fabs(theory_angle - obs_angle); + double angle_diff = fabs(theory_angle - obs_angle); - if ( angle_diff < cell->angle_tol ) { + /* within basic tolerance? */ + if ( angle_diff > cell->angle_tol ) { + continue; + } - // in the case of a brief check only - if (!her_match_idxs || !his_match_idxs) { - return 1; - } + double her_dist = rvec_length(*her_match); + double his_dist = rvec_length(*his_match); - /* If the angles are too close to 0 - or 180, one axis ill-determined */ - if (theory_angle < min_angle || - theory_angle > max_angle) { - continue; - } + score += (her_dist + his_dist) * angle_diff; - // check the third vector + /* If the angles are too close to 0 + or 180, one axis ill-determined */ + if (theory_angle < min_angle || + theory_angle > max_angle) { + continue; + } - struct rvec theory_diff = diff_vec(*his_match, *her_match); - struct rvec obs_diff = diff_vec(his_obs->obsvec, - her_obs->obsvec); + /* check that third vector adequately completes + * triangle */ - theory_angle = rvec_angle(*her_match, - theory_diff); - obs_angle = rvec_angle(her_obs->obsvec, obs_diff); + struct rvec theory_diff = diff_vec(*his_match, *her_match); + struct rvec obs_diff = diff_vec(his_obs->obsvec, + her_obs->obsvec); - if (fabs(obs_angle - theory_angle) > ANGLE_TOLERANCE) { - continue; - } + double obs_diff_dist = rvec_length(obs_diff); + theory_angle = rvec_angle(*her_match, + theory_diff); + obs_angle = rvec_angle(her_obs->obsvec, obs_diff); + angle_diff = fabs(obs_angle - theory_angle); - theory_angle = rvec_angle(*his_match, - theory_diff); - obs_angle = rvec_angle(his_obs->obsvec, obs_diff); + if (angle_diff > ANGLE_TOLERANCE) { + continue; + } + + score += (obs_diff_dist + her_dist) * angle_diff; + + 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; + } + + score += (obs_diff_dist + his_dist) * angle_diff; - if (fabs(obs_angle - theory_angle) > ANGLE_TOLERANCE) { - continue; - } + /* we add a new seed to the array */ + size_t new_size = (*match_count + 1); + new_size *= sizeof(struct Seed); - 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); + /* Reallocate the array to fit in another match */ + struct Seed *tmp_seeds = realloc(*seeds, new_size); - if ( temp_hers == NULL || temp_his == NULL ) { - apologise(); - } + if ( tmp_seeds == NULL ) { + apologise(); + } - (*her_match_idxs) = temp_hers; - (*his_match_idxs) = temp_his; + (*seeds) = tmp_seeds; - (*her_match_idxs)[*match_count] = i; - (*his_match_idxs)[*match_count] = j; - } + (*seeds)[*match_count].obs1 = her_obs; + (*seeds)[*match_count].obs2 = his_obs; + (*seeds)[*match_count].idx2 = j; + (*seeds)[*match_count].score = score; (*match_count)++; } @@ -1160,25 +1182,17 @@ static int find_seed(gsl_matrix **rotation, struct TakeTwoCell *cell) /* cell vector index matches stored in i, j and total * number stored in int matches. */ - int *i_idx = 0; - int *j_idx = 0; - int matches; + int seed_num = 0; + struct Seed *seeds = NULL; /* Check to see if any angles match from the cell vectors */ obs_vecs_match_angles(&obs_vecs[i], &obs_vecs[j], - &i_idx, &j_idx, &matches, cell); - if ( matches == 0 ) { - free(i_idx); - free(j_idx); - continue; - } + &seeds, &seed_num, cell); /* Weed out the duplicate seeds (from symmetric * reflection pairs) */ - - weed_duplicate_matches(&obs_vecs[i], &obs_vecs[j], - &i_idx, &j_idx, &matches, cell); + weed_duplicate_matches(&seeds, &matches, cell); /* We have seeds! Pass each of them through the seed-starter */ /* If a seed has the highest achieved membership, make note...*/ -- cgit v1.2.3 From 663170cf3df8168705e128359f5902145902c815 Mon Sep 17 00:00:00 2001 From: Helen Ginn Date: Sat, 21 Apr 2018 20:35:12 +0100 Subject: Temporary debug output --- libcrystfel/src/taketwo.c | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/libcrystfel/src/taketwo.c b/libcrystfel/src/taketwo.c index f92506cc..8f453308 100644 --- a/libcrystfel/src/taketwo.c +++ b/libcrystfel/src/taketwo.c @@ -1112,6 +1112,12 @@ static int grow_network(gsl_matrix *rot, int obs_idx1, int obs_idx2, *max_members = member_num; } + for (int n = 0; n < member_num; n++) + { + STATUS("*"); + } + STATUS("\n"); + /* If member_num is high enough, we want to return a yes */ if ( member_num > cell->member_thresh ) break; -- cgit v1.2.3 From 3a1a141ebf4f460665fb0d18fe063cb3393e81af Mon Sep 17 00:00:00 2001 From: Helen Ginn Date: Sat, 21 Apr 2018 20:58:15 +0100 Subject: Weed duplicates of seeds with new struct Seed. --- libcrystfel/src/taketwo.c | 52 ++++++++++++++++++++++++++--------------------- 1 file changed, 29 insertions(+), 23 deletions(-) diff --git a/libcrystfel/src/taketwo.c b/libcrystfel/src/taketwo.c index 8f453308..52602f14 100644 --- a/libcrystfel/src/taketwo.c +++ b/libcrystfel/src/taketwo.c @@ -645,7 +645,7 @@ static int obs_shares_spot_w_array(struct SpotVec *obs_vecs, int test_idx, static int obs_vecs_match_angles(struct SpotVec *her_obs, struct SpotVec *his_obs, - struct Seed *seeds, int *match_count, + struct Seed **seeds, int *match_count, struct TakeTwoCell *cell) { int i, j; @@ -734,13 +734,13 @@ static int obs_vecs_match_angles(struct SpotVec *her_obs, (*seeds)[*match_count].obs1 = her_obs; (*seeds)[*match_count].obs2 = his_obs; + (*seeds)[*match_count].idx1 = i; (*seeds)[*match_count].idx2 = j; (*seeds)[*match_count].score = score; (*match_count)++; } } - } return (*match_count > 0); } @@ -829,10 +829,8 @@ gsl_matrix *rot_mat_from_indices(struct SpotVec *her_obs, struct SpotVec *his_ob return mat; } -static int weed_duplicate_matches(struct SpotVec *her_obs, - struct SpotVec *his_obs, - int **her_match_idxs, int **his_match_idxs, - int *match_count, struct TakeTwoCell *cell) +static int weed_duplicate_matches(struct Seed **seeds, + int *match_count, struct TakeTwoCell *cell) { int num_occupied = 0; gsl_matrix **old_mats = calloc(*match_count, sizeof(gsl_matrix *)); @@ -849,11 +847,11 @@ static int weed_duplicate_matches(struct SpotVec *her_obs, /* First we weed out duplicates with previous solutions */ for (i = *match_count - 1; i >= 0; i--) { - int her_match = (*her_match_idxs)[i]; - int his_match = (*his_match_idxs)[i]; + int her_match = (*seeds)[i].idx1; + int his_match = (*seeds)[i].idx2; gsl_matrix *mat; - mat = rot_mat_from_indices(her_obs, his_obs, + mat = rot_mat_from_indices((*seeds)[i].obs1, (*seeds)[i].obs2, her_match, his_match, cell); if (mat == NULL) @@ -869,8 +867,8 @@ static int weed_duplicate_matches(struct SpotVec *her_obs, /* Found a duplicate with a previous solution */ if (sim) { - (*her_match_idxs)[i] = -1; - (*his_match_idxs)[i] = -1; + (*seeds)[i].idx1 = -1; + (*seeds)[i].idx2 = -1; break; } } @@ -881,11 +879,11 @@ static int weed_duplicate_matches(struct SpotVec *her_obs, /* Now we weed out the self-duplicates from the remaining batch */ for (i = *match_count - 1; i >= 0; i--) { - int her_match = (*her_match_idxs)[i]; - int his_match = (*his_match_idxs)[i]; + int her_match = (*seeds)[i].idx1; + int his_match = (*seeds)[i].idx2; gsl_matrix *mat; - mat = rot_mat_from_indices(her_obs, his_obs, + mat = rot_mat_from_indices((*seeds)[i].obs1, (*seeds)[i].obs2, her_match, his_match, cell); int found = 0; @@ -895,8 +893,8 @@ static int weed_duplicate_matches(struct SpotVec *her_obs, symm_rot_mats_are_similar(old_mats[j], mat, cell)) { // we have found a duplicate, so flag as bad. - (*her_match_idxs)[i] = -1; - (*his_match_idxs)[i] = -1; + (*seeds)[i].idx1 = -1; + (*seeds)[i].idx2 = -1; found = 1; duplicates++; @@ -1195,28 +1193,37 @@ static int find_seed(gsl_matrix **rotation, struct TakeTwoCell *cell) obs_vecs_match_angles(&obs_vecs[i], &obs_vecs[j], &seeds, &seed_num, cell); + if (seed_num == 0) + { + /* Nothing to clean up here */ + continue; + } + /* Weed out the duplicate seeds (from symmetric * reflection pairs) */ - weed_duplicate_matches(&seeds, &matches, cell); + weed_duplicate_matches(&seeds, &seed_num, cell); /* 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 Date: Sat, 21 Apr 2018 21:06:49 +0100 Subject: Declare integer outside for loop --- libcrystfel/src/taketwo.c | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/libcrystfel/src/taketwo.c b/libcrystfel/src/taketwo.c index 52602f14..efb5417d 100644 --- a/libcrystfel/src/taketwo.c +++ b/libcrystfel/src/taketwo.c @@ -1110,7 +1110,8 @@ static int grow_network(gsl_matrix *rot, int obs_idx1, int obs_idx2, *max_members = member_num; } - for (int n = 0; n < member_num; n++) + int n; + for (n = 0; n < member_num; n++) { STATUS("*"); } -- cgit v1.2.3 From a6de5faaf3e049a375199a6d0b92baf6a923a972 Mon Sep 17 00:00:00 2001 From: Helen Ginn Date: Sun, 22 Apr 2018 00:38:58 +0100 Subject: Change default distance to 0.12-something units for now --- libcrystfel/src/taketwo.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/libcrystfel/src/taketwo.c b/libcrystfel/src/taketwo.c index efb5417d..46033bf9 100644 --- a/libcrystfel/src/taketwo.c +++ b/libcrystfel/src/taketwo.c @@ -186,7 +186,7 @@ struct TakeTwoCell /* 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.0010*1e10) -- cgit v1.2.3 From 8d41a4d4dc0ffca5af573fa520638eab2eaa7c93 Mon Sep 17 00:00:00 2001 From: Helen Ginn Date: Sun, 22 Apr 2018 00:41:12 +0100 Subject: Now all seeds are found first, before sorting by score --- libcrystfel/src/taketwo.c | 240 +++++++++++++++++++++++++++------------------- 1 file changed, 143 insertions(+), 97 deletions(-) diff --git a/libcrystfel/src/taketwo.c b/libcrystfel/src/taketwo.c index 46033bf9..18f8365e 100644 --- a/libcrystfel/src/taketwo.c +++ b/libcrystfel/src/taketwo.c @@ -142,8 +142,8 @@ struct TheoryVec struct Seed { - struct SpotVec *obs1; - struct SpotVec *obs2; + int obs1; + int obs2; int idx1; int idx2; double score; @@ -163,9 +163,13 @@ struct TakeTwoCell UnitCell *cell; /**< Contains unit cell dimensions */ gsl_matrix **rotSymOps; unsigned int numOps; + gsl_matrix **prevSols; /**< Previous solutions to be ignored */ unsigned int numPrevs; /**< Previous solution count */ - struct SpotVec **obs_vecs; + + struct SpotVec *obs_vecs; + struct Seed *seeds; + int seed_count; int obs_vec_count; /* Options */ @@ -643,12 +647,14 @@ static int obs_shares_spot_w_array(struct SpotVec *obs_vecs, int test_idx, } -static int obs_vecs_match_angles(struct SpotVec *her_obs, - struct SpotVec *his_obs, +static int obs_vecs_match_angles(int her, int his, struct Seed **seeds, int *match_count, struct TakeTwoCell *cell) { - int i, j; + struct SpotVec *obs_vecs = cell->obs_vecs; + struct SpotVec *her_obs = &obs_vecs[her]; + struct SpotVec *his_obs = &obs_vecs[his]; + *match_count = 0; double min_angle = deg2rad(2.5); @@ -659,10 +665,17 @@ static int obs_vecs_match_angles(struct SpotVec *her_obs, /* calculate angle between all potential theoretical vectors */ + int i, j; for ( i=0; imatch_num; i++ ) { for ( j=0; jmatch_num; j++ ) { double score = 0; + if (her_obs->matches[i].asym == 0 && + his_obs->matches[j].asym == 0) + { + continue; + } + struct rvec *her_match = &her_obs->matches[i].vec; struct rvec *his_match = &his_obs->matches[j].vec; @@ -674,14 +687,15 @@ static int obs_vecs_match_angles(struct SpotVec *her_obs, double angle_diff = fabs(theory_angle - obs_angle); /* within basic tolerance? */ - if ( angle_diff > cell->angle_tol ) { + if ( angle_diff != angle_diff || + angle_diff > cell->angle_tol ) { continue; } - double her_dist = rvec_length(*her_match); - double his_dist = rvec_length(*his_match); - - score += (her_dist + his_dist) * angle_diff; + double add = angle_diff; + if (add == add) { + score += add; + } /* If the angles are too close to 0 or 180, one axis ill-determined */ @@ -697,7 +711,6 @@ static int obs_vecs_match_angles(struct SpotVec *her_obs, struct rvec obs_diff = diff_vec(his_obs->obsvec, her_obs->obsvec); - double obs_diff_dist = rvec_length(obs_diff); theory_angle = rvec_angle(*her_match, theory_diff); obs_angle = rvec_angle(her_obs->obsvec, obs_diff); @@ -707,7 +720,10 @@ static int obs_vecs_match_angles(struct SpotVec *her_obs, continue; } - score += (obs_diff_dist + her_dist) * angle_diff; + add = angle_diff; + if (add == add) { + score += add; + } theory_angle = rvec_angle(*his_match, theory_diff); @@ -717,7 +733,10 @@ static int obs_vecs_match_angles(struct SpotVec *her_obs, continue; } - score += (obs_diff_dist + his_dist) * angle_diff; + add = angle_diff; + if (add == add) { + score += add; + } /* we add a new seed to the array */ size_t new_size = (*match_count + 1); @@ -732,11 +751,11 @@ static int obs_vecs_match_angles(struct SpotVec *her_obs, (*seeds) = tmp_seeds; - (*seeds)[*match_count].obs1 = her_obs; - (*seeds)[*match_count].obs2 = his_obs; + (*seeds)[*match_count].obs1 = her; + (*seeds)[*match_count].obs2 = his; (*seeds)[*match_count].idx1 = i; (*seeds)[*match_count].idx2 = j; - (*seeds)[*match_count].score = score; + (*seeds)[*match_count].score = score * 1000; (*match_count)++; } @@ -814,10 +833,13 @@ static signed int finish_solution(gsl_matrix *rot, struct SpotVec *obs_vecs, return 1; } -gsl_matrix *rot_mat_from_indices(struct SpotVec *her_obs, struct SpotVec *his_obs, +gsl_matrix *rot_mat_from_indices(int her, int his, int her_match, int his_match, struct TakeTwoCell *cell) { + struct SpotVec *obs_vecs = cell->obs_vecs; + struct SpotVec *her_obs = &obs_vecs[her]; + struct SpotVec *his_obs = &obs_vecs[his]; struct rvec i_obsvec = her_obs->obsvec; struct rvec j_obsvec = his_obs->obsvec; struct rvec i_cellvec = her_obs->matches[her_match].vec; @@ -926,7 +948,7 @@ static signed int find_next_index(gsl_matrix *rot, int *obs_members, int *match_members, int start, int member_num, int *match_found, struct TakeTwoCell *cell) { - struct SpotVec *obs_vecs = *(cell->obs_vecs); + struct SpotVec *obs_vecs = cell->obs_vecs; int obs_vec_count = cell->obs_vec_count; gsl_matrix *sub = gsl_matrix_calloc(3, 3); gsl_matrix *mul = gsl_matrix_calloc(3, 3); @@ -1032,7 +1054,8 @@ static int grow_network(gsl_matrix *rot, int obs_idx1, int obs_idx2, int match_idx1, int match_idx2, int *max_members, struct TakeTwoCell *cell) { - struct SpotVec *obs_vecs = *(cell->obs_vecs); + + struct SpotVec *obs_vecs = cell->obs_vecs; int obs_vec_count = cell->obs_vec_count; int *obs_members; int *match_members; @@ -1098,9 +1121,6 @@ static int grow_network(gsl_matrix *rot, int obs_idx1, int obs_idx2, continue; } - /* we have elongated membership - so reset dead_ends counter */ - // dead_ends = 0; - obs_members[member_num] = next_index; match_members[member_num] = match_found; @@ -1136,7 +1156,8 @@ static int start_seed(int i, int j, int i_match, int j_match, gsl_matrix **rotation, int *max_members, struct TakeTwoCell *cell) { - struct SpotVec *obs_vecs = *(cell->obs_vecs); + STATUS("Start seed\n"); + struct SpotVec *obs_vecs = cell->obs_vecs; gsl_matrix *rot_mat; @@ -1157,22 +1178,26 @@ static int start_seed(int i, int j, int i_match, int j_match, return success; } -static int find_seed(gsl_matrix **rotation, struct TakeTwoCell *cell) +static int sort_seed_by_score(const void *av, const void *bv) { - struct SpotVec *obs_vecs = *(cell->obs_vecs); - int obs_vec_count = cell->obs_vec_count; + struct Seed *a = (struct Seed *)av; + struct Seed *b = (struct Seed *)bv; + return a->score > b->score; +} - /* META: Maximum achieved maximum network membership */ - int max_max_members = 0; - gsl_matrix *best_rotation = NULL; -// unsigned long start_time = time(NULL); +static int find_seeds(struct TakeTwoCell *cell) +{ + struct SpotVec *obs_vecs = cell->obs_vecs; + int obs_vec_count = cell->obs_vec_count; /* loop round pairs of vectors to try and find a suitable * seed to start building a self-consistent network of vectors */ int i, j; + STATUS("Total vectors: %i\n", obs_vec_count); + for ( i=0; iseed_count + seed_num; + new_size *= sizeof(struct Seed); + + struct Seed *tmp = realloc(cell->seeds, new_size); - if (seed_idx1 < 0 || seed_idx2 < 0) { + if (tmp == NULL) { + apologise(); + } + + cell->seeds = tmp; + + for (int i = 0; i < seed_num; i++) + { + if (seeds[i].idx1 < 0 || seeds[i].idx2 < 0) + { continue; } - int max_members = 0; - - int success = start_seed(i, j, - seed_idx1, seed_idx2, - rotation, &max_members, - cell); - - if (success) { - free(seeds); - gsl_matrix_free(best_rotation); - return success; - } else { - if (max_members > max_max_members) { - max_max_members = max_members; - gsl_matrix_free(best_rotation); - best_rotation = *rotation; - *rotation = NULL; - } else { - gsl_matrix_free(*rotation); - *rotation = NULL; - } - } + cell->seeds[cell->seed_count] = seeds[i]; + cell->seed_count++; } - - free(seeds); } - } /* yes this } is meant to be here */ + } - *rotation = best_rotation; + qsort(cell->seeds, cell->seed_count, sizeof(struct Seed), sort_seed_by_score); - if (best_rotation != NULL) + for (int i = 0; i < 10; i++) { - int num = cell->numPrevs + 1; - gsl_matrix **tmp = malloc(sizeof(gsl_matrix *) * num); + struct Seed seed = cell->seeds[i]; + STATUS("%i %i %i %i %.3f\n", seed.idx1, seed.idx2, + seed.obs1, seed.obs2, seed.score); + } - if (tmp) { - cell->prevSols = tmp; - } else { - return 0; + return 1; +} + +static int start_seeds(gsl_matrix **rotation, struct TakeTwoCell *cell) +{ + struct Seed *seeds = cell->seeds; + int seed_num = cell->seed_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; kprevSols[num - 1] = best_rotation; - cell->numPrevs = num; + int max_members = 0; + + int seed_obs1 = seeds[k].obs1; + int seed_obs2 = seeds[k].obs2; + + int success = start_seed(seed_obs1, seed_obs2, + seed_idx1, seed_idx2, + rotation, &max_members, + cell); + + if (success) { + free(seeds); + return success; + } } - return (best_rotation != NULL); + free(seeds); + return 0; } + 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) @@ -1362,7 +1402,7 @@ struct sortme double dist; }; -static int sort_func(const void *av, const void *bv) +static int sort_theory_distances(const void *av, const void *bv) { struct sortme *a = (struct sortme *)av; struct sortme *b = (struct sortme *)bv; @@ -1371,9 +1411,10 @@ static int sort_func(const void *av, const void *bv) static int match_obs_to_cell_vecs(struct TheoryVec *cell_vecs, int cell_vec_count, - struct SpotVec *obs_vecs, int obs_vec_count, struct TakeTwoCell *cell) { + struct SpotVec *obs_vecs = cell->obs_vecs; + int obs_vec_count = cell->obs_vec_count; int i, j; for ( i=0; iobs_vecs, count*sizeof(struct SpotVec)); if ( temp_obs_vecs == NULL ) { return 0; } else { - *obs_vecs = temp_obs_vecs; + cell->obs_vecs = temp_obs_vecs; /* initialise all SpotVec struct members */ @@ -1475,15 +1515,16 @@ static int gen_observed_vecs(struct rvec *rlps, int rlp_count, spot_vec.her_rlp = &rlps[i]; spot_vec.his_rlp = &rlps[j]; - (*obs_vecs)[count - 1] = spot_vec; + cell->obs_vecs[count - 1] = spot_vec; } } } /* Sort such that the shortest distances are searched first. */ - qsort(*obs_vecs, count, sizeof(struct SpotVec), compare_spot_vecs); + qsort(cell->obs_vecs, count, sizeof(struct SpotVec), compare_spot_vecs); - *obs_vec_count = count; + cell->obs_vec_count = count; + STATUS("Generated observed vectors.\n"); return 1; } @@ -1586,6 +1627,9 @@ static void cleanup_taketwo_cell(struct TakeTwoCell *ttCell) gsl_matrix_free(ttCell->rotSymOps[i]); } + cleanup_taketwo_obs_vecs(ttCell->obs_vecs, + ttCell->obs_vec_count); + free(ttCell->vec1Tmp); free(ttCell->vec2Tmp); free(ttCell->twiz1Tmp); @@ -1612,14 +1656,16 @@ static UnitCell *run_taketwo(UnitCell *cell, const struct taketwo_options *opts, int cell_vec_count = 0; struct TheoryVec *theory_vecs = NULL; UnitCell *result; - int obs_vec_count = 0; - struct SpotVec *obs_vecs = NULL; int success = 0; gsl_matrix *solution = NULL; + /* Initialise TakeTwoCell */ struct TakeTwoCell ttCell; ttCell.cell = cell; + ttCell.seeds = NULL; + ttCell.seed_count = 0; ttCell.rotSymOps = NULL; + ttCell.obs_vecs = NULL; ttCell.prevSols = NULL; ttCell.twiz1Tmp = gsl_matrix_calloc(3, 3); ttCell.twiz2Tmp = gsl_matrix_calloc(3, 3); @@ -1627,18 +1673,16 @@ static UnitCell *run_taketwo(UnitCell *cell, const struct taketwo_options *opts, ttCell.vec2Tmp = gsl_vector_calloc(3); ttCell.numOps = 0; ttCell.numPrevs = 0; + ttCell.obs_vec_count = 0; success = generate_rotation_sym_ops(&ttCell); success = gen_theoretical_vecs(cell, &theory_vecs, &cell_vec_count); if ( !success ) return NULL; - success = gen_observed_vecs(rlps, rlp_count, &obs_vecs, &obs_vec_count); + success = gen_observed_vecs(rlps, rlp_count, &ttCell); if ( !success ) return NULL; - ttCell.obs_vecs = &obs_vecs; - ttCell.obs_vec_count = obs_vec_count; - if ( opts->member_thresh < 0 ) { ttCell.member_thresh = NETWORK_MEMBER_THRESHOLD; } else { @@ -1664,21 +1708,23 @@ static UnitCell *run_taketwo(UnitCell *cell, const struct taketwo_options *opts, } success = match_obs_to_cell_vecs(theory_vecs, cell_vec_count, - obs_vecs, obs_vec_count, &ttCell); + &ttCell); free(theory_vecs); if ( !success ) return NULL; - find_seed(&solution, &ttCell); + STATUS("Find seeds\n"); + find_seeds(&ttCell); + start_seeds(&solution, &ttCell); if ( solution == NULL ) { - cleanup_taketwo_obs_vecs(obs_vecs, obs_vec_count); return NULL; } + STATUS("Returning something.\n"); + result = transform_cell_gsl(cell, solution); - cleanup_taketwo_obs_vecs(obs_vecs, obs_vec_count); gsl_matrix_free(solution); cleanup_taketwo_cell(&ttCell); -- cgit v1.2.3 From 8f14adf07254041401a81fbf50f5d3a2912704ab Mon Sep 17 00:00:00 2001 From: Helen Ginn Date: Sun, 22 Apr 2018 01:08:26 +0100 Subject: Search through entire vector array each time --- libcrystfel/src/taketwo.c | 18 ++++++++++++++++-- 1 file changed, 16 insertions(+), 2 deletions(-) diff --git a/libcrystfel/src/taketwo.c b/libcrystfel/src/taketwo.c index 18f8365e..8303016e 100644 --- a/libcrystfel/src/taketwo.c +++ b/libcrystfel/src/taketwo.c @@ -120,6 +120,7 @@ struct SpotVec struct rvec obsvec; struct TheoryVec *matches; int match_num; + int in_network; double distance; struct rvec *her_rlp; struct rvec *his_rlp; @@ -957,6 +958,12 @@ static signed int find_next_index(gsl_matrix *rot, int *obs_members, for ( i=start; imember_thresh+3)*sizeof(int)); match_members = malloc((cell->member_thresh+3)*sizeof(int)); @@ -1095,7 +1109,7 @@ static int grow_network(gsl_matrix *rot, int obs_idx1, int obs_idx2, signed int next_index = find_next_index(rot, obs_members, match_members, - start, member_num, + 0, member_num, &match_found, cell); if ( member_num < 2 ) { @@ -1114,13 +1128,13 @@ static int grow_network(gsl_matrix *rot, int obs_idx1, int obs_idx2, /* We have not had too many dead ends. Try removing the last member and continue. */ - start = obs_members[member_num - 1] + 1; member_num--; dead_ends++; continue; } + obs_vecs[next_index].in_network = 1; obs_members[member_num] = next_index; match_members[member_num] = match_found; -- cgit v1.2.3 From 596cd5f5199367e7634e0568abe9dede854fefb8 Mon Sep 17 00:00:00 2001 From: Helen Ginn Date: Sun, 22 Apr 2018 01:08:59 +0100 Subject: Reset default reciprocal distance to 0.15*1e10 m-1 --- libcrystfel/src/taketwo.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/libcrystfel/src/taketwo.c b/libcrystfel/src/taketwo.c index 8303016e..add0b099 100644 --- a/libcrystfel/src/taketwo.c +++ b/libcrystfel/src/taketwo.c @@ -191,7 +191,7 @@ struct TakeTwoCell /* Maximum distance between two rlp sizes to consider info for indexing */ -#define MAX_RECIP_DISTANCE (0.12*1e10) +#define MAX_RECIP_DISTANCE (0.15*1e10) /* Tolerance for two lengths in reciprocal space to be considered the same */ #define RECIP_TOLERANCE (0.0010*1e10) -- cgit v1.2.3 From 2c6b2bde551748616e0f25908bed6f790e700a5b Mon Sep 17 00:00:00 2001 From: Helen Ginn Date: Sun, 22 Apr 2018 01:09:43 +0100 Subject: Remove some shouty messages --- libcrystfel/src/taketwo.c | 16 ++-------------- 1 file changed, 2 insertions(+), 14 deletions(-) diff --git a/libcrystfel/src/taketwo.c b/libcrystfel/src/taketwo.c index add0b099..5d09a3a6 100644 --- a/libcrystfel/src/taketwo.c +++ b/libcrystfel/src/taketwo.c @@ -1144,12 +1144,14 @@ static int grow_network(gsl_matrix *rot, int obs_idx1, int obs_idx2, *max_members = member_num; } + /* int n; for (n = 0; n < member_num; n++) { STATUS("*"); } STATUS("\n"); + */ /* If member_num is high enough, we want to return a yes */ if ( member_num > cell->member_thresh ) break; @@ -1170,7 +1172,6 @@ static int start_seed(int i, int j, int i_match, int j_match, gsl_matrix **rotation, int *max_members, struct TakeTwoCell *cell) { - STATUS("Start seed\n"); struct SpotVec *obs_vecs = cell->obs_vecs; gsl_matrix *rot_mat; @@ -1210,8 +1211,6 @@ static int find_seeds(struct TakeTwoCell *cell) */ int i, j; - STATUS("Total vectors: %i\n", obs_vec_count); - for ( i=0; iseeds, cell->seed_count, sizeof(struct Seed), sort_seed_by_score); - for (int i = 0; i < 10; i++) - { - struct Seed seed = cell->seeds[i]; - STATUS("%i %i %i %i %.3f\n", seed.idx1, seed.idx2, - seed.obs1, seed.obs2, seed.score); - } - return 1; } @@ -1538,7 +1530,6 @@ static int gen_observed_vecs(struct rvec *rlps, int rlp_count, qsort(cell->obs_vecs, count, sizeof(struct SpotVec), compare_spot_vecs); cell->obs_vec_count = count; - STATUS("Generated observed vectors.\n"); return 1; } @@ -1728,7 +1719,6 @@ static UnitCell *run_taketwo(UnitCell *cell, const struct taketwo_options *opts, if ( !success ) return NULL; - STATUS("Find seeds\n"); find_seeds(&ttCell); start_seeds(&solution, &ttCell); @@ -1736,8 +1726,6 @@ static UnitCell *run_taketwo(UnitCell *cell, const struct taketwo_options *opts, return NULL; } - STATUS("Returning something.\n"); - result = transform_cell_gsl(cell, solution); gsl_matrix_free(solution); cleanup_taketwo_cell(&ttCell); -- cgit v1.2.3 From a2ba046cf1a9c15bd8af6d30ddd5f2c358ee78eb Mon Sep 17 00:00:00 2001 From: Helen Ginn Date: Sun, 22 Apr 2018 01:09:55 +0100 Subject: Documentation is wrong notice --- libcrystfel/src/taketwo.c | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/libcrystfel/src/taketwo.c b/libcrystfel/src/taketwo.c index 5d09a3a6..4b097868 100644 --- a/libcrystfel/src/taketwo.c +++ b/libcrystfel/src/taketwo.c @@ -39,7 +39,9 @@ * Match observed vectors to theoretical vectors (match_obs_to_cell_vecs()) * * --- Business bit --- - * Find starting seeds (find_seed()): + * ... n.b. rearranging to find all seeds in advance. + * + * Find starting seeds (find_seeds()): * - Loop through pairs of observed vectors * - If they share a spot, find matching pairs of theoretical vectors * - Remove all duplicate matches due to symmetry operations -- cgit v1.2.3 From 480d6349dc1c2f517e7fd653199f4952da367edf Mon Sep 17 00:00:00 2001 From: Helen Ginn Date: Sun, 22 Apr 2018 01:10:24 +0100 Subject: Return success even if rotation is not from full network membership --- libcrystfel/src/taketwo.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/libcrystfel/src/taketwo.c b/libcrystfel/src/taketwo.c index 4b097868..bc74a13e 100644 --- a/libcrystfel/src/taketwo.c +++ b/libcrystfel/src/taketwo.c @@ -1308,7 +1308,7 @@ static int start_seeds(gsl_matrix **rotation, struct TakeTwoCell *cell) } free(seeds); - return 0; + return (rotation != NULL); } -- cgit v1.2.3 From 8401b8ee23a80ad2258fac28fe7644c6f4b27976 Mon Sep 17 00:00:00 2001 From: Helen Ginn Date: Sun, 22 Apr 2018 01:39:52 +0100 Subject: Different limitations for making new seeds --- libcrystfel/src/taketwo.c | 19 ++++++++++++------- 1 file changed, 12 insertions(+), 7 deletions(-) diff --git a/libcrystfel/src/taketwo.c b/libcrystfel/src/taketwo.c index bc74a13e..93f348de 100644 --- a/libcrystfel/src/taketwo.c +++ b/libcrystfel/src/taketwo.c @@ -673,12 +673,6 @@ static int obs_vecs_match_angles(int her, int his, for ( j=0; jmatch_num; j++ ) { double score = 0; - if (her_obs->matches[i].asym == 0 && - his_obs->matches[j].asym == 0) - { - continue; - } - struct rvec *her_match = &her_obs->matches[i].vec; struct rvec *his_match = &his_obs->matches[j].vec; @@ -1214,14 +1208,21 @@ static int find_seeds(struct TakeTwoCell *cell) int i, j; for ( i=0; i + MAX_RECIP_DISTANCE) { + continue; + } + /** Check to see if there is a shared spot - opportunity * for optimisation by generating a look-up table * by spot instead of by vector. */ int shared = obs_vecs_share_spot(&obs_vecs[i], &obs_vecs[j]); - if ( !shared ) continue; /* cell vector index matches stored in i, j and total @@ -1272,6 +1273,10 @@ static int find_seeds(struct TakeTwoCell *cell) qsort(cell->seeds, cell->seed_count, sizeof(struct Seed), sort_seed_by_score); + if (cell->seed_count > 1000) { + cell->seed_count = 1000; + } + return 1; } -- cgit v1.2.3 From 8ee6527f2eea3e87970b7a8d7c2df1c4c0389afa Mon Sep 17 00:00:00 2001 From: Helen Ginn Date: Fri, 27 Apr 2018 22:20:56 +0200 Subject: Fix typo --- libcrystfel/src/taketwo.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/libcrystfel/src/taketwo.c b/libcrystfel/src/taketwo.c index 93f348de..5acc99d9 100644 --- a/libcrystfel/src/taketwo.c +++ b/libcrystfel/src/taketwo.c @@ -32,7 +32,7 @@ * \class TakeTwo * Code outline. * --- Get ready for calculation --- - * Pre-calculate symmertry operations (generate_rotation_symops()) + * Pre-calculate symmetry operations (generate_rotation_symops()) * Pre-calculate theoretical vectors from unit cell dimensions * (gen_theoretical_vecs()) * Generate observed vectors from data (gen_observed_vecs()) -- cgit v1.2.3 From f06de02c7350ebf41c91cac7309751317c4c0943 Mon Sep 17 00:00:00 2001 From: Helen Ginn Date: Fri, 27 Apr 2018 22:24:10 +0200 Subject: Taketwo private structure takes care of duplicate solutions --- libcrystfel/src/taketwo.c | 165 ++++++++++++++++++++++++++++++++++------------ 1 file changed, 122 insertions(+), 43 deletions(-) diff --git a/libcrystfel/src/taketwo.c b/libcrystfel/src/taketwo.c index 5acc99d9..0dc1ed51 100644 --- a/libcrystfel/src/taketwo.c +++ b/libcrystfel/src/taketwo.c @@ -156,6 +156,12 @@ struct taketwo_private { IndexingMethod indm; UnitCell *cell; + int serial_num; /* -1 if unassigned */ + unsigned int attempts; + + gsl_matrix **prevSols; /**< Previous solutions to be ignored */ + unsigned int numPrevs; /**< Previous solution count */ + }; /** @@ -167,9 +173,6 @@ struct TakeTwoCell gsl_matrix **rotSymOps; unsigned int numOps; - gsl_matrix **prevSols; /**< Previous solutions to be ignored */ - unsigned int numPrevs; /**< Previous solution count */ - struct SpotVec *obs_vecs; struct Seed *seeds; int seed_count; @@ -861,39 +864,8 @@ static int weed_duplicate_matches(struct Seed **seeds, } signed int i, j; - int duplicates = 0; - - /* First we weed out duplicates with previous solutions */ - - for (i = *match_count - 1; i >= 0; i--) { - int her_match = (*seeds)[i].idx1; - int his_match = (*seeds)[i].idx2; - - gsl_matrix *mat; - mat = rot_mat_from_indices((*seeds)[i].obs1, (*seeds)[i].obs2, - her_match, his_match, cell); - - if (mat == NULL) - { - continue; - } - - for (j = 0; j < cell->numPrevs; j++) - { - int sim = symm_rot_mats_are_similar(cell->prevSols[j], - mat, cell); - /* Found a duplicate with a previous solution */ - if (sim) - { - (*seeds)[i].idx1 = -1; - (*seeds)[i].idx2 = -1; - break; - } - } - - gsl_matrix_free(mat); - } + int duplicates = 0; /* Now we weed out the self-duplicates from the remaining batch */ @@ -1196,8 +1168,51 @@ static int sort_seed_by_score(const void *av, const void *bv) return a->score > b->score; } +static void remove_old_solutions(struct TakeTwoCell *cell, + struct taketwo_private *tp) +{ + int duplicates = 0; + struct Seed *seeds = cell->seeds; + unsigned int total = cell->seed_count; + + /* First we remove duplicates with previous solutions */ + + int i, j; + for (i = total - 1; i >= 0; i--) { + int her_match = seeds[i].idx1; + int his_match = seeds[i].idx2; + + gsl_matrix *mat; + mat = rot_mat_from_indices(seeds[i].obs1, seeds[i].obs2, + her_match, his_match, cell); + + if (mat == NULL) + { + continue; + } + + for (j = 0; j < tp->numPrevs; j++) + { + int sim = symm_rot_mats_are_similar(tp->prevSols[j], + mat, cell); + + /* Found a duplicate with a previous solution */ + if (sim) + { + seeds[i].idx1 = -1; + seeds[i].idx2 = -1; + duplicates++; + break; + } + } -static int find_seeds(struct TakeTwoCell *cell) + gsl_matrix_free(mat); + } + + STATUS("Removing %i duplicates due to prev solutions.\n", duplicates); +} + +static int find_seeds(struct TakeTwoCell *cell, struct taketwo_private *tp) { struct SpotVec *obs_vecs = cell->obs_vecs; int obs_vec_count = cell->obs_vec_count; @@ -1663,7 +1678,8 @@ static void cleanup_taketwo_cell(struct TakeTwoCell *ttCell) * successful. **/ static UnitCell *run_taketwo(UnitCell *cell, const struct taketwo_options *opts, - struct rvec *rlps, int rlp_count) + struct rvec *rlps, int rlp_count, + struct taketwo_private *tp) { int cell_vec_count = 0; struct TheoryVec *theory_vecs = NULL; @@ -1678,13 +1694,11 @@ static UnitCell *run_taketwo(UnitCell *cell, const struct taketwo_options *opts, ttCell.seed_count = 0; ttCell.rotSymOps = NULL; ttCell.obs_vecs = NULL; - ttCell.prevSols = NULL; ttCell.twiz1Tmp = gsl_matrix_calloc(3, 3); ttCell.twiz2Tmp = gsl_matrix_calloc(3, 3); ttCell.vec1Tmp = gsl_vector_calloc(3); ttCell.vec2Tmp = gsl_vector_calloc(3); ttCell.numOps = 0; - ttCell.numPrevs = 0; ttCell.obs_vec_count = 0; success = generate_rotation_sym_ops(&ttCell); @@ -1726,20 +1740,62 @@ static UnitCell *run_taketwo(UnitCell *cell, const struct taketwo_options *opts, if ( !success ) return NULL; - find_seeds(&ttCell); + find_seeds(&ttCell, tp); +// remove_old_solutions(&ttCell, tp); start_seeds(&solution, &ttCell); if ( solution == NULL ) { return NULL; } + for (int i = 0; i < tp->numPrevs; i++) + { + gsl_matrix *sol = tp->prevSols[i]; + + int sim = symm_rot_mats_are_similar(sol, solution, &ttCell); + if (sim) + { +// STATUS("Warning! Returning previous solution.\n"); + } + } + + int new_size = (tp->numPrevs + 1) * sizeof(gsl_matrix *); + gsl_matrix **tmp = realloc(tp->prevSols, new_size); + + if (!tmp) { + apologise(); + } + + tp->prevSols = tmp; + + tp->prevSols[tp->numPrevs] = solution; + tp->numPrevs++; + result = transform_cell_gsl(cell, solution); - gsl_matrix_free(solution); cleanup_taketwo_cell(&ttCell); return result; } +/* Cleans up the per-image information for taketwo */ + +static void partial_taketwo_cleanup(struct taketwo_private *tp) +{ + if (tp->prevSols != NULL) + { + for (int i = 0; i < tp->numPrevs; i++) + { + gsl_matrix_free(tp->prevSols[i]); + } + + free(tp->prevSols); + } + + tp->attempts = 0; + tp->numPrevs = 0; + tp->prevSols = NULL; + +} /* CrystFEL interface hooks */ @@ -1753,6 +1809,24 @@ int taketwo_index(struct image *image, const struct taketwo_options *opts, int i; struct taketwo_private *tp = (struct taketwo_private *)priv; + /* Check serial number against previous for solution tracking */ + int this_serial = image->serial; + + if (tp->serial_num == this_serial) + { + tp->attempts++; + } + else + { + partial_taketwo_cleanup(tp); + tp->serial_num = this_serial; + } + + /* + STATUS("Indexing %i with %i attempts, %i crystals\n", this_serial, tp->attempts, + image->n_crystals); + */ + rlps = malloc((image_feature_count(image->features)+1)*sizeof(struct rvec)); for ( i=0; ifeatures); i++ ) { struct imagefeature *pk = image_get_feature(image->features, i); @@ -1766,7 +1840,7 @@ int taketwo_index(struct image *image, const struct taketwo_options *opts, rlps[n_rlps].v = 0.0; rlps[n_rlps++].w = 0.0; - cell = run_taketwo(tp->cell, opts, rlps, n_rlps); + cell = run_taketwo(tp->cell, opts, rlps, n_rlps, tp); free(rlps); if ( cell == NULL ) return 0; @@ -1829,14 +1903,19 @@ void *taketwo_prepare(IndexingMethod *indm, UnitCell *cell) tp->cell = cell; tp->indm = *indm; + tp->serial_num = -1; + tp->attempts = 0; + tp->prevSols = NULL; + tp->numPrevs = 0; return tp; } - void taketwo_cleanup(IndexingPrivate *pp) { struct taketwo_private *tp = (struct taketwo_private *)pp; + + partial_taketwo_cleanup(tp); free(tp); } -- cgit v1.2.3 From d673548e1f0a21fbd677923c2ab0bcbcb487c742 Mon Sep 17 00:00:00 2001 From: Helen Ginn Date: Fri, 27 Apr 2018 22:24:42 +0200 Subject: rvec add rvec function --- libcrystfel/src/taketwo.c | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/libcrystfel/src/taketwo.c b/libcrystfel/src/taketwo.c index 0dc1ed51..082e7886 100644 --- a/libcrystfel/src/taketwo.c +++ b/libcrystfel/src/taketwo.c @@ -245,6 +245,15 @@ static struct rvec new_rvec(double new_u, double new_v, double new_w) return new_rvector; } +static struct rvec rvec_add_rvec(struct rvec first, struct rvec second) +{ + struct rvec diff = new_rvec(second.u + first.u, + second.v + first.v, + second.w + first.w); + + return diff; +} + static struct rvec diff_vec(struct rvec from, struct rvec to) { -- cgit v1.2.3 From 51269217ef086a38c1244c3f7154c1efef191449 Mon Sep 17 00:00:00 2001 From: Helen Ginn Date: Fri, 27 Apr 2018 22:24:52 +0200 Subject: decompose rotation matrix into angle and axis --- libcrystfel/src/taketwo.c | 26 ++++++++++++++++++++++++++ 1 file changed, 26 insertions(+) diff --git a/libcrystfel/src/taketwo.c b/libcrystfel/src/taketwo.c index 082e7886..a24a30c0 100644 --- a/libcrystfel/src/taketwo.c +++ b/libcrystfel/src/taketwo.c @@ -406,6 +406,32 @@ 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) { -- cgit v1.2.3 From b0c9de361ded510d1fb3facbdb21925b7f3625ab Mon Sep 17 00:00:00 2001 From: Helen Ginn Date: Fri, 27 Apr 2018 22:27:45 +0200 Subject: Change score for seed --- libcrystfel/src/taketwo.c | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/libcrystfel/src/taketwo.c b/libcrystfel/src/taketwo.c index a24a30c0..bbd3883d 100644 --- a/libcrystfel/src/taketwo.c +++ b/libcrystfel/src/taketwo.c @@ -714,6 +714,9 @@ static int obs_vecs_match_angles(int her, int his, struct rvec *her_match = &her_obs->matches[i].vec; struct rvec *his_match = &his_obs->matches[j].vec; + double her_dist = rvec_length(*her_match); + double his_dist = rvec_length(*his_match); + double theory_angle = rvec_angle(*her_match, *his_match); @@ -729,7 +732,7 @@ static int obs_vecs_match_angles(int her, int his, double add = angle_diff; if (add == add) { - score += add; + score += add * her_dist * his_dist; } /* If the angles are too close to 0 @@ -751,13 +754,15 @@ static int obs_vecs_match_angles(int her, int his, obs_angle = rvec_angle(her_obs->obsvec, obs_diff); angle_diff = fabs(obs_angle - theory_angle); + double diff_dist = rvec_length(obs_diff); + if (angle_diff > ANGLE_TOLERANCE) { continue; } add = angle_diff; if (add == add) { - score += add; + score += add * her_dist * diff_dist; } theory_angle = rvec_angle(*his_match, @@ -770,7 +775,7 @@ static int obs_vecs_match_angles(int her, int his, add = angle_diff; if (add == add) { - score += add; + score += add * his_dist * diff_dist; } /* we add a new seed to the array */ -- cgit v1.2.3 From be84751456dfd993472f64994cb26542b272dbb6 Mon Sep 17 00:00:00 2001 From: Helen Ginn Date: Fri, 27 Apr 2018 22:28:29 +0200 Subject: Fix a for loop --- libcrystfel/src/taketwo.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/libcrystfel/src/taketwo.c b/libcrystfel/src/taketwo.c index bbd3883d..293d9551 100644 --- a/libcrystfel/src/taketwo.c +++ b/libcrystfel/src/taketwo.c @@ -1262,7 +1262,7 @@ static int find_seeds(struct TakeTwoCell *cell, struct taketwo_private *tp) */ int i, j; - for ( i=0; i Date: Fri, 27 Apr 2018 22:29:03 +0200 Subject: Give up if seed count already high --- libcrystfel/src/taketwo.c | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/libcrystfel/src/taketwo.c b/libcrystfel/src/taketwo.c index 293d9551..0cf1aeb6 100644 --- a/libcrystfel/src/taketwo.c +++ b/libcrystfel/src/taketwo.c @@ -1323,14 +1323,20 @@ static int find_seeds(struct TakeTwoCell *cell, struct taketwo_private *tp) cell->seeds[cell->seed_count] = seeds[i]; cell->seed_count++; } + + if (cell->seed_count > 1000) { + break; + } } + + if (cell->seed_count > 1000) { + break; + } + } qsort(cell->seeds, cell->seed_count, sizeof(struct Seed), sort_seed_by_score); - if (cell->seed_count > 1000) { - cell->seed_count = 1000; - } return 1; } -- cgit v1.2.3 From 51eafe3fbd147a05ef77c75cbdf388f785717cd0 Mon Sep 17 00:00:00 2001 From: Helen Ginn Date: Sat, 28 Apr 2018 21:51:11 +0200 Subject: Now we refine an indexing solution to match all the observed vectors to the nearest calculated vector before returning solution --- libcrystfel/src/taketwo.c | 247 +++++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 235 insertions(+), 12 deletions(-) diff --git a/libcrystfel/src/taketwo.c b/libcrystfel/src/taketwo.c index 0cf1aeb6..9d557fd2 100644 --- a/libcrystfel/src/taketwo.c +++ b/libcrystfel/src/taketwo.c @@ -122,6 +122,7 @@ struct SpotVec struct rvec obsvec; struct TheoryVec *matches; int match_num; + int assignment; int in_network; double distance; struct rvec *her_rlp; @@ -184,6 +185,13 @@ struct TakeTwoCell double angle_tol; /**< In radians */ double trace_tol; /**< Contains sqrt(4*(1-cos(angle))) */ + /** A given solution to refine */ + gsl_matrix *solution; + + double x_ang; /**< Rotations in radians to apply to x axis of solution */ + double y_ang; /**< Rotations in radians to apply to y axis of solution */ + double z_ang; /**< Rotations in radians to apply to z axis of solution */ + /**< Temporary memory always allocated for calculations */ gsl_matrix *twiz1Tmp; /**< Temporary memory always allocated for calculations */ @@ -217,6 +225,12 @@ struct TakeTwoCell /* Tolerance for rot_mats_are_similar */ #define TRACE_TOLERANCE (deg2rad(3.0)) +/* Initial step size for refinement of solutions */ +#define ANGLE_STEP_SIZE (deg2rad(0.2)) + +/* Final required step size for refinement of solutions */ +#define ANGLE_CONVERGE_SIZE (deg2rad(0.01)) + /* TODO: Multiple lattices */ @@ -349,6 +363,35 @@ static void rotation_around_axis(struct rvec c, double th, gsl_matrix_set(res, 2, 2, cos(th) + c.w*c.w*omc); } +/** Rotate GSL matrix by three angles along x, y and z axes */ +static void rotate_gsl_by_angles(gsl_matrix *sol, double x, double y, + double z, gsl_matrix *result) +{ + gsl_matrix *x_rot = gsl_matrix_alloc(3, 3); + gsl_matrix *y_rot = gsl_matrix_alloc(3, 3); + gsl_matrix *z_rot = gsl_matrix_alloc(3, 3); + gsl_matrix *xy_rot = gsl_matrix_alloc(3, 3); + gsl_matrix *xyz_rot = gsl_matrix_alloc(3, 3); + + struct rvec x_axis = new_rvec(1, 0, 0); + struct rvec y_axis = new_rvec(1, 0, 0); + struct rvec z_axis = new_rvec(1, 0, 0); + + rotation_around_axis(x_axis, x, x_rot); + rotation_around_axis(y_axis, y, y_rot); + rotation_around_axis(z_axis, z, z_rot); + + /* Collapse the rotations in x and y onto z */ + gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, x_rot, + y_rot, 0.0, xy_rot); + gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, xy_rot, + z_rot, 0.0, xyz_rot); + + /* Apply the whole rotation offset to the solution */ + gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, xyz_rot, + sol, 0.0, result); +} + /* Rotate vector (vec1) around axis (axis) by angle theta. Find value of * theta for which the angle between (vec1) and (vec2) is minimised. */ @@ -406,6 +449,7 @@ 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); @@ -432,6 +476,7 @@ static struct rvec matrix_axis(gsl_matrix *a) struct rvec v = new_rvec(x, y, z); return v; } +*/ static double matrix_trace(gsl_matrix *a) { @@ -1249,7 +1294,7 @@ static void remove_old_solutions(struct TakeTwoCell *cell, gsl_matrix_free(mat); } - STATUS("Removing %i duplicates due to prev solutions.\n", duplicates); +// STATUS("Removing %i duplicates due to prev solutions.\n", duplicates); } static int find_seeds(struct TakeTwoCell *cell, struct taketwo_private *tp) @@ -1323,16 +1368,7 @@ static int find_seeds(struct TakeTwoCell *cell, struct taketwo_private *tp) cell->seeds[cell->seed_count] = seeds[i]; cell->seed_count++; } - - if (cell->seed_count > 1000) { - break; - } } - - if (cell->seed_count > 1000) { - break; - } - } qsort(cell->seeds, cell->seed_count, sizeof(struct Seed), sort_seed_by_score); @@ -1483,6 +1519,182 @@ static int sort_theory_distances(const void *av, const void *bv) return a->dist > b->dist; } +static double obs_to_sol_score(struct TakeTwoCell *ttCell) +{ + double total = 0; + int count = 0; + int i; + gsl_matrix *solution = ttCell->solution; + gsl_matrix *full_rot = gsl_matrix_alloc(3, 3); + rotate_gsl_by_angles(solution, ttCell->x_ang, ttCell->y_ang, + ttCell->z_ang, full_rot); + + for (i = 0; i < ttCell->obs_vec_count; i++) + { + struct rvec *obs = &ttCell->obs_vecs[i].obsvec; + rvec_to_gsl(ttCell->vec1Tmp, *obs); + + /* Rotate all the observed vectors by the modified soln */ + /* ttCell->vec2Tmp = 1.0 * full_rot * ttCell->vec1Tmp */ + gsl_blas_dgemv(CblasTrans, 1.0, full_rot, ttCell->vec1Tmp, + 0.0, ttCell->vec2Tmp); + struct rvec rotated = gsl_to_rvec(ttCell->vec2Tmp); + + int j = ttCell->obs_vecs[i].assignment; + + if (j < 0) continue; + + struct rvec *match = &ttCell->obs_vecs[i].matches[j].vec; + struct rvec diff = diff_vec(rotated, *match); + + double length = rvec_length(diff); + + double addition = exp(-(1 / RECIP_TOLERANCE) * length); + total += addition; + count++; + } + + total /= (double)-count; + + gsl_matrix_free(full_rot); + + return total; +} + +/* At the moment this is just going to do a step search of 3x3x3 until + * the angles are sufficiently tiny */ +static void refine_solution(struct TakeTwoCell *ttCell) +{ + int i, j, k; + const int total = 3 * 3 * 3; + const int middle = (total - 1) / 2; + + struct rvec steps[total]; + double start = obs_to_sol_score(ttCell); + const int max_tries = 100; + + int count = 0; + double size = ANGLE_STEP_SIZE; + + /* First we create our combinations of steps */ + for (i = -1; i <= 1; i++) { + for (j = -1; j <= 1; j++) { + for (k = -1; k <= 1; k++) { + struct rvec vec = new_rvec(i, j, k); + steps[count] = vec; + count++; + } + } + } + + struct rvec current = new_rvec(ttCell->x_ang, ttCell->y_ang, + ttCell->z_ang); + + double best = start; + count = 0; + + while (size > ANGLE_CONVERGE_SIZE && count < max_tries) + { + struct rvec sized[total]; + + int best_num = middle; + for (i = 0; i < total; i++) + { + struct rvec sized_step = steps[i]; + sized_step.u *= size; + sized_step.v *= size; + sized_step.w *= size; + + struct rvec new_angles = rvec_add_rvec(current, + sized_step); + + sized[i] = new_angles; + + ttCell->x_ang = sized[i].u; + ttCell->y_ang = sized[i].v; + ttCell->z_ang = sized[i].w; + + double score = obs_to_sol_score(ttCell); + + if (score < best) + { + best = score; + best_num = i; + } + } + + if (best_num == middle) + { + size /= 2; + } + + current = sized[best_num]; + count++; + } + + ttCell->x_ang = current.u; + ttCell->y_ang = current.v; + ttCell->z_ang = current.w; + + double end = obs_to_sol_score(ttCell); + + ttCell->x_ang = 0; + ttCell->y_ang = 0; + ttCell->z_ang = 0; + + gsl_matrix *tmp = gsl_matrix_alloc(3, 3); + rotate_gsl_by_angles(ttCell->solution, current.u, + current.v, current.w, tmp); + gsl_matrix_free(ttCell->solution); + ttCell->solution = tmp; +} + +static void match_all_obs_to_sol(struct TakeTwoCell *ttCell) +{ + int i, j; + double total = 0; + int count = 0; + gsl_matrix *solution = ttCell->solution; + + for (i = 0; i < ttCell->obs_vec_count; i++) + { + struct rvec *obs = &ttCell->obs_vecs[i].obsvec; + rvec_to_gsl(ttCell->vec1Tmp, *obs); + + /* ttCell->vec2Tmp = 1.0 * solution * ttCell->vec1Tmp */ + gsl_blas_dgemv(CblasTrans, 1.0, solution, ttCell->vec1Tmp, + 0.0, ttCell->vec2Tmp); + struct rvec rotated = gsl_to_rvec(ttCell->vec2Tmp); + + double smallest = FLT_MAX; + int assigned = -1; + + for (j = 0; j < ttCell->obs_vecs[i].match_num; j++) + { + struct rvec *match = &ttCell->obs_vecs[i].matches[j].vec; + struct rvec diff = diff_vec(rotated, *match); + + double length = rvec_length(diff); + if (length < smallest) + { + smallest = length; + assigned = j; + } + } + + ttCell->obs_vecs[i].assignment = assigned; + + if (smallest != FLT_MAX) + { + double addition = exp(-(1 / RECIP_TOLERANCE) * smallest); + total += addition; + count++; + + } + } + + total /= (double)count; +} static int match_obs_to_cell_vecs(struct TheoryVec *cell_vecs, int cell_vec_count, struct TakeTwoCell *cell) @@ -1585,6 +1797,7 @@ static int gen_observed_vecs(struct rvec *rlps, int rlp_count, spot_vec.obsvec = diff; spot_vec.distance = sqrt(sqlength); spot_vec.matches = NULL; + spot_vec.assignment = -1; spot_vec.match_num = 0; spot_vec.her_rlp = &rlps[i]; spot_vec.his_rlp = &rlps[j]; @@ -1746,6 +1959,10 @@ static UnitCell *run_taketwo(UnitCell *cell, const struct taketwo_options *opts, ttCell.vec2Tmp = gsl_vector_calloc(3); ttCell.numOps = 0; ttCell.obs_vec_count = 0; + ttCell.solution = NULL; + ttCell.x_ang = 0; + ttCell.y_ang = 0; + ttCell.z_ang = 0; success = generate_rotation_sym_ops(&ttCell); @@ -1787,14 +2004,20 @@ static UnitCell *run_taketwo(UnitCell *cell, const struct taketwo_options *opts, if ( !success ) return NULL; find_seeds(&ttCell, tp); -// remove_old_solutions(&ttCell, tp); + remove_old_solutions(&ttCell, tp); start_seeds(&solution, &ttCell); if ( solution == NULL ) { return NULL; } - for (int i = 0; i < tp->numPrevs; i++) + ttCell.solution = solution; + match_all_obs_to_sol(&ttCell); + refine_solution(&ttCell); + solution = ttCell.solution; + + int i; + for (i = 0; i < tp->numPrevs; i++) { gsl_matrix *sol = tp->prevSols[i]; -- cgit v1.2.3 From c32958e6bc6bb64d1a58bbd60ca6ce2ba801b548 Mon Sep 17 00:00:00 2001 From: Helen Ginn Date: Sun, 29 Apr 2018 18:55:04 +0200 Subject: Various comments throughout code --- libcrystfel/src/taketwo.c | 18 ++++++++++++------ 1 file changed, 12 insertions(+), 6 deletions(-) diff --git a/libcrystfel/src/taketwo.c b/libcrystfel/src/taketwo.c index 9d557fd2..7eb671ee 100644 --- a/libcrystfel/src/taketwo.c +++ b/libcrystfel/src/taketwo.c @@ -1037,6 +1037,7 @@ static signed int find_next_index(gsl_matrix *rot, int *obs_members, int all_ok = 1; int matched = -1; + /* Check all existing members are happy to let in the newcomer */ for ( j=0; j MAX_RECIP_DISTANCE) { continue; @@ -1331,9 +1333,9 @@ static int find_seeds(struct TakeTwoCell *cell, struct taketwo_private *tp) int seed_num = 0; struct Seed *seeds = NULL; - /* Check to see if any angles match from the cell vectors */ - obs_vecs_match_angles(i, j, - &seeds, &seed_num, cell); + /* Check to see if any angles match from the cell + * vectors */ + obs_vecs_match_angles(i, j, &seeds, &seed_num, cell); if (seed_num == 0) { @@ -1908,6 +1910,8 @@ static void cleanup_taketwo_obs_vecs(struct SpotVec *obs_vecs, static void cleanup_taketwo_cell(struct TakeTwoCell *ttCell) { + /* n.b. solutions in ttCell are taken care of in the + * partial taketwo cleanup. */ int i; for ( i=0; inumOps; i++ ) { gsl_matrix_free(ttCell->rotSymOps[i]); @@ -2003,6 +2007,8 @@ static UnitCell *run_taketwo(UnitCell *cell, const struct taketwo_options *opts, if ( !success ) return NULL; + /* Find all the seeds, then take each one and extend them, returning + * a solution if it exceeds the NETWORK_MEMBER_THRESHOLD. */ find_seeds(&ttCell, tp); remove_old_solutions(&ttCell, tp); start_seeds(&solution, &ttCell); @@ -2011,8 +2017,8 @@ static UnitCell *run_taketwo(UnitCell *cell, const struct taketwo_options *opts, return NULL; } + /* If we have a solution, refine against vectors in the entire image */ ttCell.solution = solution; - match_all_obs_to_sol(&ttCell); refine_solution(&ttCell); solution = ttCell.solution; -- cgit v1.2.3 From 563ca20a4c35fde586db7abd0dc72f6b074c3972 Mon Sep 17 00:00:00 2001 From: Helen Ginn Date: Sun, 29 Apr 2018 18:55:13 +0200 Subject: Change angle step size --- libcrystfel/src/taketwo.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/libcrystfel/src/taketwo.c b/libcrystfel/src/taketwo.c index 7eb671ee..9b58ef58 100644 --- a/libcrystfel/src/taketwo.c +++ b/libcrystfel/src/taketwo.c @@ -226,7 +226,7 @@ struct TakeTwoCell #define TRACE_TOLERANCE (deg2rad(3.0)) /* Initial step size for refinement of solutions */ -#define ANGLE_STEP_SIZE (deg2rad(0.2)) +#define ANGLE_STEP_SIZE (deg2rad(0.5)) /* Final required step size for refinement of solutions */ #define ANGLE_CONVERGE_SIZE (deg2rad(0.01)) -- cgit v1.2.3 From 5a212a4218af591778b1e52841d09a69bab30e1d Mon Sep 17 00:00:00 2001 From: Helen Ginn Date: Sun, 29 Apr 2018 18:55:36 +0200 Subject: Remove unnecessary checks which have been superseded by in_network --- libcrystfel/src/taketwo.c | 11 ----------- 1 file changed, 11 deletions(-) diff --git a/libcrystfel/src/taketwo.c b/libcrystfel/src/taketwo.c index 9b58ef58..d5355170 100644 --- a/libcrystfel/src/taketwo.c +++ b/libcrystfel/src/taketwo.c @@ -1023,17 +1023,6 @@ static signed int find_next_index(gsl_matrix *rot, int *obs_members, if ( !shared ) continue; - int skip = 0; - for ( j=0; j Date: Sun, 29 Apr 2018 18:55:44 +0200 Subject: Fussiness. --- libcrystfel/src/taketwo.c | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/libcrystfel/src/taketwo.c b/libcrystfel/src/taketwo.c index d5355170..75c06fbb 100644 --- a/libcrystfel/src/taketwo.c +++ b/libcrystfel/src/taketwo.c @@ -1031,7 +1031,8 @@ static signed int find_next_index(gsl_matrix *rot, int *obs_members, struct SpotVec *me = &obs_vecs[i]; struct SpotVec *you = &obs_vecs[obs_members[j]]; - struct rvec you_cell = you->matches[match_members[j]].vec; + struct rvec you_cell; + you_cell = you->matches[match_members[j]].vec; struct rvec me_obs = me->obsvec; struct rvec you_obs = you->obsvec; -- cgit v1.2.3 From 6885e4c27b0cf3eff6f7910542b328ce060e2191 Mon Sep 17 00:00:00 2001 From: Helen Ginn Date: Sun, 29 Apr 2018 18:55:56 +0200 Subject: Remove commented code --- libcrystfel/src/taketwo.c | 6 ------ 1 file changed, 6 deletions(-) diff --git a/libcrystfel/src/taketwo.c b/libcrystfel/src/taketwo.c index 75c06fbb..7d0c33f2 100644 --- a/libcrystfel/src/taketwo.c +++ b/libcrystfel/src/taketwo.c @@ -1085,12 +1085,6 @@ static signed int find_next_index(gsl_matrix *rot, int *obs_members, if (all_ok) { - - for ( j=0; j Date: Sun, 29 Apr 2018 18:56:28 +0200 Subject: Mass movement of functions to higher up --- libcrystfel/src/taketwo.c | 370 ++++++++++++++++++++++++---------------------- 1 file changed, 191 insertions(+), 179 deletions(-) diff --git a/libcrystfel/src/taketwo.c b/libcrystfel/src/taketwo.c index 7d0c33f2..d4e71705 100644 --- a/libcrystfel/src/taketwo.c +++ b/libcrystfel/src/taketwo.c @@ -1094,10 +1094,198 @@ static signed int find_next_index(gsl_matrix *rot, int *obs_members, return -1; } +/** + * Reward target function for refining solution to all vectors in a + * given image. Sum of exponentials of the negative distances, which + * means that the reward decays as the distance from the nearest + * theoretical vector reduces. */ +static double obs_to_sol_score(struct TakeTwoCell *ttCell) +{ + double total = 0; + int count = 0; + int i; + gsl_matrix *solution = ttCell->solution; + gsl_matrix *full_rot = gsl_matrix_alloc(3, 3); + rotate_gsl_by_angles(solution, ttCell->x_ang, ttCell->y_ang, + ttCell->z_ang, full_rot); + + for (i = 0; i < ttCell->obs_vec_count; i++) + { + struct rvec *obs = &ttCell->obs_vecs[i].obsvec; + rvec_to_gsl(ttCell->vec1Tmp, *obs); + + /* Rotate all the observed vectors by the modified soln */ + /* ttCell->vec2Tmp = 1.0 * full_rot * ttCell->vec1Tmp */ + gsl_blas_dgemv(CblasTrans, 1.0, full_rot, ttCell->vec1Tmp, + 0.0, ttCell->vec2Tmp); + struct rvec rotated = gsl_to_rvec(ttCell->vec2Tmp); + + int j = ttCell->obs_vecs[i].assignment; + + if (j < 0) continue; + + struct rvec *match = &ttCell->obs_vecs[i].matches[j].vec; + struct rvec diff = diff_vec(rotated, *match); + + double length = rvec_length(diff); + + double addition = exp(-(1 / RECIP_TOLERANCE) * length); + total += addition; + count++; + } + + total /= (double)-count; + + gsl_matrix_free(full_rot); + + return total; +} + +/** + * Matches every observed vector in the image to its closest theoretical + * neighbour after applying the rotation matrix, in preparation for + * refining the rotation matrix to match these. */ +static void match_all_obs_to_sol(struct TakeTwoCell *ttCell) +{ + int i, j; + double total = 0; + int count = 0; + gsl_matrix *solution = ttCell->solution; + + for (i = 0; i < ttCell->obs_vec_count; i++) + { + struct rvec *obs = &ttCell->obs_vecs[i].obsvec; + rvec_to_gsl(ttCell->vec1Tmp, *obs); + + /* ttCell->vec2Tmp = 1.0 * solution * ttCell->vec1Tmp */ + gsl_blas_dgemv(CblasTrans, 1.0, solution, ttCell->vec1Tmp, + 0.0, ttCell->vec2Tmp); + struct rvec rotated = gsl_to_rvec(ttCell->vec2Tmp); + + double smallest = FLT_MAX; + int assigned = -1; + + for (j = 0; j < ttCell->obs_vecs[i].match_num; j++) + { + struct rvec *match = &ttCell->obs_vecs[i].matches[j].vec; + struct rvec diff = diff_vec(rotated, *match); + + double length = rvec_length(diff); + if (length < smallest) + { + smallest = length; + assigned = j; + } + } + + ttCell->obs_vecs[i].assignment = assigned; + + if (smallest != FLT_MAX) + { + double addition = exp(-(1 / RECIP_TOLERANCE) * smallest); + total += addition; + count++; + + } + } + + total /= (double)count; +} + +/** + * Refines a matrix against all of the observed vectors against their + * closest theoretical neighbour, by perturbing the matrix along the principle + * axes until it maximises a reward function consisting of the sum of + * the distances of individual observed vectors to their closest + * theoretical neighbour. Reward function means that noise and alternative + * lattices do not dominate the equation! + **/ +static void refine_solution(struct TakeTwoCell *ttCell) +{ + match_all_obs_to_sol(ttCell); + + int i, j, k; + const int total = 3 * 3 * 3; + const int middle = (total - 1) / 2; + + struct rvec steps[total]; + double start = obs_to_sol_score(ttCell); + const int max_tries = 100; + + int count = 0; + double size = ANGLE_STEP_SIZE; + + /* First we create our combinations of steps */ + for (i = -1; i <= 1; i++) { + for (j = -1; j <= 1; j++) { + for (k = -1; k <= 1; k++) { + struct rvec vec = new_rvec(i, j, k); + steps[count] = vec; + count++; + } + } + } + + struct rvec current = new_rvec(ttCell->x_ang, ttCell->y_ang, + ttCell->z_ang); + + double best = start; + count = 0; + + while (size > ANGLE_CONVERGE_SIZE && count < max_tries) + { + struct rvec sized[total]; + + int best_num = middle; + for (i = 0; i < total; i++) + { + struct rvec sized_step = steps[i]; + sized_step.u *= size; + sized_step.v *= size; + sized_step.w *= size; + + struct rvec new_angles = rvec_add_rvec(current, + sized_step); + + sized[i] = new_angles; + + ttCell->x_ang = sized[i].u; + ttCell->y_ang = sized[i].v; + ttCell->z_ang = sized[i].w; + + double score = obs_to_sol_score(ttCell); + + if (score < best) + { + best = score; + best_num = i; + } + } + + if (best_num == middle) + { + size /= 2; + } + + current = sized[best_num]; + count++; + } + + ttCell->x_ang = 0; + ttCell->y_ang = 0; + ttCell->z_ang = 0; + + gsl_matrix *tmp = gsl_matrix_alloc(3, 3); + rotate_gsl_by_angles(ttCell->solution, current.u, + current.v, current.w, tmp); + gsl_matrix_free(ttCell->solution); + ttCell->solution = tmp; +} + -static int grow_network(gsl_matrix *rot, int obs_idx1, int obs_idx2, - int match_idx1, int match_idx2, int *max_members, - struct TakeTwoCell *cell) +static unsigned int grow_network(gsl_matrix *rot, int obs_idx1, int obs_idx2, + int match_idx1, int match_idx2, + struct TakeTwoCell *cell) { struct SpotVec *obs_vecs = cell->obs_vecs; @@ -1505,182 +1693,6 @@ static int sort_theory_distances(const void *av, const void *bv) return a->dist > b->dist; } -static double obs_to_sol_score(struct TakeTwoCell *ttCell) -{ - double total = 0; - int count = 0; - int i; - gsl_matrix *solution = ttCell->solution; - gsl_matrix *full_rot = gsl_matrix_alloc(3, 3); - rotate_gsl_by_angles(solution, ttCell->x_ang, ttCell->y_ang, - ttCell->z_ang, full_rot); - - for (i = 0; i < ttCell->obs_vec_count; i++) - { - struct rvec *obs = &ttCell->obs_vecs[i].obsvec; - rvec_to_gsl(ttCell->vec1Tmp, *obs); - - /* Rotate all the observed vectors by the modified soln */ - /* ttCell->vec2Tmp = 1.0 * full_rot * ttCell->vec1Tmp */ - gsl_blas_dgemv(CblasTrans, 1.0, full_rot, ttCell->vec1Tmp, - 0.0, ttCell->vec2Tmp); - struct rvec rotated = gsl_to_rvec(ttCell->vec2Tmp); - - int j = ttCell->obs_vecs[i].assignment; - - if (j < 0) continue; - - struct rvec *match = &ttCell->obs_vecs[i].matches[j].vec; - struct rvec diff = diff_vec(rotated, *match); - - double length = rvec_length(diff); - - double addition = exp(-(1 / RECIP_TOLERANCE) * length); - total += addition; - count++; - } - - total /= (double)-count; - - gsl_matrix_free(full_rot); - - return total; -} - -/* At the moment this is just going to do a step search of 3x3x3 until - * the angles are sufficiently tiny */ -static void refine_solution(struct TakeTwoCell *ttCell) -{ - int i, j, k; - const int total = 3 * 3 * 3; - const int middle = (total - 1) / 2; - - struct rvec steps[total]; - double start = obs_to_sol_score(ttCell); - const int max_tries = 100; - - int count = 0; - double size = ANGLE_STEP_SIZE; - - /* First we create our combinations of steps */ - for (i = -1; i <= 1; i++) { - for (j = -1; j <= 1; j++) { - for (k = -1; k <= 1; k++) { - struct rvec vec = new_rvec(i, j, k); - steps[count] = vec; - count++; - } - } - } - - struct rvec current = new_rvec(ttCell->x_ang, ttCell->y_ang, - ttCell->z_ang); - - double best = start; - count = 0; - - while (size > ANGLE_CONVERGE_SIZE && count < max_tries) - { - struct rvec sized[total]; - - int best_num = middle; - for (i = 0; i < total; i++) - { - struct rvec sized_step = steps[i]; - sized_step.u *= size; - sized_step.v *= size; - sized_step.w *= size; - - struct rvec new_angles = rvec_add_rvec(current, - sized_step); - - sized[i] = new_angles; - - ttCell->x_ang = sized[i].u; - ttCell->y_ang = sized[i].v; - ttCell->z_ang = sized[i].w; - - double score = obs_to_sol_score(ttCell); - - if (score < best) - { - best = score; - best_num = i; - } - } - - if (best_num == middle) - { - size /= 2; - } - - current = sized[best_num]; - count++; - } - - ttCell->x_ang = current.u; - ttCell->y_ang = current.v; - ttCell->z_ang = current.w; - - double end = obs_to_sol_score(ttCell); - - ttCell->x_ang = 0; - ttCell->y_ang = 0; - ttCell->z_ang = 0; - - gsl_matrix *tmp = gsl_matrix_alloc(3, 3); - rotate_gsl_by_angles(ttCell->solution, current.u, - current.v, current.w, tmp); - gsl_matrix_free(ttCell->solution); - ttCell->solution = tmp; -} - -static void match_all_obs_to_sol(struct TakeTwoCell *ttCell) -{ - int i, j; - double total = 0; - int count = 0; - gsl_matrix *solution = ttCell->solution; - - for (i = 0; i < ttCell->obs_vec_count; i++) - { - struct rvec *obs = &ttCell->obs_vecs[i].obsvec; - rvec_to_gsl(ttCell->vec1Tmp, *obs); - - /* ttCell->vec2Tmp = 1.0 * solution * ttCell->vec1Tmp */ - gsl_blas_dgemv(CblasTrans, 1.0, solution, ttCell->vec1Tmp, - 0.0, ttCell->vec2Tmp); - struct rvec rotated = gsl_to_rvec(ttCell->vec2Tmp); - - double smallest = FLT_MAX; - int assigned = -1; - - for (j = 0; j < ttCell->obs_vecs[i].match_num; j++) - { - struct rvec *match = &ttCell->obs_vecs[i].matches[j].vec; - struct rvec diff = diff_vec(rotated, *match); - - double length = rvec_length(diff); - if (length < smallest) - { - smallest = length; - assigned = j; - } - } - - ttCell->obs_vecs[i].assignment = assigned; - - if (smallest != FLT_MAX) - { - double addition = exp(-(1 / RECIP_TOLERANCE) * smallest); - total += addition; - count++; - - } - } - - total /= (double)count; -} static int match_obs_to_cell_vecs(struct TheoryVec *cell_vecs, int cell_vec_count, struct TakeTwoCell *cell) -- cgit v1.2.3 From d459477af06611f1b4920290f9084219aa7f3e47 Mon Sep 17 00:00:00 2001 From: Helen Ginn Date: Sun, 29 Apr 2018 18:57:17 +0200 Subject: Deal with membership number higher up, and add a minimum threshold (of 5) --- libcrystfel/src/taketwo.c | 46 ++++++++++++++++------------------------------ 1 file changed, 16 insertions(+), 30 deletions(-) diff --git a/libcrystfel/src/taketwo.c b/libcrystfel/src/taketwo.c index d4e71705..a081df52 100644 --- a/libcrystfel/src/taketwo.c +++ b/libcrystfel/src/taketwo.c @@ -212,6 +212,9 @@ struct TakeTwoCell /* Threshold for network members to consider a potential solution */ #define NETWORK_MEMBER_THRESHOLD (20) +/* Minimum for network members to consider a potential solution */ +#define MINIMUM_MEMBER_THRESHOLD (5) + /* Maximum dead ends for a single branch extension during indexing */ #define MAX_DEAD_ENDS (10) @@ -1314,7 +1317,6 @@ static unsigned int grow_network(gsl_matrix *rot, int obs_idx1, int obs_idx2, match_members[0] = match_idx1; match_members[1] = match_idx2; int member_num = 2; - *max_members = 2; /* counter for dead ends which must not exceed MAX_DEAD_ENDS * before it is reset in an additional branch */ @@ -1367,19 +1369,6 @@ static unsigned int grow_network(gsl_matrix *rot, int obs_idx1, int obs_idx2, member_num++; - if (member_num > *max_members) { - *max_members = member_num; - } - - /* - int n; - for (n = 0; n < member_num; n++) - { - STATUS("*"); - } - STATUS("\n"); - */ - /* If member_num is high enough, we want to return a yes */ if ( member_num > cell->member_thresh ) break; @@ -1391,13 +1380,12 @@ static unsigned int grow_network(gsl_matrix *rot, int obs_idx1, int obs_idx2, free(obs_members); free(match_members); - return ( member_num > cell->member_thresh ); + return ( member_num ); } -static int start_seed(int i, int j, int i_match, int j_match, - gsl_matrix **rotation, int *max_members, - struct TakeTwoCell *cell) +static unsigned int start_seed(int i, int j, int i_match, int j_match, + gsl_matrix **rotation, struct TakeTwoCell *cell) { struct SpotVec *obs_vecs = cell->obs_vecs; @@ -1411,13 +1399,13 @@ static int start_seed(int i, int j, int i_match, int j_match, /* Try to expand this rotation matrix to a larger network */ - int success = grow_network(rot_mat, i, j, i_match, j_match, max_members, - cell); + int member_num = grow_network(rot_mat, i, j, i_match, j_match, + cell); /* return this matrix and if it was immediately successful */ *rotation = rot_mat; - return success; + return member_num; } static int sort_seed_by_score(const void *av, const void *bv) @@ -1555,6 +1543,7 @@ static int start_seeds(gsl_matrix **rotation, struct TakeTwoCell *cell) { struct Seed *seeds = cell->seeds; int seed_num = cell->seed_count; + int member_num = 0; /* We have seeds! Pass each of them through the seed-starter */ /* If a seed has the highest achieved membership, make note...*/ @@ -1567,24 +1556,21 @@ static int start_seeds(gsl_matrix **rotation, struct TakeTwoCell *cell) continue; } - int max_members = 0; - int seed_obs1 = seeds[k].obs1; int seed_obs2 = seeds[k].obs2; - int success = start_seed(seed_obs1, seed_obs2, - seed_idx1, seed_idx2, - rotation, &max_members, - cell); + member_num = start_seed(seed_obs1, seed_obs2, seed_idx1, + seed_idx2, rotation, cell); - if (success) { + if (member_num >= NETWORK_MEMBER_THRESHOLD) { free(seeds); - return success; + return 1; } } free(seeds); - return (rotation != NULL); + + return (member_num > MINIMUM_MEMBER_THRESHOLD && rotation != NULL); } -- cgit v1.2.3 From 72dd667956dc3eecf88e73454c47a633d6b38c72 Mon Sep 17 00:00:00 2001 From: Helen Ginn Date: Sun, 29 Apr 2018 18:57:26 +0200 Subject: Fussiness. --- libcrystfel/src/taketwo.c | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/libcrystfel/src/taketwo.c b/libcrystfel/src/taketwo.c index a081df52..fe913901 100644 --- a/libcrystfel/src/taketwo.c +++ b/libcrystfel/src/taketwo.c @@ -1484,7 +1484,8 @@ static int find_seeds(struct TakeTwoCell *cell, struct taketwo_private *tp) * for optimisation by generating a look-up table * by spot instead of by vector. */ - int shared = obs_vecs_share_spot(&obs_vecs[i], &obs_vecs[j]); + int shared = obs_vecs_share_spot(&obs_vecs[i], + &obs_vecs[j]); if ( !shared ) continue; /* cell vector index matches stored in i, j and total @@ -1533,7 +1534,8 @@ static int find_seeds(struct TakeTwoCell *cell, struct taketwo_private *tp) } } - qsort(cell->seeds, cell->seed_count, sizeof(struct Seed), sort_seed_by_score); + qsort(cell->seeds, cell->seed_count, sizeof(struct Seed), + sort_seed_by_score); return 1; -- cgit v1.2.3 From f1261b157363ea8a19e39b11884ea4b7971c6684 Mon Sep 17 00:00:00 2001 From: Helen Ginn Date: Sun, 29 Apr 2018 19:00:55 +0200 Subject: white space --- libcrystfel/src/taketwo.c | 1 - 1 file changed, 1 deletion(-) diff --git a/libcrystfel/src/taketwo.c b/libcrystfel/src/taketwo.c index fe913901..cc68ed04 100644 --- a/libcrystfel/src/taketwo.c +++ b/libcrystfel/src/taketwo.c @@ -1681,7 +1681,6 @@ static int sort_theory_distances(const void *av, const void *bv) return a->dist > b->dist; } - static int match_obs_to_cell_vecs(struct TheoryVec *cell_vecs, int cell_vec_count, struct TakeTwoCell *cell) { -- cgit v1.2.3 From b8bca47b2b199c0205e0a6a30800f0991e477842 Mon Sep 17 00:00:00 2001 From: Helen Ginn Date: Sun, 29 Apr 2018 20:43:18 +0200 Subject: Generate theoretical vectors only once instead of every time the indexing is run --- libcrystfel/src/taketwo.c | 36 +++++++++++++++--------------------- 1 file changed, 15 insertions(+), 21 deletions(-) diff --git a/libcrystfel/src/taketwo.c b/libcrystfel/src/taketwo.c index cc68ed04..1d68ccc7 100644 --- a/libcrystfel/src/taketwo.c +++ b/libcrystfel/src/taketwo.c @@ -157,8 +157,11 @@ struct taketwo_private { IndexingMethod indm; UnitCell *cell; - int serial_num; /* -1 if unassigned */ - unsigned int attempts; + int serial_num; /**< Serial of last image, -1 if unassigned */ + unsigned int attempts; /**< Attempts at indexing current image */ + + struct TheoryVec *theory_vecs; /**< Theoretical vectors for given unit cell */ + unsigned int vec_count; /**< Number of theoretical vectors */ gsl_matrix **prevSols; /**< Previous solutions to be ignored */ unsigned int numPrevs; /**< Previous solution count */ @@ -1802,7 +1805,7 @@ static int gen_observed_vecs(struct rvec *rlps, int rlp_count, static int gen_theoretical_vecs(UnitCell *cell, struct TheoryVec **cell_vecs, - int *vec_count) + unsigned int *vec_count) { double a, b, c, alpha, beta, gamma; int h_max, k_max, l_max; @@ -1927,8 +1930,6 @@ static UnitCell *run_taketwo(UnitCell *cell, const struct taketwo_options *opts, struct rvec *rlps, int rlp_count, struct taketwo_private *tp) { - int cell_vec_count = 0; - struct TheoryVec *theory_vecs = NULL; UnitCell *result; int success = 0; gsl_matrix *solution = NULL; @@ -1953,7 +1954,6 @@ static UnitCell *run_taketwo(UnitCell *cell, const struct taketwo_options *opts, success = generate_rotation_sym_ops(&ttCell); - success = gen_theoretical_vecs(cell, &theory_vecs, &cell_vec_count); if ( !success ) return NULL; success = gen_observed_vecs(rlps, rlp_count, &ttCell); @@ -1983,11 +1983,9 @@ static UnitCell *run_taketwo(UnitCell *cell, const struct taketwo_options *opts, ttCell.trace_tol = sqrt(4.0*(1.0-cos(opts->trace_tol))); } - success = match_obs_to_cell_vecs(theory_vecs, cell_vec_count, + success = match_obs_to_cell_vecs(tp->theory_vecs, tp->vec_count, &ttCell); - free(theory_vecs); - if ( !success ) return NULL; /* Find all the seeds, then take each one and extend them, returning @@ -2005,18 +2003,7 @@ static UnitCell *run_taketwo(UnitCell *cell, const struct taketwo_options *opts, refine_solution(&ttCell); solution = ttCell.solution; - int i; - for (i = 0; i < tp->numPrevs; i++) - { - gsl_matrix *sol = tp->prevSols[i]; - - int sim = symm_rot_mats_are_similar(sol, solution, &ttCell); - if (sim) - { -// STATUS("Warning! Returning previous solution.\n"); - } - } - + /* Add the current solution to the previous solutions list */ int new_size = (tp->numPrevs + 1) * sizeof(gsl_matrix *); gsl_matrix **tmp = realloc(tp->prevSols, new_size); @@ -2029,6 +2016,7 @@ static UnitCell *run_taketwo(UnitCell *cell, const struct taketwo_options *opts, tp->prevSols[tp->numPrevs] = solution; tp->numPrevs++; + /* Prepare the solution for CrystFEL friendliness */ result = transform_cell_gsl(cell, solution); cleanup_taketwo_cell(&ttCell); @@ -2165,6 +2153,10 @@ void *taketwo_prepare(IndexingMethod *indm, UnitCell *cell) tp->attempts = 0; tp->prevSols = NULL; tp->numPrevs = 0; + tp->vec_count = 0; + tp->theory_vecs = NULL; + + gen_theoretical_vecs(cell, &tp->theory_vecs, &tp->vec_count); return tp; } @@ -2174,6 +2166,8 @@ void taketwo_cleanup(IndexingPrivate *pp) struct taketwo_private *tp = (struct taketwo_private *)pp; partial_taketwo_cleanup(tp); + free(tp->theory_vecs); + free(tp); } -- cgit v1.2.3 From 9aeb19db674f211aff6847104135798620f9f150 Mon Sep 17 00:00:00 2001 From: Helen Ginn Date: Sun, 29 Apr 2018 20:53:53 +0200 Subject: Plug leaks --- libcrystfel/src/taketwo.c | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/libcrystfel/src/taketwo.c b/libcrystfel/src/taketwo.c index 1d68ccc7..480fd76a 100644 --- a/libcrystfel/src/taketwo.c +++ b/libcrystfel/src/taketwo.c @@ -396,6 +396,12 @@ static void rotate_gsl_by_angles(gsl_matrix *sol, double x, double y, /* Apply the whole rotation offset to the solution */ gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, xyz_rot, sol, 0.0, result); + + gsl_matrix_free(x_rot); + gsl_matrix_free(y_rot); + gsl_matrix_free(z_rot); + gsl_matrix_free(xy_rot); + gsl_matrix_free(xyz_rot); } -- cgit v1.2.3 From 9a12096958cd9e47116bfac8ffa5de8f523ed81a Mon Sep 17 00:00:00 2001 From: Helen Ginn Date: Sun, 29 Apr 2018 20:54:27 +0200 Subject: Plug leaks. --- libcrystfel/src/taketwo.c | 3 --- 1 file changed, 3 deletions(-) diff --git a/libcrystfel/src/taketwo.c b/libcrystfel/src/taketwo.c index 480fd76a..0b533f0f 100644 --- a/libcrystfel/src/taketwo.c +++ b/libcrystfel/src/taketwo.c @@ -679,9 +679,6 @@ static gsl_matrix *generate_rot_mat(struct rvec obs1, struct rvec obs2, gsl_matrix *fullMat; rvec_to_gsl(cell->vec1Tmp, cell2); -// gsl_vector *cell2v = rvec_to_gsl(cell2); -// gsl_vector *cell2vr = gsl_vector_calloc(3); - normalise_rvec(&obs1); normalise_rvec(&obs2); normalise_rvec(&cell1); -- cgit v1.2.3 From 5ca551032e2cc7e7391cd41be0599277363b58b6 Mon Sep 17 00:00:00 2001 From: Helen Ginn Date: Mon, 30 Apr 2018 10:48:00 +0200 Subject: Limit distances of seeds only if cell seed count is quite high --- libcrystfel/src/taketwo.c | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/libcrystfel/src/taketwo.c b/libcrystfel/src/taketwo.c index 0b533f0f..398e6fcb 100644 --- a/libcrystfel/src/taketwo.c +++ b/libcrystfel/src/taketwo.c @@ -1480,9 +1480,9 @@ static int find_seeds(struct TakeTwoCell *cell, struct taketwo_private *tp) for ( j=0; j - MAX_RECIP_DISTANCE) { + MAX_RECIP_DISTANCE && cell->seed_count > 100) { continue; } -- cgit v1.2.3 From 71a60180b43c64f0f426a8ec181b83cbd958713e Mon Sep 17 00:00:00 2001 From: Helen Ginn Date: Mon, 30 Apr 2018 10:48:37 +0200 Subject: Correctly return the best matrix, not the last one --- libcrystfel/src/taketwo.c | 16 ++++++++++++---- 1 file changed, 12 insertions(+), 4 deletions(-) diff --git a/libcrystfel/src/taketwo.c b/libcrystfel/src/taketwo.c index 398e6fcb..1d6ef6a2 100644 --- a/libcrystfel/src/taketwo.c +++ b/libcrystfel/src/taketwo.c @@ -1547,11 +1547,13 @@ static int find_seeds(struct TakeTwoCell *cell, struct taketwo_private *tp) return 1; } -static int start_seeds(gsl_matrix **rotation, struct TakeTwoCell *cell) +static unsigned int start_seeds(gsl_matrix **rotation, struct TakeTwoCell *cell) { struct Seed *seeds = cell->seeds; int seed_num = cell->seed_count; int member_num = 0; + int max_members = 0; + gsl_matrix *rot = NULL; /* We have seeds! Pass each of them through the seed-starter */ /* If a seed has the highest achieved membership, make note...*/ @@ -1568,17 +1570,23 @@ static int start_seeds(gsl_matrix **rotation, struct TakeTwoCell *cell) int seed_obs2 = seeds[k].obs2; member_num = start_seed(seed_obs1, seed_obs2, seed_idx1, - seed_idx2, rotation, cell); + seed_idx2, &rot, cell); + + if (member_num > max_members) + { + *rotation = rot; + max_members = member_num; + } if (member_num >= NETWORK_MEMBER_THRESHOLD) { free(seeds); - return 1; + return max_members; } } free(seeds); - return (member_num > MINIMUM_MEMBER_THRESHOLD && rotation != NULL); + return max_members; } -- cgit v1.2.3 From ff0e331ec93164470bb7bdc858db148aefb4e67f Mon Sep 17 00:00:00 2001 From: Helen Ginn Date: Mon, 30 Apr 2018 10:49:30 +0200 Subject: Track solution network membership and scores for diagnostics --- libcrystfel/src/taketwo.c | 38 +++++++++++++++++++++++++++++++++----- 1 file changed, 33 insertions(+), 5 deletions(-) diff --git a/libcrystfel/src/taketwo.c b/libcrystfel/src/taketwo.c index 1d6ef6a2..234b00c3 100644 --- a/libcrystfel/src/taketwo.c +++ b/libcrystfel/src/taketwo.c @@ -158,13 +158,15 @@ struct taketwo_private IndexingMethod indm; UnitCell *cell; int serial_num; /**< Serial of last image, -1 if unassigned */ - unsigned int attempts; /**< Attempts at indexing current image */ + unsigned int xtal_num; /**< last number of crystals recorded */ struct TheoryVec *theory_vecs; /**< Theoretical vectors for given unit cell */ unsigned int vec_count; /**< Number of theoretical vectors */ gsl_matrix **prevSols; /**< Previous solutions to be ignored */ unsigned int numPrevs; /**< Previous solution count */ + double *prevScores; /**< previous solution scores */ + unsigned int *membership; /**< previous solution was success or failure */ }; @@ -2003,7 +2005,7 @@ static UnitCell *run_taketwo(UnitCell *cell, const struct taketwo_options *opts, * a solution if it exceeds the NETWORK_MEMBER_THRESHOLD. */ find_seeds(&ttCell, tp); remove_old_solutions(&ttCell, tp); - start_seeds(&solution, &ttCell); + unsigned int members = start_seeds(&solution, &ttCell); if ( solution == NULL ) { return NULL; @@ -2013,18 +2015,28 @@ static UnitCell *run_taketwo(UnitCell *cell, const struct taketwo_options *opts, ttCell.solution = solution; refine_solution(&ttCell); solution = ttCell.solution; + double score = obs_to_sol_score(&ttCell); /* Add the current solution to the previous solutions list */ int new_size = (tp->numPrevs + 1) * sizeof(gsl_matrix *); gsl_matrix **tmp = realloc(tp->prevSols, new_size); + double *tmpScores = realloc(tp->prevScores, + (tp->numPrevs + 1) * sizeof(double)); + unsigned int *tmpSuccesses; + tmpSuccesses = realloc(tp->membership, + (tp->numPrevs + 1) * sizeof(unsigned int)); if (!tmp) { apologise(); } tp->prevSols = tmp; + tp->prevScores = tmpScores; + tp->membership = tmpSuccesses; tp->prevSols[tp->numPrevs] = solution; + tp->prevScores[tp->numPrevs] = score; + tp->membership[tp->numPrevs] = members; tp->numPrevs++; /* Prepare the solution for CrystFEL friendliness */ @@ -2048,7 +2060,11 @@ static void partial_taketwo_cleanup(struct taketwo_private *tp) free(tp->prevSols); } - tp->attempts = 0; + free(tp->prevScores); + free(tp->membership); + tp->prevScores = NULL; + tp->membership = NULL; + tp->xtal_num = 0; tp->numPrevs = 0; tp->prevSols = NULL; @@ -2071,12 +2087,22 @@ int taketwo_index(struct image *image, const struct taketwo_options *opts, if (tp->serial_num == this_serial) { - tp->attempts++; + tp->xtal_num = image->n_crystals; } else { + /* + for (i = 0; i < tp->numPrevs; i++) + { + STATUS("score, %i, %.5f, %i\n", + this_serial, tp->prevScores[i], + tp->membership[i]); + } + */ + partial_taketwo_cleanup(tp); tp->serial_num = this_serial; + tp->xtal_num = image->n_crystals; } /* @@ -2161,9 +2187,11 @@ void *taketwo_prepare(IndexingMethod *indm, UnitCell *cell) tp->cell = cell; tp->indm = *indm; tp->serial_num = -1; - tp->attempts = 0; + tp->xtal_num = 0; tp->prevSols = NULL; tp->numPrevs = 0; + tp->prevScores = NULL; + tp->membership = NULL; tp->vec_count = 0; tp->theory_vecs = NULL; -- cgit v1.2.3