aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--libcrystfel/src/taketwo.c247
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];