Now we refine an indexing solution to match all the observed vectors to the nearest calculated vector before returning solution

1 files changed, 235 insertions, 12 deletions

diff --git a/libcrystfel/src/taketwo.c b/libcrystfel/src/taketwo.c
index 0cf1aeb6..9d557fd2 100644
--- a/libcrystfel/src/taketwo.c
+++ b/libcrystfel/src/taketwo.c
@@ -122,6 +122,7 @@ struct SpotVec
 	struct rvec obsvec;
 	struct TheoryVec *matches;
 	int match_num;
+	int assignment;
 	int in_network;
 	double distance;
 	struct rvec *her_rlp;
@@ -184,6 +185,13 @@ struct TakeTwoCell
 	double angle_tol; /**< In radians */
 	double trace_tol; /**< Contains sqrt(4*(1-cos(angle))) */
 
+	/** A given solution to refine */
+	gsl_matrix *solution;
+
+	double x_ang; /**< Rotations in radians to apply to x axis of solution */
+	double y_ang; /**< Rotations in radians to apply to y axis of solution */
+	double z_ang; /**< Rotations in radians to apply to z axis of solution */
+
 	/**< Temporary memory always allocated for calculations */
 	gsl_matrix *twiz1Tmp;
 	/**< Temporary memory always allocated for calculations */
@@ -217,6 +225,12 @@ struct TakeTwoCell
 /* Tolerance for rot_mats_are_similar */
 #define TRACE_TOLERANCE (deg2rad(3.0))
 
+/* Initial step size for refinement of solutions */
+#define ANGLE_STEP_SIZE (deg2rad(0.2))
+
+/* Final required step size for refinement of solutions */
+#define ANGLE_CONVERGE_SIZE (deg2rad(0.01))
+
 /* TODO: Multiple lattices */
 
 
@@ -349,6 +363,35 @@ static void rotation_around_axis(struct rvec c, double th,
 	gsl_matrix_set(res, 2, 2, cos(th) + c.w*c.w*omc);
 }
 
+/** Rotate GSL matrix by three angles along x, y and z axes */
+static void rotate_gsl_by_angles(gsl_matrix *sol, double x, double y,
+                                 double z, gsl_matrix *result)
+{
+	gsl_matrix *x_rot = gsl_matrix_alloc(3, 3);
+	gsl_matrix *y_rot = gsl_matrix_alloc(3, 3);
+	gsl_matrix *z_rot = gsl_matrix_alloc(3, 3);
+	gsl_matrix *xy_rot = gsl_matrix_alloc(3, 3);
+	gsl_matrix *xyz_rot = gsl_matrix_alloc(3, 3);
+
+	struct rvec x_axis = new_rvec(1, 0, 0);
+	struct rvec y_axis = new_rvec(1, 0, 0);
+	struct rvec z_axis = new_rvec(1, 0, 0);
+
+	rotation_around_axis(x_axis, x, x_rot);
+	rotation_around_axis(y_axis, y, y_rot);
+	rotation_around_axis(z_axis, z, z_rot);
+
+	/* Collapse the rotations in x and y onto z */
+	gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, x_rot,
+	               y_rot, 0.0, xy_rot);
+	gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, xy_rot,
+	               z_rot, 0.0, xyz_rot);
+
+	/* Apply the whole rotation offset to the solution */
+	gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, xyz_rot,
+	               sol, 0.0, result);
+}
+
 
 /* Rotate vector (vec1) around axis (axis) by angle theta. Find value of
  * theta for which the angle between (vec1) and (vec2) is minimised. */
@@ -406,6 +449,7 @@ static void closest_rot_mat(struct rvec vec1, struct rvec vec2,
 	rotation_around_axis(axis, bestAngle, twizzle);
 }
 
+/*
 static double matrix_angle(gsl_matrix *m)
 {
 	double a = gsl_matrix_get(m, 0, 0);
@@ -432,6 +476,7 @@ static struct rvec matrix_axis(gsl_matrix *a)
 	struct rvec v = new_rvec(x, y, z);
 	return v;
 }
+*/
 
 static double matrix_trace(gsl_matrix *a)
 {
@@ -1249,7 +1294,7 @@ static void remove_old_solutions(struct TakeTwoCell *cell,
 		gsl_matrix_free(mat);
 	}
 
-	STATUS("Removing %i duplicates due to prev solutions.\n", duplicates);
+//	STATUS("Removing %i duplicates due to prev solutions.\n", duplicates);
 }
 
 static int find_seeds(struct TakeTwoCell *cell, struct taketwo_private *tp)
@@ -1323,16 +1368,7 @@ static int find_seeds(struct TakeTwoCell *cell, struct taketwo_private *tp)
 				cell->seeds[cell->seed_count] = seeds[i];
 				cell->seed_count++;
 			}
-
-			if (cell->seed_count > 1000) {
-				break;
-			}
 		}
-
-		if (cell->seed_count > 1000) {
-			break;
-		}
-
 	}
 
 	qsort(cell->seeds, cell->seed_count, sizeof(struct Seed), sort_seed_by_score);
@@ -1483,6 +1519,182 @@ static int sort_theory_distances(const void *av, const void *bv)
 	return a->dist > b->dist;
 }
 
+static double obs_to_sol_score(struct TakeTwoCell *ttCell)
+{
+	double total = 0;
+	int count = 0;
+	int i;
+	gsl_matrix *solution = ttCell->solution;
+	gsl_matrix *full_rot = gsl_matrix_alloc(3, 3);
+	rotate_gsl_by_angles(solution, ttCell->x_ang, ttCell->y_ang,
+	                     ttCell->z_ang, full_rot);
+
+	for (i = 0; i < ttCell->obs_vec_count; i++)
+	{
+		struct rvec *obs = &ttCell->obs_vecs[i].obsvec;
+		rvec_to_gsl(ttCell->vec1Tmp, *obs);
+
+		/* Rotate all the observed vectors by the modified soln */
+		/* ttCell->vec2Tmp = 1.0 * full_rot * ttCell->vec1Tmp */
+		gsl_blas_dgemv(CblasTrans, 1.0, full_rot, ttCell->vec1Tmp,
+		               0.0, ttCell->vec2Tmp);
+		struct rvec rotated = gsl_to_rvec(ttCell->vec2Tmp);
+
+		int j = ttCell->obs_vecs[i].assignment;
+
+		if (j < 0) continue;
+
+		struct rvec *match = &ttCell->obs_vecs[i].matches[j].vec;
+		struct rvec diff = diff_vec(rotated, *match);
+
+		double length = rvec_length(diff);
+
+		double addition = exp(-(1 / RECIP_TOLERANCE) * length);
+		total += addition;
+		count++;
+	}
+
+	total /= (double)-count;
+
+	gsl_matrix_free(full_rot);
+
+	return total;
+}
+
+/* At the moment this is just going to do a step search of 3x3x3 until
+ * the angles are sufficiently tiny */
+static void refine_solution(struct TakeTwoCell *ttCell)
+{
+	int i, j, k;
+	const int total = 3 * 3 * 3;
+	const int middle = (total - 1) / 2;
+
+	struct rvec steps[total];
+	double start = obs_to_sol_score(ttCell);
+	const int max_tries = 100;
+
+	int count = 0;
+	double size = ANGLE_STEP_SIZE;
+
+	/* First we create our combinations of steps */
+	for (i = -1; i <= 1; i++) {
+	for (j = -1; j <= 1; j++) {
+	for (k = -1; k <= 1; k++) {
+		struct rvec vec = new_rvec(i, j, k);
+		steps[count] = vec;
+		count++;
+	}
+	}
+	}
+
+	struct rvec current = new_rvec(ttCell->x_ang, ttCell->y_ang,
+	                               ttCell->z_ang);
+
+	double best = start;
+	count = 0;
+
+	while (size > ANGLE_CONVERGE_SIZE && count < max_tries)
+	{
+		struct rvec sized[total];
+
+		int best_num = middle;
+		for (i = 0; i < total; i++)
+		{
+			struct rvec sized_step = steps[i];
+			sized_step.u *= size;
+			sized_step.v *= size;
+			sized_step.w *= size;
+
+			struct rvec new_angles = rvec_add_rvec(current,
+			                                       sized_step);
+
+			sized[i] = new_angles;
+
+			ttCell->x_ang = sized[i].u;
+			ttCell->y_ang = sized[i].v;
+			ttCell->z_ang = sized[i].w;
+
+			double score = obs_to_sol_score(ttCell);
+
+			if (score < best)
+			{
+				best = score;
+				best_num = i;
+			}
+		}
+
+		if (best_num == middle)
+		{
+			size /= 2;
+		}
+
+		current = sized[best_num];
+		count++;
+	}
+
+	ttCell->x_ang = current.u;
+	ttCell->y_ang = current.v;
+	ttCell->z_ang = current.w;
+
+	double end = obs_to_sol_score(ttCell);
+
+	ttCell->x_ang = 0;
+	ttCell->y_ang = 0;
+	ttCell->z_ang = 0;
+
+	gsl_matrix *tmp = gsl_matrix_alloc(3, 3);
+	rotate_gsl_by_angles(ttCell->solution, current.u,
+	                     current.v, current.w, tmp);
+	gsl_matrix_free(ttCell->solution);
+	ttCell->solution = tmp;
+}
+
+static void match_all_obs_to_sol(struct TakeTwoCell *ttCell)
+{
+	int i, j;
+	double total = 0;
+	int count = 0;
+	gsl_matrix *solution = ttCell->solution;
+
+	for (i = 0; i < ttCell->obs_vec_count; i++)
+	{
+		struct rvec *obs = &ttCell->obs_vecs[i].obsvec;
+		rvec_to_gsl(ttCell->vec1Tmp, *obs);
+
+		/* ttCell->vec2Tmp = 1.0 * solution * ttCell->vec1Tmp */
+		gsl_blas_dgemv(CblasTrans, 1.0, solution, ttCell->vec1Tmp,
+		               0.0, ttCell->vec2Tmp);
+		struct rvec rotated = gsl_to_rvec(ttCell->vec2Tmp);
+
+		double smallest = FLT_MAX;
+		int assigned = -1;
+
+		for (j = 0; j < ttCell->obs_vecs[i].match_num; j++)
+		{
+			struct rvec *match = &ttCell->obs_vecs[i].matches[j].vec;
+			struct rvec diff = diff_vec(rotated, *match);
+
+			double length = rvec_length(diff);
+			if (length < smallest)
+			{
+				smallest = length;
+				assigned = j;
+			}
+		}
+
+		ttCell->obs_vecs[i].assignment = assigned;
+
+		if (smallest != FLT_MAX)
+		{
+			double addition = exp(-(1 / RECIP_TOLERANCE) * smallest);
+			total += addition;
+			count++;
+
+		}
+	}
+
+	total /= (double)count;
+}
 
 static int match_obs_to_cell_vecs(struct TheoryVec *cell_vecs, int cell_vec_count,
 				  struct TakeTwoCell *cell)
@@ -1585,6 +1797,7 @@ static int gen_observed_vecs(struct rvec *rlps, int rlp_count,
 				spot_vec.obsvec = diff;
 				spot_vec.distance = sqrt(sqlength);
 				spot_vec.matches = NULL;
+				spot_vec.assignment = -1;
 				spot_vec.match_num = 0;
 				spot_vec.her_rlp = &rlps[i];
 				spot_vec.his_rlp = &rlps[j];
@@ -1746,6 +1959,10 @@ static UnitCell *run_taketwo(UnitCell *cell, const struct taketwo_options *opts,
 	ttCell.vec2Tmp = gsl_vector_calloc(3);
 	ttCell.numOps = 0;
 	ttCell.obs_vec_count = 0;
+	ttCell.solution = NULL;
+	ttCell.x_ang = 0;
+	ttCell.y_ang = 0;
+	ttCell.z_ang = 0;
 
 	success = generate_rotation_sym_ops(&ttCell);
 
@@ -1787,14 +2004,20 @@ static UnitCell *run_taketwo(UnitCell *cell, const struct taketwo_options *opts,
 	if ( !success ) return NULL;
 
 	find_seeds(&ttCell, tp);
-//	remove_old_solutions(&ttCell, tp);
+	remove_old_solutions(&ttCell, tp);
 	start_seeds(&solution, &ttCell);
 
 	if ( solution == NULL ) {
 		return NULL;
 	}
 
-	for (int i = 0; i < tp->numPrevs; i++)
+	ttCell.solution = solution;
+	match_all_obs_to_sol(&ttCell);
+	refine_solution(&ttCell);
+	solution = ttCell.solution;
+
+	int i;
+	for (i = 0; i < tp->numPrevs; i++)
 	{
 		gsl_matrix *sol = tp->prevSols[i];