From 185484f94e6b5dbabad8d72c0ade8f2d6cd570ab Mon Sep 17 00:00:00 2001 From: Helen Ginn Date: Wed, 14 Jun 2017 14:29:52 +0100 Subject: Added symmetry checking - better indexing rate but something's holding back a number of crystals --- libcrystfel/src/taketwo.c | 572 ++++++++++++++++++++++++++++++---------------- 1 file changed, 374 insertions(+), 198 deletions(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/taketwo.c b/libcrystfel/src/taketwo.c index 8751d6d1..791af586 100644 --- a/libcrystfel/src/taketwo.c +++ b/libcrystfel/src/taketwo.c @@ -71,6 +71,8 @@ struct SpotVec struct rvec obsvec; struct rvec *matches; int match_num; + struct rvec *asym_matches; + int asym_match_num; double distance; struct rvec *her_rlp; struct rvec *his_rlp; @@ -100,19 +102,19 @@ struct TakeTwoCell #define RECIP_TOLERANCE (0.001*1e10) /* Threshold for network members to consider a potential solution */ -#define NETWORK_MEMBER_THRESHOLD (20) +#define NETWORK_MEMBER_THRESHOLD (30) /* Maximum network members (obviously a solution so should stop) */ -#define MAX_NETWORK_MEMBERS (NETWORK_MEMBER_THRESHOLD + 5) +#define MAX_NETWORK_MEMBERS (NETWORK_MEMBER_THRESHOLD + 3) /* Maximum dead ends for a single branch extension during indexing */ -#define MAX_DEAD_ENDS (5) +#define MAX_DEAD_ENDS (10) /* Tolerance for two angles to be considered the same */ -#define ANGLE_TOLERANCE (deg2rad(3.0)) +#define ANGLE_TOLERANCE (deg2rad(0.7)) /* Tolerance for rot_mats_are_similar */ -#define TRACE_TOLERANCE (deg2rad(3.0)) +#define TRACE_TOLERANCE (deg2rad(5.0)) /** TODO: * @@ -316,10 +318,94 @@ static double matrix_trace(gsl_matrix *a) return tr; } -static int symm_rot_mats_are_similar(gsl_matrix *rot1, gsl_matrix *rot2, - struct TakeTwoCell *cell) +static char *add_ua(const char *inp, char ua) +{ + char *pg = malloc(64); + if ( pg == NULL ) return NULL; + snprintf(pg, 63, "%s_ua%c", inp, ua); + return pg; +} + + +static char *get_chiral_holohedry(UnitCell *cell) +{ + LatticeType lattice = cell_get_lattice_type(cell); + char *pg = malloc(64); + char *pgout = 0; + + if ( pg == NULL ) return NULL; + + switch (lattice) + { + case L_TRICLINIC: + pg = "1"; + break; + + case L_MONOCLINIC: + pg = "2"; + break; + + case L_ORTHORHOMBIC: + pg = "222"; + break; + + case L_TETRAGONAL: + pg = "422"; + break; + + case L_RHOMBOHEDRAL: + pg = "3_R"; + break; + + case L_HEXAGONAL: + if ( cell_get_centering(cell) == 'H' ) { + pg = "3_H"; + } else { + pg = "622"; + } + break; + + case L_CUBIC: + pg = "432"; + break; + + default: + pg = "error"; + break; + } + + switch (lattice) + { + case L_TRICLINIC: + case L_ORTHORHOMBIC: + case L_RHOMBOHEDRAL: + case L_CUBIC: + pgout = strdup(pg); + break; + + case L_MONOCLINIC: + case L_TETRAGONAL: + case L_HEXAGONAL: + pgout = add_ua(pg, cell_get_unique_axis(cell)); + break; + + default: + break; + } + + return pgout; +} + + +static SymOpList *sym_ops_for_cell(UnitCell *cell) { + SymOpList *rawList; + char *pg = get_chiral_holohedry(cell); + rawList = get_pointgroup(pg); + free(pg); + + return rawList; } static int rot_mats_are_similar(gsl_matrix *rot1, gsl_matrix *rot2, @@ -346,6 +432,27 @@ static int rot_mats_are_similar(gsl_matrix *rot1, gsl_matrix *rot2, return (tr < max); } +static int symm_rot_mats_are_similar(gsl_matrix *rot1, gsl_matrix *rot2, + struct TakeTwoCell *cell) +{ + int i; + + for (i = 0; i < cell->numOps; i++) { + gsl_matrix *testRot = gsl_matrix_alloc(3, 3); + gsl_matrix *symOp = cell->rotSymOps[i]; + + gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, rot1, symOp, + 0.0, testRot); + + if (rot_mats_are_similar(testRot, rot2, NULL)) { + return 1; + } + + gsl_matrix_free(testRot); + } + + return 0; +} static gsl_matrix *rotation_between_vectors(struct rvec a, struct rvec b) { @@ -469,6 +576,9 @@ static int obs_vecs_match_angles(struct SpotVec *her_obs, int i, j; *match_count = 0; + double min_angle = deg2rad(2.5); + double max_angle = deg2rad(187.5); + // *her_match_idx = -1; // *his_match_idx = -1; @@ -478,49 +588,57 @@ static int obs_vecs_match_angles(struct SpotVec *her_obs, /* calculate angle between all potential theoretical vectors */ for ( i=0; imatch_num; i++ ) { - for ( j=0; jmatch_num; j++ ) { - - struct rvec *her_match = &her_obs->matches[i]; - struct rvec *his_match = &his_obs->matches[j]; - - double theory_angle = rvec_angle(*her_match, - *his_match); - - /* is this angle a match? */ + for ( j=0; jmatch_num; j++ ) { - double angle_diff = fabs(theory_angle - obs_angle); + struct rvec *her_match = &her_obs->matches[i]; + struct rvec *his_match = &his_obs->matches[j]; - if ( angle_diff < ANGLE_TOLERANCE ) { - size_t new_size = (*match_count + 1) * - sizeof(int); - if (her_match_idxs && his_match_idxs) - { - /* Reallocate the array to fit in another match */ - int *temp_hers; - temp_hers = realloc(*her_match_idxs, new_size); - int *temp_his; - temp_his = realloc(*his_match_idxs, new_size); + double theory_angle = rvec_angle(*her_match, + *his_match); - if ( temp_hers == NULL || temp_his == NULL ) { - apologise(); - } + /* is this angle a match? */ + + /* If the angles are too close to 0 + or 180, one axis ill-determined */ + if (theory_angle < min_angle || + theory_angle > max_angle) { + continue; + } - (*her_match_idxs) = temp_hers; - (*his_match_idxs) = temp_his; + double angle_diff = fabs(theory_angle - obs_angle); - (*her_match_idxs)[*match_count] = i; - (*his_match_idxs)[*match_count] = j; + if ( angle_diff < ANGLE_TOLERANCE ) { + + size_t new_size = (*match_count + 1) * + sizeof(int); + if (her_match_idxs && his_match_idxs) + { + /* Reallocate the array to fit in another match */ + int *temp_hers; + temp_hers = realloc(*her_match_idxs, new_size); + int *temp_his; + temp_his = realloc(*his_match_idxs, new_size); + + if ( temp_hers == NULL || temp_his == NULL ) { + apologise(); } - (*match_count)++; + (*her_match_idxs) = temp_hers; + (*his_match_idxs) = temp_his; + + (*her_match_idxs)[*match_count] = i; + (*his_match_idxs)[*match_count] = j; } + + (*match_count)++; } } + } return (*match_count > 0); } - +// DELETE ME IF UNUSED static int obs_angles_match_array(struct SpotVec *obs_vecs, int test_idx, int *obs_members, int *match_members, int num) { @@ -530,7 +648,7 @@ static int obs_angles_match_array(struct SpotVec *obs_vecs, int test_idx, * is identical - too much faff. **/ - int i, j; + int i = 0; struct SpotVec *her_obs = &obs_vecs[test_idx]; for ( i=0; i *max_members) { *max_members = member_num; } @@ -777,10 +896,23 @@ static int grow_network(gsl_matrix *rot, struct SpotVec *obs_vecs, /* If member_num is high enough, we want to return a yes */ if ( member_num > NETWORK_MEMBER_THRESHOLD ) break; } + + /* + int i; + for (i = 0; i < member_num; i++) + { + STATUS("."); + } - finalise_solution(rot, obs_vecs, obs_members, - match_members, member_num); + if ( member_num > NETWORK_MEMBER_THRESHOLD ) { + STATUS(" yes (%i)\n", member_num); + } else { + STATUS(" nope (%i)\n", member_num); + }*/ + finish_solution(rot, obs_vecs, obs_members, + match_members, member_num); + return ( member_num > NETWORK_MEMBER_THRESHOLD ); } @@ -808,18 +940,26 @@ static int start_seed(struct SpotVec *obs_vecs, int obs_vec_count, int i, } 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) + 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 *)); + if (old_mats == NULL) + { + apologise(); + return 0; + } + signed int i, j; - for (i = *match_count - 1; i >= 0; i++) { + 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]; @@ -828,6 +968,8 @@ static int weed_duplicate_matches(struct SpotVec *her_obs, gsl_matrix *mat = generate_rot_mat(i_obsvec, j_obsvec, i_cellvec, j_cellvec); + int found = 0; + for (j = 0; j < num_occupied; j++) { if (old_mats[j] && symm_rot_mats_are_similar(old_mats[j], mat, cell)) @@ -835,16 +977,21 @@ static int weed_duplicate_matches(struct SpotVec *her_obs, // 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); - } else { - // we have not found a duplicate, add to list. - old_mats[num_occupied] = mat; - num_occupied++; } } + + 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]); @@ -852,6 +999,8 @@ static int weed_duplicate_matches(struct SpotVec *her_obs, } free(old_mats); + + return 1; } static int find_seed(struct SpotVec *obs_vecs, int obs_vec_count, @@ -861,7 +1010,7 @@ static int find_seed(struct SpotVec *obs_vecs, int obs_vec_count, int max_max_members = 0; gsl_matrix *best_rotation = NULL; - unsigned long start_time = time(NULL); +// unsigned long start_time = time(NULL); /* loop round pairs of vectors to try and find a suitable * seed to start building a self-consistent network of vectors @@ -885,7 +1034,10 @@ static int find_seed(struct SpotVec *obs_vecs, int obs_vec_count, int *i_idx = 0; int *j_idx = 0; int matches; - + + // double dist1Pc = 100 * obs_vecs[i].distance / MAX_RECIP_DISTANCE; +// double dist2Pc = 100 * obs_vecs[j].distance / MAX_RECIP_DISTANCE; + /* Check to see if any angles match from the cell vectors */ obs_vecs_match_angles(&obs_vecs[i], &obs_vecs[j], &i_idx, &j_idx, &matches); @@ -932,8 +1084,8 @@ static int find_seed(struct SpotVec *obs_vecs, int obs_vec_count, } } - unsigned long now_time = time(NULL); - unsigned int seconds = now_time - start_time; +// unsigned long now_time = time(NULL); +// unsigned int seconds = now_time - start_time; // Commented out for the time being. /* @@ -957,121 +1109,96 @@ static int find_seed(struct SpotVec *obs_vecs, int obs_vec_count, return (best_rotation != NULL); } - -static char *add_ua(const char *inp, char ua) -{ - char *pg = malloc(64); - if ( pg == NULL ) return NULL; - snprintf(pg, 63, "%s_ua%c", inp, ua); - return pg; -} - - -static char *get_chiral_holohedry(UnitCell *cell) +static void set_gsl_matrix(gsl_matrix *mat, double asx, double asy, double asz, + double bsx, double bsy, double bsz, + double csx, double csy, double csz) { - LatticeType lattice = cell_get_lattice_type(cell); - char *pg = malloc(64); - char *pgout; - - if ( pg == NULL ) return NULL; - - switch (lattice) - { - case L_TRICLINIC: - pg = "1"; - break; - - case L_MONOCLINIC: - pg = "2"; - break; - - case L_ORTHORHOMBIC: - pg = "222"; - break; - - case L_TETRAGONAL: - pg = "422"; - break; - - case L_RHOMBOHEDRAL: - pg = "3_R"; - break; - - case L_HEXAGONAL: - if ( cell_get_centering(cell) == 'H' ) { - pg = "3_H"; - } else { - pg = "622"; - } - break; - - case L_CUBIC: - pg = "432"; - break; - - default: - pg = "error"; - break; - } - - switch (lattice) - { - case L_TRICLINIC: - case L_ORTHORHOMBIC: - case L_RHOMBOHEDRAL: - case L_CUBIC: - pgout = strdup(pg); - break; - - case L_MONOCLINIC: - case L_TETRAGONAL: - case L_HEXAGONAL: - pgout = add_ua(pg, cell_get_unique_axis(cell)); - break; - - default: - break; - } - - return pgout; + gsl_matrix_set(mat, 0, 0, asx); + gsl_matrix_set(mat, 0, 1, asy); + gsl_matrix_set(mat, 0, 2, asz); + gsl_matrix_set(mat, 1, 0, bsx); + gsl_matrix_set(mat, 1, 1, bsy); + gsl_matrix_set(mat, 1, 2, bsz); + gsl_matrix_set(mat, 2, 0, csx); + gsl_matrix_set(mat, 2, 1, csy); + gsl_matrix_set(mat, 2, 2, csz); } - static int generate_rotation_sym_ops(struct TakeTwoCell *ttCell) { - SymOpList *rawList; - - char *pg = get_chiral_holohedry(ttCell->cell); - rawList = get_pointgroup(pg); - free(pg); + SymOpList *rawList = sym_ops_for_cell(ttCell->cell); /* Now we must convert these into rotation matrices rather than hkl * transformations (affects triclinic, monoclinic, rhombohedral and * hexagonal space groups only) */ + + double asx, asy, asz; + double bsx, bsy, bsz; + double csx, csy, csz; + gsl_matrix *recip = gsl_matrix_alloc(3, 3); + gsl_matrix *cart = gsl_matrix_alloc(3, 3); cell_get_reciprocal(ttCell->cell, &asx, &asy, &asz, &bsx, &bsy, &bsz, &csx, &csy, &csz); - inv = gsl_matrix_alloc(3, 3); - gsl_matrix_set(inv, 0, 0, asx); - gsl_matrix_set(inv, 0, 1, asy); - gsl_matrix_set(inv, 0, 2, asz); - gsl_matrix_set(inv, 1, 0, bsx); - gsl_matrix_set(inv, 1, 1, bsy); - gsl_matrix_set(inv, 1, 2, bsz); - gsl_matrix_set(inv, 2, 0, csx); - gsl_matrix_set(inv, 2, 1, csy); - gsl_matrix_set(inv, 2, 2, csz); + set_gsl_matrix(recip, asx, asy, asz, bsx, bsy, bsz, csx, csy, csz); - // FIXME: Finish me! Along the lines of: - /* - MatrixPtr newMat = matFromCCP4(&spaceGroup->symop[j]); - MatrixPtr ortho = newMat; - ortho->preMultiply(*unitCell); - ortho->multiply(*reverse); - */ + cell_get_cartesian(ttCell->cell, &asx, &asy, &asz, &bsx, &bsy, + &bsz, &csx, &csy, &csz); + + set_gsl_matrix(cart, asx, asy, asz, bsx, bsy, bsz, csx, csy, csz); + + int i, j, k; + int numOps = num_equivs(rawList, NULL); + + ttCell->rotSymOps = malloc(numOps * sizeof(gsl_matrix *)); + ttCell->numOps = numOps; + + if (ttCell->rotSymOps == NULL) { + apologise(); + return 0; + } + + for (i = 0; i < numOps; i++) + { + gsl_matrix *symOp = gsl_matrix_alloc(3, 3); + IntegerMatrix *op = get_symop(rawList, NULL, i); + + for (j = 0; j < 3; j++) { + for (k = 0; k < 3; k++) { + gsl_matrix_set(symOp, j, k, intmat_get(op, j, k)); + } + } + + gsl_matrix *first = gsl_matrix_calloc(3, 3); + gsl_matrix *second = gsl_matrix_calloc(3, 3); + + /* Each equivalence is of the form: + cartesian * symOp * reciprocal. + First multiplication: symOp * reciprocal */ + + gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, + 1.0, symOp, recip, + 0.0, first); + + /* Second multiplication: cartesian * first */ + + gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, + 1.0, cart, first, + 0.0, second); + + ttCell->rotSymOps[i] = second; + + gsl_matrix_free(symOp); + gsl_matrix_free(first); + } - gsl_matrix_free(inv); + // FIXME: Finish me! Along the lines of: + + gsl_matrix_free(cart); + gsl_matrix_free(recip); + + return 1; } @@ -1090,7 +1217,8 @@ static int sort_func(const void *av, const void *bv) static int match_obs_to_cell_vecs(struct rvec *cell_vecs, int cell_vec_count, - struct SpotVec *obs_vecs, int obs_vec_count) + struct SpotVec *obs_vecs, int obs_vec_count, + int is_asymmetric) { int i, j; @@ -1098,9 +1226,8 @@ static int match_obs_to_cell_vecs(struct rvec *cell_vecs, int cell_vec_count, int count = 0; struct sortme *for_sort = NULL; - + for ( j=0; j