From 12641fdf021f45fe7a4b19a384b80f3cfe8f9277 Mon Sep 17 00:00:00 2001 From: cppxfel Date: Wed, 5 Jul 2017 10:12:32 +0200 Subject: Shoved a bunch of the stuff being passed from function to function into TakeTwoCell, and shoved gsl_vector allocation into it as well to reduce alloc events --- libcrystfel/src/taketwo.c | 133 ++++++++++++++++++++++------------------------ 1 file changed, 64 insertions(+), 69 deletions(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/taketwo.c b/libcrystfel/src/taketwo.c index 4cd90cd4..c1011a96 100644 --- a/libcrystfel/src/taketwo.c +++ b/libcrystfel/src/taketwo.c @@ -92,6 +92,12 @@ struct TakeTwoCell UnitCell *cell; gsl_matrix **rotSymOps; unsigned int numOps; + struct SpotVec **obs_vecs; // Pointer to an array + int obs_vec_count; + gsl_matrix *twiz1Tmp; + gsl_matrix *twiz2Tmp; + gsl_vector *vec1Tmp; + gsl_vector *vec2Tmp; }; @@ -118,7 +124,7 @@ struct TakeTwoCell #define ANGLE_TOLERANCE (deg2rad(0.6)) /* Tolerance for rot_mats_are_similar */ -#define TRACE_TOLERANCE (deg2rad(2.5)) +#define TRACE_TOLERANCE (deg2rad(3.0)) /** TODO: * @@ -467,13 +473,11 @@ static void rotation_between_vectors(struct rvec a, struct rvec b, } -static gsl_vector *rvec_to_gsl(struct rvec v) +static void rvec_to_gsl(gsl_vector *newVec, struct rvec v) { - gsl_vector *a = gsl_vector_alloc(3); - gsl_vector_set(a, 0, v.u); - gsl_vector_set(a, 1, v.v); - gsl_vector_set(a, 2, v.w); - return a; + gsl_vector_set(newVec, 0, v.u); + gsl_vector_set(newVec, 1, v.v); + gsl_vector_set(newVec, 2, v.w); } @@ -493,11 +497,13 @@ struct rvec gsl_to_rvec(gsl_vector *a) * to 'pre-scan' vectors beforehand. */ static gsl_matrix *generate_rot_mat(struct rvec obs1, struct rvec obs2, struct rvec cell1, struct rvec cell2, - gsl_matrix *twiz1, gsl_matrix *twiz2) + struct TakeTwoCell *cell) { gsl_matrix *fullMat; - gsl_vector *cell2v = rvec_to_gsl(cell2); - gsl_vector *cell2vr = gsl_vector_calloc(3); + rvec_to_gsl(cell->vec1Tmp, cell2); + +// gsl_vector *cell2v = rvec_to_gsl(cell2); +// gsl_vector *cell2vr = gsl_vector_calloc(3); normalise_rvec(&obs1); normalise_rvec(&obs2); @@ -506,29 +512,26 @@ 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. */ - rotation_between_vectors(cell1, obs1, twiz1); + rotation_between_vectors(cell1, obs1, cell->twiz1Tmp); normalise_rvec(&obs1); /* Multiply cell2 by rotateSpotDiffMatrix --> cell2vr */ - gsl_blas_dgemv(CblasNoTrans, 1.0, twiz1, cell2v, - 0.0, cell2vr); + gsl_blas_dgemv(CblasNoTrans, 1.0, cell->twiz1Tmp, cell->vec1Tmp, + 0.0, cell->vec2Tmp); /* Now we twirl around the firstAxisUnit until the rotated simulated * vector matches the second observed vector as closely as possible. */ - closest_rot_mat(gsl_to_rvec(cell2vr), obs2, obs1, twiz2); + closest_rot_mat(gsl_to_rvec(cell->vec2Tmp), obs2, obs1, cell->twiz2Tmp); /* 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, - twiz1, twiz2, 0.0, fullMat); + cell->twiz1Tmp, cell->twiz2Tmp, 0.0, fullMat); gsl_matrix_transpose(fullMat); - gsl_vector_free(cell2v); - gsl_vector_free(cell2vr); - return fullMat; } @@ -664,15 +667,13 @@ static int obs_vecs_match_angles(struct SpotVec *her_obs, static signed int finish_solution(gsl_matrix *rot, struct SpotVec *obs_vecs, int *obs_members, int *match_members, - int member_num) + int member_num, struct TakeTwoCell *cell) { 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); + gsl_matrix **rotations = malloc(sizeof(*rotations)* pow(member_num, 2) + - member_num); int i, j, count; @@ -690,7 +691,7 @@ static signed int finish_solution(gsl_matrix *rot, struct SpotVec *obs_vecs, rotations[count] = generate_rot_mat(i_obsvec, j_obsvec, i_cellvec, j_cellvec, - twiz1, twiz2); + cell); count++; } @@ -725,8 +726,6 @@ 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; } @@ -738,8 +737,6 @@ 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) { @@ -760,8 +757,7 @@ 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, - twiz1, twiz2); + i_cellvec, j_cellvec, cell); int found = 0; @@ -799,14 +795,16 @@ static int weed_duplicate_matches(struct SpotVec *her_obs, return 1; } -static signed int find_next_index(gsl_matrix *rot, struct SpotVec *obs_vecs, - int obs_vec_count, int *obs_members, +static signed int find_next_index(gsl_matrix *rot, int *obs_members, int *match_members, int start, int member_num, int *match_found, struct TakeTwoCell *cell) { + struct SpotVec *obs_vecs = *(cell->obs_vecs); + int obs_vec_count = cell->obs_vec_count; + gsl_matrix *sub = gsl_matrix_calloc(3, 3); + gsl_matrix *mul = gsl_matrix_calloc(3, 3); + int i, j, k; - gsl_matrix *twiz1 = gsl_matrix_calloc(3, 3); - gsl_matrix *twiz2 = gsl_matrix_calloc(3, 3); for ( i=start; iobs_vecs); + int obs_vec_count = cell->obs_vec_count; + /* indices of members of the self-consistent network of vectors */ int obs_members[MAX_NETWORK_MEMBERS]; int match_members[MAX_NETWORK_MEMBERS]; @@ -927,8 +927,7 @@ static int grow_network(gsl_matrix *rot, struct SpotVec *obs_vecs, int match_found = -1; - signed int next_index = find_next_index(rot, obs_vecs, - obs_vec_count, obs_members, + signed int next_index = find_next_index(rot, obs_members, match_members, start, member_num, &match_found, cell); @@ -971,30 +970,29 @@ static int grow_network(gsl_matrix *rot, struct SpotVec *obs_vecs, } finish_solution(rot, obs_vecs, obs_members, - match_members, member_num); + match_members, member_num, cell); return ( member_num > NETWORK_MEMBER_THRESHOLD ); } -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, struct TakeTwoCell *cell) +static int start_seed(int i, int j, int i_match, int j_match, + gsl_matrix **rotation, int *max_members, + struct TakeTwoCell *cell) { + struct SpotVec *obs_vecs = *(cell->obs_vecs); + gsl_matrix *rot_mat; - gsl_matrix *twiz1 = gsl_matrix_calloc(3, 3); - gsl_matrix *twiz2 = gsl_matrix_calloc(3, 3); 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], - twiz1, twiz2); + cell); /* Try to expand this rotation matrix to a larger network */ - int success = grow_network(rot_mat, obs_vecs, obs_vec_count, - i, j, i_match, j_match, max_members, + int success = grow_network(rot_mat, i, j, i_match, j_match, max_members, cell); /* return this matrix and if it was immediately successful */ @@ -1003,9 +1001,11 @@ static int start_seed(struct SpotVec *obs_vecs, int obs_vec_count, int i, return success; } -static int find_seed(struct SpotVec *obs_vecs, int obs_vec_count, - gsl_matrix **rotation, struct TakeTwoCell *cell) +static int find_seed(gsl_matrix **rotation, struct TakeTwoCell *cell) { + struct SpotVec *obs_vecs = *(cell->obs_vecs); + int obs_vec_count = cell->obs_vec_count; + /* META: Maximum achieved maximum network membership */ int max_max_members = 0; gsl_matrix *best_rotation = NULL; @@ -1064,8 +1064,7 @@ static int find_seed(struct SpotVec *obs_vecs, int obs_vec_count, int max_members = 0; - int success = start_seed(obs_vecs, obs_vec_count, - i, j, + int success = start_seed(i, j, i_idx[k], j_idx[k], rotation, &max_members, cell); @@ -1085,21 +1084,6 @@ static int find_seed(struct SpotVec *obs_vecs, int obs_vec_count, *rotation = NULL; } } - -// unsigned long now_time = time(NULL); -// unsigned int seconds = now_time - start_time; - - // Commented out for the time being. - /* - if (seconds > 30) { - // Heading towards CrystFEL cutoff so - // return your best guess and run - free(i_idx); free(j_idx); - - *rotation = best_rotation; - return (best_rotation != NULL); - } - */ } free(i_idx); @@ -1457,6 +1441,10 @@ static void cleanup_taketwo_cell(struct TakeTwoCell *ttCell) gsl_matrix_free(ttCell->rotSymOps[i]); } + free(ttCell->vec1Tmp); + free(ttCell->vec2Tmp); + free(ttCell->twiz1Tmp); + free(ttCell->twiz2Tmp); free(ttCell->rotSymOps); } @@ -1488,6 +1476,10 @@ static UnitCell *run_taketwo(UnitCell *cell, struct rvec *rlps, int rlp_count) struct TakeTwoCell ttCell; ttCell.cell = cell; ttCell.rotSymOps = NULL; + ttCell.twiz1Tmp = gsl_matrix_calloc(3, 3); + ttCell.twiz2Tmp = gsl_matrix_calloc(3, 3); + ttCell.vec1Tmp = gsl_vector_calloc(3); + ttCell.vec2Tmp = gsl_vector_calloc(3); ttCell.numOps = 0; success = generate_rotation_sym_ops(&ttCell); @@ -1499,6 +1491,9 @@ static UnitCell *run_taketwo(UnitCell *cell, struct rvec *rlps, int rlp_count) success = gen_observed_vecs(rlps, rlp_count, &obs_vecs, &obs_vec_count); if ( !success ) return NULL; + ttCell.obs_vecs = &obs_vecs; + ttCell.obs_vec_count = obs_vec_count; + success = match_obs_to_cell_vecs(asym_vecs, asym_vec_count, obs_vecs, obs_vec_count, 1); @@ -1510,7 +1505,7 @@ static UnitCell *run_taketwo(UnitCell *cell, struct rvec *rlps, int rlp_count) if ( !success ) return NULL; - find_seed(obs_vecs, obs_vec_count, &solution, &ttCell); + find_seed(&solution, &ttCell); if ( solution == NULL ) { cleanup_taketwo_obs_vecs(obs_vecs, obs_vec_count); -- cgit v1.2.3