Taketwo private structure takes care of duplicate solutions

1 files changed, 122 insertions, 43 deletions

diff --git a/libcrystfel/src/taketwo.c b/libcrystfel/src/taketwo.c
index 5acc99d9..0dc1ed51 100644
--- a/libcrystfel/src/taketwo.c
+++ b/libcrystfel/src/taketwo.c
@@ -156,6 +156,12 @@ struct taketwo_private
 {
 	IndexingMethod indm;
 	UnitCell       *cell;
+	int   serial_num; /* -1 if unassigned */
+	unsigned int attempts;
+
+	gsl_matrix     **prevSols; /**< Previous solutions to be ignored */
+	unsigned int   numPrevs; /**< Previous solution count */
+
 };
 
 /**
@@ -167,9 +173,6 @@ struct TakeTwoCell
 	gsl_matrix     **rotSymOps;
 	unsigned int   numOps;
 
-	gsl_matrix     **prevSols; /**< Previous solutions to be ignored */
-	unsigned int   numPrevs; /**< Previous solution count */
-
 	struct SpotVec *obs_vecs;
 	struct Seed *seeds;
 	int seed_count;
@@ -861,39 +864,8 @@ static int weed_duplicate_matches(struct Seed **seeds,
 	}
 
 	signed int i, j;
-	int duplicates = 0;
-
-	/* First we weed out duplicates with previous solutions */
-
-	for (i = *match_count - 1; i >= 0; i--) {
-		int her_match = (*seeds)[i].idx1;
-		int his_match = (*seeds)[i].idx2;
-
-		gsl_matrix *mat;
-		mat = rot_mat_from_indices((*seeds)[i].obs1, (*seeds)[i].obs2,
-		                           her_match, his_match, cell);
-
-		if (mat == NULL)
-		{
-			continue;
-		}
-
-		for (j = 0; j < cell->numPrevs; j++)
-		{
-			int sim = symm_rot_mats_are_similar(cell->prevSols[j],
-			                                    mat, cell);
 
-			/* Found a duplicate with a previous solution */
-			if (sim)
-			{
-				(*seeds)[i].idx1 = -1;
-				(*seeds)[i].idx2 = -1;
-				break;
-			}
-		}
-
-		gsl_matrix_free(mat);
-	}
+	int duplicates = 0;
 
 	/* Now we weed out the self-duplicates from the remaining batch */
 
@@ -1196,8 +1168,51 @@ static int sort_seed_by_score(const void *av, const void *bv)
 	return a->score > b->score;
 }
 
+static void remove_old_solutions(struct TakeTwoCell *cell,
+                                 struct taketwo_private *tp)
+{
+	int duplicates = 0;
+	struct Seed *seeds = cell->seeds;
+	unsigned int total = cell->seed_count;
+
+	/* First we remove duplicates with previous solutions */
+
+	int i, j;
+	for (i = total - 1; i >= 0; i--) {
+		int her_match = seeds[i].idx1;
+		int his_match = seeds[i].idx2;
+
+		gsl_matrix *mat;
+		mat = rot_mat_from_indices(seeds[i].obs1, seeds[i].obs2,
+		                           her_match, his_match, cell);
+
+		if (mat == NULL)
+		{
+			continue;
+		}
+
+		for (j = 0; j < tp->numPrevs; j++)
+		{
+			int sim = symm_rot_mats_are_similar(tp->prevSols[j],
+			                                    mat, cell);
+
+			/* Found a duplicate with a previous solution */
+			if (sim)
+			{
+				seeds[i].idx1 = -1;
+				seeds[i].idx2 = -1;
+				duplicates++;
+				break;
+			}
+		}
 
-static int find_seeds(struct TakeTwoCell *cell)
+		gsl_matrix_free(mat);
+	}
+
+	STATUS("Removing %i duplicates due to prev solutions.\n", duplicates);
+}
+
+static int find_seeds(struct TakeTwoCell *cell, struct taketwo_private *tp)
 {
 	struct SpotVec *obs_vecs = cell->obs_vecs;
 	int obs_vec_count = cell->obs_vec_count;
@@ -1663,7 +1678,8 @@ static void cleanup_taketwo_cell(struct TakeTwoCell *ttCell)
  * successful.
  **/
 static UnitCell *run_taketwo(UnitCell *cell, const struct taketwo_options *opts,
-                             struct rvec *rlps, int rlp_count)
+                             struct rvec *rlps, int rlp_count,
+                             struct taketwo_private *tp)
 {
 	int cell_vec_count = 0;
 	struct TheoryVec *theory_vecs = NULL;
@@ -1678,13 +1694,11 @@ static UnitCell *run_taketwo(UnitCell *cell, const struct taketwo_options *opts,
 	ttCell.seed_count = 0;
 	ttCell.rotSymOps = NULL;
 	ttCell.obs_vecs = NULL;
-	ttCell.prevSols = NULL;
 	ttCell.twiz1Tmp = gsl_matrix_calloc(3, 3);
 	ttCell.twiz2Tmp = gsl_matrix_calloc(3, 3);
 	ttCell.vec1Tmp = gsl_vector_calloc(3);
 	ttCell.vec2Tmp = gsl_vector_calloc(3);
 	ttCell.numOps = 0;
-	ttCell.numPrevs = 0;
 	ttCell.obs_vec_count = 0;
 
 	success = generate_rotation_sym_ops(&ttCell);
@@ -1726,20 +1740,62 @@ static UnitCell *run_taketwo(UnitCell *cell, const struct taketwo_options *opts,
 
 	if ( !success ) return NULL;
 
-	find_seeds(&ttCell);
+	find_seeds(&ttCell, tp);
+//	remove_old_solutions(&ttCell, tp);
 	start_seeds(&solution, &ttCell);
 
 	if ( solution == NULL ) {
 		return NULL;
 	}
 
+	for (int i = 0; i < tp->numPrevs; i++)
+	{
+		gsl_matrix *sol = tp->prevSols[i];
+
+		int sim = symm_rot_mats_are_similar(sol, solution, &ttCell);
+		if (sim)
+		{
+//			STATUS("Warning! Returning previous solution.\n");
+		}
+	}
+
+	int new_size = (tp->numPrevs + 1) * sizeof(gsl_matrix *);
+	gsl_matrix **tmp = realloc(tp->prevSols, new_size);
+
+	if (!tmp) {
+		apologise();
+	}
+
+	tp->prevSols = tmp;
+
+	tp->prevSols[tp->numPrevs] = solution;
+	tp->numPrevs++;
+
 	result = transform_cell_gsl(cell, solution);
-	gsl_matrix_free(solution);
 	cleanup_taketwo_cell(&ttCell);
 
 	return result;
 }
 
+/* Cleans up the per-image information for taketwo */
+
+static void partial_taketwo_cleanup(struct taketwo_private *tp)
+{
+	if (tp->prevSols != NULL)
+	{
+		for (int i = 0; i < tp->numPrevs; i++)
+		{
+			gsl_matrix_free(tp->prevSols[i]);
+		}
+
+		free(tp->prevSols);
+	}
+
+	tp->attempts = 0;
+	tp->numPrevs = 0;
+	tp->prevSols = NULL;
+
+}
 
 /* CrystFEL interface hooks */
 
@@ -1753,6 +1809,24 @@ int taketwo_index(struct image *image, const struct taketwo_options *opts,
 	int i;
 	struct taketwo_private *tp = (struct taketwo_private *)priv;
 
+	/* Check serial number against previous for solution tracking */
+	int this_serial = image->serial;
+
+	if (tp->serial_num == this_serial)
+	{
+		tp->attempts++;
+	}
+	else
+	{
+		partial_taketwo_cleanup(tp);
+		tp->serial_num = this_serial;
+	}
+
+	/*
+	STATUS("Indexing %i with %i attempts, %i crystals\n", this_serial, tp->attempts,
+	       image->n_crystals);
+	*/
+
 	rlps = malloc((image_feature_count(image->features)+1)*sizeof(struct rvec));
 	for ( i=0; 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);
 }