aboutsummaryrefslogtreecommitdiff
path: root/libcrystfel/src/taketwo.c
diff options
context:
space:
mode:
authorcppxfel <helenginn@Helen'sMacBookPro>2017-06-19 17:47:00 +0200
committercppxfel <helenginn@Helen'sMacBookPro>2017-06-19 17:47:00 +0200
commit9b110b965ce64bcbf24af9a2336533522394a887 (patch)
tree70d0d193bd940fea79594c1c4199da36185ca645 /libcrystfel/src/taketwo.c
parent242e756dcb83321121af8cb434ebf57777e68470 (diff)
Reduced amount of gsl_matrix reallocation
Diffstat (limited to 'libcrystfel/src/taketwo.c')
-rw-r--r--libcrystfel/src/taketwo.c74
1 files changed, 47 insertions, 27 deletions
diff --git a/libcrystfel/src/taketwo.c b/libcrystfel/src/taketwo.c
index cc3d83f0..c9314395 100644
--- a/libcrystfel/src/taketwo.c
+++ b/libcrystfel/src/taketwo.c
@@ -108,7 +108,7 @@ struct TakeTwoCell
#define MAX_NETWORK_MEMBERS (NETWORK_MEMBER_THRESHOLD + 3)
/* Maximum dead ends for a single branch extension during indexing */
-#define MAX_DEAD_ENDS (0)
+#define MAX_DEAD_ENDS (1)
/* Tolerance for two angles to be considered the same */
#define ANGLE_TOLERANCE (deg2rad(0.5))
@@ -227,12 +227,11 @@ static void show_rvec(struct rvec r2)
* functions called under the core functions, still specialised (Level 3)
* ------------------------------------------------------------------------*/
-static gsl_matrix *rotation_around_axis(struct rvec c, double th)
+static gsl_matrix *rotation_around_axis(struct rvec c, double th,
+ gsl_matrix *res)
{
double omc = 1.0 - cos(th);
double s = sin(th);
- gsl_matrix *res = gsl_matrix_alloc(3, 3);
-
gsl_matrix_set(res, 0, 0, cos(th) + c.u*c.u*omc);
gsl_matrix_set(res, 0, 1, c.u*c.v*omc - c.w*s);
gsl_matrix_set(res, 0, 2, c.u*c.w*omc + c.v*s);
@@ -253,7 +252,7 @@ static gsl_matrix *rotation_around_axis(struct rvec c, double th)
* that @result has already been allocated. Will upload the maths to the
* shared Google drive. */
static gsl_matrix *closest_rot_mat(struct rvec vec1, struct rvec vec2,
- struct rvec axis)
+ struct rvec axis, gsl_matrix *twizzle)
{
/* Let's have unit vectors */
normalise_rvec(&vec1);
@@ -303,7 +302,7 @@ static gsl_matrix *closest_rot_mat(struct rvec vec1, struct rvec vec2,
/* Return an identity matrix which has been rotated by
* theta around "axis" */
- return rotation_around_axis(axis, bestAngle);
+ rotation_around_axis(axis, bestAngle, twizzle);
}
@@ -461,12 +460,13 @@ static int symm_rot_mats_are_similar(gsl_matrix *rot1, gsl_matrix *rot2,
return 0;
}
-static gsl_matrix *rotation_between_vectors(struct rvec a, struct rvec b)
+static gsl_matrix *rotation_between_vectors(struct rvec a, struct rvec b,
+ gsl_matrix *twizzle)
{
double th = rvec_angle(a, b);
struct rvec c = rvec_cross(a, b);
normalise_rvec(&c);
- return rotation_around_axis(c, th);
+ rotation_around_axis(c, th, twizzle);
}
@@ -495,10 +495,9 @@ struct rvec gsl_to_rvec(gsl_vector *a)
* intensive on the number crunching side so simple angle checks are used
* to 'pre-scan' vectors beforehand. */
static gsl_matrix *generate_rot_mat(struct rvec obs1, struct rvec obs2,
- struct rvec cell1, struct rvec cell2)
+ struct rvec cell1, struct rvec cell2,
+ gsl_matrix *twiz1, gsl_matrix *twiz2)
{
- gsl_matrix *rotateSpotDiffMatrix;
- gsl_matrix *secondTwizzleMatrix;
gsl_matrix *fullMat;
gsl_vector *cell2v = rvec_to_gsl(cell2);
gsl_vector *cell2vr = gsl_vector_calloc(3);
@@ -510,30 +509,28 @@ static gsl_matrix *generate_rot_mat(struct rvec obs1, struct rvec obs2,
/* Rotate reciprocal space so that the first simulated vector lines up
* with the observed vector. */
- rotateSpotDiffMatrix = rotation_between_vectors(cell1, obs1);
+ rotation_between_vectors(cell1, obs1, twiz1);
normalise_rvec(&obs1);
/* Multiply cell2 by rotateSpotDiffMatrix --> cell2vr */
- gsl_blas_dgemv(CblasNoTrans, 1.0, rotateSpotDiffMatrix, cell2v,
+ gsl_blas_dgemv(CblasNoTrans, 1.0, twiz1, cell2v,
0.0, cell2vr);
/* Now we twirl around the firstAxisUnit until the rotated simulated
* vector matches the second observed vector as closely as possible. */
- secondTwizzleMatrix = closest_rot_mat(gsl_to_rvec(cell2vr), obs2, obs1);
+ closest_rot_mat(gsl_to_rvec(cell2vr), obs2, obs1, twiz2);
/* We want to apply the first matrix and then the second matrix,
* so we multiply these. */
fullMat = gsl_matrix_calloc(3, 3);
gsl_blas_dgemm(CblasTrans, CblasTrans, 1.0,
- rotateSpotDiffMatrix, secondTwizzleMatrix, 0.0, fullMat);
+ twiz1, twiz2, 0.0, fullMat);
gsl_matrix_transpose(fullMat);
gsl_vector_free(cell2v);
gsl_vector_free(cell2vr);
- gsl_matrix_free(secondTwizzleMatrix);
- gsl_matrix_free(rotateSpotDiffMatrix);
return fullMat;
}
@@ -711,11 +708,12 @@ static signed int finish_solution(gsl_matrix *rot, struct SpotVec *obs_vecs,
int *obs_members, int *match_members,
int member_num)
{
- gsl_matrix *sub;
- gsl_matrix *mul;
- sub = gsl_matrix_calloc(3, 3);
- mul = gsl_matrix_calloc(3, 3);
+ gsl_matrix *sub = gsl_matrix_calloc(3, 3);
+ gsl_matrix *mul = gsl_matrix_calloc(3, 3);
+ gsl_matrix *twiz1 = gsl_matrix_calloc(3, 3);
+ gsl_matrix *twiz2 = gsl_matrix_calloc(3, 3);
+
gsl_matrix **rotations = malloc(sizeof(*rotations)* pow(member_num, 2) - member_num);
int i, j, count;
@@ -733,7 +731,8 @@ static signed int finish_solution(gsl_matrix *rot, struct SpotVec *obs_vecs,
struct rvec j_cellvec = j_vec.matches[match_members[j]];
rotations[count] = generate_rot_mat(i_obsvec, j_obsvec,
- i_cellvec, j_cellvec);
+ i_cellvec, j_cellvec,
+ twiz1, twiz2);
count++;
}
@@ -768,6 +767,8 @@ static signed int finish_solution(gsl_matrix *rot, struct SpotVec *obs_vecs,
free(rotations);
gsl_matrix_free(sub);
gsl_matrix_free(mul);
+ gsl_matrix_free(twiz1);
+ gsl_matrix_free(twiz2);
return 1;
}
@@ -783,7 +784,9 @@ static signed int find_next_index(gsl_matrix *rot, struct SpotVec *obs_vecs,
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 < start + 1000; i++ ) {
@@ -836,7 +839,8 @@ static signed int find_next_index(gsl_matrix *rot, struct SpotVec *obs_vecs,
test_rot = generate_rot_mat(obs_vecs[idx_k].obsvec,
obs_vecs[i].obsvec,
member_match,
- a_match);
+ a_match,
+ twiz1, twiz2);
double trace = 0;
int ok = rot_mats_are_similar(rot, test_rot,
@@ -857,6 +861,8 @@ static signed int find_next_index(gsl_matrix *rot, struct SpotVec *obs_vecs,
free(test_idx);
gsl_matrix_free(sub);
gsl_matrix_free(mul);
+ gsl_matrix_free(twiz1);
+ gsl_matrix_free(twiz2);
return i;
}
@@ -869,7 +875,8 @@ static signed int find_next_index(gsl_matrix *rot, struct SpotVec *obs_vecs,
gsl_matrix_free(sub);
gsl_matrix_free(mul);
-
+ gsl_matrix_free(twiz1);
+ gsl_matrix_free(twiz2);
/* give up. */
@@ -977,12 +984,16 @@ 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)
{
+ 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],
- obs_vecs[j].matches[j_match]);
+ obs_vecs[j].matches[j_match],
+ twiz1, twiz2);
/* Try to expand this rotation matrix to a larger network */
@@ -991,6 +1002,8 @@ static int start_seed(struct SpotVec *obs_vecs, int obs_vec_count, int i,
/* return this matrix and if it was immediately successful */
*rotation = rot_mat;
+ gsl_matrix_free(twiz1);
+ gsl_matrix_free(twiz2);
return success;
}
@@ -1003,6 +1016,9 @@ static int weed_duplicate_matches(struct SpotVec *her_obs,
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();
@@ -1022,7 +1038,8 @@ static int weed_duplicate_matches(struct SpotVec *her_obs,
struct rvec j_cellvec = his_obs->matches[his_match];
gsl_matrix *mat = generate_rot_mat(i_obsvec, j_obsvec,
- i_cellvec, j_cellvec);
+ i_cellvec, j_cellvec,
+ twiz1, twiz2);
int found = 0;
@@ -1056,6 +1073,9 @@ static int weed_duplicate_matches(struct SpotVec *her_obs,
free(old_mats);
+ gsl_matrix_free(twiz1);
+ gsl_matrix_free(twiz2);
+
return 1;
}