From 22951c51e383e8a2f1a01d8bca60e4da9732633c Mon Sep 17 00:00:00 2001 From: Thomas White Date: Mon, 17 Oct 2016 17:29:16 +0200 Subject: Initial TakeTwo import Only changes from Helen's code so far: 1. Stripping trailing spaces 2. Tweaking includes (<> -> "") 3. Adding initial CrystFEL hooks at the bottom of taketwo.c 4. Moving definition of struct SpotVec to taketwo.c 5. Removing prototype for run_taketwo from taketwo.h (comment moved to .c) 6. Authorship/copyright boilerplate --- libcrystfel/Makefile.am | 6 +- libcrystfel/src/index.c | 24 +- libcrystfel/src/index.h | 6 + libcrystfel/src/taketwo.c | 819 ++++++++++++++++++++++++++++++++++++++++++++++ libcrystfel/src/taketwo.h | 41 +++ 5 files changed, 892 insertions(+), 4 deletions(-) create mode 100644 libcrystfel/src/taketwo.c create mode 100644 libcrystfel/src/taketwo.h diff --git a/libcrystfel/Makefile.am b/libcrystfel/Makefile.am index d7142f6a..7bb08685 100644 --- a/libcrystfel/Makefile.am +++ b/libcrystfel/Makefile.am @@ -10,7 +10,8 @@ libcrystfel_la_SOURCES = src/reflist.c src/utils.c src/cell.c src/detector.c \ src/render.c src/index.c src/dirax.c src/mosflm.c \ src/cell-utils.c src/integer_matrix.c src/crystal.c \ src/xds.c src/integration.c src/predict-refine.c \ - src/histogram.c src/events.c src/felix.c + src/histogram.c src/events.c src/felix.c \ + src/taketwo.c if HAVE_FFTW libcrystfel_la_SOURCES += src/asdf.c @@ -30,7 +31,8 @@ libcrystfel_la_include_HEADERS = ${top_srcdir}/version.h \ src/integer_matrix.h src/crystal.h \ src/xds.h src/predict-refine.h \ src/integration.h src/histogram.h \ - src/events.h src/asdf.h src/felix.h + src/events.h src/asdf.h src/felix.h \ + src/taketwo.h AM_CPPFLAGS = -DDATADIR=\""$(datadir)"\" -I$(top_builddir)/lib -Wall AM_CPPFLAGS += -I$(top_srcdir)/lib @LIBCRYSTFEL_CFLAGS@ diff --git a/libcrystfel/src/index.c b/libcrystfel/src/index.c index c8c86414..e0bebfa4 100644 --- a/libcrystfel/src/index.c +++ b/libcrystfel/src/index.c @@ -3,12 +3,12 @@ * * Perform indexing (somehow) * - * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2012-2016 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * Copyright © 2012 Lorenzo Galli * * Authors: - * 2010-2014 Thomas White + * 2010-2016 Thomas White * 2010-2011 Richard Kirian * 2012 Lorenzo Galli * 2013 Cornelius Gati @@ -56,6 +56,7 @@ #include "cell-utils.h" #include "felix.h" #include "predict-refine.h" +#include "taketwo.h" static int debug_index(struct image *image) @@ -117,6 +118,10 @@ IndexingPrivate **prepare_indexing(IndexingMethod *indm, UnitCell *cell, options); break; + case INDEXING_TAKETWO : + iprivs[n] = taketwo_prepare(&indm[n], cell, det, ltl); + break; + default : ERROR("Don't know how to prepare indexing method %i\n", indm[n]); @@ -193,6 +198,10 @@ void cleanup_indexing(IndexingMethod *indms, IndexingPrivate **privs) free(privs[n]); break; + case INDEXING_TAKETWO : + taketwo_cleanup(privs[n]); + break; + default : ERROR("Don't know how to clean up indexing method %i\n", indms[n]); @@ -268,6 +277,10 @@ static int try_indexer(struct image *image, IndexingMethod indm, r = felix_index(image, ipriv); break; + case INDEXING_TAKETWO : + r = taketwo_index(image, ipriv); + break; + default : ERROR("Unrecognised indexing method: %i\n", indm); return 0; @@ -568,6 +581,10 @@ char *indexer_str(IndexingMethod indm) strcpy(str, "xds"); break; + case INDEXING_TAKETWO : + strcpy(str, "taketwo"); + break; + case INDEXING_SIMULATION : strcpy(str, "simulation"); break; @@ -660,6 +677,9 @@ IndexingMethod *build_indexer_list(const char *str) } else if ( strcmp(methods[i], "felix") == 0) { list[++nmeth] = INDEXING_DEFAULTS_FELIX; + } else if ( strcmp(methods[i], "taketwo") == 0) { + list[++nmeth] = INDEXING_DEFAULTS_TAKETWO; + } else if ( strcmp(methods[i], "none") == 0) { list[++nmeth] = INDEXING_NONE; diff --git a/libcrystfel/src/index.h b/libcrystfel/src/index.h index 2fb5a13d..6ac9792a 100644 --- a/libcrystfel/src/index.h +++ b/libcrystfel/src/index.h @@ -58,6 +58,10 @@ | INDEXING_USE_CELL_PARAMETERS \ | INDEXING_RETRY | INDEXING_REFINE) +#define INDEXING_DEFAULTS_TAKETWO (INDEXING_TAKETWO \ + | INDEXING_USE_CELL_PARAMETERS \ + | INDEXING_RETRY | INDEXING_REFINE) + /* Axis check is needed for XDS, because it likes to permute the axes */ #define INDEXING_DEFAULTS_XDS (INDEXING_XDS | INDEXING_USE_LATTICE_TYPE \ | INDEXING_USE_CELL_PARAMETERS \ @@ -75,6 +79,7 @@ * @INDEXING_SIMULATION: Dummy value * @INDEXING_DEBUG: Results injector for debugging * @INDEXING_ASDF: Use in-built "asdf" indexer + * @INDEXING_TAKETWO: Use in-built "taketwo" indexer * @INDEXING_CHECK_CELL_COMBINATIONS: Check linear combinations of unit cell * axes for agreement with given cell. * @INDEXING_CHECK_CELL_AXES: Check unit cell axes for agreement with given @@ -107,6 +112,7 @@ typedef enum { INDEXING_SIMULATION = 6, INDEXING_DEBUG = 7, INDEXING_ASDF = 8, + INDEXING_TAKETWO = 9, /* Bits at the top of the IndexingMethod are flags which modify the * behaviour of the indexer. */ diff --git a/libcrystfel/src/taketwo.c b/libcrystfel/src/taketwo.c new file mode 100644 index 00000000..32300113 --- /dev/null +++ b/libcrystfel/src/taketwo.c @@ -0,0 +1,819 @@ +/* + * taketwo.c + * + * Rewrite of TakeTwo algorithm (Acta D72 (8) 956-965) for CrystFEL + * + * Copyright © 2016 Helen Ginn + * Copyright © 2016 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. + * + * Authors: + * 2016 Helen Ginn + * 2016 Thomas White + * + * This file is part of CrystFEL. + * + * CrystFEL is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * CrystFEL is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with CrystFEL. If not, see . + * + */ + +#include +#include + +#include "cell-utils.h" +#include "index.h" +#include "taketwo.h" + +/** + * spotvec + * @obsvec: an observed vector between two spots + * @matches: array of matching theoretical vectors from unit cell + * @match_num: number of matches + * @distance: length of obsvec (do I need this?) + * @her_rlp: pointer to first rlp position for difference vec + * @his_rlp: pointer to second rlp position for difference vec + * + * Structure representing 3D vector between two potential Bragg peaks + * in reciprocal space, and an array of potential matching theoretical + * vectors from unit cell/centering considerations. + **/ +struct SpotVec +{ + struct rvec obsvec; + struct rvec *matches; + int match_num; + double distance; + struct rvec *her_rlp; + struct rvec *his_rlp; +}; + +/* Maximum distance between two rlp sizes to consider info for indexing */ +#define MAX_RECIP_DISTANCE 0.15 + +/* Tolerance for two lengths in reciprocal space to be considered the same */ +#define RECIP_TOLERANCE 0.001 + +/* Threshold for network members to consider a potential solution */ +#define NETWORK_MEMBER_THRESHOLD 20 + +/* Maximum network members (obviously a solution so should stop) */ +#define MAX_NETWORK_MEMBERS 100 + +/* Maximum dead ends for a single branch extension during indexing */ +#define MAX_DEAD_ENDS 5 + +/* FIXME: Tolerance for two angles to be considered the same - is there some + * pre-existing degrees-to-radians in CrystFEL to use here or not? + */ +#define ANGLE_TOLERANCE 1.0 * M_PI / 180 + +/** TODO: + * + * - May need to be capable of playing with the tolerances/#defined stuff. + */ + + +/* ------------------------------------------------------------------------ + * apologetic function + * ------------------------------------------------------------------------*/ + +void apologise() +{ + printf("Error - could not allocate memory for indexing.\n"); +} + +/* ------------------------------------------------------------------------ + * functions concerning aspects of rvec which are very likely to be + * incorporated somewhere else in CrystFEL and therefore may need to be + * deleted and references connected to a pre-existing function. (Lowest level) + * ------------------------------------------------------------------------*/ + +static struct rvec new_rvec(double new_u, double new_v, double new_w) +{ + struct rvec new_rvector; + new_rvector.u = new_u; + new_rvector.v = new_v; + new_rvector.w = new_w; + + return new_rvector; +} + +static struct rvec diff_vec(struct rvec from, struct rvec to) +{ + struct rvec diff = new_rvec(to.u - from.u, + to.v - from.v, + to.w - from.w); + + return diff; +} + +static double sq_length(struct rvec vec) +{ + double sqlength = (vec.u * vec.u + vec.v * vec.v + vec.w * vec.w); + + return sqlength; +} + +static double rvec_length(struct rvec vec) +{ + return sqrt(sq_length(vec)); +} + +static void normalise_rvec(struct rvec *vec) +{ + double length = rvec_length(*vec); + vec->u /= length; + vec->v /= length; + vec->w /= length; +} + +static double rvec_cosine(struct rvec v1, struct rvec v2) +{ + double dot_prod = v1.u * v2.u + v1.v * v2.v + v1.w * v2.w; + double v1_length = rvec_length(v1); + double v2_length = rvec_length(v2); + + double cos_theta = dot_prod / (v1_length * v2_length); + + return cos_theta; +} + +static double rvec_angle(struct rvec v1, struct rvec v2) +{ + double cos_theta = rvec_cosine(v1, v2); + double angle = acos(cos_theta); + + return angle; +} + +/* ------------------------------------------------------------------------ + * functions called under the core functions, still specialised (Level 3) + * ------------------------------------------------------------------------*/ + +/* Rotate vector (vec1) around axis (axis) by angle theta. Find value of + * theta for which the angle between (vec1) and (vec2) is minimised. + * Behold! Finally an analytical solution for this one. Assuming + * that @result has already been allocated. Will upload the maths to the + * shared Google drive. */ +static int closest_rot_mat(struct rvec vec1, struct rvec vec2, + struct rvec axis, gsl_matrix *result) +{ + /* Let's have unit vectors */ + normalise_rvec(&vec1); + normalise_rvec(&vec2); + normalise_rvec(&axis); + + /* Redeclaring these to try and maintain readability and + * check-ability against the maths I wrote down */ + double a = vec2.u; double b = vec2.w; double c = vec2.v; + double p = vec1.u; double q = vec1.w; double r = vec1.v; + double x = axis.u; double y = axis.w; double z = axis.v; + + /* Components in handwritten maths online when I upload it */ + double A = a*(p*x*x - p + x*y*q + x*z*r) + + b*(p*x*y + q*y*y - q + r*y*z) + + c*(p*x*z + q*y*z + r*z*z - r); + + double B = a*(y*r - z*q) + b*(p*z - r*x) + c*(q*x - p*y); + + double tan_theta = - B / A; + + /* Not sure why I always have to take the + M_PI solution. Work + * this one out. But it always works!? */ + double theta = atan(tan_theta) + M_PI; + + /* FIXME: fill in the gsl_matrix *result with an identity matrix + * which has subsequently been rotated by theta around the axis. + * Help from Tom, probably. */ + + result = result; + + return 1; +} + +static int rot_mats_are_similar(gsl_matrix *rot1, gsl_matrix *rot2) +{ + /* FIXME: write me! + * Write code for that fancy algorithm here from XPLOR */ + + return 0; +} + +/* Code me in gsl_matrix language to copy the contents of the function + * in cppxfel (IndexingSolution::createSolution). This function is quite + * intensive on the number crunching side so simple angle checks are used + * to 'pre-scan' vectors beforehand. */ +static int generate_rot_mat(struct rvec obs1, struct rvec obs2, + struct rvec cell1, struct rvec cell2, + gsl_matrix *mat) +{ + /* FIXME: Write this code - assume *mat has already been allocated + * and insert a CrystFEL-friendly rotation matrix. Help from Tom. */ + return 1; +} + +static int obs_vecs_share_spot(struct SpotVec *her_obs, struct SpotVec *his_obs) +{ + /* FIXME: Disgusting... can I tone this down a bit? */ + if ( (her_obs->her_rlp == his_obs->her_rlp) || + (her_obs->her_rlp == his_obs->his_rlp) || + (her_obs->his_rlp == his_obs->her_rlp) || + (her_obs->his_rlp == his_obs->his_rlp) ) { + return 1; + } + + return 0; +} + +static int obs_shares_spot_w_array(struct SpotVec *obs_vecs, int test_idx, + int *members[MAX_NETWORK_MEMBERS], int num) +{ + int i; + struct SpotVec *her_obs = &obs_vecs[test_idx]; + + for ( i=0; iobsvec, his_obs->obsvec); + + /* 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? */ + + double angle_diff = fabs(theory_angle - obs_angle); + + if ( angle_diff < ANGLE_TOLERANCE ) { + *her_match_idx = i; + *his_match_idx = j; + + return 1; + } + } + } + + return 0; +} + +static int obs_angles_match_array(struct SpotVec *obs_vecs, int test_idx, + int *members[MAX_NETWORK_MEMBERS], int num) +{ + /* note: this is just a preliminary check to reduce unnecessary + * computation later down the line, but is not entirely accurate. + * For example, I have not checked that the 'matching cell vector' + * is identical - too much faff. + **/ + + int i; + struct SpotVec *her_obs = &obs_vecs[test_idx]; + + for ( i=0; i MAX_DEAD_ENDS ) { + break; + } + + /* We have not had too many dead ends. Try removing + the last member and continue. */ + int failed = members[member_num - 1]; + start = failed + 1; + member_num--; + dead_ends++; + + continue; + } + + /* we have elongated membership - so reset dead_ends counter */ + dead_ends = 0; + + members[member_num] = next_index; + start = next_index + 1; + member_num++; + + /* If member_num is high enough, we want to return a yes */ + + if ( member_num > NETWORK_MEMBER_THRESHOLD ) { + break; + } + } + + /* Deal with this shit after coffee */ + + return ( member_num > NETWORK_MEMBER_THRESHOLD ); +} + +static int find_seed_and_network(struct SpotVec *obs_vecs, int obs_vec_count, + gsl_matrix **rotation) +{ + /* 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 RECIP_TOLERANCE ) { + continue; + } + + /* we have a match, add to array! */ + + count++; + int new_size = count*sizeof(struct rvec); + struct rvec *temp_matches; + temp_matches = realloc(obs_vecs[i].matches, new_size); + + if ( temp_matches == NULL ) { + return 0; + } else { + obs_vecs[i].matches = temp_matches; + temp_matches[count - 1] = cell_vecs[j]; + obs_vecs[i].match_num = count; + } + } + } + + return 1; +} + +static int gen_observed_vecs(struct rvec *rlps, int rlp_count, + struct SpotVec **obs_vecs, int *obs_vec_count) +{ + int i, j; + int count = 0; + + /* maximum distance squared for comparisons */ + double max_sq_length = pow(MAX_RECIP_DISTANCE, 2); + + /* Indentation... bending the rules a bit? */ + for ( i=0; i max_sq_length ) { + continue; + } + + count++; + + struct SpotVec *temp_obs_vecs; + temp_obs_vecs = realloc(*obs_vecs, + count*sizeof(struct SpotVec)); + + if ( temp_obs_vecs == NULL ) { + return 0; + } else { + *obs_vecs = temp_obs_vecs; + + /* initialise all SpotVec struct members */ + + struct SpotVec spot_vec; + spot_vec.obsvec = diff; + spot_vec.distance = sqrt(sqlength); + spot_vec.matches = NULL; + spot_vec.match_num = 0; + spot_vec.her_rlp = &rlps[i]; + spot_vec.his_rlp = &rlps[j]; + + (*obs_vecs)[count - 1] = spot_vec; + } + } + } + + *obs_vec_count = count; + + return 1; +} + +static int gen_theoretical_vecs(UnitCell *cell, struct rvec **cell_vecs, + int *vec_count) +{ + double a, b, c, alpha, beta, gamma; + int h_max, k_max, l_max; + + cell_get_parameters(cell, &a, &b, &c, &alpha, &beta, &gamma); + + /* find maximum Miller (h, k, l) indices for a given resolution */ + h_max = MAX_RECIP_DISTANCE / a + 1; + k_max = MAX_RECIP_DISTANCE / b + 1; + l_max = MAX_RECIP_DISTANCE / c + 1; + + int h, k, l; + int count = 0; + + for ( h=-h_max; h<=+h_max; h++ ) { + for ( k=-k_max; k<=+k_max; k++ ) { + for ( l=-l_max; l<=+l_max; l++ ) { + + /* Exclude systematic absences from centering concerns */ + if ( forbidden_reflection(cell, h, k, l) ) { + continue; + } + + struct rvec cell_vec = new_rvec(h, k, l); + + /* FIXME: transform int (h, k, l) to reciprocal coordinates. + * Don't want to do this manually if there is already a + * function in libcrystfel to do this. Would like to map + * Miller index (5, 13, 2) onto (0.05, 0.13, 0.02) for example, + * if I had a 100 x 100 x 100 Å cubic cell. */ + + cell_vec = cell_vec; + + /* Assumption that "cell_vec" now has transformed coords */ + + /* add this to our array - which may require expanding */ + count++; + + struct rvec *temp_cell_vecs; + temp_cell_vecs = realloc(*cell_vecs, count*sizeof(struct rvec)); + + if ( temp_cell_vecs == NULL ) { + return 0; + } else { + *cell_vecs = temp_cell_vecs; + (*cell_vecs)[count - 1] = cell_vec; + } + } + } + } + + *vec_count = count; + + return 1; +} + +static void generate_basis_vectors(UnitCell *cell, gsl_matrix *rot, + struct rvec *a_star, struct rvec *b_star, + struct rvec *c_star) +{ + /* FIXME: more matrix stuff - multiply cell matrix by rotation matrix + * and extract the reciprocal axes from the definition of the matrix. + */ +} + +/* ------------------------------------------------------------------------ + * cleanup functions - called from run_taketwo(). + * ------------------------------------------------------------------------*/ + +static int cleanup_taketwo_cell_vecs(struct rvec *cell_vecs) +{ + free(cell_vecs); +} + +static int cleanup_taketwo_obs_vecs(struct SpotVec *obs_vecs, + int obs_vec_count) +{ + for ( int i=0; iltl = ltl; + tp->template = cell; + tp->indm = *indm; + + return (IndexingPrivate *)tp; +} + + +void taketwo_cleanup(IndexingPrivate *pp) +{ + struct taketwo_private *tp = (struct taketwo_private *)pp; + free(tp); +} diff --git a/libcrystfel/src/taketwo.h b/libcrystfel/src/taketwo.h new file mode 100644 index 00000000..35f7cbe3 --- /dev/null +++ b/libcrystfel/src/taketwo.h @@ -0,0 +1,41 @@ +/* + * taketwo.h + * + * Rewrite of TakeTwo algorithm (Acta D72 (8) 956-965) for CrystFEL + * + * Copyright © 2016 Helen Ginn + * Copyright © 2016 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. + * + * Authors: + * 2016 Helen Ginn + * 2016 Thomas White + * + * This file is part of CrystFEL. + * + * CrystFEL is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * CrystFEL is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with CrystFEL. If not, see . + * + */ + +#ifndef TAKETWO_H +#define TAKETWO_H + +#include "cell.h" + +extern IndexingPrivate *taketwo_prepare(IndexingMethod *indm, UnitCell *cell, + struct detector *det, float *ltl); +extern int taketwo_index(struct image *image, IndexingPrivate *ipriv); +extern void taketwo_cleanup(IndexingPrivate *pp); + +#endif /* TAKETWO_H */ -- cgit v1.2.3