From fb1832f80cebe48f8bf066f37973058b7ce3e334 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Tue, 4 Mar 2014 15:09:45 +0100 Subject: Introduce "ambigator" --- .gitignore | 1 + Makefile.am | 5 +- src/ambigator.c | 406 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 411 insertions(+), 1 deletion(-) create mode 100644 src/ambigator.c diff --git a/.gitignore b/.gitignore index 40000069..5091a898 100644 --- a/.gitignore +++ b/.gitignore @@ -34,6 +34,7 @@ src/partialator src/check_hkl src/partial_sim src/cell_explorer +src/ambigator src/.dirstamp *~ doc/reference/libcrystfel/* diff --git a/Makefile.am b/Makefile.am index ad39c1ee..0a187d94 100644 --- a/Makefile.am +++ b/Makefile.am @@ -4,7 +4,8 @@ SUBDIRS = lib doc/reference/libcrystfel libcrystfel ACLOCAL_AMFLAGS = -I m4 bin_PROGRAMS = src/pattern_sim src/process_hkl src/get_hkl src/indexamajig \ - src/compare_hkl src/partialator src/check_hkl src/partial_sim + src/compare_hkl src/partialator src/check_hkl src/partial_sim \ + src/ambigator noinst_PROGRAMS = tests/list_check tests/integration_check \ tests/pr_p_gradient_check tests/symmetry_check \ @@ -85,6 +86,8 @@ endif src_partialator_SOURCES = src/partialator.c src/post-refinement.c \ src/hrs-scaling.c +src_ambigator_SOURCES = src/ambigator.c + if HAVE_CAIRO if HAVE_PANGO if HAVE_PANGOCAIRO diff --git a/src/ambigator.c b/src/ambigator.c new file mode 100644 index 00000000..1e417ce2 --- /dev/null +++ b/src/ambigator.c @@ -0,0 +1,406 @@ +/* + * ambigator.c + * + * Resolve indexing ambiguities + * + * Copyright © 2014 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. + * Copyright © 2014 Wolfgang Brehm + * + * Authors: + * 2014 Thomas White + * 2014 Wolfgang Brehm + * + * 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 . + * + */ + + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "post-refinement.h" +#include "hrs-scaling.h" +#include "scaling-report.h" + + +static void show_help(const char *s) +{ + printf("Syntax: %s [options]\n\n", s); + printf( +"Resolve indexing ambiguities.\n" +"\n" +" -h, --help Display this help message.\n" +"\n" +" -i, --input= Input stream.\n" +" -o, --output= Output stream.\n" +" -y, --symmetry= Apparent (\"source\") symmetry.\n" +" -e Actual (\"target\") symmetry.\n" +" -n, --iterations= Iterate times.\n" +); +} + + +static RefList *asymm_and_merge(RefList *in, const SymOpList *sym) +{ + Reflection *refl; + RefListIterator *iter; + RefList *asym; + + asym = reflist_new(); + if ( asym == NULL ) return NULL; + + for ( refl = first_refl(in, &iter); + refl != NULL; + refl = next_refl(refl, iter) ) + { + signed int h, k, l; + signed int ha, ka, la; + Reflection *cr; + + get_indices(refl, &h, &k, &l); + + get_asymm(sym, h, k, l, &ha, &ka, &la); + + cr = find_refl(asym, ha, ka, la); + if ( cr == NULL ) { + cr = add_refl(asym, ha, ka, la); + assert(cr != NULL); + copy_data(cr, refl); + } else { + const double i = get_intensity(cr); + const int r = get_redundancy(cr); + set_intensity(cr, (r*i + get_intensity(refl))/(r+1)); + set_redundancy(cr, r+1); + } + } + + return asym; +} + + +static float corr(Crystal *a, Crystal *b, int *pn) +{ + Reflection *aref; + RefListIterator *iter; + float s_xy = 0.0; + float s_x = 0.0; + float s_y = 0.0; + float s_x2 = 0.0; + float s_y2 = 0.0; + int n = 0; + float t1, t2; + + for ( aref = first_refl(crystal_get_reflections(a), &iter); + aref != NULL; + aref = next_refl(aref, iter) ) + { + signed int h, k, l; + Reflection *bref; + float aint, bint; + + get_indices(aref, &h, &k, &l); + + bref = find_refl(crystal_get_reflections(b), h, k, l); + if ( bref == NULL ) continue; + + aint = get_intensity(aref); + bint = get_intensity(bref); + + s_xy += aint*bint; + s_x += aint; + s_y += bint; + s_x2 += aint*aint; + s_y2 += bint*bint; + n++; + + } + + *pn = n; + + t1 = s_x2 - s_x*s_x / n; + t2 = s_y2 - s_y*s_y / n; + + if ( (t1 < 0.0) || (t2 <= 0.0) ) return 0.0; + + return (s_xy - s_x*s_y)/n/sqrt(t1*t2); +} + + +static void detwin(Crystal **crystals, int n_crystals, SymOpList *amb, + int *assignments) +{ + int i; + int nch = 0; + + for ( i=0; i g ) { + assignments[i] = 1 - assignments[i]; + nch++; + } + + progress_bar(i, n_crystals-1, "Calculating"); + + } + + STATUS("Changed %i assignments this time.\n", nch); +} + + +int main(int argc, char *argv[]) +{ + int c; + char *infile = NULL; + char *outfile = NULL; + char *s_sym_str = NULL; + SymOpList *s_sym; + char *e_sym_str = NULL; + SymOpList *e_sym; + SymOpList *amb; + int n_iter = 1; + int n_crystals, n_chunks, max_crystals; + Crystal **crystals; + Stream *st; + int i; + int *assignments; + + /* Long options */ + const struct option longopts[] = { + {"help", 0, NULL, 'h'}, + {"input", 1, NULL, 'i'}, + {"output", 1, NULL, 'o'}, + {"symmetry", 1, NULL, 'y'}, + {"iterations", 1, NULL, 'n'}, + + {0, 0, NULL, 0} + }; + + /* Short options */ + while ((c = getopt_long(argc, argv, "hi:o:y:n:e:", + longopts, NULL)) != -1) + { + + switch (c) { + + case 'h' : + show_help(argv[0]); + return 0; + + case 'i' : + infile = strdup(optarg); + break; + + case 'o' : + outfile = strdup(optarg); + break; + + case 'y' : + s_sym_str = strdup(optarg); + break; + + case 'e' : + e_sym_str = strdup(optarg); + break; + + case 'n' : + n_iter = atoi(optarg); + break; + + case 0 : + break; + + case '?' : + break; + + default : + ERROR("Unhandled option '%c'\n", c); + break; + + } + + } + + if ( infile == NULL ) { + infile = strdup("-"); + } + st = open_stream_for_read(infile); + if ( st == NULL ) { + ERROR("Failed to open input stream '%s'\n", infile); + return 1; + } + free(infile); + + /* Sanitise output filename */ + if ( outfile == NULL ) { + outfile = strdup("partialator.hkl"); + } + + if ( s_sym_str == NULL ) { + ERROR("You must specify the input symmetry (with -y)\n"); + return 1; + } + s_sym = get_pointgroup(s_sym_str); + free(s_sym_str); + + if ( e_sym_str == NULL ) { + e_sym = NULL; + amb = NULL; + } else { + e_sym = get_pointgroup(e_sym_str); + free(e_sym_str); + if ( e_sym == NULL ) return 1; + amb = get_ambiguities(s_sym, e_sym); + if ( amb == NULL ) return 1; + STATUS("Ambiguity operations:\n"); + describe_symmetry(amb); + } + + crystals = NULL; + n_crystals = 0; + max_crystals = 0; + n_chunks = 0; + do { + + struct image cur; + int i; + + cur.det = NULL; + + if ( read_chunk(st, &cur) != 0 ) { + break; + } + + image_feature_list_free(cur.features); + + for ( i=0; i Date: Tue, 4 Mar 2014 15:25:59 +0100 Subject: Fussiness --- src/ambigator.c | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/ambigator.c b/src/ambigator.c index 1e417ce2..d86dc63c 100644 --- a/src/ambigator.c +++ b/src/ambigator.c @@ -146,7 +146,6 @@ static float corr(Crystal *a, Crystal *b, int *pn) s_x2 += aint*aint; s_y2 += bint*bint; n++; - } *pn = n; @@ -156,7 +155,7 @@ static float corr(Crystal *a, Crystal *b, int *pn) if ( (t1 < 0.0) || (t2 <= 0.0) ) return 0.0; - return (s_xy - s_x*s_y)/n/sqrt(t1*t2); + return ((s_xy - s_x*s_y)/n)/sqrt(t1*t2); } @@ -181,8 +180,8 @@ static void detwin(Crystal **crystals, int n_crystals, SymOpList *amb, cc = corr(crystals[i], crystals[j], &n); - if ( i == j ) continue; if ( n < 3 ) continue; + if ( i == j ) continue; if ( assignments[i] == assignments[j] ) { f += cc; -- cgit v1.2.3 From 9bb0d608b82f731356a58c24637522a61105ba29 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Tue, 4 Mar 2014 15:26:03 +0100 Subject: Report the number of assignments changed since the start --- src/ambigator.c | 23 +++++++++++++++++++++-- 1 file changed, 21 insertions(+), 2 deletions(-) diff --git a/src/ambigator.c b/src/ambigator.c index d86dc63c..8a3879e3 100644 --- a/src/ambigator.c +++ b/src/ambigator.c @@ -221,10 +221,13 @@ int main(int argc, char *argv[]) SymOpList *amb; int n_iter = 1; int n_crystals, n_chunks, max_crystals; + int n_dif; Crystal **crystals; Stream *st; int i; int *assignments; + int *orig_assignments; + gsl_rng *rng; /* Long options */ const struct option longopts[] = { @@ -384,22 +387,38 @@ int main(int argc, char *argv[]) } - assignments = calloc(n_crystals, sizeof(int)); + assignments = malloc(n_crystals*sizeof(int)); if ( assignments == NULL ) { ERROR("Couldn't allocate memory for assignments.\n"); return 1; } + orig_assignments = malloc(n_crystals*sizeof(int)); + if ( orig_assignments == NULL ) { + ERROR("Couldn't allocate memory for original assignments.\n"); + return 1; + } + + rng = gsl_rng_alloc(gsl_rng_mt19937); for ( i=0; i 0.5); + orig_assignments[i] = assignments[i]; } for ( i=0; i Date: Tue, 4 Mar 2014 15:42:51 +0100 Subject: D'oh --- src/ambigator.c | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/ambigator.c b/src/ambigator.c index 8a3879e3..c169d63c 100644 --- a/src/ambigator.c +++ b/src/ambigator.c @@ -400,7 +400,7 @@ int main(int argc, char *argv[]) } rng = gsl_rng_alloc(gsl_rng_mt19937); - for ( i=0; i 0.5); orig_assignments[i] = assignments[i]; } @@ -410,7 +410,7 @@ int main(int argc, char *argv[]) } n_dif = 0; - for ( i=0; i Date: Tue, 4 Mar 2014 16:37:44 +0100 Subject: Use faster data structure for list during correlation --- src/ambigator.c | 130 +++++++++++++++++++++++++++++++++++++++----------------- 1 file changed, 90 insertions(+), 40 deletions(-) diff --git a/src/ambigator.c b/src/ambigator.c index c169d63c..34f68742 100644 --- a/src/ambigator.c +++ b/src/ambigator.c @@ -73,12 +73,26 @@ static void show_help(const char *s) ); } +#define SERIAL(h, k, l) ((((h)+512)<<20) + (((k)+512)<<10) + ((l)+512)) +#define GET_H(serial) ((((serial) & 0x3ff00000)>>20)-512) +#define GET_K(serial) ((((serial) & 0x000ffc00)>>10)-512) +#define GET_L(serial) (((serial) & 0x000003ff)-512) -static RefList *asymm_and_merge(RefList *in, const SymOpList *sym) +struct flist +{ + int n; + unsigned int *s; + float *i; +}; + + +static struct flist *asymm_and_merge(RefList *in, const SymOpList *sym) { Reflection *refl; RefListIterator *iter; RefList *asym; + struct flist *f; + int n; asym = reflist_new(); if ( asym == NULL ) return NULL; @@ -108,14 +122,40 @@ static RefList *asymm_and_merge(RefList *in, const SymOpList *sym) } } - return asym; + f = malloc(sizeof(struct flist)); + if ( f == NULL ) { + ERROR("Failed to allocate flist\n"); + return NULL; + } + + n = num_reflections(asym); + f->s = malloc(n*sizeof(unsigned int)); + f->i = malloc(n*sizeof(float)); + if ( (f->s == NULL) || (f->i == NULL) ) { + ERROR("Failed to allocate flist\n"); + return NULL; + } + + f->n = 0; + for ( refl = first_refl(asym, &iter); + refl != NULL; + refl = next_refl(refl, iter) ) + { + signed int h, k, l; + + get_indices(refl, &h, &k, &l); + f->s[f->n] = SERIAL(h, k, l); + f->i[f->n] = get_intensity(refl); + f->n++; + } + reflist_free(asym); + + return f; } -static float corr(Crystal *a, Crystal *b, int *pn) +static float corr(struct flist *a, struct flist *b, int *pn) { - Reflection *aref; - RefListIterator *iter; float s_xy = 0.0; float s_x = 0.0; float s_y = 0.0; @@ -123,31 +163,48 @@ static float corr(Crystal *a, Crystal *b, int *pn) float s_y2 = 0.0; int n = 0; float t1, t2; + int ap = 0; + int bp = 0; + int done = 0; - for ( aref = first_refl(crystal_get_reflections(a), &iter); - aref != NULL; - aref = next_refl(aref, iter) ) - { - signed int h, k, l; - Reflection *bref; - float aint, bint; + while ( 1 ) { - get_indices(aref, &h, &k, &l); + while ( a->s[ap] > b->s[bp] ) { + if ( ++bp == b->n ) { + done = 1; + break; + } + } + if ( done ) break; - bref = find_refl(crystal_get_reflections(b), h, k, l); - if ( bref == NULL ) continue; + while ( a->s[ap] < b->s[bp] ) { + if ( ++ap == a->n ) { + done = 1; + break; + } + } + if ( done ) break; - aint = get_intensity(aref); - bint = get_intensity(bref); + if ( a->s[ap] == b->s[bp] ) { - s_xy += aint*bint; - s_x += aint; - s_y += bint; - s_x2 += aint*aint; - s_y2 += bint*bint; - n++; - } + float aint, bint; + + aint = a->i[ap]; + bint = b->i[bp]; + + s_xy += aint*bint; + s_x += aint; + s_y += bint; + s_x2 += aint*aint; + s_y2 += bint*bint; + n++; + + if ( ++ap == a->n ) break; + if ( ++bp == b->n ) break; + + } + } *pn = n; t1 = s_x2 - s_x*s_x / n; @@ -159,7 +216,7 @@ static float corr(Crystal *a, Crystal *b, int *pn) } -static void detwin(Crystal **crystals, int n_crystals, SymOpList *amb, +static void detwin(struct flist **crystals, int n_crystals, SymOpList *amb, int *assignments) { int i; @@ -222,7 +279,7 @@ int main(int argc, char *argv[]) int n_iter = 1; int n_crystals, n_chunks, max_crystals; int n_dif; - Crystal **crystals; + struct flist **crystals; Stream *st; int i; int *assignments; @@ -340,6 +397,7 @@ int main(int argc, char *argv[]) for ( i=0; i Date: Tue, 4 Mar 2014 17:26:52 +0100 Subject: Add --highres and --lowres --- src/ambigator.c | 42 ++++++++++++++++++++++++++++++++++++++---- 1 file changed, 38 insertions(+), 4 deletions(-) diff --git a/src/ambigator.c b/src/ambigator.c index 34f68742..6580b97f 100644 --- a/src/ambigator.c +++ b/src/ambigator.c @@ -51,6 +51,7 @@ #include #include #include +#include #include "post-refinement.h" #include "hrs-scaling.h" @@ -70,6 +71,8 @@ static void show_help(const char *s) " -y, --symmetry= Apparent (\"source\") symmetry.\n" " -e Actual (\"target\") symmetry.\n" " -n, --iterations= Iterate times.\n" +" --highres= High resolution cutoff in A.\n" +" --lowres= Low resolution cutoff in A.\n" ); } @@ -86,7 +89,8 @@ struct flist }; -static struct flist *asymm_and_merge(RefList *in, const SymOpList *sym) +static struct flist *asymm_and_merge(RefList *in, const SymOpList *sym, + UnitCell *cell, double rmin, double rmax) { Reflection *refl; RefListIterator *iter; @@ -104,9 +108,14 @@ static struct flist *asymm_and_merge(RefList *in, const SymOpList *sym) signed int h, k, l; signed int ha, ka, la; Reflection *cr; + double res; get_indices(refl, &h, &k, &l); + res = 2.0*resolution(cell, h, k, l); + if ( res < rmin ) continue; + if ( res > rmax ) continue; + get_asymm(sym, h, k, l, &ha, &ka, &la); cr = find_refl(asym, ha, ka, la); @@ -285,6 +294,9 @@ int main(int argc, char *argv[]) int *assignments; int *orig_assignments; gsl_rng *rng; + float highres, lowres; + double rmin = 0.0; /* m^-1 */ + double rmax = INFINITY; /* m^-1 */ /* Long options */ const struct option longopts[] = { @@ -294,6 +306,9 @@ int main(int argc, char *argv[]) {"symmetry", 1, NULL, 'y'}, {"iterations", 1, NULL, 'n'}, + {"highres", 1, NULL, 2}, + {"lowres", 1, NULL, 3}, + {0, 0, NULL, 0} }; @@ -328,6 +343,22 @@ int main(int argc, char *argv[]) n_iter = atoi(optarg); break; + case 2 : + if ( sscanf(optarg, "%e", &highres) != 1 ) { + ERROR("Invalid value for --highres\n"); + return 1; + } + rmax = 1.0 / (highres/1e10); + break; + + case 3 : + if ( sscanf(optarg, "%e", &lowres) != 1 ) { + ERROR("Invalid value for --lowres\n"); + return 1; + } + rmin = 1.0 / (lowres/1e10); + break; + case 0 : break; @@ -398,10 +429,10 @@ int main(int argc, char *argv[]) Crystal *cr; RefList *list; + UnitCell *cell; cr = cur.crystals[i]; - - cell_free(crystal_get_cell(cr)); + cell = crystal_get_cell(cr); if ( n_crystals == max_crystals ) { @@ -422,7 +453,10 @@ int main(int argc, char *argv[]) } list = crystal_get_reflections(cr); - crystals[n_crystals] = asymm_and_merge(list, s_sym); + crystals[n_crystals] = asymm_and_merge(list, s_sym, + cell, + rmin, rmax); + cell_free(cell); n_crystals++; reflist_free(list); -- cgit v1.2.3 From b1f4e521baf644069f02dd9efe869e63aed9225e Mon Sep 17 00:00:00 2001 From: Thomas White Date: Tue, 4 Mar 2014 18:35:04 +0100 Subject: Fix the algorithm --- src/ambigator.c | 17 ++++++++++++----- 1 file changed, 12 insertions(+), 5 deletions(-) diff --git a/src/ambigator.c b/src/ambigator.c index 6580b97f..e18636ff 100644 --- a/src/ambigator.c +++ b/src/ambigator.c @@ -219,9 +219,9 @@ static float corr(struct flist *a, struct flist *b, int *pn) t1 = s_x2 - s_x*s_x / n; t2 = s_y2 - s_y*s_y / n; - if ( (t1 < 0.0) || (t2 <= 0.0) ) return 0.0; + if ( (t1 <= 0.0) || (t2 <= 0.0) ) return 0.0; - return ((s_xy - s_x*s_y)/n)/sqrt(t1*t2); + return (s_xy - s_x*s_y/n) / sqrt(t1*t2); } @@ -230,6 +230,8 @@ static void detwin(struct flist **crystals, int n_crystals, SymOpList *amb, { int i; int nch = 0; + float mf = 0.0; + int nmf = 0; for ( i=0; i g ) { + if ( f < g ) { assignments[i] = 1 - assignments[i]; nch++; } @@ -271,7 +278,7 @@ static void detwin(struct flist **crystals, int n_crystals, SymOpList *amb, } - STATUS("Changed %i assignments this time.\n", nch); + STATUS("Mean f = %f, changed %i assignments this time.\n", mf/nmf, nch); } -- cgit v1.2.3 From 9127ee2452cf5251e257ff338b026bf633ba16f7 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Tue, 4 Mar 2014 21:19:38 +0100 Subject: Don't use "-i" --- src/ambigator.c | 20 ++++++++------------ 1 file changed, 8 insertions(+), 12 deletions(-) diff --git a/src/ambigator.c b/src/ambigator.c index e18636ff..2af39e8a 100644 --- a/src/ambigator.c +++ b/src/ambigator.c @@ -60,13 +60,12 @@ static void show_help(const char *s) { - printf("Syntax: %s [options]\n\n", s); + printf("Syntax: %s [options] input.stream\n\n", s); printf( "Resolve indexing ambiguities.\n" "\n" " -h, --help Display this help message.\n" "\n" -" -i, --input= Input stream.\n" " -o, --output= Output stream.\n" " -y, --symmetry= Apparent (\"source\") symmetry.\n" " -e Actual (\"target\") symmetry.\n" @@ -285,7 +284,7 @@ static void detwin(struct flist **crystals, int n_crystals, SymOpList *amb, int main(int argc, char *argv[]) { int c; - char *infile = NULL; + const char *infile; char *outfile = NULL; char *s_sym_str = NULL; SymOpList *s_sym; @@ -308,7 +307,6 @@ int main(int argc, char *argv[]) /* Long options */ const struct option longopts[] = { {"help", 0, NULL, 'h'}, - {"input", 1, NULL, 'i'}, {"output", 1, NULL, 'o'}, {"symmetry", 1, NULL, 'y'}, {"iterations", 1, NULL, 'n'}, @@ -320,7 +318,7 @@ int main(int argc, char *argv[]) }; /* Short options */ - while ((c = getopt_long(argc, argv, "hi:o:y:n:e:", + while ((c = getopt_long(argc, argv, "ho:y:n:e:", longopts, NULL)) != -1) { @@ -330,10 +328,6 @@ int main(int argc, char *argv[]) show_help(argv[0]); return 0; - case 'i' : - infile = strdup(optarg); - break; - case 'o' : outfile = strdup(optarg); break; @@ -380,15 +374,17 @@ int main(int argc, char *argv[]) } - if ( infile == NULL ) { - infile = strdup("-"); + if ( argc != (optind+1) ) { + ERROR("Please provide exactly one stream filename.\n"); + return 1; } + + infile = argv[optind++]; st = open_stream_for_read(infile); if ( st == NULL ) { ERROR("Failed to open input stream '%s'\n", infile); return 1; } - free(infile); /* Sanitise output filename */ if ( outfile == NULL ) { -- cgit v1.2.3 From 82e03c74e1f054d885b1572b73a719583c364985 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Wed, 5 Mar 2014 07:57:08 +0100 Subject: Fix progress bar for detwinning --- src/ambigator.c | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/ambigator.c b/src/ambigator.c index 2af39e8a..e790eba0 100644 --- a/src/ambigator.c +++ b/src/ambigator.c @@ -240,6 +240,8 @@ static void detwin(struct flist **crystals, int n_crystals, SymOpList *amb, int p = 0; int q = 0; + progress_bar(i, n_crystals-1, "Calculating"); + for ( j=0; j Date: Wed, 5 Mar 2014 13:54:52 +0100 Subject: Symmetry for detwinning should be the actual symmetry, not the apparent symmetry --- src/ambigator.c | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/ambigator.c b/src/ambigator.c index e790eba0..d60c484b 100644 --- a/src/ambigator.c +++ b/src/ambigator.c @@ -67,8 +67,8 @@ static void show_help(const char *s) " -h, --help Display this help message.\n" "\n" " -o, --output= Output stream.\n" -" -y, --symmetry= Apparent (\"source\") symmetry.\n" -" -e Actual (\"target\") symmetry.\n" +" -y, --symmetry= Actual (\"target\") symmetry.\n" +" -w Apparent (\"source\" or \"twinned\") symmetry.\n" " -n, --iterations= Iterate times.\n" " --highres= High resolution cutoff in A.\n" " --lowres= Low resolution cutoff in A.\n" @@ -318,7 +318,7 @@ int main(int argc, char *argv[]) }; /* Short options */ - while ((c = getopt_long(argc, argv, "ho:y:n:e:", + while ((c = getopt_long(argc, argv, "ho:y:n:w:", longopts, NULL)) != -1) { @@ -336,7 +336,7 @@ int main(int argc, char *argv[]) s_sym_str = strdup(optarg); break; - case 'e' : + case 'w' : e_sym_str = strdup(optarg); break; -- cgit v1.2.3 From 4ce94e129805d131441e3cfc0e78b6ef36439cef Mon Sep 17 00:00:00 2001 From: Thomas White Date: Wed, 5 Mar 2014 13:56:01 +0100 Subject: Print warning if any crystals have no correlations at all --- src/ambigator.c | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/src/ambigator.c b/src/ambigator.c index d60c484b..e685066d 100644 --- a/src/ambigator.c +++ b/src/ambigator.c @@ -231,6 +231,7 @@ static void detwin(struct flist **crystals, int n_crystals, SymOpList *amb, int nch = 0; float mf = 0.0; int nmf = 0; + int ndud = 0; for ( i=0; i 0 ) { + STATUS("Warning: %i crystals had no correlation\n", ndud); + } + STATUS("Mean f = %f, changed %i assignments this time.\n", mf/nmf, nch); } -- cgit v1.2.3 From 9c59236459d2b34176d06065f6244979ec90f467 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Wed, 5 Mar 2014 13:59:55 +0100 Subject: Fix includes --- src/ambigator.c | 11 ++--------- 1 file changed, 2 insertions(+), 9 deletions(-) diff --git a/src/ambigator.c b/src/ambigator.c index e685066d..95d0b097 100644 --- a/src/ambigator.c +++ b/src/ambigator.c @@ -41,22 +41,15 @@ #include #include +#include #include -#include #include #include -#include -#include -#include -#include #include #include +#include #include -#include "post-refinement.h" -#include "hrs-scaling.h" -#include "scaling-report.h" - static void show_help(const char *s) { -- cgit v1.2.3 From 7e93b4229db92f84be0888a680c35c4eb2c9a371 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Wed, 5 Mar 2014 14:16:01 +0100 Subject: Add an assertion --- src/ambigator.c | 1 + 1 file changed, 1 insertion(+) diff --git a/src/ambigator.c b/src/ambigator.c index 95d0b097..14361578 100644 --- a/src/ambigator.c +++ b/src/ambigator.c @@ -150,6 +150,7 @@ static struct flist *asymm_and_merge(RefList *in, const SymOpList *sym, f->n++; } reflist_free(asym); + assert(f->n == n); return f; } -- cgit v1.2.3 From a72c3c13881071dabaf43abbf25b7669505573ee Mon Sep 17 00:00:00 2001 From: Thomas White Date: Wed, 5 Mar 2014 14:16:55 +0100 Subject: This line moved further up an earlier commit --- src/ambigator.c | 1 - 1 file changed, 1 deletion(-) diff --git a/src/ambigator.c b/src/ambigator.c index 14361578..d7ae18fe 100644 --- a/src/ambigator.c +++ b/src/ambigator.c @@ -265,7 +265,6 @@ static void detwin(struct flist **crystals, int n_crystals, SymOpList *amb, f /= p; g /= q; - if ( (p==0) || (q==0) ) continue; mf += f; nmf++; -- cgit v1.2.3 From 388edd2e2d8ca6e457ecd653d2d5e65bd3eb39f3 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Wed, 5 Mar 2014 14:24:39 +0100 Subject: e->w --- src/ambigator.c | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/src/ambigator.c b/src/ambigator.c index d7ae18fe..5bc8e7d1 100644 --- a/src/ambigator.c +++ b/src/ambigator.c @@ -291,8 +291,8 @@ int main(int argc, char *argv[]) char *outfile = NULL; char *s_sym_str = NULL; SymOpList *s_sym; - char *e_sym_str = NULL; - SymOpList *e_sym; + char *w_sym_str = NULL; + SymOpList *w_sym; SymOpList *amb; int n_iter = 1; int n_crystals, n_chunks, max_crystals; @@ -340,7 +340,7 @@ int main(int argc, char *argv[]) break; case 'w' : - e_sym_str = strdup(optarg); + w_sym_str = strdup(optarg); break; case 'n' : @@ -401,14 +401,14 @@ int main(int argc, char *argv[]) s_sym = get_pointgroup(s_sym_str); free(s_sym_str); - if ( e_sym_str == NULL ) { - e_sym = NULL; + if ( w_sym_str == NULL ) { + w_sym = NULL; amb = NULL; } else { - e_sym = get_pointgroup(e_sym_str); - free(e_sym_str); - if ( e_sym == NULL ) return 1; - amb = get_ambiguities(s_sym, e_sym); + w_sym = get_pointgroup(w_sym_str); + free(w_sym_str); + if ( w_sym == NULL ) return 1; + amb = get_ambiguities(w_sym, s_sym); if ( amb == NULL ) return 1; STATUS("Ambiguity operations:\n"); describe_symmetry(amb); -- cgit v1.2.3 From 21828ab4aedcde9885b8627d5d7e7b499081a60b Mon Sep 17 00:00:00 2001 From: Thomas White Date: Wed, 5 Mar 2014 15:00:50 +0100 Subject: Add reindexing --- src/ambigator.c | 75 ++++++++++++++++++++++++++++++++++++++++++--------------- 1 file changed, 55 insertions(+), 20 deletions(-) diff --git a/src/ambigator.c b/src/ambigator.c index 5bc8e7d1..0a5415cf 100644 --- a/src/ambigator.c +++ b/src/ambigator.c @@ -77,12 +77,14 @@ struct flist { int n; unsigned int *s; + unsigned int *s_reidx; float *i; }; static struct flist *asymm_and_merge(RefList *in, const SymOpList *sym, - UnitCell *cell, double rmin, double rmax) + UnitCell *cell, double rmin, double rmax, + SymOpList *amb) { Reflection *refl; RefListIterator *iter; @@ -131,8 +133,9 @@ static struct flist *asymm_and_merge(RefList *in, const SymOpList *sym, n = num_reflections(asym); f->s = malloc(n*sizeof(unsigned int)); + f->s_reidx = malloc(n*sizeof(unsigned int)); f->i = malloc(n*sizeof(float)); - if ( (f->s == NULL) || (f->i == NULL) ) { + if ( (f->s == NULL) || (f->i == NULL) || (f->s_reidx == NULL) ) { ERROR("Failed to allocate flist\n"); return NULL; } @@ -143,9 +146,16 @@ static struct flist *asymm_and_merge(RefList *in, const SymOpList *sym, refl = next_refl(refl, iter) ) { signed int h, k, l; + signed int hr, kr, lr; + signed int hra, kra, lra; get_indices(refl, &h, &k, &l); f->s[f->n] = SERIAL(h, k, l); + + get_equiv(amb, NULL, 1, h, k, l, &hr, &kr, &lr); + get_asymm(sym, hr, kr, lr, &hra, &kra, &lra); + f->s_reidx[f->n] = SERIAL(hra, kra, lra); + f->i[f->n] = get_intensity(refl); f->n++; } @@ -156,7 +166,7 @@ static struct flist *asymm_and_merge(RefList *in, const SymOpList *sym, } -static float corr(struct flist *a, struct flist *b, int *pn) +static float corr(struct flist *a, struct flist *b, int *pn, int a_reidx) { float s_xy = 0.0; float s_x = 0.0; @@ -168,10 +178,20 @@ static float corr(struct flist *a, struct flist *b, int *pn) int ap = 0; int bp = 0; int done = 0; + unsigned int *sa; + unsigned int *sb; + + if ( a_reidx ) { + sa = a->s_reidx; + } else { + sa = a->s; + } + + sb = b->s; while ( 1 ) { - while ( a->s[ap] > b->s[bp] ) { + while ( sa[ap] > sb[bp] ) { if ( ++bp == b->n ) { done = 1; break; @@ -179,7 +199,7 @@ static float corr(struct flist *a, struct flist *b, int *pn) } if ( done ) break; - while ( a->s[ap] < b->s[bp] ) { + while ( sa[ap] < sb[bp] ) { if ( ++ap == a->n ) { done = 1; break; @@ -187,7 +207,7 @@ static float corr(struct flist *a, struct flist *b, int *pn) } if ( done ) break; - if ( a->s[ap] == b->s[bp] ) { + if ( sa[ap] == sb[bp] ) { float aint, bint; @@ -218,8 +238,7 @@ static float corr(struct flist *a, struct flist *b, int *pn) } -static void detwin(struct flist **crystals, int n_crystals, SymOpList *amb, - int *assignments) +static void detwin(struct flist **crystals, int n_crystals, int *assignments) { int i; int nch = 0; @@ -239,21 +258,36 @@ static void detwin(struct flist **crystals, int n_crystals, SymOpList *amb, for ( j=0; j 2 ) { - if ( n < 3 ) continue; + if ( assignments[i] == assignments[j] ) { + f += cc; + p++; + } else { + g += cc; + q++; + } + + } + + if ( n_reidx > 2 ) { + + if ( assignments[i] == assignments[j] ) { + g += cc_reidx; + q++; + } else { + f += cc_reidx; + p++; + } - if ( assignments[i] == assignments[j] ) { - f += cc; - p++; - } else { - g += cc; - q++; } } @@ -461,7 +495,8 @@ int main(int argc, char *argv[]) list = crystal_get_reflections(cr); crystals[n_crystals] = asymm_and_merge(list, s_sym, cell, - rmin, rmax); + rmin, rmax, + amb); cell_free(cell); n_crystals++; @@ -496,7 +531,7 @@ int main(int argc, char *argv[]) } for ( i=0; i Date: Wed, 5 Mar 2014 15:10:11 +0100 Subject: Work with or without -w --- src/ambigator.c | 21 ++++++++++++++------- 1 file changed, 14 insertions(+), 7 deletions(-) diff --git a/src/ambigator.c b/src/ambigator.c index 0a5415cf..f864bb1e 100644 --- a/src/ambigator.c +++ b/src/ambigator.c @@ -152,9 +152,11 @@ static struct flist *asymm_and_merge(RefList *in, const SymOpList *sym, get_indices(refl, &h, &k, &l); f->s[f->n] = SERIAL(h, k, l); - get_equiv(amb, NULL, 1, h, k, l, &hr, &kr, &lr); - get_asymm(sym, hr, kr, lr, &hra, &kra, &lra); - f->s_reidx[f->n] = SERIAL(hra, kra, lra); + if ( amb != NULL ) { + get_equiv(amb, NULL, 1, h, k, l, &hr, &kr, &lr); + get_asymm(sym, hr, kr, lr, &hra, &kra, &lra); + f->s_reidx[f->n] = SERIAL(hra, kra, lra); + } f->i[f->n] = get_intensity(refl); f->n++; @@ -238,7 +240,8 @@ static float corr(struct flist *a, struct flist *b, int *pn, int a_reidx) } -static void detwin(struct flist **crystals, int n_crystals, int *assignments) +static void detwin(struct flist **crystals, int n_crystals, int *assignments, + SymOpList *amb) { int i; int nch = 0; @@ -264,7 +267,11 @@ static void detwin(struct flist **crystals, int n_crystals, int *assignments) if ( i == j ) continue; cc = corr(crystals[i], crystals[j], &n, 0); - cc_reidx = corr(crystals[i], crystals[j], &n_reidx, 1); + + if ( amb != NULL ) { + cc_reidx = corr(crystals[i], crystals[j], + &n_reidx, 1); + } if ( n > 2 ) { @@ -278,7 +285,7 @@ static void detwin(struct flist **crystals, int n_crystals, int *assignments) } - if ( n_reidx > 2 ) { + if ( (amb != NULL) && (n_reidx > 2) ) { if ( assignments[i] == assignments[j] ) { g += cc_reidx; @@ -531,7 +538,7 @@ int main(int argc, char *argv[]) } for ( i=0; i Date: Wed, 5 Mar 2014 16:13:23 +0100 Subject: Formatting --- libcrystfel/src/symmetry.c | 2 -- 1 file changed, 2 deletions(-) diff --git a/libcrystfel/src/symmetry.c b/libcrystfel/src/symmetry.c index 81d87b25..dec8a520 100644 --- a/libcrystfel/src/symmetry.c +++ b/libcrystfel/src/symmetry.c @@ -1110,8 +1110,6 @@ void get_equiv(const SymOpList *ops, const SymOpMask *m, int idx, } - - if ( idx >= n ) { ERROR("Index %i out of range for point group '%s'\n", idx, -- cgit v1.2.3 From fba54b335b3c9dcdbf31d7ae92b6ed307b2efd21 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Wed, 5 Mar 2014 16:13:45 +0100 Subject: Fix reindexing --- src/ambigator.c | 65 +++++++++++++++++++++++++++++++++++++++++++-------------- 1 file changed, 49 insertions(+), 16 deletions(-) diff --git a/src/ambigator.c b/src/ambigator.c index f864bb1e..bd85ea8e 100644 --- a/src/ambigator.c +++ b/src/ambigator.c @@ -76,9 +76,12 @@ static void show_help(const char *s) struct flist { int n; + unsigned int *s; - unsigned int *s_reidx; float *i; + + unsigned int *s_reidx; + float *i_reidx; }; @@ -135,7 +138,9 @@ static struct flist *asymm_and_merge(RefList *in, const SymOpList *sym, f->s = malloc(n*sizeof(unsigned int)); f->s_reidx = malloc(n*sizeof(unsigned int)); f->i = malloc(n*sizeof(float)); - if ( (f->s == NULL) || (f->i == NULL) || (f->s_reidx == NULL) ) { + f->i_reidx = malloc(n*sizeof(float)); + if ( (f->s == NULL) || (f->i == NULL) + || (f->s_reidx == NULL) || (f->i_reidx == NULL) ) { ERROR("Failed to allocate flist\n"); return NULL; } @@ -146,23 +151,51 @@ static struct flist *asymm_and_merge(RefList *in, const SymOpList *sym, refl = next_refl(refl, iter) ) { signed int h, k, l; - signed int hr, kr, lr; - signed int hra, kra, lra; get_indices(refl, &h, &k, &l); f->s[f->n] = SERIAL(h, k, l); - if ( amb != NULL ) { + f->i[f->n] = get_intensity(refl); + f->n++; + } + assert(f->n == n); + + if ( amb != NULL ) { + + RefList *reidx = reflist_new(); + if ( reidx == NULL ) return NULL; + + for ( refl = first_refl(asym, &iter); + refl != NULL; + refl = next_refl(refl, iter) ) + { + signed int h, k, l; + signed int hr, kr, lr; + signed int hra, kra, lra; + Reflection *cr; + + get_indices(refl, &h, &k, &l); get_equiv(amb, NULL, 1, h, k, l, &hr, &kr, &lr); get_asymm(sym, hr, kr, lr, &hra, &kra, &lra); - f->s_reidx[f->n] = SERIAL(hra, kra, lra); + cr = add_refl(reidx, hra, kra, lra); + copy_data(cr, refl); } - f->i[f->n] = get_intensity(refl); - f->n++; + n = 0; + for ( refl = first_refl(reidx, &iter); + refl != NULL; + refl = next_refl(refl, iter) ) + { + signed int h, k, l; + get_indices(refl, &h, &k, &l); + f->s_reidx[n] = SERIAL(h, k, l); + f->i_reidx[n++] = get_intensity(refl); + } + + reflist_free(reidx); } + reflist_free(asym); - assert(f->n == n); return f; } @@ -181,19 +214,19 @@ static float corr(struct flist *a, struct flist *b, int *pn, int a_reidx) int bp = 0; int done = 0; unsigned int *sa; - unsigned int *sb; + float *ia; if ( a_reidx ) { sa = a->s_reidx; + ia = a->i_reidx; } else { sa = a->s; + ia = a->i; } - sb = b->s; - while ( 1 ) { - while ( sa[ap] > sb[bp] ) { + while ( sa[ap] > b->s[bp] ) { if ( ++bp == b->n ) { done = 1; break; @@ -201,7 +234,7 @@ static float corr(struct flist *a, struct flist *b, int *pn, int a_reidx) } if ( done ) break; - while ( sa[ap] < sb[bp] ) { + while ( sa[ap] < b->s[bp] ) { if ( ++ap == a->n ) { done = 1; break; @@ -209,11 +242,11 @@ static float corr(struct flist *a, struct flist *b, int *pn, int a_reidx) } if ( done ) break; - if ( sa[ap] == sb[bp] ) { + if ( sa[ap] == b->s[bp] ) { float aint, bint; - aint = a->i[ap]; + aint = ia[ap]; bint = b->i[bp]; s_xy += aint*bint; -- cgit v1.2.3 From 533c8aad7b256cceeb1ae9cae346ce00827cc3b1 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Wed, 5 Mar 2014 16:29:41 +0100 Subject: Add --end-assignments --- src/ambigator.c | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/src/ambigator.c b/src/ambigator.c index bd85ea8e..7aa9a2ae 100644 --- a/src/ambigator.c +++ b/src/ambigator.c @@ -363,6 +363,7 @@ int main(int argc, char *argv[]) int c; const char *infile; char *outfile = NULL; + char *end_ass_fn = NULL; char *s_sym_str = NULL; SymOpList *s_sym; char *w_sym_str = NULL; @@ -390,6 +391,7 @@ int main(int argc, char *argv[]) {"highres", 1, NULL, 2}, {"lowres", 1, NULL, 3}, + {"end-assignments", 1, NULL, 4}, {0, 0, NULL, 0} }; @@ -437,6 +439,10 @@ int main(int argc, char *argv[]) rmin = 1.0 / (lowres/1e10); break; + case 4 : + end_ass_fn = strdup(optarg); + break; + case 0 : break; @@ -574,6 +580,18 @@ int main(int argc, char *argv[]) detwin(crystals, n_crystals, assignments, amb); } + if ( end_ass_fn != NULL ) { + FILE *fh = fopen(end_ass_fn, "w"); + if ( fh == NULL ) { + ERROR("Failed to open '%s'\n", end_ass_fn); + } else { + for ( i=0; i Date: Thu, 6 Mar 2014 12:07:45 +0100 Subject: Add scripts/fg-graph --- Makefile.am | 2 +- scripts/fg-graph | 61 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 62 insertions(+), 1 deletion(-) create mode 100755 scripts/fg-graph diff --git a/Makefile.am b/Makefile.am index 0a187d94..bdbbb829 100644 --- a/Makefile.am +++ b/Makefile.am @@ -176,7 +176,7 @@ script_DATA = scripts/alternate-stream scripts/cell-please \ scripts/random-image scripts/README scripts/sequence-image \ scripts/split-indexed scripts/stream_grep scripts/zone-axes \ scripts/Rsplit_surface scripts/Rsplit_surface.py \ - scripts/clean-stream.py + scripts/clean-stream.py scripts/fg-graph EXTRA_DIST += $(script_DATA) diff --git a/scripts/fg-graph b/scripts/fg-graph new file mode 100755 index 00000000..701277c3 --- /dev/null +++ b/scripts/fg-graph @@ -0,0 +1,61 @@ +#!/bin/bash + +# Copyright © 2014 Wolfgang Brehm +# Copyright © 2014 Deutsches Elektronen-Synchrotron DESY, +# a research centre of the Helmholtz Association. +# +# Authors: +# 2013-2014 Wolfgang Brehm +# 2014 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 . +# + +FILE=fg.dat +NPATT=4000 +NITER=4 +CMIN=-0.1 +CMAX=0.3 + +gnuplot -p << EOF +set terminal pngcairo size 1600,1000 enhanced monochrome font 'Helvetica,20' linewidth 2 +set output "correlation.png" +set xlabel "Iteration i" +set ylabel "Correlation" +rnd(x) = x - floor(x) < 0.5 ? floor(x) : ceil(x) +round(x1,x2)=rnd(10**x2*x1)/10.0**x2 # funktion mit zwei Argumenten, x1 wird gerundet x2 anzahl der Stellen nach dem Komma +set xzeroaxis lc rgb "black" lt 1 +set key top left +set xrange [0:${NITER}*${NPATT}] +set yrange [${CMIN}:${CMAX}] +set ytics nomirror +set xtics nomirror +set x2tics nomirror + +set x2range [0:${NITER}] +set x2label "Iteration i gemessen in n" +set arrow from ${NPATT},${CMIN} to ${NPATT},${CMAX} nohead lc rgb "black" + +plot\ + 1 lc rgb "black" lt 1 notitle,\ + "${FILE}" u 0:1 ps 0.2 pt 7 lc rgb "orange" notitle,\ + "${FILE}" u 0:2 ps 0.2 pt 7 lc rgb "#0088FF" notitle,\ + "${FILE}" u 0:(rand(0)>0.5?\$1:1/0) ps 0.2 pt 7 lc rgb "orange" notitle ,\ + "${FILE}" u (round(\$0, -3)):(\$1>\$2?\$1:\$2) lw 3 lc rgb "black" smooth unique notitle,\ + "${FILE}" u (round(\$0, -3)):(\$1<\$2?\$1:\$2) lw 3 lc rgb "black" smooth unique notitle,\ + "${FILE}" u (round(\$0, -3)):1 lw 3 lc rgb "red" smooth unique notitle,\ + "${FILE}" u (round(\$0, -3)):2 lw 3 lc rgb "blue" smooth unique notitle +EOF -- cgit v1.2.3 From c4649fe727361c43482b0508f1a5d3e82d8cb57a Mon Sep 17 00:00:00 2001 From: Thomas White Date: Thu, 6 Mar 2014 12:08:08 +0100 Subject: Update --help --- src/ambigator.c | 1 + 1 file changed, 1 insertion(+) diff --git a/src/ambigator.c b/src/ambigator.c index 7aa9a2ae..289fc36f 100644 --- a/src/ambigator.c +++ b/src/ambigator.c @@ -65,6 +65,7 @@ static void show_help(const char *s) " -n, --iterations= Iterate times.\n" " --highres= High resolution cutoff in A.\n" " --lowres= Low resolution cutoff in A.\n" +" --end-assignments= Save end assignments to file .\n" ); } -- cgit v1.2.3 From ad347841d21095fe414e09f42eccd08883a2eb51 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Thu, 6 Mar 2014 12:08:51 +0100 Subject: Add --fg-graph --- src/ambigator.c | 25 +++++++++++++++++++++++-- 1 file changed, 23 insertions(+), 2 deletions(-) diff --git a/src/ambigator.c b/src/ambigator.c index 289fc36f..6c4e0ff2 100644 --- a/src/ambigator.c +++ b/src/ambigator.c @@ -66,6 +66,7 @@ static void show_help(const char *s) " --highres= High resolution cutoff in A.\n" " --lowres= Low resolution cutoff in A.\n" " --end-assignments= Save end assignments to file .\n" +" --fg-graph= Save f and g correlation values to file .\n" ); } @@ -275,7 +276,7 @@ static float corr(struct flist *a, struct flist *b, int *pn, int a_reidx) static void detwin(struct flist **crystals, int n_crystals, int *assignments, - SymOpList *amb) + SymOpList *amb, FILE *fh) { int i; int nch = 0; @@ -341,6 +342,8 @@ static void detwin(struct flist **crystals, int n_crystals, int *assignments, f /= p; g /= q; + fprintf(fh, "%5.3f %5.3f\n", f, g); + mf += f; nmf++; @@ -365,6 +368,7 @@ int main(int argc, char *argv[]) const char *infile; char *outfile = NULL; char *end_ass_fn = NULL; + char *fg_graph_fn = NULL; char *s_sym_str = NULL; SymOpList *s_sym; char *w_sym_str = NULL; @@ -382,6 +386,7 @@ int main(int argc, char *argv[]) float highres, lowres; double rmin = 0.0; /* m^-1 */ double rmax = INFINITY; /* m^-1 */ + FILE *fgfh = NULL; /* Long options */ const struct option longopts[] = { @@ -393,6 +398,7 @@ int main(int argc, char *argv[]) {"highres", 1, NULL, 2}, {"lowres", 1, NULL, 3}, {"end-assignments", 1, NULL, 4}, + {"fg-graph", 1, NULL, 5}, {0, 0, NULL, 0} }; @@ -444,6 +450,10 @@ int main(int argc, char *argv[]) end_ass_fn = strdup(optarg); break; + case 5 : + fg_graph_fn = strdup(optarg); + break; + case 0 : break; @@ -577,8 +587,19 @@ int main(int argc, char *argv[]) orig_assignments[i] = assignments[i]; } + if ( fg_graph_fn != NULL ) { + fgfh = fopen(fg_graph_fn, "w"); + if ( fgfh == NULL ) { + ERROR("Failed to open '%s'\n", fg_graph_fn); + } + } + for ( i=0; i Date: Thu, 6 Mar 2014 12:10:13 +0100 Subject: Update AUTHORS and README --- AUTHORS | 3 +++ README | 1 + 2 files changed, 4 insertions(+) diff --git a/AUTHORS b/AUTHORS index 91876b85..4b24d960 100644 --- a/AUTHORS +++ b/AUTHORS @@ -42,3 +42,6 @@ * Alexandra Tolstikova SASE spectrum simulation + +* Wolfgang Brehm + "Detwinning" algorithm (ambigator) diff --git a/README b/README index e5d3b5db..83f44b38 100644 --- a/README +++ b/README @@ -19,6 +19,7 @@ Authors: Cornelius Gati Fedor Chervinskii Alexandra Tolstikova + Wolfgang Brehm Please read the AUTHORS file for a full list of contributions and contributors. See the COPYING file for licence conditions. Summary: GPLv3+. -- cgit v1.2.3 From 8591e4c373f9b9921569b2d8c8e7b57bb555d563 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Thu, 6 Mar 2014 14:45:44 +0100 Subject: Ignore "twin-proof" reflections --- src/ambigator.c | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/src/ambigator.c b/src/ambigator.c index 6c4e0ff2..5be72053 100644 --- a/src/ambigator.c +++ b/src/ambigator.c @@ -117,6 +117,22 @@ static struct flist *asymm_and_merge(RefList *in, const SymOpList *sym, get_asymm(sym, h, k, l, &ha, &ka, &la); + if ( amb != NULL ) { + + signed int hr, kr, lr; + signed int hra, kra, lra; + + get_equiv(amb, NULL, 1, ha, ka, la, &hr, &kr, &lr); + get_asymm(sym, hr, kr, lr, &hra, &kra, &lra); + + /* Skip twin-proof reflections */ + if ( (ha==hra) && (ka==kra) && (la==lra) ) { + //STATUS("%i %i %i is twin proof\n", h, k, l); + continue; + } + + } + cr = find_refl(asym, ha, ka, la); if ( cr == NULL ) { cr = add_refl(asym, ha, ka, la); @@ -179,6 +195,7 @@ static struct flist *asymm_and_merge(RefList *in, const SymOpList *sym, get_indices(refl, &h, &k, &l); get_equiv(amb, NULL, 1, h, k, l, &hr, &kr, &lr); get_asymm(sym, hr, kr, lr, &hra, &kra, &lra); + cr = add_refl(reidx, hra, kra, lra); copy_data(cr, refl); } -- cgit v1.2.3 From b361d7c9bd36a831c712aac55e60f733aa15f7b3 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Thu, 6 Mar 2014 16:13:20 +0100 Subject: Add --ncorr --- src/ambigator.c | 162 ++++++++++++++++++++++++++++++++++++++++++++------------ 1 file changed, 128 insertions(+), 34 deletions(-) diff --git a/src/ambigator.c b/src/ambigator.c index 5be72053..8d7cd4f2 100644 --- a/src/ambigator.c +++ b/src/ambigator.c @@ -67,6 +67,7 @@ static void show_help(const char *s) " --lowres= Low resolution cutoff in A.\n" " --end-assignments= Save end assignments to file .\n" " --fg-graph= Save f and g correlation values to file .\n" +" --ncorr= Use correlations per crystal. Default 1000\n" ); } @@ -292,8 +293,93 @@ static float corr(struct flist *a, struct flist *b, int *pn, int a_reidx) } -static void detwin(struct flist **crystals, int n_crystals, int *assignments, - SymOpList *amb, FILE *fh) +struct cc_list +{ + signed int *ind; + float *cc; + + signed int *ind_reidx; + float *cc_reidx; +}; + + +static struct cc_list *calc_ccs(struct flist **crystals, int n_crystals, + int ncorr, SymOpList *amb) +{ + struct cc_list *ccs; + int i; + + assert(n_crystals >= ncorr); + ncorr++; /* Extra value at end for sentinel */ + + ccs = malloc(n_crystals*sizeof(struct cc_list)); + if ( ccs == NULL ) return NULL; + + for ( i=0; i 2 ) { - - if ( assignments[i] == assignments[j] ) { - f += cc; - p++; - } else { - g += cc; - q++; - } - - } + } - if ( (amb != NULL) && (n_reidx > 2) ) { + for ( k=0; (ccs[i].ind_reidx[k] != 0); k++ ) { - if ( assignments[i] == assignments[j] ) { - g += cc_reidx; - q++; - } else { - f += cc_reidx; - p++; - } + int j = ccs[i].ind_reidx[k]-1; + float cc = ccs[i].cc_reidx[k]; + if ( assignments[i] == assignments[j] ) { + g += cc; + q++; + } else { + f += cc; + p++; } } @@ -404,6 +480,8 @@ int main(int argc, char *argv[]) double rmin = 0.0; /* m^-1 */ double rmax = INFINITY; /* m^-1 */ FILE *fgfh = NULL; + struct cc_list *ccs; + int ncorr = 1000; /* Long options */ const struct option longopts[] = { @@ -416,6 +494,7 @@ int main(int argc, char *argv[]) {"lowres", 1, NULL, 3}, {"end-assignments", 1, NULL, 4}, {"fg-graph", 1, NULL, 5}, + {"ncorr", 1, NULL, 6}, {0, 0, NULL, 0} }; @@ -471,6 +550,13 @@ int main(int argc, char *argv[]) fg_graph_fn = strdup(optarg); break; + case 6 : + if ( sscanf(optarg, "%i", &ncorr) != 1 ) { + ERROR("Invalid value for --ncorr\n"); + return 1; + } + break; + case 0 : break; @@ -611,8 +697,16 @@ int main(int argc, char *argv[]) } } + ccs = calc_ccs(crystals, n_crystals, ncorr, amb); + if ( ccs == NULL ) { + ERROR("Failed to allocate CCs\n"); + return 1; + } + + /* FIXME: Free crystals */ + for ( i=0; i Date: Thu, 6 Mar 2014 17:55:51 +0100 Subject: Add --stop-after --- src/ambigator.c | 15 ++++++++++++++- 1 file changed, 14 insertions(+), 1 deletion(-) diff --git a/src/ambigator.c b/src/ambigator.c index 8d7cd4f2..fe295cf7 100644 --- a/src/ambigator.c +++ b/src/ambigator.c @@ -68,6 +68,7 @@ static void show_help(const char *s) " --end-assignments= Save end assignments to file .\n" " --fg-graph= Save f and g correlation values to file .\n" " --ncorr= Use correlations per crystal. Default 1000\n" +" --stop-after= Use at most the first crystals.\n" ); } @@ -482,6 +483,7 @@ int main(int argc, char *argv[]) FILE *fgfh = NULL; struct cc_list *ccs; int ncorr = 1000; + int stop_after = 0; /* Long options */ const struct option longopts[] = { @@ -495,6 +497,7 @@ int main(int argc, char *argv[]) {"end-assignments", 1, NULL, 4}, {"fg-graph", 1, NULL, 5}, {"ncorr", 1, NULL, 6}, + {"stop-after", 1, NULL, 7}, {0, 0, NULL, 0} }; @@ -557,6 +560,13 @@ int main(int argc, char *argv[]) } break; + case 7 : + if ( sscanf(optarg, "%i", &stop_after) != 1 ) { + ERROR("Invalid value for --stop-after\n"); + return 1; + } + break; + case 0 : break; @@ -659,14 +669,17 @@ int main(int argc, char *argv[]) amb); cell_free(cell); n_crystals++; - reflist_free(list); + if ( stop_after && (n_crystals == stop_after) ) break; + } fprintf(stderr, "Loaded %i crystals from %i chunks\r", n_crystals, ++n_chunks); + if ( stop_after && (n_crystals == stop_after) ) break; + } while ( 1 ); fprintf(stderr, "\n"); -- cgit v1.2.3 From a2637fde88a99e222c98a15ed214f673a083821d Mon Sep 17 00:00:00 2001 From: Thomas White Date: Thu, 6 Mar 2014 17:55:58 +0100 Subject: Don't break if there are no reflections in the list --- src/ambigator.c | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/ambigator.c b/src/ambigator.c index fe295cf7..018ed186 100644 --- a/src/ambigator.c +++ b/src/ambigator.c @@ -245,6 +245,11 @@ static float corr(struct flist *a, struct flist *b, int *pn, int a_reidx) ia = a->i; } + if ( (a->n == 0) || (b->n == 0) ) { + *pn = 0; + return 0.0; + } + while ( 1 ) { while ( sa[ap] > b->s[bp] ) { -- cgit v1.2.3 From ea0907eaa8fb63f5a930ec286bfd02ae8a4c7d9f Mon Sep 17 00:00:00 2001 From: Thomas White Date: Thu, 6 Mar 2014 17:56:08 +0100 Subject: Remove unused line --- src/ambigator.c | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/ambigator.c b/src/ambigator.c index 018ed186..96309a8d 100644 --- a/src/ambigator.c +++ b/src/ambigator.c @@ -401,8 +401,6 @@ static void detwin(struct cc_list *ccs, int n_crystals, int *assignments, int p = 0; int q = 0; - //progress_bar(i, n_crystals-1, "Calculating"); - for ( k=0; (ccs[i].ind[k] != 0); k++ ) { int j = ccs[i].ind[k]-1; -- cgit v1.2.3 From 0b4786372baf696b2b8f34fc4ae455187c4f6f54 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Thu, 6 Mar 2014 18:23:39 +0100 Subject: Correlate against randomly chosen crystals instead of the next few --- src/ambigator.c | 29 ++++++++++++++++++++++++----- 1 file changed, 24 insertions(+), 5 deletions(-) diff --git a/src/ambigator.c b/src/ambigator.c index 96309a8d..9aed3b3f 100644 --- a/src/ambigator.c +++ b/src/ambigator.c @@ -40,6 +40,9 @@ #include #include #include +#include +#include +#include #include #include @@ -310,10 +313,11 @@ struct cc_list static struct cc_list *calc_ccs(struct flist **crystals, int n_crystals, - int ncorr, SymOpList *amb) + int ncorr, SymOpList *amb, gsl_rng *rng) { struct cc_list *ccs; int i; + gsl_permutation *p; assert(n_crystals >= ncorr); ncorr++; /* Extra value at end for sentinel */ @@ -321,9 +325,12 @@ static struct cc_list *calc_ccs(struct flist **crystals, int n_crystals, ccs = malloc(n_crystals*sizeof(struct cc_list)); if ( ccs == NULL ) return NULL; + p = gsl_permutation_alloc(n_crystals); + gsl_permutation_init(p); + for ( i=0; idata, n_crystals, sizeof(size_t)); + + for ( l=0; l Date: Thu, 6 Mar 2014 19:57:39 +0100 Subject: Show mean number of correlations per crystal --- src/ambigator.c | 15 +++++++++++++-- 1 file changed, 13 insertions(+), 2 deletions(-) diff --git a/src/ambigator.c b/src/ambigator.c index 9aed3b3f..03c0c476 100644 --- a/src/ambigator.c +++ b/src/ambigator.c @@ -313,11 +313,14 @@ struct cc_list static struct cc_list *calc_ccs(struct flist **crystals, int n_crystals, - int ncorr, SymOpList *amb, gsl_rng *rng) + int ncorr, SymOpList *amb, gsl_rng *rng, + float *pmean_nac) { struct cc_list *ccs; int i; gsl_permutation *p; + long int mean_nac = 0; + int nmean_nac = 0; assert(n_crystals >= ncorr); ncorr++; /* Extra value at end for sentinel */ @@ -365,6 +368,8 @@ static struct cc_list *calc_ccs(struct flist **crystals, int n_crystals, } ccs[i].ind[k] = 0; + mean_nac += k; + nmean_nac++; if ( amb != NULL ) { @@ -390,6 +395,8 @@ static struct cc_list *calc_ccs(struct flist **crystals, int n_crystals, } ccs[i].ind_reidx[k] = 0; + mean_nac += k; + nmean_nac++; } @@ -399,6 +406,8 @@ static struct cc_list *calc_ccs(struct flist **crystals, int n_crystals, gsl_permutation_free(p); + *pmean_nac = (float)mean_nac/nmean_nac; + return ccs; } @@ -506,6 +515,7 @@ int main(int argc, char *argv[]) struct cc_list *ccs; int ncorr = 1000; int stop_after = 0; + float mean_nac; /* Long options */ const struct option longopts[] = { @@ -732,11 +742,12 @@ int main(int argc, char *argv[]) } } - ccs = calc_ccs(crystals, n_crystals, ncorr, amb, rng); + ccs = calc_ccs(crystals, n_crystals, ncorr, amb, rng, &mean_nac); if ( ccs == NULL ) { ERROR("Failed to allocate CCs\n"); return 1; } + STATUS("Mean number of correlations per crystal: %.1f\n", mean_nac); /* FIXME: Free crystals */ -- cgit v1.2.3 From 788ceb929968341e714f11d980944d9ff42ffd8e Mon Sep 17 00:00:00 2001 From: Thomas White Date: Thu, 6 Mar 2014 19:57:49 +0100 Subject: Fussiness --- src/ambigator.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ambigator.c b/src/ambigator.c index 03c0c476..5f7f921b 100644 --- a/src/ambigator.c +++ b/src/ambigator.c @@ -344,9 +344,9 @@ static struct cc_list *calc_ccs(struct flist **crystals, int n_crystals, return NULL; } - k = 0; gsl_ran_shuffle(rng, p->data, n_crystals, sizeof(size_t)); + k = 0; for ( l=0; l Date: Thu, 6 Mar 2014 21:33:44 +0100 Subject: Fix segfault --- src/ambigator.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ambigator.c b/src/ambigator.c index 5f7f921b..028ed4bb 100644 --- a/src/ambigator.c +++ b/src/ambigator.c @@ -467,7 +467,7 @@ static void detwin(struct cc_list *ccs, int n_crystals, int *assignments, f /= p; g /= q; - fprintf(fh, "%5.3f %5.3f\n", f, g); + if ( fh != NULL ) fprintf(fh, "%5.3f %5.3f\n", f, g); mf += f; nmf++; -- cgit v1.2.3 From 017e55112a161349f4859c00d63eb0ea1be563a7 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Sun, 9 Mar 2014 21:34:18 +0100 Subject: Parallelise the CC calculation --- src/ambigator.c | 201 ++++++++++++++++++++++++++++++++++++++++++-------------- 1 file changed, 151 insertions(+), 50 deletions(-) diff --git a/src/ambigator.c b/src/ambigator.c index 028ed4bb..fbe29091 100644 --- a/src/ambigator.c +++ b/src/ambigator.c @@ -52,6 +52,7 @@ #include #include #include +#include static void show_help(const char *s) @@ -312,39 +313,132 @@ struct cc_list }; -static struct cc_list *calc_ccs(struct flist **crystals, int n_crystals, - int ncorr, SymOpList *amb, gsl_rng *rng, - float *pmean_nac) +struct queue_args +{ + int n_started; + int n_finished; + int n_to_do; + int mean_nac; + int nmean_nac; + gsl_permutation *p; + + struct cc_list *ccs; + struct flist **crystals; + int n_crystals; + int ncorr; + SymOpList *amb; + gsl_rng *rng; +}; + + +struct cc_job { struct cc_list *ccs; int i; + int mean_nac; + int nmean_nac; + + struct flist **crystals; + int n_crystals; + int ncorr; + SymOpList *amb; gsl_permutation *p; - long int mean_nac = 0; +}; + + +static void *get_task(void *vp) +{ + struct queue_args *qargs = vp; + struct cc_job *job; + + if ( qargs->n_started == qargs->n_to_do ) return NULL; + + job = malloc(sizeof(struct cc_job)); + if ( job == NULL ) return NULL; + + job->ccs = qargs->ccs; + job->i = qargs->n_started++; + + job->crystals = qargs->crystals; + job->n_crystals = qargs->n_crystals; + job->ncorr = qargs->ncorr; + job->amb = qargs->amb; + + gsl_ran_shuffle(qargs->rng, qargs->p->data, qargs->n_crystals, + sizeof(size_t)); + job->p = gsl_permutation_alloc(qargs->n_crystals); + gsl_permutation_memcpy(job->p, qargs->p); + + return job; +} + + + +static void final(void *qp, void *wp) +{ + struct queue_args *qargs = qp; + struct cc_job *job = wp; + + qargs->mean_nac += job->mean_nac; + qargs->nmean_nac += job->nmean_nac; + + gsl_permutation_free(job->p); + free(job); + + qargs->n_finished++; + progress_bar(qargs->n_finished, qargs->n_to_do, "Calculating CCs"); +} + + +static void work(void *wp, int cookie) +{ + struct cc_job *job = wp; + int i = job->i; + int k, l; + struct cc_list *ccs = job->ccs; + gsl_permutation *p = job->p; + struct flist **crystals = job->crystals; + int n_crystals = job->n_crystals; + int ncorr = job->ncorr; + SymOpList *amb = job->amb; + int mean_nac = 0; int nmean_nac = 0; - assert(n_crystals >= ncorr); - ncorr++; /* Extra value at end for sentinel */ + ccs[i].ind = malloc(ncorr*sizeof(int)); + ccs[i].cc = malloc(ncorr*sizeof(float)); + ccs[i].ind_reidx = calloc(ncorr, sizeof(int)); + ccs[i].cc_reidx = calloc(ncorr, sizeof(float)); + if ( (ccs[i].ind==NULL) || (ccs[i].cc==NULL) || + (ccs[i].ind_reidx==NULL) || (ccs[i].cc_reidx==NULL) ) { + return; + } - ccs = malloc(n_crystals*sizeof(struct cc_list)); - if ( ccs == NULL ) return NULL; + k = 0; + for ( l=0; ldata, n_crystals, sizeof(size_t)); + } + ccs[i].ind[k] = 0; + mean_nac += k; + nmean_nac++; + + if ( amb != NULL ) { k = 0; for ( l=0; lmean_nac = mean_nac; + job->nmean_nac = nmean_nac; +} - cc = corr(crystals[i], crystals[j], &n, 1); - if ( n < 4 ) continue; +static struct cc_list *calc_ccs(struct flist **crystals, int n_crystals, + int ncorr, SymOpList *amb, gsl_rng *rng, + float *pmean_nac) +{ + struct cc_list *ccs; + int nthreads = 8; + struct queue_args qargs; - ccs[i].ind_reidx[k] = j+1; - ccs[i].cc_reidx[k] = cc; - k++; + assert(n_crystals >= ncorr); + ncorr++; /* Extra value at end for sentinel */ - if ( k == ncorr-1 ) break; + ccs = malloc(n_crystals*sizeof(struct cc_list)); + if ( ccs == NULL ) return NULL; - } - ccs[i].ind_reidx[k] = 0; - mean_nac += k; - nmean_nac++; + qargs.p = gsl_permutation_alloc(n_crystals); + gsl_permutation_init(qargs.p); - } + qargs.n_started = 0; + qargs.n_finished = 0; + qargs.n_to_do = n_crystals; + qargs.ccs = ccs; + qargs.mean_nac = 0; + qargs.nmean_nac = 0; - progress_bar(i, n_crystals-1, "Calculating CCs"); + qargs.rng = rng; + qargs.crystals = crystals; + qargs.n_crystals = n_crystals; + qargs.ncorr = ncorr; + qargs.amb = amb; - } + run_threads(nthreads, work, get_task, final, &qargs, n_crystals, + 0, 0, 0); - gsl_permutation_free(p); + gsl_permutation_free(qargs.p); - *pmean_nac = (float)mean_nac/nmean_nac; + *pmean_nac = (float)qargs.mean_nac/qargs.nmean_nac; return ccs; } -- cgit v1.2.3 From 37deaefe474f698d2785189ffc9480f68c01007f Mon Sep 17 00:00:00 2001 From: Thomas White Date: Sun, 9 Mar 2014 21:36:50 +0100 Subject: Add -j --- src/ambigator.c | 17 +++++++++++++---- 1 file changed, 13 insertions(+), 4 deletions(-) diff --git a/src/ambigator.c b/src/ambigator.c index fbe29091..72894eb9 100644 --- a/src/ambigator.c +++ b/src/ambigator.c @@ -73,6 +73,7 @@ static void show_help(const char *s) " --fg-graph= Save f and g correlation values to file .\n" " --ncorr= Use correlations per crystal. Default 1000\n" " --stop-after= Use at most the first crystals.\n" +" -j Use threads for CC calculation.\n" ); } @@ -474,10 +475,9 @@ static void work(void *wp, int cookie) static struct cc_list *calc_ccs(struct flist **crystals, int n_crystals, int ncorr, SymOpList *amb, gsl_rng *rng, - float *pmean_nac) + float *pmean_nac, int nthreads) { struct cc_list *ccs; - int nthreads = 8; struct queue_args qargs; assert(n_crystals >= ncorr); @@ -617,6 +617,7 @@ int main(int argc, char *argv[]) int ncorr = 1000; int stop_after = 0; float mean_nac; + int n_threads = 1; /* Long options */ const struct option longopts[] = { @@ -636,7 +637,7 @@ int main(int argc, char *argv[]) }; /* Short options */ - while ((c = getopt_long(argc, argv, "ho:y:n:w:", + while ((c = getopt_long(argc, argv, "ho:y:n:w:j:", longopts, NULL)) != -1) { @@ -662,6 +663,13 @@ int main(int argc, char *argv[]) n_iter = atoi(optarg); break; + case 'j' : + if ( sscanf(optarg, "%i", &n_threads) != 1 ) { + ERROR("Invalid value for -j\n"); + return 1; + } + break; + case 2 : if ( sscanf(optarg, "%e", &highres) != 1 ) { ERROR("Invalid value for --highres\n"); @@ -843,7 +851,8 @@ int main(int argc, char *argv[]) } } - ccs = calc_ccs(crystals, n_crystals, ncorr, amb, rng, &mean_nac); + ccs = calc_ccs(crystals, n_crystals, ncorr, amb, rng, &mean_nac, + n_threads); if ( ccs == NULL ) { ERROR("Failed to allocate CCs\n"); return 1; -- cgit v1.2.3 From 90cd484d9ccaace998a27857b268722eb3f96722 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Sun, 9 Mar 2014 21:50:10 +0100 Subject: Threads should shuffle their own permutations --- src/ambigator.c | 46 +++++++++++++++++++++++++++++++++------------- 1 file changed, 33 insertions(+), 13 deletions(-) diff --git a/src/ambigator.c b/src/ambigator.c index 72894eb9..438af0d9 100644 --- a/src/ambigator.c +++ b/src/ambigator.c @@ -328,7 +328,7 @@ struct queue_args int n_crystals; int ncorr; SymOpList *amb; - gsl_rng *rng; + gsl_rng **rngs; }; @@ -343,7 +343,7 @@ struct cc_job int n_crystals; int ncorr; SymOpList *amb; - gsl_permutation *p; + gsl_rng **rngs; }; @@ -364,11 +364,7 @@ static void *get_task(void *vp) job->n_crystals = qargs->n_crystals; job->ncorr = qargs->ncorr; job->amb = qargs->amb; - - gsl_ran_shuffle(qargs->rng, qargs->p->data, qargs->n_crystals, - sizeof(size_t)); - job->p = gsl_permutation_alloc(qargs->n_crystals); - gsl_permutation_memcpy(job->p, qargs->p); + job->rngs = qargs->rngs; return job; } @@ -383,7 +379,6 @@ static void final(void *qp, void *wp) qargs->mean_nac += job->mean_nac; qargs->nmean_nac += job->nmean_nac; - gsl_permutation_free(job->p); free(job); qargs->n_finished++; @@ -397,13 +392,17 @@ static void work(void *wp, int cookie) int i = job->i; int k, l; struct cc_list *ccs = job->ccs; - gsl_permutation *p = job->p; struct flist **crystals = job->crystals; int n_crystals = job->n_crystals; int ncorr = job->ncorr; SymOpList *amb = job->amb; int mean_nac = 0; int nmean_nac = 0; + gsl_permutation *p; + + p = gsl_permutation_alloc(n_crystals); + gsl_permutation_init(p); + gsl_ran_shuffle(job->rngs[cookie], p->data, n_crystals, sizeof(size_t)); ccs[i].ind = malloc(ncorr*sizeof(int)); ccs[i].cc = malloc(ncorr*sizeof(float)); @@ -468,27 +467,46 @@ static void work(void *wp, int cookie) } + gsl_permutation_free(p); + job->mean_nac = mean_nac; job->nmean_nac = nmean_nac; } +static gsl_rng **setup_random(gsl_rng *rng, int n) +{ + gsl_rng **rngs; + int i; + + rngs = malloc(n * sizeof(gsl_rng *)); + if ( rngs == NULL ) return NULL; + + for ( i=0; i= ncorr); ncorr++; /* Extra value at end for sentinel */ + qargs.rngs = setup_random(rng, nthreads); + ccs = malloc(n_crystals*sizeof(struct cc_list)); if ( ccs == NULL ) return NULL; - qargs.p = gsl_permutation_alloc(n_crystals); - gsl_permutation_init(qargs.p); - qargs.n_started = 0; qargs.n_finished = 0; qargs.n_to_do = n_crystals; @@ -496,7 +514,6 @@ static struct cc_list *calc_ccs(struct flist **crystals, int n_crystals, qargs.mean_nac = 0; qargs.nmean_nac = 0; - qargs.rng = rng; qargs.crystals = crystals; qargs.n_crystals = n_crystals; qargs.ncorr = ncorr; @@ -505,6 +522,9 @@ static struct cc_list *calc_ccs(struct flist **crystals, int n_crystals, run_threads(nthreads, work, get_task, final, &qargs, n_crystals, 0, 0, 0); + for ( i=0; i Date: Sun, 9 Mar 2014 21:51:08 +0100 Subject: Fussiness --- src/ambigator.c | 1 - 1 file changed, 1 deletion(-) diff --git a/src/ambigator.c b/src/ambigator.c index 438af0d9..4387d252 100644 --- a/src/ambigator.c +++ b/src/ambigator.c @@ -370,7 +370,6 @@ static void *get_task(void *vp) } - static void final(void *qp, void *wp) { struct queue_args *qargs = qp; -- cgit v1.2.3 From 18f3a3a08794c933ed422a51580cf39e3c324cf7 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Mon, 10 Mar 2014 10:40:49 +0100 Subject: ncorr limit and default --- src/ambigator.c | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/src/ambigator.c b/src/ambigator.c index 4387d252..42838b41 100644 --- a/src/ambigator.c +++ b/src/ambigator.c @@ -633,7 +633,8 @@ int main(int argc, char *argv[]) double rmax = INFINITY; /* m^-1 */ FILE *fgfh = NULL; struct cc_list *ccs; - int ncorr = 1000; + int ncorr; + int ncorr_set = 0; int stop_after = 0; float mean_nac; int n_threads = 1; @@ -717,6 +718,8 @@ int main(int argc, char *argv[]) if ( sscanf(optarg, "%i", &ncorr) != 1 ) { ERROR("Invalid value for --ncorr\n"); return 1; + } else { + ncorr_set = 1; } break; @@ -870,6 +873,10 @@ int main(int argc, char *argv[]) } } + if ( !ncorr_set || (ncorr > n_crystals) ) { + ncorr = n_crystals; + } + ccs = calc_ccs(crystals, n_crystals, ncorr, amb, rng, &mean_nac, n_threads); if ( ccs == NULL ) { -- cgit v1.2.3 From 2544135897b6e920d74f824c26bbf9a75df6d275 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Mon, 10 Mar 2014 10:44:03 +0100 Subject: Free the crystals --- src/ambigator.c | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/src/ambigator.c b/src/ambigator.c index 42838b41..3b83d749 100644 --- a/src/ambigator.c +++ b/src/ambigator.c @@ -885,7 +885,14 @@ int main(int argc, char *argv[]) } STATUS("Mean number of correlations per crystal: %.1f\n", mean_nac); - /* FIXME: Free crystals */ + for ( i=0; is); + free(crystals[i]->i); + free(crystals[i]->s_reidx); + free(crystals[i]->i_reidx); + free(crystals[i]); + } + free(crystals); for ( i=0; i Date: Mon, 10 Mar 2014 10:59:12 +0100 Subject: Remove --stop-after --- src/ambigator.c | 13 ------------- 1 file changed, 13 deletions(-) diff --git a/src/ambigator.c b/src/ambigator.c index 3b83d749..a7aedd39 100644 --- a/src/ambigator.c +++ b/src/ambigator.c @@ -635,7 +635,6 @@ int main(int argc, char *argv[]) struct cc_list *ccs; int ncorr; int ncorr_set = 0; - int stop_after = 0; float mean_nac; int n_threads = 1; @@ -651,7 +650,6 @@ int main(int argc, char *argv[]) {"end-assignments", 1, NULL, 4}, {"fg-graph", 1, NULL, 5}, {"ncorr", 1, NULL, 6}, - {"stop-after", 1, NULL, 7}, {0, 0, NULL, 0} }; @@ -723,13 +721,6 @@ int main(int argc, char *argv[]) } break; - case 7 : - if ( sscanf(optarg, "%i", &stop_after) != 1 ) { - ERROR("Invalid value for --stop-after\n"); - return 1; - } - break; - case 0 : break; @@ -834,15 +825,11 @@ int main(int argc, char *argv[]) n_crystals++; reflist_free(list); - if ( stop_after && (n_crystals == stop_after) ) break; - } fprintf(stderr, "Loaded %i crystals from %i chunks\r", n_crystals, ++n_chunks); - if ( stop_after && (n_crystals == stop_after) ) break; - } while ( 1 ); fprintf(stderr, "\n"); -- cgit v1.2.3 From 1601fc3597a153a0cb5f868321697ebed7cf0b99 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Mon, 10 Mar 2014 11:54:31 +0100 Subject: Add man page --- Makefile.am | 2 +- doc/man/ambigator.1 | 103 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 104 insertions(+), 1 deletion(-) create mode 100644 doc/man/ambigator.1 diff --git a/Makefile.am b/Makefile.am index bdbbb829..cc5c6485 100644 --- a/Makefile.am +++ b/Makefile.am @@ -140,7 +140,7 @@ man_MANS = doc/man/crystfel.7 doc/man/check_hkl.1 doc/man/compare_hkl.1 \ doc/man/crystfel_geometry.5 doc/man/get_hkl.1 doc/man/hdfsee.1 \ doc/man/indexamajig.1 doc/man/partialator.1 doc/man/partial_sim.1 \ doc/man/pattern_sim.1 doc/man/process_hkl.1 doc/man/render_hkl.1 \ - doc/man/cell_explorer.1 + doc/man/cell_explorer.1 doc/man/ambigator.1 EXTRA_DIST += $(man_MANS) diff --git a/doc/man/ambigator.1 b/doc/man/ambigator.1 new file mode 100644 index 00000000..b0978fb7 --- /dev/null +++ b/doc/man/ambigator.1 @@ -0,0 +1,103 @@ +.\" +.\" ambigator man page +.\" +.\" Copyright © 2014 Deutsches Elektronen-Synchrotron DESY, +.\" a research centre of the Helmholtz Association. +.\" +.\" Part of CrystFEL - crystallography with a FEL +.\" + +.TH AMBIGATOR 1 +.SH NAME +ambigator \- Resolve indexing ambiguities +.SH SYNOPSIS +.PP +.B ambigator \fIinput.stream\fR \fB[-o\fR \fIoutput.stream\fR\fB] [options] + +.B ambigator --help + +.SH DESCRIPTION +This program resolves indexing ambiguities using a simplified variant of the clustering algorithm described by (FIXME: Reference). + +The algorithm works by calculating the individual correlation coefficients between the intensities from one crystal and those from each of the other crystals in turn. The mean \fIf\fR is then taken over all crystals with the same indexing assignment, and separately the mean \fIg\fR is taken over all crystals with the opposite indexing assignment. The indexing assignment for a crystal is changed if \fIg\fR > \fIf\fR. Every crystal is visited once in turn, and the pass over all the crystals repeated several times. + +Only one indexing ambiguity can be resolved at once. In other words, each crystal is considered to be indexable in one of only two ways. If the true indexing ambiguity has more possibilities than this, the resolution must be performed by running \fBambigator\fR multiple times with a different ambiguity operator each time. + +If the ambiguity operator is known (or, equivalently, the actual and apparent symmetries are both known), then the algorithm can be enhanced by including in \fIf\fR the correlation coefficients of all the patterns with the opposite indexing assignment to the current one, but after reindexing the other pattern first. Likewise, \fg\fR includes the correlation coefficients of the patterns with the same indexing assignment after reindexing. This enhances the algorithm to an extent roughly equivalent to doubling the number of crystals. + +The default behaviour is to compare each crystal to every other crystal. This leads to a computational complexity proportional to the square of the number of crystals. If the number of crystals is large, the number of comparisons can be limited without compromising the algorithm much. In this case, the crystals to correlate against will be selected randomly. + + +.SH OPTIONS +.PD 0 +.IP "\fB-o\fR \fIfilename\fR" +.IP \fB--output=\fR\fIfilename\fR +.PD +Write a re-indexed version of the input stream to \fIfilename\fR. This stream can then be merged normally using \fBprocess_hkl\fR or \fBpartialator\fR, but using the actual symmetry instead of the apparent one. +.IP +\fBWARNING\fR: There is no default filename. The default behaviour is not to output any reindexed stream! + +.PD 0 +.IP "\fB-y\fR \fIpg\fR" +.IP \fB--symmetry=\fR\fIpg\fR +.PD +Set the actual symmetry of the crystals. If you're not sure, set this to the highest symmetry which you want to assume, which might be \fB-1\fR to assume Friedel's Law alone or \fB1\fR (the default) for no symmetry at all. The algorithm will work significantly better if you can use a higher symmetry here. + +.PD 0 +.IP "\fB-w\fR \fIpg\fR" +.PD +Set the apparent symmetry of the crystals. The ambiguity operator will be determined by comparing this to the actual symmetry. If there is more than one operator, only the first will be used. +.IP +Using this option improves the algorithm to an extent roughly equivalent to doubling the number of crystals. + +.PD 0 +.IP "\fB-n\fR \fIn\fR" +.IP \fB--iterations=\fR\fIn\fR +The number of passes through the data to make. Extra iterations are not expensive once the initial correlation calculation has been performed, so set this value quite high. Two or three iterations are normally sufficient unless the number of correlations (see \fB--ncorr\fR) is small compared to the number of crystals. The default is \fB--iterations=6\fR. + +.PD 0 +.IP "\fB-j\fR \fIn\fR" +Number of threads to use for the CC calculation. + +.PD 0 +.IP \fB--highres=\fR\fId\fR +High resolution cutoff in Angstroms. + +.PD 0 +.IP \fB--lowres=\fR\fId\fR +Low resolution cutoff in Angstroms. + +.PD 0 +.IP \fB--end-assignments=\fR\fIfilename\fR +Write the end assignments to \fIfilename\fR. The file will be a list of 0 or 1, one value per line, in the same order as the crystals appear in the input stream. 1 means that the pattern should be reindexed according to the ambiguity operator. + +.PD 0 +.IP \fB--fg-graph=\fR\fIfilename\fR +Write f and g values to \fIfilename\fR, one line per crystal, repeating all crystals as they are visited by the algorithm. Plot these using \fBfg-graph\fR from the CrystFEL script folder to evaluate the ambiguity resolution. + +.PD 0 +.IP \fB--ncorr=\fR\fIn\fR +Use \fIn\fR correlations per crystal. The default is to correlate against every crystal. If the CC calculation is too slow, try \fB--ncorr=1000\fR. Note that this option sets the maximum number of correlations, and some crystals might not have enough common reflections to correlate to the number requested. The mean number of actual correlations per crystal will be output by the program after the CC calculation, and if this number is much smaller than \fIn\fR then this option will not have a significant effect. + + +.SH AUTHOR +This page was written by Thomas White. + +.SH REPORTING BUGS +Report bugs to , or visit . + +.SH COPYRIGHT AND DISCLAIMER +Copyright © 2014 Deutsches Elektronen-Synchrotron DESY, a research centre of the Helmholtz Association. + +ambigator, and this manual, are 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 . + +.SH SEE ALSO +.BR crystfel (7), +.BR indexamajig (1). -- cgit v1.2.3 From 48e14bee0fb56d9a5b4f462105be7dc27cd79347 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Mon, 10 Mar 2014 11:54:37 +0100 Subject: Set default number of iterations to 6 --- src/ambigator.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ambigator.c b/src/ambigator.c index a7aedd39..851c2d02 100644 --- a/src/ambigator.c +++ b/src/ambigator.c @@ -619,7 +619,7 @@ int main(int argc, char *argv[]) char *w_sym_str = NULL; SymOpList *w_sym; SymOpList *amb; - int n_iter = 1; + int n_iter = 6; int n_crystals, n_chunks, max_crystals; int n_dif; struct flist **crystals; -- cgit v1.2.3 From 9f768cc81e9b09a153c0ca77ccaeb91c5a0766d3 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Mon, 10 Mar 2014 13:47:39 +0100 Subject: Fix reading and writing of reflections in 2.1 stream --- libcrystfel/src/stream.c | 104 ++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 102 insertions(+), 2 deletions(-) diff --git a/libcrystfel/src/stream.c b/libcrystfel/src/stream.c index db8a4f72..80256829 100644 --- a/libcrystfel/src/stream.c +++ b/libcrystfel/src/stream.c @@ -133,6 +133,61 @@ static void write_peaks(struct image *image, FILE *ofh) } +static RefList *read_stream_reflections_2_1(FILE *fh) +{ + char *rval = NULL; + int first = 1; + RefList *out; + + out = reflist_new(); + + do { + + char line[1024]; + signed int h, k, l; + float intensity, sigma, fs, ss; + char phs[1024]; + int cts; + int r; + Reflection *refl; + + rval = fgets(line, 1023, fh); + if ( rval == NULL ) continue; + chomp(line); + + if ( strcmp(line, REFLECTION_END_MARKER) == 0 ) return out; + + r = sscanf(line, "%i %i %i %f %s %f %i %f %f", + &h, &k, &l, &intensity, phs, &sigma, &cts, &fs, &ss); + if ( (r != 9) && (!first) ) { + reflist_free(out); + return NULL; + } + + first = 0; + if ( r == 9 ) { + + double ph; + char *v; + + refl = add_refl(out, h, k, l); + set_intensity(refl, intensity); + set_detector_pos(refl, 0.0, fs, ss); + set_esd_intensity(refl, sigma); + set_redundancy(refl, cts); + + ph = strtod(phs, &v); + if ( v != phs ) set_phase(refl, deg2rad(ph)); + + } + + } while ( rval != NULL ); + + /* Got read error of some kind before finding PEAK_LIST_END_MARKER */ + return NULL; +} + + static RefList *read_stream_reflections(FILE *fh) { char *rval = NULL; @@ -214,7 +269,50 @@ static void write_stream_reflections(FILE *fh, RefList *list) h, k, l, intensity, esd_i, pk, bg, fs, ss); } +} + + +static void write_stream_reflections_2_1(FILE *fh, RefList *list) +{ + Reflection *refl; + RefListIterator *iter; + + fprintf(fh, " h k l I phase sigma(I) " + " counts fs/px ss/px\n"); + + for ( refl = first_refl(list, &iter); + refl != NULL; + refl = next_refl(refl, iter) ) + { + + signed int h, k, l; + double intensity, esd_i, ph; + int red; + double fs, ss; + char phs[16]; + int have_phase; + + get_indices(refl, &h, &k, &l); + get_detector_pos(refl, &fs, &ss); + intensity = get_intensity(refl); + esd_i = get_esd_intensity(refl); + red = get_redundancy(refl); + ph = get_phase(refl, &have_phase); + + /* Reflections with redundancy = 0 are not written */ + if ( red == 0 ) continue; + if ( have_phase ) { + snprintf(phs, 16, "%8.2f", rad2deg(ph)); + } else { + strncpy(phs, " -", 15); + } + + fprintf(fh, + "%3i %3i %3i %10.2f %s %10.2f %7i %6.1f %6.1f\n", + h, k, l, intensity, phs, esd_i, red, fs, ss); + + } } @@ -276,7 +374,9 @@ static void write_crystal(Stream *st, Crystal *cr, int include_reflections) if ( AT_LEAST_VERSION(st, 2, 2) ) { write_stream_reflections(st->fh, reflist); } else { - write_reflections_to_file(st->fh, reflist); + /* This function writes like a normal reflection + * list was written in stream 2.1 */ + write_stream_reflections_2_1(st->fh, reflist); } fprintf(st->fh, REFLECTION_END_MARKER"\n"); @@ -465,7 +565,7 @@ void read_crystal(Stream *st, struct image *image) if ( AT_LEAST_VERSION(st, 2, 2) ) { reflist = read_stream_reflections(st->fh); } else { - reflist = read_reflections_from_file(st->fh); + reflist = read_stream_reflections_2_1(st->fh); } if ( reflist == NULL ) { ERROR("Failed while reading reflections\n"); -- cgit v1.2.3 From afebe167c384a280412353d805b7c97a249e2a6d Mon Sep 17 00:00:00 2001 From: Thomas White Date: Mon, 10 Mar 2014 13:55:16 +0100 Subject: Remove old stuff --- src/ambigator.c | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/ambigator.c b/src/ambigator.c index 851c2d02..03561f6c 100644 --- a/src/ambigator.c +++ b/src/ambigator.c @@ -321,7 +321,6 @@ struct queue_args int n_to_do; int mean_nac; int nmean_nac; - gsl_permutation *p; struct cc_list *ccs; struct flist **crystals; @@ -524,7 +523,6 @@ static struct cc_list *calc_ccs(struct flist **crystals, int n_crystals, for ( i=0; i Date: Mon, 10 Mar 2014 14:52:57 +0100 Subject: Output the reindexed stream --- src/ambigator.c | 89 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 89 insertions(+) diff --git a/src/ambigator.c b/src/ambigator.c index 03561f6c..69439a24 100644 --- a/src/ambigator.c +++ b/src/ambigator.c @@ -605,6 +605,93 @@ static void detwin(struct cc_list *ccs, int n_crystals, int *assignments, } +static void reindex_reflections(FILE *fh, FILE *ofh, int assignment, + SymOpList *amb) +{ + int first = 1; + + do { + + char *rval; + char line[1024]; + int n; + signed int h, k, l; + int r; + + rval = fgets(line, 1023, fh); + if ( rval == NULL ) break; + + if ( strcmp(line, REFLECTION_END_MARKER) == 0 ) { + fputs(line, ofh); + return; + } + + if ( first ) { + fputs(line, ofh); + first = 0; + continue; + } + + r = sscanf(line, "%i %i %i%n", &h, &k, &l, &n); + + /* See scanf() manual page about %n to see why <3 is used */ + if ( (r < 3) && !first ) return; + + get_equiv(amb, NULL, assignment, h, k, l, &h, &k, &l); + + fprintf(ofh, "%4i %4i %4i%s", h, k, l, line+n); + + } while ( 1 ); +} + + +/* This is nasty, but means the output includes absolutely everything in the + * input, even stuff ignored by read_chunk() */ +static void write_reindexed_stream(const char *infile, const char *outfile, + int *assignments, SymOpList *amb) +{ + FILE *fh; + FILE *ofh; + int i; + + fh = fopen(infile, "r"); + if ( fh == NULL ) { + ERROR("Failed to open '%s'\n", infile); + return; + } + + ofh = fopen(outfile, "w"); + if ( ofh == NULL ) { + ERROR("Failed to open '%s'\n", outfile); + return; + } + + i = 0; + do { + + char *rval; + char line[1024]; + + rval = fgets(line, 1023, fh); + if ( rval == NULL ) break; + + if ( strcmp(line, REFLECTION_START_MARKER"\n") == 0 ) { + reindex_reflections(fh, ofh, assignments[i++], amb); + } else { + fputs(line, ofh); + } + + } while ( 1 ); + + if ( !feof(fh) ) { + ERROR("Error reading stream.\n"); + } + + fclose(fh); + fclose(ofh); +} + + int main(int argc, char *argv[]) { int c; @@ -906,6 +993,8 @@ int main(int argc, char *argv[]) STATUS("%i assignments are different from their starting values.\n", n_dif); + write_reindexed_stream(infile, outfile, assignments, amb); + free(assignments); gsl_rng_free(rng); -- cgit v1.2.3 From bd4b10315c16eb8e52362fb345b51182a9ee1dbc Mon Sep 17 00:00:00 2001 From: Thomas White Date: Mon, 10 Mar 2014 15:25:18 +0100 Subject: Tweak scripts/fg-graph --- scripts/fg-graph | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/scripts/fg-graph b/scripts/fg-graph index 701277c3..30ca522a 100755 --- a/scripts/fg-graph +++ b/scripts/fg-graph @@ -33,10 +33,10 @@ CMAX=0.3 gnuplot -p << EOF set terminal pngcairo size 1600,1000 enhanced monochrome font 'Helvetica,20' linewidth 2 set output "correlation.png" -set xlabel "Iteration i" +set xlabel "Number of crystals" set ylabel "Correlation" rnd(x) = x - floor(x) < 0.5 ? floor(x) : ceil(x) -round(x1,x2)=rnd(10**x2*x1)/10.0**x2 # funktion mit zwei Argumenten, x1 wird gerundet x2 anzahl der Stellen nach dem Komma +round(x1,x2)=rnd(10**x2*x1)/10.0**x2 set xzeroaxis lc rgb "black" lt 1 set key top left set xrange [0:${NITER}*${NPATT}] @@ -46,7 +46,7 @@ set xtics nomirror set x2tics nomirror set x2range [0:${NITER}] -set x2label "Iteration i gemessen in n" +set x2label "Number of passes over all crystals" set arrow from ${NPATT},${CMIN} to ${NPATT},${CMAX} nohead lc rgb "black" plot\ -- cgit v1.2.3 From cc2502e38dc93057cadf271614eea1ca3794ad13 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Mon, 10 Mar 2014 17:03:03 +0100 Subject: Update merge tests for new stream format --- tests/first_merge_check | 6 ++++-- tests/fourth_merge_check | 4 +++- tests/partialator_merge_check_1 | 6 ++++-- tests/partialator_merge_check_2 | 6 ++++-- tests/partialator_merge_check_3 | 8 +++++--- tests/second_merge_check | 6 ++++-- tests/third_merge_check | 6 ++++-- 7 files changed, 28 insertions(+), 14 deletions(-) diff --git a/tests/first_merge_check b/tests/first_merge_check index d649e284..5cf7cc70 100755 --- a/tests/first_merge_check +++ b/tests/first_merge_check @@ -35,8 +35,10 @@ End of reflections EOF cat > first_merge_check_ans.hkl << EOF - h k l I phase sigma(I) counts fs/px ss/px - 1 0 0 150.03 - 35.36 2 0.0 0.0 +CrystFEL reflection list version 2.0 +Symmetry: 1 + h k l I phase sigma(I) nmeas + 1 0 0 150.03 - 35.36 2 End of reflections EOF diff --git a/tests/fourth_merge_check b/tests/fourth_merge_check index e2d2487f..bbe32b9b 100755 --- a/tests/fourth_merge_check +++ b/tests/fourth_merge_check @@ -20,7 +20,9 @@ End of reflections EOF cat > fourth_merge_check_ans.hkl << EOF - h k l I phase sigma(I) counts fs/px ss/px +CrystFEL reflection list version 2.0 +Symmetry: 1 + h k l I phase sigma(I) nmeas End of reflections EOF diff --git a/tests/partialator_merge_check_1 b/tests/partialator_merge_check_1 index 05b74e97..f1a6adf5 100755 --- a/tests/partialator_merge_check_1 +++ b/tests/partialator_merge_check_1 @@ -32,8 +32,10 @@ EOF # We merge two patterns, without scaling or partiality, the result should just # be an average. cat > partialator_merge_check_1_ans.hkl << EOF - h k l I phase sigma(I) counts fs/px ss/px - 1 0 0 150.00 - 35.36 2 0.0 0.0 +CrystFEL reflection list version 2.0 +Symmetry: unknown + h k l I phase sigma(I) nmeas + 1 0 0 150.00 - 35.36 2 End of reflections EOF diff --git a/tests/partialator_merge_check_2 b/tests/partialator_merge_check_2 index be48249c..5513f5ec 100755 --- a/tests/partialator_merge_check_2 +++ b/tests/partialator_merge_check_2 @@ -33,8 +33,10 @@ EOF # be the mean but with the standard deviation should be zero because the scaling # factor can absorb the difference. cat > partialator_merge_check_2_ans.hkl << EOF - h k l I phase sigma(I) counts fs/px ss/px - 1 0 0 150.00 - 0.00 2 0.0 0.0 +CrystFEL reflection list version 2.0 +Symmetry: unknown + h k l I phase sigma(I) nmeas + 1 0 0 150.00 - 0.00 2 End of reflections EOF diff --git a/tests/partialator_merge_check_3 b/tests/partialator_merge_check_3 index 9a2293fd..274d7c2d 100755 --- a/tests/partialator_merge_check_3 +++ b/tests/partialator_merge_check_3 @@ -34,9 +34,11 @@ EOF # W cat > partialator_merge_check_3_ans.hkl << EOF - h k l I phase sigma(I) counts fs/px ss/px - 1 0 0 78.40 - 3.15 2 0.0 0.0 - 2 0 0 145.51 - 1.70 2 0.0 0.0 +CrystFEL reflection list version 2.0 +Symmetry: unknown + h k l I phase sigma(I) nmeas + 1 0 0 78.40 - 3.15 2 + 2 0 0 145.51 - 1.70 2 End of reflections EOF diff --git a/tests/second_merge_check b/tests/second_merge_check index aec5b120..d34cc75e 100755 --- a/tests/second_merge_check +++ b/tests/second_merge_check @@ -35,8 +35,10 @@ End of reflections EOF cat > second_merge_check_ans.hkl << EOF - h k l I phase sigma(I) counts fs/px ss/px - 1 0 0 150.03 - 35.36 2 0.0 0.0 +CrystFEL reflection list version 2.0 +Symmetry: -1 + h k l I phase sigma(I) nmeas + 1 0 0 150.03 - 35.36 2 End of reflections EOF diff --git a/tests/third_merge_check b/tests/third_merge_check index 569d5ddf..81f9ad78 100755 --- a/tests/third_merge_check +++ b/tests/third_merge_check @@ -50,8 +50,10 @@ End of reflections EOF cat > third_merge_check_ans.hkl << EOF - h k l I phase sigma(I) counts fs/px ss/px - 1 0 0 133.36 - 27.22 3 0.0 0.0 +CrystFEL reflection list version 2.0 +Symmetry: 1 + h k l I phase sigma(I) nmeas + 1 0 0 133.36 - 27.22 3 End of reflections EOF -- cgit v1.2.3 From 65d9edef292567e4ed31859bd191af5cee274d92 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Mon, 10 Mar 2014 20:09:51 +0100 Subject: Fix writing of reindexed stream --- src/ambigator.c | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/ambigator.c b/src/ambigator.c index 69439a24..77f23e8c 100644 --- a/src/ambigator.c +++ b/src/ambigator.c @@ -621,7 +621,7 @@ static void reindex_reflections(FILE *fh, FILE *ofh, int assignment, rval = fgets(line, 1023, fh); if ( rval == NULL ) break; - if ( strcmp(line, REFLECTION_END_MARKER) == 0 ) { + if ( strcmp(line, REFLECTION_END_MARKER"\n") == 0 ) { fputs(line, ofh); return; } @@ -675,10 +675,10 @@ static void write_reindexed_stream(const char *infile, const char *outfile, rval = fgets(line, 1023, fh); if ( rval == NULL ) break; + fputs(line, ofh); + if ( strcmp(line, REFLECTION_START_MARKER"\n") == 0 ) { reindex_reflections(fh, ofh, assignments[i++], amb); - } else { - fputs(line, ofh); } } while ( 1 ); -- cgit v1.2.3 From 7c19d12fd954bc08708f710c3da7cbe9769ed615 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Mon, 10 Mar 2014 20:09:59 +0100 Subject: Don't try to write reindexed stream without knowing the operator --- src/ambigator.c | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/src/ambigator.c b/src/ambigator.c index 77f23e8c..4f008633 100644 --- a/src/ambigator.c +++ b/src/ambigator.c @@ -993,7 +993,12 @@ int main(int argc, char *argv[]) STATUS("%i assignments are different from their starting values.\n", n_dif); - write_reindexed_stream(infile, outfile, assignments, amb); + if ( amb != NULL ) { + write_reindexed_stream(infile, outfile, assignments, amb); + } else { + ERROR("Can only write stream with known ambiguity operator.\n"); + ERROR("Try again with -w\n"); + } free(assignments); gsl_rng_free(rng); -- cgit v1.2.3 From 41b7558be6bf25250cc4ad6e1b720837f7a6549e Mon Sep 17 00:00:00 2001 From: Thomas White Date: Mon, 10 Mar 2014 20:13:57 +0100 Subject: Check symmetry options --- src/ambigator.c | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) diff --git a/src/ambigator.c b/src/ambigator.c index 4f008633..da7e0e11 100644 --- a/src/ambigator.c +++ b/src/ambigator.c @@ -842,6 +842,7 @@ int main(int argc, char *argv[]) return 1; } s_sym = get_pointgroup(s_sym_str); + if ( s_sym == NULL ) return 1; free(s_sym_str); if ( w_sym_str == NULL ) { @@ -852,9 +853,19 @@ int main(int argc, char *argv[]) free(w_sym_str); if ( w_sym == NULL ) return 1; amb = get_ambiguities(w_sym, s_sym); - if ( amb == NULL ) return 1; + if ( amb == NULL ) { + ERROR("Couldn't find ambiguity operator.\n"); + ERROR("Check that your values for -y and -w are " + "correct.\n"); + return 1; + } STATUS("Ambiguity operations:\n"); describe_symmetry(amb); + if ( num_equivs(amb, NULL) != 2 ) { + ERROR("There must be only one ambiguity operator.\n"); + ERROR("Try again with a different value for -w.\n"); + return 1; + } } crystals = NULL; -- cgit v1.2.3