From c24bf8f1f75c18beec6a90df4e7fa620b5ce8d82 Mon Sep 17 00:00:00 2001
From: Helen Ginn <helen@strubi.ox.ac.uk>
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(-)

(limited to 'libcrystfel/src')

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 <gsl/gsl_matrix.h>
 #include <gsl/gsl_blas.h>
 #include <float.h>
@@ -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 <helen@strubi.ox.ac.uk>
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(-)

(limited to 'libcrystfel/src')

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; i<her_obs->match_num; i++ ) {
         for ( j=0; j<his_obs->match_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; 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,
@@ -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; 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 */
@@ -1307,20 +1316,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));
+		*match_array = malloc(count*sizeof(struct TheoryVec));
 		*match_count = count;
 		for ( j=0; j<count; j++ ) {
 			(*match_array)[j] = for_sort[j].v;
@@ -1395,9 +1399,8 @@ static int gen_observed_vecs(struct rvec *rlps, int rlp_count,
 }
 
 
-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,
+				int *vec_count)
 {
 	double a, b, c, alpha, beta, gamma;
 	int h_max, k_max, l_max;
@@ -1421,7 +1424,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++ ) {
@@ -1437,50 +1439,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;
 }
 
@@ -1495,7 +1484,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);
@@ -1532,9 +1520,7 @@ static UnitCell *run_taketwo(UnitCell *cell, const struct taketwo_options *opts,
                              struct rvec *rlps, int rlp_count)
 {
 	int cell_vec_count = 0;
-	int asym_vec_count = 0;
-	struct rvec *cell_vecs = NULL;
-	struct rvec *asym_vecs = NULL;
+	struct TheoryVec *theory_vecs = NULL;
 	UnitCell *result;
 	int obs_vec_count = 0;
 	struct SpotVec *obs_vecs = NULL;
@@ -1552,8 +1538,7 @@ static UnitCell *run_taketwo(UnitCell *cell, const struct taketwo_options *opts,
 
 	success = generate_rotation_sym_ops(&ttCell);
 
-	success = gen_theoretical_vecs(cell, &cell_vecs, &asym_vecs,
-	                               &cell_vec_count, &asym_vec_count);
+	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);
@@ -1586,14 +1571,10 @@ 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);
+	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 <helen@strubi.ox.ac.uk>
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(-)

(limited to 'libcrystfel/src')

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 <helen@strubi.ox.ac.uk>
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(-)

(limited to 'libcrystfel/src')

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; i<her_obs->match_num; i++ ) {
-        for ( j=0; j<his_obs->match_num; j++ ) {
+	for ( i=0; i<her_obs->match_num; i++ ) {
+		for ( j=0; j<his_obs->match_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 <helen@strubi.ox.ac.uk>
Date: Sat, 21 Apr 2018 20:35:12 +0100
Subject: Temporary debug output

---
 libcrystfel/src/taketwo.c | 6 ++++++
 1 file changed, 6 insertions(+)

(limited to 'libcrystfel/src')

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 <helen@strubi.ox.ac.uk>
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(-)

(limited to 'libcrystfel/src')

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<matches; k++ ) {
-				if (i_idx[k] < 0 || j_idx[k] < 0) {
+			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 max_members = 0;
 
 				int success = start_seed(i, j,
-							 i_idx[k], j_idx[k],
+							 seed_idx1, seed_idx2,
 							 rotation, &max_members,
 							 cell);
 
 				if (success) {
-					free(i_idx); free(j_idx);
+					free(seeds);
 					gsl_matrix_free(best_rotation);
 					return success;
 				} else {
@@ -1232,8 +1239,7 @@ static int find_seed(gsl_matrix **rotation, struct TakeTwoCell *cell)
 				}
 			}
 
-			free(i_idx);
-			free(j_idx);
+			free(seeds);
 		}
 	} /* yes this } is meant to be here */
 
-- 
cgit v1.2.3


From b23b43b669a9a467bd5ed2f523dca943a75415e5 Mon Sep 17 00:00:00 2001
From: Helen Ginn <helen@strubi.ox.ac.uk>
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(-)

(limited to 'libcrystfel/src')

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 <helen@strubi.ox.ac.uk>
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(-)

(limited to 'libcrystfel/src')

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 <helen@strubi.ox.ac.uk>
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(-)

(limited to 'libcrystfel/src')

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; i<her_obs->match_num; i++ ) {
 		for ( j=0; j<his_obs->match_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; i<obs_vec_count; i++ ) {
 		for ( j=0; j<i; j++ ) {
 
@@ -1191,7 +1216,7 @@ static int find_seed(gsl_matrix **rotation, struct TakeTwoCell *cell)
 			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],
+			obs_vecs_match_angles(i, j,
 			                      &seeds, &seed_num, cell);
 
 			if (seed_num == 0)
@@ -1201,69 +1226,84 @@ static int find_seed(gsl_matrix **rotation, struct TakeTwoCell *cell)
 			}
 
 			/* Weed out the duplicate seeds (from symmetric
-			 * reflection pairs)
-			 */
+			 * reflection pairs) */
 			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<seed_num; k++ ) {
-				int seed_idx1 = seeds[k].idx1;
-				int seed_idx2 = seeds[k].idx2;
+			/* 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);
 
-				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; k<seed_num; k++ ) {
+		int seed_idx1 = seeds[k].idx1;
+		int seed_idx2 = seeds[k].idx2;
+
+		if (seed_idx1 < 0 || seed_idx2 < 0) {
+			continue;
 		}
 
-		cell->prevSols[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; i<obs_vec_count; i++ ) {
@@ -1413,7 +1454,7 @@ static int match_obs_to_cell_vecs(struct TheoryVec *cell_vecs, int cell_vec_coun
 		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);
+		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++ ) {
@@ -1434,7 +1475,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;
@@ -1445,7 +1486,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]);
 
@@ -1457,13 +1497,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 */
 
@@ -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 <helen@strubi.ox.ac.uk>
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(-)

(limited to 'libcrystfel/src')

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; 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);
@@ -1060,6 +1067,13 @@ static int grow_network(gsl_matrix *rot, int obs_idx1, int obs_idx2,
 	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));
@@ -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 <helen@strubi.ox.ac.uk>
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(-)

(limited to 'libcrystfel/src')

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 <helen@strubi.ox.ac.uk>
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(-)

(limited to 'libcrystfel/src')

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; i<obs_vec_count; i++ ) {
 		for ( j=0; j<i; j++ ) {
 
@@ -1271,13 +1270,6 @@ static int find_seeds(struct TakeTwoCell *cell)
 
 	qsort(cell->seeds, 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 <helen@strubi.ox.ac.uk>
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(-)

(limited to 'libcrystfel/src')

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 <helen@strubi.ox.ac.uk>
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(-)

(limited to 'libcrystfel/src')

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 <helen@strubi.ox.ac.uk>
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(-)

(limited to 'libcrystfel/src')

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; j<his_obs->match_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<obs_vec_count; i++ ) {
+
 		for ( j=0; j<i; j++ ) {
 
+			/** Only check distances which are accumulatively less than
+ * 			the limit */
+			if (obs_vecs[j].distance + obs_vecs[i].distance >
+			    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 <helen@strubi.ox.ac.uk>
Date: Fri, 27 Apr 2018 22:20:56 +0200
Subject: Fix typo

---
 libcrystfel/src/taketwo.c | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

(limited to 'libcrystfel/src')

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 <helen@strubi.ox.ac.uk>
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(-)

(limited to 'libcrystfel/src')

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; i<image_feature_count(image->features); 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 <helen@strubi.ox.ac.uk>
Date: Fri, 27 Apr 2018 22:24:42 +0200
Subject: rvec add rvec function

---
 libcrystfel/src/taketwo.c | 9 +++++++++
 1 file changed, 9 insertions(+)

(limited to 'libcrystfel/src')

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 <helen@strubi.ox.ac.uk>
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(+)

(limited to 'libcrystfel/src')

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 <helen@strubi.ox.ac.uk>
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(-)

(limited to 'libcrystfel/src')

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 <helen@strubi.ox.ac.uk>
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(-)

(limited to 'libcrystfel/src')

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<obs_vec_count; i++ ) {
+	for ( i=1; i<obs_vec_count; i++ ) {
 
 		for ( j=0; j<i; j++ ) {
 
-- 
cgit v1.2.3


From c7d9278a28efb9c0f7166167c1848f339864778b Mon Sep 17 00:00:00 2001
From: Helen Ginn <helen@strubi.ox.ac.uk>
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(-)

(limited to 'libcrystfel/src')

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 <helen@strubi.ox.ac.uk>
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(-)

(limited to 'libcrystfel/src')

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 <helen@strubi.ox.ac.uk>
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(-)

(limited to 'libcrystfel/src')

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<member_num && all_ok; j++ ) {
 
 			struct SpotVec *me = &obs_vecs[i];
@@ -1187,6 +1188,7 @@ static int grow_network(gsl_matrix *rot, int obs_idx1, int obs_idx2,
 			continue;
 		}
 
+		/* 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;
@@ -1311,8 +1313,8 @@ static int find_seeds(struct TakeTwoCell *cell, struct taketwo_private *tp)
 
 		for ( j=0; j<i; j++ ) {
 
-			/** Only check distances which are accumulatively less than
- * 			the limit */
+			/** Only check distances which are accumulatively less
+			* than the limit */
 			if (obs_vecs[j].distance + obs_vecs[i].distance >
 			    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; i<ttCell->numOps; 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 <helen@strubi.ox.ac.uk>
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(-)

(limited to 'libcrystfel/src')

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 <helen@strubi.ox.ac.uk>
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(-)

(limited to 'libcrystfel/src')

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<member_num && skip == 0; j++ ) {
-			if (i == obs_members[j]) {
-				skip = 1;
-			}
-		}
-
-		if (skip) {
-			continue;
-		}
-
 		int all_ok = 1;
 		int matched = -1;
 
-- 
cgit v1.2.3


From 5bdf81f67abe5824b7635ce9855b14c41033682a Mon Sep 17 00:00:00 2001
From: Helen Ginn <helen@strubi.ox.ac.uk>
Date: Sun, 29 Apr 2018 18:55:44 +0200
Subject: Fussiness.

---
 libcrystfel/src/taketwo.c | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

(limited to 'libcrystfel/src')

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 <helen@strubi.ox.ac.uk>
Date: Sun, 29 Apr 2018 18:55:56 +0200
Subject: Remove commented code

---
 libcrystfel/src/taketwo.c | 6 ------
 1 file changed, 6 deletions(-)

(limited to 'libcrystfel/src')

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<member_num; j++ ) {
-			//	STATUS("%i ", obs_members[j]);
-			}
-			//STATUS("%i\n", i);
-
 			return i;
 		}
 	}
-- 
cgit v1.2.3


From 67fb3922a3b0535f8dea5a46a90f5b3e4226fca1 Mon Sep 17 00:00:00 2001
From: Helen Ginn <helen@strubi.ox.ac.uk>
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(-)

(limited to 'libcrystfel/src')

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 <helen@strubi.ox.ac.uk>
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(-)

(limited to 'libcrystfel/src')

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 <helen@strubi.ox.ac.uk>
Date: Sun, 29 Apr 2018 18:57:26 +0200
Subject: Fussiness.

---
 libcrystfel/src/taketwo.c | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

(limited to 'libcrystfel/src')

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 <helen@strubi.ox.ac.uk>
Date: Sun, 29 Apr 2018 19:00:55 +0200
Subject: white space

---
 libcrystfel/src/taketwo.c | 1 -
 1 file changed, 1 deletion(-)

(limited to 'libcrystfel/src')

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 <helen@strubi.ox.ac.uk>
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(-)

(limited to 'libcrystfel/src')

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 <helen@strubi.ox.ac.uk>
Date: Sun, 29 Apr 2018 20:53:53 +0200
Subject: Plug leaks

---
 libcrystfel/src/taketwo.c | 6 ++++++
 1 file changed, 6 insertions(+)

(limited to 'libcrystfel/src')

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 <helen@strubi.ox.ac.uk>
Date: Sun, 29 Apr 2018 20:54:27 +0200
Subject: Plug leaks.

---
 libcrystfel/src/taketwo.c | 3 ---
 1 file changed, 3 deletions(-)

(limited to 'libcrystfel/src')

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 <helen@strubi.ox.ac.uk>
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(-)

(limited to 'libcrystfel/src')

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<i; j++ ) {
 
 			/** Only check distances which are accumulatively less
-			* than the limit */
+			* than the limit if we can easily generate seeds */
 			if (obs_vecs[j].distance + obs_vecs[i].distance >
-			    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 <helen@strubi.ox.ac.uk>
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(-)

(limited to 'libcrystfel/src')

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 <helen@strubi.ox.ac.uk>
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(-)

(limited to 'libcrystfel/src')

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