aboutsummaryrefslogtreecommitdiff
path: root/libcrystfel/src
diff options
context:
space:
mode:
authorHelen Ginn <ginn@rescomp1.(none)>2017-01-17 15:15:29 +0000
committerHelen Ginn <ginn@rescomp1.(none)>2017-01-17 15:15:29 +0000
commit1db7d84152ecd8cea4922cabb9db725943f83151 (patch)
treec9bf960697e0c55e3ea661f5eef4f08f0be03d23 /libcrystfel/src
parent743f8019286e7f68bd7c88580673ef7ccd63e3a8 (diff)
General beginning of tidying and removing some debug statements and return best solution if they're all too rubbish to hit the threshold - not tested
Diffstat (limited to 'libcrystfel/src')
-rw-r--r--libcrystfel/src/taketwo.c163
1 files changed, 64 insertions, 99 deletions
diff --git a/libcrystfel/src/taketwo.c b/libcrystfel/src/taketwo.c
index 52540d15..b5604d1a 100644
--- a/libcrystfel/src/taketwo.c
+++ b/libcrystfel/src/taketwo.c
@@ -352,7 +352,7 @@ 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, int printDebug)
+ struct rvec cell1, struct rvec cell2)
{
gsl_matrix *rotateSpotDiffMatrix;
gsl_matrix *secondTwizzleMatrix;
@@ -360,46 +360,25 @@ static gsl_matrix *generate_rot_mat(struct rvec obs1, struct rvec obs2,
gsl_vector *cell2v = rvec_to_gsl(cell2);
gsl_vector *cell2vr = gsl_vector_calloc(3);
- normalise_rvec(&obs1);
- normalise_rvec(&obs2);
- normalise_rvec(&cell1);
- normalise_rvec(&cell2);
-
- if (0 && printDebug)
- {
- STATUS("--\n");
- show_rvec(obs1);
- show_rvec(obs2);
- show_rvec(cell1);
- show_rvec(cell2);
- STATUS(".....\n");
- }
+ normalise_rvec(&obs1);
+ normalise_rvec(&obs2);
+ normalise_rvec(&cell1);
+ normalise_rvec(&cell2);
- /* Rotate reciprocal space so that the first simulated vector lines up
- * with the observed vector. */
- rotateSpotDiffMatrix = rotation_between_vectors(cell1, obs1);
- if (printDebug)
- {
- // STATUS("rotateSpotDiffMatrix:\n");
- // show_matrix(rotateSpotDiffMatrix);
- }
+ /* Rotate reciprocal space so that the first simulated vector lines up
+ * with the observed vector. */
+ rotateSpotDiffMatrix = rotation_between_vectors(cell1, obs1);
- normalise_rvec(&obs1);
+ normalise_rvec(&obs1);
- /* Multiply cell2 by rotateSpotDiffMatrix --> cell2vr */
- gsl_blas_dgemv(CblasNoTrans, 1.0, rotateSpotDiffMatrix, cell2v,
- 0.0, cell2vr);
+ /* Multiply cell2 by rotateSpotDiffMatrix --> cell2vr */
+ gsl_blas_dgemv(CblasNoTrans, 1.0, rotateSpotDiffMatrix, cell2v,
+ 0.0, cell2vr);
- /* Now we twirl around the firstAxisUnit until the rotated simulated
- * vector matches the second observed vector as closely as possible. */
+ /* 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);
-
- if (printDebug)
- {
- // STATUS("secondTwizzleMatrix:\n");
- // show_matrix(secondTwizzleMatrix);
- }
+ secondTwizzleMatrix = closest_rot_mat(gsl_to_rvec(cell2vr), obs2, obs1);
/* We want to apply the first matrix and then the second matrix,
* so we multiply these. */
@@ -407,12 +386,6 @@ static gsl_matrix *generate_rot_mat(struct rvec obs1, struct rvec obs2,
gsl_blas_dgemm(CblasTrans, CblasTrans, 1.0,
rotateSpotDiffMatrix, secondTwizzleMatrix, 0.0, fullMat);
gsl_matrix_transpose(fullMat);
-
- if (printDebug)
- {
- STATUS("final matrix:\n");
- show_matrix(fullMat);
- }
return fullMat;
}
@@ -611,16 +584,12 @@ static signed int find_next_index(gsl_matrix *rot, struct SpotVec *obs_vecs,
test_rot = generate_rot_mat(obs_vecs[idx0].obsvec,
obs_vecs[i].obsvec,
member_match,
- a_match, 0);
+ a_match);
int ok = rot_mats_are_similar(rot, test_rot);
if (ok)
{
- /* generate_rot_mat(obs_vecs[idx0].obsvec,
- obs_vecs[i].obsvec,
- member_match,
- a_match, 0);*/
*match_found = j;
return i;
}
@@ -635,11 +604,12 @@ 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 match_idx1, int match_idx2, int *max_members,
+ int accept_anything)
{
/* indices of members of the self-consistent network of vectors */
int obs_members[MAX_NETWORK_MEMBERS];
- int match_members[MAX_NETWORK_MEMBERS];
+ int match_members[MAX_NETWORK_MEMBERS];
/* initialise the ones we know already */
obs_members[0] = obs_idx1;
@@ -703,6 +673,11 @@ static int grow_network(gsl_matrix *rot, struct SpotVec *obs_vecs,
start = next_index + 1;
member_num++;
+
+ if (member_num > *max_members) {
+ *max_members = member_num;
+ }
+
STATUS("elongating seed, now %i members. They are: ", member_num);
int c;
for (c = 0; c < member_num; c++) {
@@ -732,7 +707,8 @@ static signed int spot_idx(struct rvec *rlp)
}
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 j, int i_match, int j_match, gsl_matrix **rotation,
+ int *max_members, int accept_anything)
{
STATUS("starting new seed for obsvs %i and %i..\n", i, j);
gsl_matrix *rot_mat = gsl_matrix_calloc(3, 3);
@@ -741,7 +717,7 @@ static int start_seed(struct SpotVec *obs_vecs, int obs_vec_count, int i,
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], 0);
+ obs_vecs[j].matches[j_match]);
/* Try to expand this rotation matrix to a larger network */
STATUS("idx: %i %i %i %i spots %i %i and %i %i\n",
@@ -753,7 +729,8 @@ static int start_seed(struct SpotVec *obs_vecs, int obs_vec_count, int i,
show_matrix(rot_mat);
int success = grow_network(rot_mat, obs_vecs, obs_vec_count,
- i, j, i_match, j_match);
+ i, j, i_match, j_match, max_members,
+ accept_anything);
/* return this matrix or free it and try again */
if ( success ) {
@@ -766,13 +743,17 @@ static int start_seed(struct SpotVec *obs_vecs, int obs_vec_count, int i,
}
static int find_seed(struct SpotVec *obs_vecs, int obs_vec_count,
- gsl_matrix **rotation)
+ gsl_matrix **rotation)
{
+ /* META: Maximum achieved maximum network membership */
+ int max_max_members = 0;
+ int best_i, best_j, best_i_idx, best_j_idx;
+
/* loop round pairs of vectors to try and find a suitable
* seed to start building a self-consistent network of vectors
*/
int i, j;
-
+
for ( i=0; i<obs_vec_count-1; i++ ) {
for ( j=i+1; j<obs_vec_count; j++ ) {
@@ -796,20 +777,45 @@ static int find_seed(struct SpotVec *obs_vecs, int obs_vec_count,
if ( matches == 0 ) continue;
STATUS("...ok, %i matches\n", matches);
- /* We have seeds! Pass each of them through the seed-starter */
+ /* We have seeds! Pass each of them through the seed-starter */
+ /* If a seed has the highest achieved membership, make note...*/
int k;
for ( k=0; k<1; k++ ) {
+ int max_members = 0;
+
int success = start_seed(obs_vecs, obs_vec_count, i, j,
i_idx[k], j_idx[k],
- rotation);
+ rotation, &max_members,
+ 0);
- if (success) { return success; }
+ if (success) {
+ return success;
+ } else {
+ if (max_members > max_max_members) {
+ max_max_members = max_members;
+ best_i = i;
+ best_j = j;
+ best_i_idx = i_idx[k];
+ best_j_idx = j_idx[k];
+ }
+ }
}
}
} /* yes this } is meant to be here */
- return 0;
+ int max_members = 0;
+ int success = start_seed(obs_vecs, obs_vec_count, best_i, best_j,
+ best_i_idx, best_j_idx,
+ rotation, &max_members,
+ 1);
+
+ if (success) {
+ return success;
+ } else {
+ /* Well, shit... something went wrong. */
+ return 0;
+ }
}
@@ -873,12 +879,6 @@ static int match_obs_to_cell_vecs(struct rvec *cell_vecs, int cell_vec_count,
obs_vecs[i].match_num = count;
for ( j=0; j<count; j++ ) {
obs_vecs[i].matches[j] = for_sort[j].v;
- if ( (i==0) || (i==1) || (i==2) ) {
- /* STATUS("fluff %i %14f %14f : ", j,
- rvec_length(for_sort[j].v)/1e10,
- for_sort[j].dist/1e10);
- show_rvec(for_sort[j].v);*/
- }
}
free(for_sort);
@@ -1066,15 +1066,6 @@ global_nrlps = rlp_count;
find_seed(obs_vecs, obs_vec_count, &solution);
if ( solution == NULL ) return NULL;
-/*
- double x0 = gsl_matrix_get(solution, 0, 0);
- gsl_matrix_set(solution, 0, 0, -x0);
-
- double x1 = gsl_matrix_get(solution, 0, 1);
- gsl_matrix_set(solution, 0, 1, -x1);
-
- double x2 = gsl_matrix_get(solution, 0, 2);
- gsl_matrix_set(solution, 0, 2, -x2);*/
result = transform_cell_gsl(cell, solution);
@@ -1107,32 +1098,6 @@ int taketwo_index(struct image *image, IndexingPrivate *ipriv)
rlps[n_rlps].u = 0.0;
rlps[n_rlps].v = 0.0;
rlps[n_rlps++].w = 0.0;
-
-/*
- FILE *fh = fopen("../../spots.csv", "r");
- rlps = malloc(500*sizeof(struct rvec));
- while ( (fh != NULL) && !feof(fh) ) {
- float x, y, z;
- float dd;
- int di;
- char line[1024];
- fgets(line, 1024, fh);
-
- if ( sscanf(line, "%f,%f,%f,%f,%i,%i,%f,%f,%i",
- &x, &y, &z, &dd, &di, &di, &dd, &dd, &di) == 9 ) {
- if ( isnan(x) ) continue;
- STATUS("got %f %f %f\n", x, y, z);
- rlps[n_rlps].u = x*1e10;
- rlps[n_rlps].v = y*1e10;
- rlps[n_rlps].w = z*1e10;
- n_rlps++;
- }
- }
- rlps[n_rlps].u = 0.0;
- rlps[n_rlps].v = 0.0;
- rlps[n_rlps++].w = 0.0;
-
- if ( fh != NULL ) fclose(fh);*/
cell = run_taketwo(tp->cell, rlps, n_rlps);
if ( cell == NULL ) return 0;