WARNING - profanity - uploading to my server so I can gvalgrind

1 files changed, 166 insertions, 136 deletions

diff --git a/libcrystfel/src/taketwo.c b/libcrystfel/src/taketwo.c
index d1f2b9b8..b37fc8f1 100644
--- a/libcrystfel/src/taketwo.c
+++ b/libcrystfel/src/taketwo.c
@@ -774,21 +774,84 @@ static signed int finish_solution(gsl_matrix *rot, struct SpotVec *obs_vecs,
 	return 1;
 }
 
+static int weed_duplicate_matches(struct SpotVec *her_obs,
+                                   struct SpotVec *his_obs,
+                                   int **her_match_idxs, int **his_match_idxs,
+                                   int *match_count, struct TakeTwoCell *cell)
+{
+	int num_occupied = 0;
+	gsl_matrix **old_mats = calloc(*match_count, sizeof(gsl_matrix *));
+	gsl_matrix *twiz1 = gsl_matrix_calloc(3, 3);
+	gsl_matrix *twiz2 = gsl_matrix_calloc(3, 3);
+	
+	if (old_mats == NULL)
+	{
+		apologise();
+		return 0;
+	}
+
+	signed int i, j;
+	int duplicates = 0;
+	
+	for (i = *match_count - 1; i >= 0; i--) {
+		int her_match = (*her_match_idxs)[i];
+		int his_match = (*his_match_idxs)[i];
+	
+		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];
+
+		gsl_matrix *mat = generate_rot_mat(i_obsvec, j_obsvec,
+		                                   i_cellvec, j_cellvec,
+		                                   twiz1, twiz2);
+
+		int found = 0;
+		
+		for (j = 0; j < num_occupied; j++) {
+			if (old_mats[j] &&
+			    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;
+				found = 1;
+	
+				duplicates++;
+
+				gsl_matrix_free(mat);
+			}
+		}
+		
+		if (!found) {
+			// we have not found a duplicate, add to list.
+			old_mats[num_occupied] = mat;
+			num_occupied++;
+		}
+	}
+	
+	for (i = 0; i < num_occupied; i++) {
+		if (old_mats[i]) {
+			gsl_matrix_free(old_mats[i]);
+		}
+	}
+
+	free(old_mats);
+	
+	STATUS("Found %i\n", num_occupied);
+	
+	return 1;
+}
+
 static signed int find_next_index(gsl_matrix *rot, struct SpotVec *obs_vecs,
 				  int obs_vec_count, int *obs_members,
 				  int *match_members, int start, int member_num,
-				  int *match_found)
+				  int *match_found, struct TakeTwoCell *cell)
 {
 	int i;
-	
-	gsl_matrix *sub;
-	gsl_matrix *mul;
-	sub = gsl_matrix_calloc(3, 3);
-	mul = gsl_matrix_calloc(3, 3);
 	gsl_matrix *twiz1 = gsl_matrix_calloc(3, 3);
 	gsl_matrix *twiz2 = gsl_matrix_calloc(3, 3);
 	
-
 	for ( i=start; i<obs_vec_count; i++ ) {
 
 		/* first we check for a shared spot - harshest condition */
@@ -798,11 +861,11 @@ static signed int find_next_index(gsl_matrix *rot, struct SpotVec *obs_vecs,
 		if ( !shared ) continue;
 
 		/* now we check that angles between all vectors match */
-	//	int matches = obs_angles_match_array(obs_vecs, i, obs_members,
-	//	match_members, member_num);
+		/*	int matches = obs_angles_match_array(obs_vecs, i, obs_members,
+		 match_members, member_num);
 
-	//	 if ( !matches ) continue;
-		 
+		 if ( !matches ) continue;
+		 */
 		/* final test: does the corresponding rotation matrix
 		 * match the others? NOTE: have not tested to see if
 		 * every combination of test/existing member has to be checked
@@ -818,66 +881,100 @@ static signed int find_next_index(gsl_matrix *rot, struct SpotVec *obs_vecs,
 		obs_vecs_match_angles(&obs_vecs[obs_members[0]], &obs_vecs[i],
 				      &member_idx, &test_idx, &match_num);
 
-		/* if ok is set to 0, give up on this vector before
-		 * checking the next value of j */
+		weed_duplicate_matches(&obs_vecs[obs_members[0]], &obs_vecs[i],
+				       &member_idx, &test_idx, &match_num, cell);
+
+		free(member_idx);
+		free(test_idx);
+
+		// If we have no matches for the first vector,
+		// then this is fruitless so give up now */
+		if (match_num == 0) {
+			continue;
+		}
+
+		/* this observed vec (obs_vecs[i]) must agree with all
+		the other observed vecs in the network (obs_members[member_num])
+		*/
+		 
 		int j, k;
+			
+		int all_ok = 1;
+		
+		for ( j=0; j<member_num && all_ok; j++ ) {
+			// me: member in question (obs_vecs[i])
+			// you: obs_vecs[obs_members[j]]
+			
+			struct SpotVec *me = &obs_vecs[i];
+			struct SpotVec *you = &obs_vecs[obs_members[j]];
+			
+			struct rvec me_obs = me->obsvec;
+			struct rvec you_obs = you->obsvec;
 
-		for ( j=0; j<match_num; j++ ) {
-			int all_ok = 1;
-			for ( k=0; k<member_num; k++ )
-			{
+			int *your_idxs = 0;
+			int *my_idxs = 0;
+
+			obs_vecs_match_angles(me, you, &my_idxs,
+					      &your_idxs, &match_num);
+			
+			STATUS("2\n");
+
+			if (match_num == 0) {
+				all_ok = 0;
+				STATUS("shit\n");
+				continue;
+			}
+			
+			weed_duplicate_matches(me, you,
+					       &my_idxs, &your_idxs,
+				               &match_num, cell);
+
+			STATUS("3\n");
+
+			for ( k=0; k<match_num; k++ ) {
 				gsl_matrix *test_rot;
-				struct rvec member_match;
 
-				int idx_k = obs_members[k];
+				if (my_idxs[k] < 0 || your_idxs[k] < 0) {
+					continue;
+				}
 
-				/* First observed vector and matching theoretical */
-				member_match = obs_vecs[idx_k].matches[match_members[k]];
+				STATUS("4\n");
 
-				/* Potential match with another vector */
-				struct rvec a_match = obs_vecs[i].matches[test_idx[j]];
+				struct rvec me_cell = me->matches[my_idxs[k]];								
+				struct rvec you_cell;
 
-				test_rot = generate_rot_mat(obs_vecs[idx_k].obsvec,
-							    obs_vecs[i].obsvec,
-							    member_match,
-							    a_match,
-							    twiz1, twiz2);
+				you_cell = you->matches[your_idxs[k]];
+					
+				test_rot = generate_rot_mat(me_obs,
+					    you_obs, me_cell, you_cell,
+					    twiz1, twiz2);			
 
 				double trace = 0;
 				int ok = rot_mats_are_similar(rot, test_rot,
-							      sub, mul,
-							      &trace);
+						       twiz1, twiz2, &trace);
+				
 				gsl_matrix_free(test_rot);
 
+				STATUS("5\n");
+
 				if (!ok) {
 					all_ok = 0;
+					STATUS("fuck\n");
+
 					break;
 				}
-
-			}
-
-			if (all_ok) {
-				*match_found = test_idx[j];
-				free(member_idx);
-				free(test_idx);
-				gsl_matrix_free(sub);
-				gsl_matrix_free(mul);
-				gsl_matrix_free(twiz1);
-				gsl_matrix_free(twiz2);
-
-				return i;
 			}
+			
+			free(my_idxs);
+			free(your_idxs);
 		}
-
-		free(member_idx); member_idx = NULL;
-		free(test_idx); member_idx = NULL;
-
-	}
 	
-	gsl_matrix_free(sub);
-	gsl_matrix_free(mul);
-	gsl_matrix_free(twiz1);
-	gsl_matrix_free(twiz2);
+		if (all_ok) {
+			STATUS("phew\n");
+			*match_found = i;
+			return i;
+		}
+	}
 
 	/* give up. */
 
@@ -887,7 +984,8 @@ static signed int find_next_index(gsl_matrix *rot, struct SpotVec *obs_vecs,
 
 static int grow_network(gsl_matrix *rot, struct SpotVec *obs_vecs,
 			int obs_vec_count, int obs_idx1, int obs_idx2,
-			int match_idx1, int match_idx2, int *max_members)
+			int match_idx1, int match_idx2, int *max_members,
+			struct TakeTwoCell *cell)
 {
 	/* indices of members of the self-consistent network of vectors */
 	int obs_members[MAX_NETWORK_MEMBERS];
@@ -920,7 +1018,7 @@ static int grow_network(gsl_matrix *rot, struct SpotVec *obs_vecs,
 							obs_vec_count, obs_members,
 							match_members,
 							start, member_num,
-							&match_found);
+							&match_found, cell);
 
 		if ( member_num < 2 ) {
 			return 0;
@@ -939,6 +1037,7 @@ static int grow_network(gsl_matrix *rot, struct SpotVec *obs_vecs,
 			start++;
 			member_num--;
 			dead_ends++;
+			STATUS("\n");
 
 			continue;
 		}
@@ -983,13 +1082,12 @@ static int grow_network(gsl_matrix *rot, struct SpotVec *obs_vecs,
 
 static int start_seed(struct SpotVec *obs_vecs, int obs_vec_count, int i,
 		      int j, int i_match, int j_match, gsl_matrix **rotation,
-		      int *max_members)
+		      int *max_members, struct TakeTwoCell *cell)
 {
+	gsl_matrix *rot_mat;
 	gsl_matrix *twiz1 = gsl_matrix_calloc(3, 3);
 	gsl_matrix *twiz2 = gsl_matrix_calloc(3, 3);
 	
-	gsl_matrix *rot_mat;
-
 	rot_mat = generate_rot_mat(obs_vecs[i].obsvec,
 				   obs_vecs[j].obsvec,
 				   obs_vecs[i].matches[i_match],
@@ -999,87 +1097,15 @@ static int start_seed(struct SpotVec *obs_vecs, int obs_vec_count, int i,
 	/* Try to expand this rotation matrix to a larger network */
 
 	int success = grow_network(rot_mat, obs_vecs, obs_vec_count,
-				   i, j, i_match, j_match, max_members);
+				   i, j, i_match, j_match, max_members,
+				   cell);
 
 	/* return this matrix and if it was immediately successful */
 	*rotation = rot_mat;
-	gsl_matrix_free(twiz1);
-	gsl_matrix_free(twiz2);
 
 	return success;
 }
 
-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)
-{
-	int num_occupied = 0;
-	gsl_matrix **old_mats = calloc(*match_count, sizeof(gsl_matrix *));
-
-	gsl_matrix *twiz1 = gsl_matrix_calloc(3, 3);
-	gsl_matrix *twiz2 = gsl_matrix_calloc(3, 3);
-	
-	if (old_mats == NULL)
-	{
-		apologise();
-		return 0;
-	}
-
-	signed int i, j;
-	int duplicates = 0;
-	
-	for (i = *match_count - 1; i >= 0; i--) {
-		int her_match = (*her_match_idxs)[i];
-		int his_match = (*his_match_idxs)[i];
-	
-		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];
-
-		gsl_matrix *mat = generate_rot_mat(i_obsvec, j_obsvec,
-		                                   i_cellvec, j_cellvec,
-		                                   twiz1, twiz2);
-
-		int found = 0;
-		
-		for (j = 0; j < num_occupied; j++) {
-			if (old_mats[j] && mat &&
-			    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;
-				found = 1;
-	
-				duplicates++;
-			}
-		}
-
-		if (!found) {
-			// we have not found a duplicate, add to list.
-			old_mats[num_occupied] = mat;
-			num_occupied++;
-		} else {
-			gsl_matrix_free(mat);
-		}
-	}
-	
-	for (i = 0; i < num_occupied; i++) {
-		if (old_mats[i]) {
-			gsl_matrix_free(old_mats[i]);
-		}
-	}
-
-	free(old_mats);
-	
-	gsl_matrix_free(twiz1);
-	gsl_matrix_free(twiz2);
-
-	return 1;
-}
-
 static int find_seed(struct SpotVec *obs_vecs, int obs_vec_count,
 		     gsl_matrix **rotation, struct TakeTwoCell *cell)
 {
@@ -1129,7 +1155,7 @@ static int find_seed(struct SpotVec *obs_vecs, int obs_vec_count,
 			 */
 
 			weed_duplicate_matches(&obs_vecs[i], &obs_vecs[j],
-							&i_idx, &j_idx, &matches, cell);
+					       &i_idx, &j_idx, &matches, cell);
 
 			/* We have seeds! Pass each of them through the seed-starter  */
 			/* If a seed has the highest achieved membership, make note...*/
@@ -1141,9 +1167,11 @@ static int find_seed(struct SpotVec *obs_vecs, int obs_vec_count,
 
 				int max_members = 0;
 
-				int success = start_seed(obs_vecs, obs_vec_count, i, j,
+				int success = start_seed(obs_vecs, obs_vec_count,
+							 i, j,
 							 i_idx[k], j_idx[k],
-							 rotation, &max_members);
+							 rotation, &max_members,
+							 cell);
 
 				if (success) {
 					free(i_idx); free(j_idx);
@@ -1593,12 +1621,14 @@ static UnitCell *run_taketwo(UnitCell *cell, struct rvec *rlps, int rlp_count)
 	}
 
 //	STATUS("Returning result.\n");
-	
+
 	result = transform_cell_gsl(cell, solution);
 	gsl_matrix_free(solution);
 	cleanup_taketwo_obs_vecs(obs_vecs, obs_vec_count);
 	cleanup_taketwo_cell(&ttCell);
-
+	
+	STATUS("good.\n");
+	
 	return result;
 }