Correlate against randomly chosen crystals instead of the next few

1 files 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 <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;