From fb1832f80cebe48f8bf066f37973058b7ce3e334 Mon Sep 17 00:00:00 2001
From: Thomas White <taw@physics.org>
Date: Tue, 4 Mar 2014 15:09:45 +0100
Subject: Introduce "ambigator"

---
 src/ambigator.c | 406 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 406 insertions(+)
 create mode 100644 src/ambigator.c

(limited to 'src/ambigator.c')

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 <taw@physics.org>
+ *   2014 Wolfgang Brehm <wolfgang.brehm@gmail.com>
+ *
+ * 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 <http://www.gnu.org/licenses/>.
+ *
+ */
+
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <stdarg.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include <unistd.h>
+#include <getopt.h>
+#include <assert.h>
+
+#include <utils.h>
+#include <hdf5-file.h>
+#include <symmetry.h>
+#include <stream.h>
+#include <geometry.h>
+#include <peaks.h>
+#include <thread-pool.h>
+#include <beam-parameters.h>
+#include <reflist.h>
+#include <reflist-utils.h>
+
+#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=<filename>     Input stream.\n"
+"  -o, --output=<filename>    Output stream.\n"
+"  -y, --symmetry=<sym>       Apparent (\"source\") symmetry.\n"
+"  -e <sym>                   Actual (\"target\") symmetry.\n"
+"  -n, --iterations=<n>       Iterate <n> 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<n_crystals; i++ ) {
+
+		int j;
+		float f = 0.0;
+		float g = 0.0;;
+		int p = 0;
+		int q = 0;
+
+		for ( j=0; j<n_crystals; j++ ) {
+
+			float cc;
+			int n;
+
+			cc = corr(crystals[i], crystals[j], &n);
+
+			if ( i == j ) continue;
+			if ( n < 3 ) continue;
+
+			if ( assignments[i] == assignments[j] ) {
+				f += cc;
+				p++;
+			} else {
+				g += cc;
+				q++;
+			}
+
+		}
+
+		f /= p;
+		g /= q;
+
+		if ( f > 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<cur.n_crystals; i++ ) {
+
+			Crystal *cr;
+
+			cr = cur.crystals[i];
+
+			cell_free(crystal_get_cell(cr));
+
+			if ( n_crystals == max_crystals ) {
+
+				Crystal **crystals_new;
+				size_t nsz;
+
+				nsz = (max_crystals+1024)*sizeof(Crystal *);
+				crystals_new = realloc(crystals, nsz);
+				if ( crystals_new == NULL ) {
+					fprintf(stderr, "Failed to allocate "
+					        "memory for crystals.\n");
+					break;
+				}
+
+				max_crystals += 1024;
+				crystals = crystals_new;
+
+			}
+
+			crystals[n_crystals] = cr;
+			n_crystals++;
+
+		}
+
+		fprintf(stderr, "Loaded %i crystals from %i chunks\r",
+		        n_crystals, ++n_chunks);
+
+	} while ( 1 );
+	fprintf(stderr, "\n");
+
+	close_stream(st);
+
+	for ( i=0; i<n_crystals; i++ ) {
+
+		RefList *list, *merged;
+
+		list = crystal_get_reflections(crystals[i]);
+		merged = asymm_and_merge(list, s_sym);
+		reflist_free(list);
+		crystal_set_reflections(crystals[i], merged);
+
+	}
+
+	assignments = calloc(n_crystals, sizeof(int));
+	if ( assignments == NULL ) {
+		ERROR("Couldn't allocate memory for assignments.\n");
+		return 1;
+	}
+
+	for ( i=0; i<n_crystals/2; i++ ) {
+		assignments[i] = 1;
+	}
+
+	for ( i=0; i<n_iter; i++ ) {
+		detwin(crystals, n_crystals, amb, assignments);
+	}
+
+	free(assignments);
+	free(crystals);
+
+	return 0;
+}
-- 
cgit v1.2.3


From 3bcd0cb4cb5b328d54b4e8342215f47168216ace Mon Sep 17 00:00:00 2001
From: Thomas White <taw@physics.org>
Date: Tue, 4 Mar 2014 15:25:59 +0100
Subject: Fussiness

---
 src/ambigator.c | 5 ++---
 1 file changed, 2 insertions(+), 3 deletions(-)

(limited to 'src/ambigator.c')

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 <taw@physics.org>
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(-)

(limited to 'src/ambigator.c')

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<n_crystals/2; i++ ) {
-		assignments[i] = 1;
+		assignments[i] = (random_flat(rng, 1.0) > 0.5);
+		orig_assignments[i] = assignments[i];
 	}
 
 	for ( i=0; i<n_iter; i++ ) {
 		detwin(crystals, n_crystals, amb, assignments);
 	}
 
+	n_dif = 0;
+	for ( i=0; i<n_crystals/2; i++ ) {
+		if ( orig_assignments[i] != assignments[i] ) n_dif++;
+	}
+	STATUS("%i assignments are different from their starting values.\n",
+	       n_dif);
+
 	free(assignments);
 	free(crystals);
+	gsl_rng_free(rng);
 
 	return 0;
 }
-- 
cgit v1.2.3


From c0356bb0ea8627eba0752b63c7fb0fa86088183f Mon Sep 17 00:00:00 2001
From: Thomas White <taw@physics.org>
Date: Tue, 4 Mar 2014 15:42:51 +0100
Subject: D'oh

---
 src/ambigator.c | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

(limited to 'src/ambigator.c')

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<n_crystals/2; i++ ) {
+	for ( i=0; i<n_crystals; i++ ) {
 		assignments[i] = (random_flat(rng, 1.0) > 0.5);
 		orig_assignments[i] = assignments[i];
 	}
@@ -410,7 +410,7 @@ int main(int argc, char *argv[])
 	}
 
 	n_dif = 0;
-	for ( i=0; i<n_crystals/2; i++ ) {
+	for ( i=0; i<n_crystals; i++ ) {
 		if ( orig_assignments[i] != assignments[i] ) n_dif++;
 	}
 	STATUS("%i assignments are different from their starting values.\n",
-- 
cgit v1.2.3


From 84253d30dc723d0a48776490c9d8b4f8164fa077 Mon Sep 17 00:00:00 2001
From: Thomas White <taw@physics.org>
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(-)

(limited to 'src/ambigator.c')

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<cur.n_crystals; i++ ) {
 
 			Crystal *cr;
+			RefList *list;
 
 			cr = cur.crystals[i];
 
@@ -347,10 +405,10 @@ int main(int argc, char *argv[])
 
 			if ( n_crystals == max_crystals ) {
 
-				Crystal **crystals_new;
+				struct flist **crystals_new;
 				size_t nsz;
 
-				nsz = (max_crystals+1024)*sizeof(Crystal *);
+				nsz = (max_crystals+1024)*sizeof(struct flist *);
 				crystals_new = realloc(crystals, nsz);
 				if ( crystals_new == NULL ) {
 					fprintf(stderr, "Failed to allocate "
@@ -363,9 +421,12 @@ int main(int argc, char *argv[])
 
 			}
 
-			crystals[n_crystals] = cr;
+			list = crystal_get_reflections(cr);
+			crystals[n_crystals] = asymm_and_merge(list, s_sym);
 			n_crystals++;
 
+			reflist_free(list);
+
 		}
 
 		fprintf(stderr, "Loaded %i crystals from %i chunks\r",
@@ -376,17 +437,6 @@ int main(int argc, char *argv[])
 
 	close_stream(st);
 
-	for ( i=0; i<n_crystals; i++ ) {
-
-		RefList *list, *merged;
-
-		list = crystal_get_reflections(crystals[i]);
-		merged = asymm_and_merge(list, s_sym);
-		reflist_free(list);
-		crystal_set_reflections(crystals[i], merged);
-
-	}
-
 	assignments = malloc(n_crystals*sizeof(int));
 	if ( assignments == NULL ) {
 		ERROR("Couldn't allocate memory for assignments.\n");
-- 
cgit v1.2.3


From 079aceb278183a61f29a31b42dca094ea08c6a67 Mon Sep 17 00:00:00 2001
From: Thomas White <taw@physics.org>
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(-)

(limited to 'src/ambigator.c')

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 <beam-parameters.h>
 #include <reflist.h>
 #include <reflist-utils.h>
+#include <cell-utils.h>
 
 #include "post-refinement.h"
 #include "hrs-scaling.h"
@@ -70,6 +71,8 @@ static void show_help(const char *s)
 "  -y, --symmetry=<sym>       Apparent (\"source\") symmetry.\n"
 "  -e <sym>                   Actual (\"target\") symmetry.\n"
 "  -n, --iterations=<n>       Iterate <n> times.\n"
+"      --highres=<n>          High resolution cutoff in A.\n"
+"      --lowres=<n>           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 <taw@physics.org>
Date: Tue, 4 Mar 2014 18:35:04 +0100
Subject: Fix the algorithm

---
 src/ambigator.c | 17 ++++++++++++-----
 1 file changed, 12 insertions(+), 5 deletions(-)

(limited to 'src/ambigator.c')

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<n_crystals; i++ ) {
 
@@ -244,10 +246,11 @@ static void detwin(struct flist **crystals, int n_crystals, SymOpList *amb,
 			float cc;
 			int n;
 
+			if ( i == j ) continue;
+
 			cc = corr(crystals[i], crystals[j], &n);
 
 			if ( n < 3 ) continue;
-			if ( i == j ) continue;
 
 			if ( assignments[i] == assignments[j] ) {
 				f += cc;
@@ -261,8 +264,12 @@ static void detwin(struct flist **crystals, int n_crystals, SymOpList *amb,
 
 		f /= p;
 		g /= q;
+		if ( (p==0) || (q==0) ) continue;
+
+		mf += f;
+		nmf++;
 
-		if ( f > 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 <taw@physics.org>
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(-)

(limited to 'src/ambigator.c')

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=<filename>     Input stream.\n"
 "  -o, --output=<filename>    Output stream.\n"
 "  -y, --symmetry=<sym>       Apparent (\"source\") symmetry.\n"
 "  -e <sym>                   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 <taw@physics.org>
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(-)

(limited to 'src/ambigator.c')

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<n_crystals; j++ ) {
 
 			float cc;
@@ -273,8 +275,6 @@ static void detwin(struct flist **crystals, int n_crystals, SymOpList *amb,
 			nch++;
 		}
 
-		progress_bar(i, n_crystals-1, "Calculating");
-
 	}
 
 	STATUS("Mean f = %f, changed %i assignments this time.\n", mf/nmf, nch);
-- 
cgit v1.2.3


From bd889c3bb7653880068498d8386a81c65045167f Mon Sep 17 00:00:00 2001
From: Thomas White <taw@physics.org>
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(-)

(limited to 'src/ambigator.c')

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=<filename>    Output stream.\n"
-"  -y, --symmetry=<sym>       Apparent (\"source\") symmetry.\n"
-"  -e <sym>                   Actual (\"target\") symmetry.\n"
+"  -y, --symmetry=<sym>       Actual (\"target\") symmetry.\n"
+"  -w <sym>                   Apparent (\"source\" or \"twinned\") symmetry.\n"
 "  -n, --iterations=<n>       Iterate <n> times.\n"
 "      --highres=<n>          High resolution cutoff in A.\n"
 "      --lowres=<n>           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 <taw@physics.org>
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(+)

(limited to 'src/ambigator.c')

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<n_crystals; i++ ) {
 
@@ -263,6 +264,11 @@ static void detwin(struct flist **crystals, int n_crystals, SymOpList *amb,
 
 		}
 
+		if ( (p==0) || (q==0) ) {
+			ndud++;
+			continue;
+		}
+
 		f /= p;
 		g /= q;
 		if ( (p==0) || (q==0) ) continue;
@@ -277,6 +283,10 @@ static void detwin(struct flist **crystals, int n_crystals, SymOpList *amb,
 
 	}
 
+	if ( ndud > 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 <taw@physics.org>
Date: Wed, 5 Mar 2014 13:59:55 +0100
Subject: Fix includes

---
 src/ambigator.c | 11 ++---------
 1 file changed, 2 insertions(+), 9 deletions(-)

(limited to 'src/ambigator.c')

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 <getopt.h>
 #include <assert.h>
 
+#include <image.h>
 #include <utils.h>
-#include <hdf5-file.h>
 #include <symmetry.h>
 #include <stream.h>
-#include <geometry.h>
-#include <peaks.h>
-#include <thread-pool.h>
-#include <beam-parameters.h>
 #include <reflist.h>
 #include <reflist-utils.h>
+#include <cell.h>
 #include <cell-utils.h>
 
-#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 <taw@physics.org>
Date: Wed, 5 Mar 2014 14:16:01 +0100
Subject: Add an assertion

---
 src/ambigator.c | 1 +
 1 file changed, 1 insertion(+)

(limited to 'src/ambigator.c')

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 <taw@physics.org>
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(-)

(limited to 'src/ambigator.c')

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 <taw@physics.org>
Date: Wed, 5 Mar 2014 14:24:39 +0100
Subject: e->w

---
 src/ambigator.c | 18 +++++++++---------
 1 file changed, 9 insertions(+), 9 deletions(-)

(limited to 'src/ambigator.c')

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 <taw@physics.org>
Date: Wed, 5 Mar 2014 15:00:50 +0100
Subject: Add reindexing

---
 src/ambigator.c | 75 ++++++++++++++++++++++++++++++++++++++++++---------------
 1 file changed, 55 insertions(+), 20 deletions(-)

(limited to 'src/ambigator.c')

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<n_crystals; j++ ) {
 
-			float cc;
-			int n;
+			int n, n_reidx;
+			float cc, cc_reidx;
 
 			if ( i == j ) continue;
 
-			cc = corr(crystals[i], crystals[j], &n);
+			cc = corr(crystals[i], crystals[j], &n, 0);
+			cc_reidx = corr(crystals[i], crystals[j], &n_reidx, 1);
+
+			if ( n > 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<n_iter; i++ ) {
-		detwin(crystals, n_crystals, amb, assignments);
+		detwin(crystals, n_crystals, assignments);
 	}
 
 	n_dif = 0;
-- 
cgit v1.2.3


From e93d7bca330bb8718196815c3e3e917c3ae97128 Mon Sep 17 00:00:00 2001
From: Thomas White <taw@physics.org>
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(-)

(limited to 'src/ambigator.c')

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<n_iter; i++ ) {
-		detwin(crystals, n_crystals, assignments);
+		detwin(crystals, n_crystals, assignments, amb);
 	}
 
 	n_dif = 0;
-- 
cgit v1.2.3


From fba54b335b3c9dcdbf31d7ae92b6ed307b2efd21 Mon Sep 17 00:00:00 2001
From: Thomas White <taw@physics.org>
Date: Wed, 5 Mar 2014 16:13:45 +0100
Subject: Fix reindexing

---
 src/ambigator.c | 65 +++++++++++++++++++++++++++++++++++++++++++--------------
 1 file changed, 49 insertions(+), 16 deletions(-)

(limited to 'src/ambigator.c')

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 <taw@physics.org>
Date: Wed, 5 Mar 2014 16:29:41 +0100
Subject: Add --end-assignments

---
 src/ambigator.c | 18 ++++++++++++++++++
 1 file changed, 18 insertions(+)

(limited to 'src/ambigator.c')

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<n_crystals; i++ ) {
+				fprintf(fh, "%i\n", assignments[i]);
+			}
+		}
+		fclose(fh);
+	}
+
 	n_dif = 0;
 	for ( i=0; i<n_crystals; i++ ) {
 		if ( orig_assignments[i] != assignments[i] ) n_dif++;
-- 
cgit v1.2.3


From c4649fe727361c43482b0508f1a5d3e82d8cb57a Mon Sep 17 00:00:00 2001
From: Thomas White <taw@physics.org>
Date: Thu, 6 Mar 2014 12:08:08 +0100
Subject: Update --help

---
 src/ambigator.c | 1 +
 1 file changed, 1 insertion(+)

(limited to 'src/ambigator.c')

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=<n>       Iterate <n> times.\n"
 "      --highres=<n>          High resolution cutoff in A.\n"
 "      --lowres=<n>           Low resolution cutoff in A.\n"
+"      --end-assignments=<fn> Save end assignments to file <fn>.\n"
 );
 }
 
-- 
cgit v1.2.3


From ad347841d21095fe414e09f42eccd08883a2eb51 Mon Sep 17 00:00:00 2001
From: Thomas White <taw@physics.org>
Date: Thu, 6 Mar 2014 12:08:51 +0100
Subject: Add --fg-graph

---
 src/ambigator.c | 25 +++++++++++++++++++++++--
 1 file changed, 23 insertions(+), 2 deletions(-)

(limited to 'src/ambigator.c')

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=<n>          High resolution cutoff in A.\n"
 "      --lowres=<n>           Low resolution cutoff in A.\n"
 "      --end-assignments=<fn> Save end assignments to file <fn>.\n"
+"      --fg-graph=<fn>        Save f and g correlation values to file <fn>.\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<n_iter; i++ ) {
-		detwin(crystals, n_crystals, assignments, amb);
+		detwin(crystals, n_crystals, assignments, amb, fgfh);
+	}
+
+	if ( fgfh != NULL ) {
+		fclose(fgfh);
 	}
 
 	if ( end_ass_fn != NULL ) {
-- 
cgit v1.2.3


From 8591e4c373f9b9921569b2d8c8e7b57bb555d563 Mon Sep 17 00:00:00 2001
From: Thomas White <taw@physics.org>
Date: Thu, 6 Mar 2014 14:45:44 +0100
Subject: Ignore "twin-proof" reflections

---
 src/ambigator.c | 17 +++++++++++++++++
 1 file changed, 17 insertions(+)

(limited to 'src/ambigator.c')

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 <taw@physics.org>
Date: Thu, 6 Mar 2014 16:13:20 +0100
Subject: Add --ncorr

---
 src/ambigator.c | 162 ++++++++++++++++++++++++++++++++++++++++++++------------
 1 file changed, 128 insertions(+), 34 deletions(-)

(limited to 'src/ambigator.c')

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=<n>           Low resolution cutoff in A.\n"
 "      --end-assignments=<fn> Save end assignments to file <fn>.\n"
 "      --fg-graph=<fn>        Save f and g correlation values to file <fn>.\n"
+"      --ncorr=<n>            Use <n> 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<n_crystals; i++ ) {
+
+		int j, k;
+
+		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 NULL;
+		}
+
+		k = 0;
+		for ( j=(i+1)%n_crystals; j!=i; j=(j+1)%n_crystals ) {
+
+			int n;
+			float cc;
+
+			cc = corr(crystals[i], crystals[j], &n, 0);
+
+			if ( n < 4 ) continue;
+
+			ccs[i].ind[k] = j+1;
+			ccs[i].cc[k] = cc;
+			k++;
+
+			if ( k == ncorr-1 ) break;
+
+		}
+		ccs[i].ind[k] = 0;
+
+		if ( amb != NULL ) {
+
+			k = 0;
+			for ( j=(i+1)%n_crystals; j!=i; j=(j+1)%n_crystals ) {
+
+				int n;
+				float cc;
+
+				cc = corr(crystals[i], crystals[j], &n, 1);
+
+				if ( n < 4 ) continue;
+
+				ccs[i].ind_reidx[k] = j+1;
+				ccs[i].cc_reidx[k] = cc;
+				k++;
+
+				if ( k == ncorr-1 ) break;
+
+			}
+			ccs[i].ind_reidx[k] = 0;
+
+		}
+
+		progress_bar(i, n_crystals-1, "Calculating CCs");
+
+	}
+
+	return ccs;
+}
+
+
+static void detwin(struct cc_list *ccs, int n_crystals, int *assignments,
+                   FILE *fh, struct flist **crystals)
 {
 	int i;
 	int nch = 0;
@@ -303,50 +389,40 @@ static void detwin(struct flist **crystals, int n_crystals, int *assignments,
 
 	for ( i=0; i<n_crystals; i++ ) {
 
-		int j;
+		int k;
 		float f = 0.0;
 		float g = 0.0;;
 		int p = 0;
 		int q = 0;
 
-		progress_bar(i, n_crystals-1, "Calculating");
+		//progress_bar(i, n_crystals-1, "Calculating");
 
-		for ( j=0; j<n_crystals; j++ ) {
+		for ( k=0; (ccs[i].ind[k] != 0); k++ ) {
 
-			int n, n_reidx;
-			float cc, cc_reidx;
-
-			if ( i == j ) continue;
-
-			cc = corr(crystals[i], crystals[j], &n, 0);
+			int j = ccs[i].ind[k]-1;
+			float cc = ccs[i].cc[k];
 
-			if ( amb != NULL ) {
-				cc_reidx = corr(crystals[i], crystals[j],
-				                &n_reidx, 1);
+			if ( assignments[i] == assignments[j] ) {
+				f += cc;
+				p++;
+			} else {
+				g += cc;
+				q++;
 			}
 
-			if ( n > 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<n_iter; i++ ) {
-		detwin(crystals, n_crystals, assignments, amb, fgfh);
+		detwin(ccs, n_crystals, assignments, fgfh, crystals);
 	}
 
 	if ( fgfh != NULL ) {
-- 
cgit v1.2.3


From fc647d846dc8fbec9ef13c41382be260ff5e8e9a Mon Sep 17 00:00:00 2001
From: Thomas White <taw@physics.org>
Date: Thu, 6 Mar 2014 17:55:51 +0100
Subject: Add --stop-after

---
 src/ambigator.c | 15 ++++++++++++++-
 1 file changed, 14 insertions(+), 1 deletion(-)

(limited to 'src/ambigator.c')

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=<fn> Save end assignments to file <fn>.\n"
 "      --fg-graph=<fn>        Save f and g correlation values to file <fn>.\n"
 "      --ncorr=<n>            Use <n> correlations per crystal.  Default 1000\n"
+"      --stop-after=<n>       Use at most the first <n> 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 <taw@physics.org>
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(+)

(limited to 'src/ambigator.c')

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 <taw@physics.org>
Date: Thu, 6 Mar 2014 17:56:08 +0100
Subject: Remove unused line

---
 src/ambigator.c | 2 --
 1 file changed, 2 deletions(-)

(limited to 'src/ambigator.c')

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 <taw@physics.org>
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(-)

(limited to 'src/ambigator.c')

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 <unistd.h>
 #include <getopt.h>
 #include <assert.h>
+#include <gsl/gsl_rng.h>
+#include <gsl/gsl_permutation.h>
+#include <gsl/gsl_randist.h>
 
 #include <image.h>
 #include <utils.h>
@@ -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; i<n_crystals; i++ ) {
 
-		int j, k;
+		int k, l;
 
 		ccs[i].ind = malloc(ncorr*sizeof(int));
 		ccs[i].cc = malloc(ncorr*sizeof(float));
@@ -335,11 +342,17 @@ static struct cc_list *calc_ccs(struct flist **crystals, int n_crystals,
 		}
 
 		k = 0;
-		for ( j=(i+1)%n_crystals; j!=i; j=(j+1)%n_crystals ) {
+		gsl_ran_shuffle(rng, p->data, n_crystals, sizeof(size_t));
+
+		for ( l=0; l<n_crystals; l++ ) {
 
 			int n;
+			int j;
 			float cc;
 
+			j = gsl_permutation_get(p, l);
+			if ( i == j ) continue;
+
 			cc = corr(crystals[i], crystals[j], &n, 0);
 
 			if ( n < 4 ) continue;
@@ -356,11 +369,15 @@ static struct cc_list *calc_ccs(struct flist **crystals, int n_crystals,
 		if ( amb != NULL ) {
 
 			k = 0;
-			for ( j=(i+1)%n_crystals; j!=i; j=(j+1)%n_crystals ) {
+			for ( l=0; l<n_crystals; l++ ) {
 
 				int n;
+				int j;
 				float cc;
 
+				j = gsl_permutation_get(p, l);
+				if ( i == j ) continue;
+
 				cc = corr(crystals[i], crystals[j], &n, 1);
 
 				if ( n < 4 ) continue;
@@ -380,6 +397,8 @@ static struct cc_list *calc_ccs(struct flist **crystals, int n_crystals,
 
 	}
 
+	gsl_permutation_free(p);
+
 	return ccs;
 }
 
@@ -713,7 +732,7 @@ int main(int argc, char *argv[])
 		}
 	}
 
-	ccs = calc_ccs(crystals, n_crystals, ncorr, amb);
+	ccs = calc_ccs(crystals, n_crystals, ncorr, amb, rng);
 	if ( ccs == NULL ) {
 		ERROR("Failed to allocate CCs\n");
 		return 1;
-- 
cgit v1.2.3


From 017ab36327202c4548460b7a12ae21db066fc0d2 Mon Sep 17 00:00:00 2001
From: Thomas White <taw@physics.org>
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(-)

(limited to 'src/ambigator.c')

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 <taw@physics.org>
Date: Thu, 6 Mar 2014 19:57:49 +0100
Subject: Fussiness

---
 src/ambigator.c | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

(limited to 'src/ambigator.c')

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<n_crystals; l++ ) {
 
 			int n;
-- 
cgit v1.2.3


From 6563d65f5efd836c7152c0ab75ea9a3afa995879 Mon Sep 17 00:00:00 2001
From: Thomas White <taw@physics.org>
Date: Thu, 6 Mar 2014 21:33:44 +0100
Subject: Fix segfault

---
 src/ambigator.c | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

(limited to 'src/ambigator.c')

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 <taw@physics.org>
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(-)

(limited to 'src/ambigator.c')

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 <reflist-utils.h>
 #include <cell.h>
 #include <cell-utils.h>
+#include <thread-pool.h>
 
 
 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; l<n_crystals; l++ ) {
 
-	p = gsl_permutation_alloc(n_crystals);
-	gsl_permutation_init(p);
+		int n;
+		int j;
+		float cc;
 
-	for ( i=0; i<n_crystals; i++ ) {
+		j = gsl_permutation_get(p, l);
+		if ( i == j ) continue;
 
-		int k, l;
+		cc = corr(crystals[i], crystals[j], &n, 0);
 
-		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 NULL;
-		}
+		if ( n < 4 ) continue;
+
+		ccs[i].ind[k] = j+1;
+		ccs[i].cc[k] = cc;
+		k++;
+
+		if ( k == ncorr-1 ) break;
 
-		gsl_ran_shuffle(rng, p->data, n_crystals, sizeof(size_t));
+	}
+	ccs[i].ind[k] = 0;
+	mean_nac += k;
+	nmean_nac++;
+
+	if ( amb != NULL ) {
 
 		k = 0;
 		for ( l=0; l<n_crystals; l++ ) {
@@ -356,57 +450,64 @@ static struct cc_list *calc_ccs(struct flist **crystals, int n_crystals,
 			j = gsl_permutation_get(p, l);
 			if ( i == j ) continue;
 
-			cc = corr(crystals[i], crystals[j], &n, 0);
+			cc = corr(crystals[i], crystals[j], &n, 1);
 
 			if ( n < 4 ) continue;
 
-			ccs[i].ind[k] = j+1;
-			ccs[i].cc[k] = cc;
+			ccs[i].ind_reidx[k] = j+1;
+			ccs[i].cc_reidx[k] = cc;
 			k++;
 
 			if ( k == ncorr-1 ) break;
 
 		}
-		ccs[i].ind[k] = 0;
+		ccs[i].ind_reidx[k] = 0;
 		mean_nac += k;
 		nmean_nac++;
 
-		if ( amb != NULL ) {
-
-			k = 0;
-			for ( l=0; l<n_crystals; l++ ) {
-
-				int n;
-				int j;
-				float cc;
+	}
 
-				j = gsl_permutation_get(p, l);
-				if ( i == j ) continue;
+	job->mean_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 <taw@physics.org>
Date: Sun, 9 Mar 2014 21:36:50 +0100
Subject: Add -j

---
 src/ambigator.c | 17 +++++++++++++----
 1 file changed, 13 insertions(+), 4 deletions(-)

(limited to 'src/ambigator.c')

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=<fn>        Save f and g correlation values to file <fn>.\n"
 "      --ncorr=<n>            Use <n> correlations per crystal.  Default 1000\n"
 "      --stop-after=<n>       Use at most the first <n> crystals.\n"
+"  -j <n>                     Use <n> 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 <taw@physics.org>
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(-)

(limited to 'src/ambigator.c')

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<n; i++ ) {
+		rngs[i] = gsl_rng_alloc(gsl_rng_mt19937);
+		gsl_rng_set(rngs[i], gsl_rng_get(rng));
+	}
+
+	return rngs;
+}
+
+
 static struct cc_list *calc_ccs(struct flist **crystals, int n_crystals,
                                 int ncorr, SymOpList *amb, gsl_rng *rng,
                                 float *pmean_nac, int nthreads)
 {
 	struct cc_list *ccs;
 	struct queue_args qargs;
+	int i;
 
 	assert(n_crystals >= 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<nthreads; i++ ) {
+		gsl_rng_free(qargs.rngs[i]);
+	}
 	gsl_permutation_free(qargs.p);
 
 	*pmean_nac = (float)qargs.mean_nac/qargs.nmean_nac;
-- 
cgit v1.2.3


From f9f8eac82c1b5a947965626fa3d471e4a1f33988 Mon Sep 17 00:00:00 2001
From: Thomas White <taw@physics.org>
Date: Sun, 9 Mar 2014 21:51:08 +0100
Subject: Fussiness

---
 src/ambigator.c | 1 -
 1 file changed, 1 deletion(-)

(limited to 'src/ambigator.c')

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 <taw@bitwiz.org.uk>
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(-)

(limited to 'src/ambigator.c')

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 <taw@bitwiz.org.uk>
Date: Mon, 10 Mar 2014 10:44:03 +0100
Subject: Free the crystals

---
 src/ambigator.c | 10 ++++++++--
 1 file changed, 8 insertions(+), 2 deletions(-)

(limited to 'src/ambigator.c')

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; i<n_crystals; i++ ) {
+		free(crystals[i]->s);
+		free(crystals[i]->i);
+		free(crystals[i]->s_reidx);
+		free(crystals[i]->i_reidx);
+		free(crystals[i]);
+	}
+	free(crystals);
 
 	for ( i=0; i<n_iter; i++ ) {
 		detwin(ccs, n_crystals, assignments, fgfh, crystals);
@@ -915,7 +922,6 @@ int main(int argc, char *argv[])
 	       n_dif);
 
 	free(assignments);
-	free(crystals);
 	gsl_rng_free(rng);
 
 	return 0;
-- 
cgit v1.2.3


From 8839635714fe802f088fc2ce4d788a57224482f8 Mon Sep 17 00:00:00 2001
From: Thomas White <taw@bitwiz.org.uk>
Date: Mon, 10 Mar 2014 10:59:12 +0100
Subject: Remove --stop-after

---
 src/ambigator.c | 13 -------------
 1 file changed, 13 deletions(-)

(limited to 'src/ambigator.c')

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 48e14bee0fb56d9a5b4f462105be7dc27cd79347 Mon Sep 17 00:00:00 2001
From: Thomas White <taw@bitwiz.org.uk>
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(-)

(limited to 'src/ambigator.c')

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 afebe167c384a280412353d805b7c97a249e2a6d Mon Sep 17 00:00:00 2001
From: Thomas White <taw@physics.org>
Date: Mon, 10 Mar 2014 13:55:16 +0100
Subject: Remove old stuff

---
 src/ambigator.c | 2 --
 1 file changed, 2 deletions(-)

(limited to 'src/ambigator.c')

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<nthreads; i++ ) {
 		gsl_rng_free(qargs.rngs[i]);
 	}
-	gsl_permutation_free(qargs.p);
 
 	*pmean_nac = (float)qargs.mean_nac/qargs.nmean_nac;
 
-- 
cgit v1.2.3


From 9a193a4b42a335a9367f8aefed7fa1029aff7f9a Mon Sep 17 00:00:00 2001
From: Thomas White <taw@physics.org>
Date: Mon, 10 Mar 2014 14:52:57 +0100
Subject: Output the reindexed stream

---
 src/ambigator.c | 89 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 89 insertions(+)

(limited to 'src/ambigator.c')

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 65d9edef292567e4ed31859bd191af5cee274d92 Mon Sep 17 00:00:00 2001
From: Thomas White <taw@bitwiz.org.uk>
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(-)

(limited to 'src/ambigator.c')

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 <taw@bitwiz.org.uk>
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(-)

(limited to 'src/ambigator.c')

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 <taw@bitwiz.org.uk>
Date: Mon, 10 Mar 2014 20:13:57 +0100
Subject: Check symmetry options

---
 src/ambigator.c | 13 ++++++++++++-
 1 file changed, 12 insertions(+), 1 deletion(-)

(limited to 'src/ambigator.c')

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