aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2018-05-04 15:21:59 +0200
committerThomas White <taw@physics.org>2018-05-04 15:21:59 +0200
commit61c158de4db8fa30ad51cbab9ee44950a182b091 (patch)
tree1496d63a644d346c7470d16fcfc6c14c4dc6a471
parented8e1c9751d7839aa6cd74422f93da513840dcf4 (diff)
parentff0e331ec93164470bb7bdc858db148aefb4e67f (diff)
Merge branch 'taketwo'
-rw-r--r--libcrystfel/src/taketwo.c1116
1 files changed, 842 insertions, 274 deletions
diff --git a/libcrystfel/src/taketwo.c b/libcrystfel/src/taketwo.c
index 0857f813..234b00c3 100644
--- a/libcrystfel/src/taketwo.c
+++ b/libcrystfel/src/taketwo.c
@@ -28,6 +28,69 @@
*
*/
+/**
+ * \class TakeTwo
+ * Code outline.
+ * --- Get ready for calculation ---
+ * 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())
+ * Match observed vectors to theoretical vectors (match_obs_to_cell_vecs())
+ *
+ * --- Business bit ---
+ * ... 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
+ * - 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 <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>
#include <float.h>
@@ -57,40 +120,90 @@
struct SpotVec
{
struct rvec obsvec;
- struct rvec *matches;
+ struct TheoryVec *matches;
int match_num;
- struct rvec *asym_matches;
- int asym_match_num;
+ int assignment;
+ int in_network;
double distance;
struct rvec *her_rlp;
struct rvec *his_rlp;
};
+/**
+ * theoryvec
+ */
+
+struct TheoryVec
+{
+ struct rvec vec;
+ int asym;
+};
+
+
+/**
+ * seed
+ */
+
+struct Seed
+{
+ int obs1;
+ int obs2;
+ int idx1;
+ int idx2;
+ double score;
+};
struct taketwo_private
{
IndexingMethod indm;
UnitCell *cell;
+ int serial_num; /**< Serial of last image, -1 if unassigned */
+ 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 */
+
};
-// 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;
+ struct Seed *seeds;
+ int seed_count;
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))) */
+
+ /** 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 */
gsl_matrix *twiz2Tmp;
+ /**< Temporary memory always allocated for calculations */
gsl_vector *vec1Tmp;
+ /**< Temporary memory always allocated for calculations */
gsl_vector *vec2Tmp;
};
@@ -104,6 +217,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)
@@ -117,6 +233,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.5))
+
+/* Final required step size for refinement of solutions */
+#define ANGLE_CONVERGE_SIZE (deg2rad(0.01))
+
/* TODO: Multiple lattices */
@@ -145,6 +267,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)
{
@@ -240,6 +371,41 @@ 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);
+
+ 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);
+}
+
/* Rotate vector (vec1) around axis (axis) by angle theta. Find value of
* theta for which the angle between (vec1) and (vec2) is minimised. */
@@ -297,6 +463,34 @@ 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)
{
@@ -487,9 +681,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);
@@ -553,96 +744,118 @@ 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)
+static int obs_vecs_match_angles(int her, int his,
+ struct Seed **seeds, int *match_count,
+ struct TakeTwoCell *cell)
{
- int i, j;
- *match_count = 0;
+ 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);
+ 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 */
+
+ int i, j;
+ for ( i=0; i<her_obs->match_num; i++ ) {
+ for ( j=0; j<his_obs->match_num; j++ ) {
+ double score = 0;
- /* calculate angle between all potential theoretical vectors */
+ struct rvec *her_match = &her_obs->matches[i].vec;
+ struct rvec *his_match = &his_obs->matches[j].vec;
- for ( i=0; i<her_obs->match_num; i++ ) {
- for ( j=0; j<his_obs->match_num; j++ ) {
+ double her_dist = rvec_length(*her_match);
+ double his_dist = rvec_length(*his_match);
- struct rvec *her_match = &her_obs->matches[i];
- struct rvec *his_match = &his_obs->matches[j];
+ double theory_angle = rvec_angle(*her_match,
+ *his_match);
- double theory_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);
+ /* within basic tolerance? */
+ if ( angle_diff != angle_diff ||
+ angle_diff > cell->angle_tol ) {
+ continue;
+ }
+
+ double add = angle_diff;
+ if (add == add) {
+ score += add * her_dist * his_dist;
+ }
- if ( angle_diff < cell->angle_tol ) {
+ /* If the angles are too close to 0
+ or 180, one axis ill-determined */
+ if (theory_angle < min_angle ||
+ theory_angle > max_angle) {
+ continue;
+ }
- // in the case of a brief check only
- if (!her_match_idxs || !his_match_idxs) {
- return 1;
- }
+ /* check that third vector adequately completes
+ * triangle */
- /* 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 the third vector
+ theory_angle = rvec_angle(*her_match,
+ theory_diff);
+ obs_angle = rvec_angle(her_obs->obsvec, obs_diff);
+ angle_diff = fabs(obs_angle - theory_angle);
- struct rvec theory_diff = diff_vec(*his_match, *her_match);
- struct rvec obs_diff = diff_vec(his_obs->obsvec,
- her_obs->obsvec);
+ double diff_dist = rvec_length(obs_diff);
- theory_angle = rvec_angle(*her_match,
- theory_diff);
- obs_angle = rvec_angle(her_obs->obsvec, obs_diff);
+ if (angle_diff > ANGLE_TOLERANCE) {
+ continue;
+ }
- if (fabs(obs_angle - theory_angle) > ANGLE_TOLERANCE) {
- continue;
- }
+ add = angle_diff;
+ if (add == add) {
+ score += add * her_dist * diff_dist;
+ }
- theory_angle = rvec_angle(*his_match,
- theory_diff);
- obs_angle = rvec_angle(his_obs->obsvec, obs_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;
- }
+ if (fabs(obs_angle - theory_angle) > ANGLE_TOLERANCE) {
+ continue;
+ }
- size_t new_size = (*match_count + 1) *
- sizeof(int);
- if (her_match_idxs && his_match_idxs)
- {
- /* Reallocate the array to fit in another match */
- int *temp_hers;
- int *temp_his;
- temp_hers = realloc(*her_match_idxs, new_size);
- temp_his = realloc(*his_match_idxs, new_size);
+ add = angle_diff;
+ if (add == add) {
+ score += add * his_dist * diff_dist;
+ }
- if ( temp_hers == NULL || temp_his == NULL ) {
- apologise();
- }
+ /* we add a new seed to the array */
+ size_t new_size = (*match_count + 1);
+ new_size *= sizeof(struct Seed);
- (*her_match_idxs) = temp_hers;
- (*his_match_idxs) = temp_his;
+ /* Reallocate the array to fit in another match */
+ struct Seed *tmp_seeds = realloc(*seeds, new_size);
- (*her_match_idxs)[*match_count] = i;
- (*his_match_idxs)[*match_count] = j;
- }
+ if ( tmp_seeds == NULL ) {
+ apologise();
+ }
+
+ (*seeds) = tmp_seeds;
+
+ (*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 * 1000;
(*match_count)++;
}
}
- }
return (*match_count > 0);
}
@@ -672,8 +885,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,
@@ -716,10 +929,26 @@ static signed int finish_solution(gsl_matrix *rot, struct SpotVec *obs_vecs,
return 1;
}
-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)
+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;
+ 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 Seed **seeds,
+ int *match_count, struct TakeTwoCell *cell)
{
int num_occupied = 0;
gsl_matrix **old_mats = calloc(*match_count, sizeof(gsl_matrix *));
@@ -731,19 +960,18 @@ static int weed_duplicate_matches(struct SpotVec *her_obs,
}
signed int i, j;
+
int duplicates = 0;
- for (i = *match_count - 1; i >= 0; i--) {
- int her_match = (*her_match_idxs)[i];
- int his_match = (*his_match_idxs)[i];
+ /* Now we weed out the self-duplicates from the remaining batch */
- 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];
+ for (i = *match_count - 1; i >= 0; i--) {
+ int her_match = (*seeds)[i].idx1;
+ int his_match = (*seeds)[i].idx2;
- gsl_matrix *mat = generate_rot_mat(i_obsvec, j_obsvec,
- i_cellvec, j_cellvec, cell);
+ gsl_matrix *mat;
+ mat = rot_mat_from_indices((*seeds)[i].obs1, (*seeds)[i].obs2,
+ her_match, his_match, cell);
int found = 0;
@@ -752,8 +980,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++;
@@ -785,7 +1013,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);
@@ -794,42 +1022,42 @@ static signed int find_next_index(gsl_matrix *rot, int *obs_members,
for ( i=start; i<obs_vec_count; i++ ) {
+ /* If we've considered this vector before, ignore it */
+ if (obs_vecs[i].in_network == 1)
+ {
+ continue;
+ }
+
/* first we check for a shared spot - harshest condition */
int shared = obs_shares_spot_w_array(obs_vecs, i, obs_members,
member_num);
if ( !shared ) continue;
- int skip = 0;
- for ( j=0; j<member_num && skip == 0; j++ ) {
- if (i == obs_members[j]) {
- skip = 1;
- }
- }
-
- if (skip) {
- continue;
- }
-
int all_ok = 1;
int matched = -1;
+ /* Check all existing members are happy to let in the newcomer */
for ( j=0; j<member_num && all_ok; j++ ) {
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_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; k<me->match_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,
@@ -844,6 +1072,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) {
@@ -863,12 +1096,6 @@ static signed int find_next_index(gsl_matrix *rot, int *obs_members,
if (all_ok) {
-
- for ( j=0; j<member_num; j++ ) {
- // STATUS("%i ", obs_members[j]);
- }
- //STATUS("%i\n", i);
-
return i;
}
}
@@ -878,16 +1105,212 @@ 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++;
-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)
+ }
+ }
+
+ 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)
{
- struct SpotVec *obs_vecs = *(cell->obs_vecs);
+ 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 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;
int obs_vec_count = cell->obs_vec_count;
int *obs_members;
int *match_members;
+ /* Clear the in_network status of all vectors to start */
+ int i;
+ for (i = 0; i < obs_vec_count; i++)
+ {
+ obs_vecs[i].in_network = 0;
+ }
+
/* indices of members of the self-consistent network of vectors */
obs_members = malloc((cell->member_thresh+3)*sizeof(int));
match_members = malloc((cell->member_thresh+3)*sizeof(int));
@@ -902,7 +1325,6 @@ static 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 */
@@ -923,7 +1345,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 ) {
@@ -942,25 +1364,19 @@ 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;
}
- /* we have elongated membership - so reset dead_ends counter */
- // dead_ends = 0;
-
+ /* Elongation of the network was successful */
+ obs_vecs[next_index].in_network = 1;
obs_members[member_num] = next_index;
match_members[member_num] = match_found;
member_num++;
- if (member_num > *max_members) {
- *max_members = member_num;
- }
-
/* If member_num is high enough, we want to return a yes */
if ( member_num > cell->member_thresh ) break;
@@ -972,126 +1388,210 @@ static 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);
+ struct SpotVec *obs_vecs = cell->obs_vecs;
gsl_matrix *rot_mat;
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 */
- 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 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;
+}
+
+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;
+ }
+ }
- /* META: Maximum achieved maximum network membership */
- int max_max_members = 0;
- gsl_matrix *best_rotation = NULL;
+ gsl_matrix_free(mat);
+ }
+
+// STATUS("Removing %i duplicates due to prev solutions.\n", duplicates);
+}
-// unsigned long start_time = time(NULL);
+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;
/* loop round pairs of vectors to try and find a suitable
* seed to start building a self-consistent network of vectors
*/
int i, j;
- for ( i=0; i<obs_vec_count; i++ ) {
+ for ( i=1; i<obs_vec_count; i++ ) {
+
for ( j=0; j<i; j++ ) {
+ /** Only check distances which are accumulatively less
+ * than the limit if we can easily generate seeds */
+ if (obs_vecs[j].distance + obs_vecs[i].distance >
+ MAX_RECIP_DISTANCE && cell->seed_count > 100) {
+ 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]);
-
+ 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
* number stored in int matches.
*/
- int *i_idx = 0;
- int *j_idx = 0;
- int matches;
-
- /* 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);
+ 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);
+
+ if (seed_num == 0)
+ {
+ /* Nothing to clean up here */
continue;
}
/* Weed out the duplicate seeds (from symmetric
- * reflection pairs)
- */
+ * reflection pairs) */
+ weed_duplicate_matches(&seeds, &seed_num, cell);
+
+ /* Add all the new seeds to the full list */
+
+ size_t new_size = cell->seed_count + seed_num;
+ new_size *= sizeof(struct Seed);
+
+ struct Seed *tmp = realloc(cell->seeds, new_size);
- weed_duplicate_matches(&obs_vecs[i], &obs_vecs[j],
- &i_idx, &j_idx, &matches, cell);
+ if (tmp == NULL) {
+ apologise();
+ }
+
+ cell->seeds = tmp;
- /* 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<matches; k++ ) {
- if (i_idx[k] < 0 || j_idx[k] < 0) {
+ 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,
- i_idx[k], j_idx[k],
- rotation, &max_members,
- cell);
-
- if (success) {
- free(i_idx); free(j_idx);
- 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++;
}
+ }
+ }
+
+ qsort(cell->seeds, cell->seed_count, sizeof(struct Seed),
+ sort_seed_by_score);
+
+
+ return 1;
+}
+
+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...*/
+ int k;
+ for ( k=0; k<seed_num; k++ ) {
+ int seed_idx1 = seeds[k].idx1;
+ int seed_idx2 = seeds[k].idx2;
+
+ if (seed_idx1 < 0 || seed_idx2 < 0) {
+ continue;
+ }
+
+ int seed_obs1 = seeds[k].obs1;
+ int seed_obs2 = seeds[k].obs2;
- free(i_idx);
- free(j_idx);
+ member_num = start_seed(seed_obs1, seed_obs2, seed_idx1,
+ seed_idx2, &rot, cell);
+
+ if (member_num > max_members)
+ {
+ *rotation = rot;
+ max_members = member_num;
+ }
+
+ if (member_num >= NETWORK_MEMBER_THRESHOLD) {
+ free(seeds);
+ return max_members;
}
- } /* yes this } is meant to be here */
+ }
+
+ free(seeds);
- *rotation = best_rotation;
- return (best_rotation != NULL);
+ return max_members;
}
+
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)
@@ -1184,25 +1684,24 @@ static int generate_rotation_sym_ops(struct TakeTwoCell *ttCell)
return 1;
}
-
struct sortme
{
- struct rvec v;
+ struct TheoryVec v;
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;
return a->dist > b->dist;
}
-
-static int match_obs_to_cell_vecs(struct rvec *cell_vecs, int cell_vec_count,
- struct SpotVec *obs_vecs, int obs_vec_count,
- int is_asymmetric, struct TakeTwoCell *cell)
+static int match_obs_to_cell_vecs(struct TheoryVec *cell_vecs, int cell_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; i<obs_vec_count; i++ ) {
@@ -1212,7 +1711,7 @@ static int match_obs_to_cell_vecs(struct rvec *cell_vecs, int cell_vec_count,
for ( j=0; j<cell_vec_count; j++ ) {
/* get distance for unit cell vector */
- double cell_length = rvec_length(cell_vecs[j]);
+ double cell_length = rvec_length(cell_vecs[j].vec);
double obs_length = obs_vecs[i].distance;
/* check if this matches the observed length */
@@ -1235,20 +1734,15 @@ static int match_obs_to_cell_vecs(struct rvec *cell_vecs, int cell_vec_count,
/* Pointers to relevant things */
- struct rvec **match_array;
+ struct TheoryVec **match_array;
int *match_count;
- if (!is_asymmetric) {
- match_array = &(obs_vecs[i].matches);
- match_count = &(obs_vecs[i].match_num);
- } else {
- match_array = &(obs_vecs[i].asym_matches);
- match_count = &(obs_vecs[i].asym_match_num);
- }
+ match_array = &(obs_vecs[i].matches);
+ match_count = &(obs_vecs[i].match_num);
/* Sort in order to get most agreeable matches first */
- qsort(for_sort, count, sizeof(struct sortme), sort_func);
- *match_array = malloc(count*sizeof(struct rvec));
+ qsort(for_sort, count, sizeof(struct sortme), sort_theory_distances);
+ *match_array = malloc(count*sizeof(struct TheoryVec));
*match_count = count;
for ( j=0; j<count; j++ ) {
(*match_array)[j] = for_sort[j].v;
@@ -1268,7 +1762,7 @@ static int compare_spot_vecs(const void *av, const void *bv)
}
static int gen_observed_vecs(struct rvec *rlps, int rlp_count,
- struct SpotVec **obs_vecs, int *obs_vec_count)
+ struct TakeTwoCell *cell)
{
int i, j;
int count = 0;
@@ -1279,7 +1773,6 @@ static int gen_observed_vecs(struct rvec *rlps, int rlp_count,
for ( i=0; i<rlp_count-1 && count < MAX_OBS_VECTORS; i++ ) {
for ( j=i+1; j<rlp_count; j++ ) {
-
/* calculate difference vector between rlps */
struct rvec diff = diff_vec(rlps[i], rlps[j]);
@@ -1291,13 +1784,13 @@ static int gen_observed_vecs(struct rvec *rlps, int rlp_count,
count++;
struct SpotVec *temp_obs_vecs;
- temp_obs_vecs = realloc(*obs_vecs,
+ temp_obs_vecs = realloc(cell->obs_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 */
@@ -1305,29 +1798,27 @@ 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];
- (*obs_vecs)[count - 1] = spot_vec;
+ cell->obs_vecs[count - 1] = spot_vec;
}
}
}
- /* Sort such that the shortest and least error-prone distances
- are searched first.
- */
- qsort(*obs_vecs, count, sizeof(struct SpotVec), compare_spot_vecs);
+ /* Sort such that the shortest distances are searched first. */
+ qsort(cell->obs_vecs, count, sizeof(struct SpotVec), compare_spot_vecs);
- *obs_vec_count = count;
+ cell->obs_vec_count = count;
return 1;
}
-static int gen_theoretical_vecs(UnitCell *cell, struct rvec **cell_vecs,
- struct rvec **asym_vecs, int *vec_count,
- int *asym_vec_count)
+static int gen_theoretical_vecs(UnitCell *cell, struct TheoryVec **cell_vecs,
+ unsigned int *vec_count)
{
double a, b, c, alpha, beta, gamma;
int h_max, k_max, l_max;
@@ -1351,7 +1842,6 @@ static int gen_theoretical_vecs(UnitCell *cell, struct rvec **cell_vecs,
int h, k, l;
int _h, _k, _l;
int count = 0;
- int asym_count = 0;
for ( h=-h_max; h<=+h_max; h++ ) {
for ( k=-k_max; k<=+k_max; k++ ) {
@@ -1367,50 +1857,37 @@ static int gen_theoretical_vecs(UnitCell *cell, struct rvec **cell_vecs,
if (h == _h && k == _k && l == _l) {
asymmetric = 1;
- asym_count++;
}
cell_vec.u = h*asx + k*bsx + l*csx;
cell_vec.v = h*asy + k*bsy + l*csy;
cell_vec.w = h*asz + k*bsz + l*csz;
+ struct TheoryVec theory;
+ theory.vec = cell_vec;
+ theory.asym = asymmetric;
+
/* add this to our array - which may require expanding */
count++;
- struct rvec *temp_cell_vecs;
- temp_cell_vecs = realloc(*cell_vecs, count*sizeof(struct rvec));
- struct rvec *temp_asym_vecs = NULL;
-
- if (asymmetric) {
- temp_asym_vecs = realloc(*asym_vecs,
- count*sizeof(struct rvec));
- }
+ struct TheoryVec *temp_cell_vecs;
+ temp_cell_vecs = realloc(*cell_vecs,
+ count*sizeof(struct TheoryVec));
if ( temp_cell_vecs == NULL ) {
return 0;
- } else if (asymmetric && temp_asym_vecs == NULL) {
- return 0;
} else {
*cell_vecs = temp_cell_vecs;
- (*cell_vecs)[count - 1] = cell_vec;
-
- if (!asymmetric) {
- continue;
- }
-
- *asym_vecs = temp_asym_vecs;
- (*asym_vecs)[asym_count - 1] = cell_vec;
+ (*cell_vecs)[count - 1] = theory;
}
}
}
}
*vec_count = count;
- *asym_vec_count = asym_count;
free_symoplist(rawList);
-
return 1;
}
@@ -1425,7 +1902,6 @@ static void cleanup_taketwo_obs_vecs(struct SpotVec *obs_vecs,
int i;
for ( i=0; i<obs_vec_count; i++ ) {
free(obs_vecs[i].matches);
- free(obs_vecs[i].asym_matches);
}
free(obs_vecs);
@@ -1433,11 +1909,16 @@ 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; i<ttCell->numOps; i++ ) {
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);
@@ -1459,39 +1940,38 @@ 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;
- int asym_vec_count = 0;
- struct rvec *cell_vecs = NULL;
- struct rvec *asym_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.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.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);
- success = gen_theoretical_vecs(cell, &cell_vecs, &asym_vecs,
- &cell_vec_count, &asym_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 {
@@ -1516,32 +1996,79 @@ 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(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);
-
- free(cell_vecs);
- free(asym_vecs);
+ success = match_obs_to_cell_vecs(tp->theory_vecs, tp->vec_count,
+ &ttCell);
if ( !success ) return NULL;
- find_seed(&solution, &ttCell);
+ /* 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);
+ unsigned int members = start_seeds(&solution, &ttCell);
if ( solution == NULL ) {
- cleanup_taketwo_obs_vecs(obs_vecs, obs_vec_count);
return NULL;
}
+ /* If we have a solution, refine against vectors in the entire image */
+ 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 */
result = transform_cell_gsl(cell, solution);
- gsl_matrix_free(solution);
- cleanup_taketwo_obs_vecs(obs_vecs, obs_vec_count);
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);
+ }
+
+ free(tp->prevScores);
+ free(tp->membership);
+ tp->prevScores = NULL;
+ tp->membership = NULL;
+ tp->xtal_num = 0;
+ tp->numPrevs = 0;
+ tp->prevSols = NULL;
+
+}
/* CrystFEL interface hooks */
@@ -1555,6 +2082,34 @@ 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->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;
+ }
+
+ /*
+ 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; i<image_feature_count(image->features); i++ ) {
struct imagefeature *pk = image_get_feature(image->features, i);
@@ -1568,7 +2123,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;
@@ -1631,14 +2186,27 @@ void *taketwo_prepare(IndexingMethod *indm, UnitCell *cell)
tp->cell = cell;
tp->indm = *indm;
+ tp->serial_num = -1;
+ tp->xtal_num = 0;
+ tp->prevSols = NULL;
+ tp->numPrevs = 0;
+ tp->prevScores = NULL;
+ tp->membership = NULL;
+ tp->vec_count = 0;
+ tp->theory_vecs = NULL;
+
+ gen_theoretical_vecs(cell, &tp->theory_vecs, &tp->vec_count);
return tp;
}
-
void taketwo_cleanup(IndexingPrivate *pp)
{
struct taketwo_private *tp = (struct taketwo_private *)pp;
+
+ partial_taketwo_cleanup(tp);
+ free(tp->theory_vecs);
+
free(tp);
}