Work on making iteration work

9 files changed, 212 insertions, 52 deletions

diff --git a/src/geometry.c b/src/geometry.c
index 9f9b953c..1cf839c2 100644
--- a/src/geometry.c
+++ b/src/geometry.c
@@ -263,10 +263,11 @@ double integrate_all(struct image *image, RefList *reflections)
 {
 	double itot = 0.0;
 	Reflection *refl;
+	RefListIterator *iter;
 
-	for ( refl = first_refl(reflections);
+	for ( refl = first_refl(reflections, &iter);
 	      refl != NULL;
-	      refl = next_refl(refl) ) {
+	      refl = next_refl(refl, iter) ) {
 
 		float x, y, intensity;
 		double xp, yp;
diff --git a/src/hrs-scaling.c b/src/hrs-scaling.c
index d5d96ffc..9f1ada75 100644
--- a/src/hrs-scaling.c
+++ b/src/hrs-scaling.c
@@ -52,6 +52,7 @@ static void s_uhavha(signed int hat, signed int kat, signed int lat,
 		if ( !get_scalable(refl) ) continue;
 
 		ic = get_intensity(refl) / get_partiality(refl);
+		STATUS("%f / %f = %f\n", get_intensity(refl), get_partiality(refl), ic);
 		sigi = sqrt(fabs(ic));
 
 		uha_val += 1.0 / pow(sigi, 2.0);
@@ -117,6 +118,7 @@ static double iterate_scale(struct image *images, int n,
 
 	uh_arr = new_list_intensity();
 	vh_arr = new_list_intensity();
+	STATUS("%i\n", n_ref);
 	for ( h=0; h<n_ref; h++ ) {
 
 		double uh, vh;
@@ -124,6 +126,7 @@ static double iterate_scale(struct image *images, int n,
 
 		uh = s_uh(images, n, it->h, it->k, it->l, sym);
 		vh = s_vh(images, n, it->h, it->k, it->l, sym);
+		STATUS("%f %f\n", uh, vh);
 
 		set_intensity(uh_arr, it->h, it->k, it->l, uh);
 		set_intensity(vh_arr, it->h, it->k, it->l, vh);
diff --git a/src/partialator.c b/src/partialator.c
index 37c8fbd6..3251f7e8 100644
--- a/src/partialator.c
+++ b/src/partialator.c
@@ -116,10 +116,11 @@ static void select_scalable_reflections(struct image *images, int n)
 	for ( m=0; m<n; m++ ) {
 
 		Reflection *refl;
+		RefListIterator *iter;
 
-		for ( refl = first_refl(images[m].reflections);
+		for ( refl = first_refl(images[m].reflections, &iter);
 		      refl != NULL;
-		      refl = next_refl(refl) ) {
+		      refl = next_refl(refl, iter) ) {
 
 			int scalable = 1;
 			double v;
@@ -275,6 +276,7 @@ int main(int argc, char *argv[])
 		RefList *peaks;
 		RefList *transfer;
 		Reflection *refl;
+		RefListIterator *iter;
 
 		if ( find_chunk(fh, &cell, &filename) == 1 ) {
 			ERROR("Couldn't get all of the filenames and cells"
@@ -327,9 +329,9 @@ int main(int argc, char *argv[])
 		transfer = find_intersections(&images[i], cell, 0);
 		images[i].reflections = reflist_new();
 
-		for ( refl = first_refl(transfer);
+		for ( refl = first_refl(transfer, &iter);
 		      refl != NULL;
-		      refl = next_refl(refl) ) {
+		      refl = next_refl(refl, iter) ) {
 
 			Reflection *peak;
 			Reflection *new;
@@ -338,6 +340,8 @@ int main(int argc, char *argv[])
 			int clamp1, clamp2;
 
 			get_indices(refl, &h, &k, &l);
+			STATUS("%3i %3i %3i\n", h, k, l);
+
 			peak = find_refl(peaks, h, k, l);
 			if ( peak == NULL ) {
 				if ( (h==0) && (k==0) && (l==0) ) continue;
diff --git a/src/peaks.c b/src/peaks.c
index 06c5ee35..4ae0a70b 100644
--- a/src/peaks.c
+++ b/src/peaks.c
@@ -669,6 +669,7 @@ void output_intensities(struct image *image, UnitCell *cell,
 	double bsx, bsy, bsz;
 	double csx, csy, csz;
 	Reflection *refl;
+	RefListIterator *iter;
 
 	if ( image->reflections == NULL ) {
 		find_projected_peaks(image, cell, circular_domain, domain_r);
@@ -683,9 +684,9 @@ void output_intensities(struct image *image, UnitCell *cell,
 	                          &bsx, &bsy, &bsz,
 	                          &csx, &csy, &csz);
 
-	for ( refl = first_refl(image->reflections);
+	for ( refl = first_refl(image->reflections, &iter);
 	      refl != NULL;
-	      refl = next_refl(refl) ) {
+	      refl = next_refl(refl, iter) ) {
 
 		float x, y, intensity;
 		double d;
diff --git a/src/post-refinement.c b/src/post-refinement.c
index 4fa7ad47..0bd26f76 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -257,6 +257,7 @@ static double pr_iterate(struct image *image, const double *i_full,
 	gsl_vector *shifts;
 	int param;
 	Reflection *refl;
+	RefListIterator *iter;
 	RefList *reflections;
 	double max_shift;
 
@@ -266,9 +267,9 @@ static double pr_iterate(struct image *image, const double *i_full,
 	v = gsl_vector_calloc(NUM_PARAMS);
 
 	/* Construct the equations, one per reflection in this image */
-	for ( refl = first_refl(reflections);
+	for ( refl = first_refl(reflections, &iter);
 	      refl != NULL;
-	      refl = next_refl(refl) ) {
+	      refl = next_refl(refl, iter) ) {
 
 		signed int hind, kind, lind;
 		signed int ha, ka, la;
diff --git a/src/reflist.c b/src/reflist.c
index e543c20d..4b08b3af 100644
--- a/src/reflist.c
+++ b/src/reflist.c
@@ -24,6 +24,7 @@ struct _reflection {
 	struct _reflection *parent;   /* Parent node */
 	struct _reflection *next;     /* Another reflection with the same
 	                               * indices, or NULL */
+	struct _reflection *prev;
 
 	signed int h;
 	signed int k;
@@ -71,6 +72,7 @@ static Reflection *new_node(unsigned int serial)
 	new = calloc(1, sizeof(struct _reflection));
 	new->serial = serial;
 	new->next = NULL;
+	new->prev = NULL;
 	new->child[0] = NULL;
 	new->child[1] = NULL;
 
@@ -110,24 +112,36 @@ void reflist_free(RefList *list)
 
 /********************************** Search ************************************/
 
-static Reflection *recursive_search(Reflection *refl, unsigned int search)
+/* Return the first reflection in 'list' with the given indices, or NULL */
+Reflection *find_refl(RefList *list, INDICES)
 {
-	int i;
+	unsigned int search = SERIAL(h, k, l);
+	Reflection *refl = list->head->child[0];
 
-	/* Hit the bottom of the tree? */
-	if ( refl == NULL ) return NULL;
+	while ( refl != NULL ) {
 
-	/* Is this the correct reflection? */
-	if ( refl->serial == search ) return refl;
+		if ( search < refl->serial ) {
 
-	for ( i=0; i<2; i++ ) {
+			if ( refl->child[0] != NULL ) {
+				refl = refl->child[0];
+			} else {
+				/* Hit the bottom of the tree */
+				return NULL;
+			}
 
-		if ( search < refl->serial ) {
-			return recursive_search(refl->child[0], search);
-		}
+		} else if ( search > refl->serial ) {
+
+			if ( refl->child[1] != NULL ) {
+				refl = refl->child[1];
+			} else {
+				/* Hit the bottom of the tree */
+				return NULL;
+			}
+
+		} else {
+
+			return refl;
 
-		if ( search >= refl->serial ) {
-			return recursive_search(refl->child[1], search);
 		}
 
 	}
@@ -136,14 +150,6 @@ static Reflection *recursive_search(Reflection *refl, unsigned int search)
 }
 
 
-/* Return the first reflection in 'list' with the given indices, or NULL */
-Reflection *find_refl(RefList *list, INDICES)
-{
-	unsigned int search = SERIAL(h, k, l);
-	return recursive_search(list->head, search);
-}
-
-
 /* Find the next reflection in 'refl's list with the same indices, or NULL */
 Reflection *next_found_refl(Reflection *refl)
 {
@@ -247,6 +253,7 @@ static void insert_node(Reflection *head, Reflection *new)
 	while ( refl != NULL ) {
 
 		if ( new->serial < refl->serial ) {
+
 			if ( refl->child[0] != NULL ) {
 				refl = refl->child[0];
 			} else {
@@ -254,7 +261,9 @@ static void insert_node(Reflection *head, Reflection *new)
 				new->parent = refl;
 				return;
 			}
-		} else {
+
+		} else if ( new->serial > refl->serial ) {
+
 			if ( refl->child[1] != NULL ) {
 				refl = refl->child[1];
 			} else {
@@ -262,8 +271,20 @@ static void insert_node(Reflection *head, Reflection *new)
 				new->parent = refl;
 				return;
 			}
-		}
 
+		} else {
+
+			/* New reflection is identical to a previous one */
+			do {
+				if ( refl->next == NULL ) {
+					refl->next = new;
+					new->prev = refl;
+					return;
+				}
+				refl = refl->next;
+			} while ( 1 );
+
+		}
 
 	}
 }
@@ -294,6 +315,33 @@ void delete_refl(Reflection *refl)
 	int i;
 	Reflection **parent_pos = NULL;
 
+	/* Is this a duplicate, and not the first one? */
+	if ( refl->prev != NULL ) {
+		refl->prev->next = refl->next;
+		free(refl);
+		return;
+	}
+
+	/* Is this the first duplicate of many? */
+	if ( refl->next != NULL ) {
+
+		refl->next->parent = refl->parent;
+		refl->next->child[0] = refl->child[0];
+		refl->next->child[1] = refl->child[1];
+		refl->next->prev = NULL;
+
+		for ( i=0; i<2; i++ ) {
+			if ( refl->parent->child[i] == refl ) {
+				refl->parent->child[i] = refl->next;
+			}
+		}
+
+		free(refl);
+
+		return;
+
+	}
+
 	/* Remove parent's reference */
 	for ( i=0; i<2; i++ ) {
 		if ( refl->parent->child[i] == refl ) {
@@ -344,24 +392,80 @@ void delete_refl(Reflection *refl)
 
 /********************************* Iteration **********************************/
 
-Reflection *first_refl(RefList *list)
+
+struct _reflistiterator {
+
+	int stack_size;
+	int stack_ptr;
+	Reflection **stack;
+
+};
+
+Reflection *first_refl(RefList *list, RefListIterator **piter)
 {
-	return list->head->child[0];
+	RefListIterator *iter;
+
+	iter = malloc(sizeof(struct _reflistiterator));
+	iter->stack_size = 32;
+	iter->stack = malloc(iter->stack_size*sizeof(Reflection *));
+	iter->stack_ptr = 0;
+	*piter = iter;
+
+	Reflection *refl = list->head->child[0];
+
+	do {
+
+		if ( refl != NULL ) {
+			iter->stack[iter->stack_ptr++] = refl;
+			if ( iter->stack_ptr == iter->stack_size ) {
+				iter->stack_size += 32;
+				iter->stack = realloc(iter->stack,
+				         iter->stack_size*sizeof(Reflection *));
+			}
+			refl = refl->child[0];
+			continue;
+		}
+
+		if ( iter->stack_ptr == 0 ) {
+			free(iter->stack);
+			free(iter);
+			return NULL;
+		}
+
+		refl = iter->stack[iter->stack_ptr--];
+
+		return refl;
+
+	} while ( 1 );
 }
 
 
-Reflection *next_refl(Reflection *refl)
+Reflection *next_refl(Reflection *refl, RefListIterator *iter)
 {
-	/* Does a left child exist? */
-	if ( refl->child[0] != NULL ) return refl->child[0];
+	do {
 
-	/* Otherwise move up the tree to find the next right child */
-	while ( refl->child[1] != NULL ) {
-		refl = refl->parent;
-		if ( refl == NULL ) return NULL;
-	}
+		refl = refl->child[1];;
+
+		if ( refl != NULL ) {
+			iter->stack[iter->stack_ptr++] = refl;
+			if ( iter->stack_ptr == iter->stack_size ) {
+				iter->stack_size += 32;
+				iter->stack = realloc(iter->stack,
+				         iter->stack_size*sizeof(Reflection *));
+			}
+			refl = refl->child[0];
+			continue;
+		}
+		if ( iter->stack_ptr == 0 ) {
+			free(iter->stack);
+			free(iter);
+			return NULL;
+		}
+		refl = iter->stack[iter->stack_ptr--];
+
+		return refl;
 
-	return refl->child[1];
+	} while ( 1 );
 }
 
 
diff --git a/src/reflist.h b/src/reflist.h
index 78aadf2a..e25a16ad 100644
--- a/src/reflist.h
+++ b/src/reflist.h
@@ -19,6 +19,7 @@
 
 typedef struct _reflist RefList;
 typedef struct _reflection Reflection;
+typedef struct _reflistiterator RefListIterator;
 
 #define INDICES signed int h, signed int k, signed int l
 
@@ -56,8 +57,8 @@ extern Reflection *add_refl(RefList *list, INDICES);
 extern void delete_refl(Reflection *refl);
 
 /* Iteration */
-extern Reflection *first_refl(RefList *list);
-extern Reflection *next_refl(Reflection *refl);
+extern Reflection *first_refl(RefList *list, RefListIterator **iterator);
+extern Reflection *next_refl(Reflection *refl, RefListIterator *iter);
 
 /* Voodoo */
 extern void optimise_reflist(RefList *list);
diff --git a/src/templates.c b/src/templates.c
index 9abbcb17..79e60aec 100644
--- a/src/templates.c
+++ b/src/templates.c
@@ -232,13 +232,14 @@ static double integrate_all_rot(struct image *image, RefList *reflections,
 	double cosr, sinr;
 	int num_int = 0;
 	Reflection *refl;
+	RefListIterator *iter;
 
 	cosr = cos(rot);
 	sinr = sin(rot);
 
-	for ( refl = first_refl(reflections);
+	for ( refl = first_refl(reflections, &iter);
 	      refl != NULL;
-	      refl = next_refl(refl) ) {
+	      refl = next_refl(refl, iter) ) {
 
 		float xp, yp;
 		double x, y;
@@ -266,14 +267,15 @@ static double mean_distance(struct image *image, RefList *reflections,
 	double cosr, sinr;
 	int num_dist = 0;
 	Reflection *refl;
+	RefListIterator *iter;
 
 	cosr = cos(rot);
 	sinr = sin(rot);
 
 	/* For each template peak */
-	for ( refl = first_refl(reflections);
+	for ( refl = first_refl(reflections, &iter);
 	      refl != NULL;
-	      refl = next_refl(refl) ) {
+	      refl = next_refl(refl, iter) ) {
 
 		float xp, yp;
 		int j;
diff --git a/tests/list_check.c b/tests/list_check.c
index 00e27e00..ffae8d12 100644
--- a/tests/list_check.c
+++ b/tests/list_check.c
@@ -27,6 +27,7 @@ struct refltemp {
 	signed int l;
 	int del;
 	int dup;
+	int found;
 };
 
 #define RANDOM_INDEX (128*random()/RAND_MAX - 256*random()/RAND_MAX)
@@ -38,6 +39,8 @@ static int test_lists(int num_items)
 	RefList *list;
 	int i;
 	signed int h, k, l;
+	Reflection *refl;
+	RefListIterator *iter;
 
 	check = malloc(num_items * sizeof(struct refltemp));
 	list = reflist_new();
@@ -53,11 +56,11 @@ static int test_lists(int num_items)
 		int j;
 		int duplicate = 0;
 
-		if ( random() > RAND_MAX/2 ) {
+		//if ( random() > RAND_MAX/2 ) {
 			h = RANDOM_INDEX;
 			k = RANDOM_INDEX;
 			l = RANDOM_INDEX;
-		} /* else use the same as last time */
+		//} /* else use the same as last time */
 
 		/* Count the number of times this reflection appeared before */
 		for ( j=0; j<i; j++ ) {
@@ -83,9 +86,49 @@ static int test_lists(int num_items)
 		check[i].l = l;
 		check[i].del = 0;
 		check[i].dup = duplicate;
+		check[i].found = 0;
 
 	}
 
+	/* Iterate over the list and check we find everything */
+	for ( refl = first_refl(list, &iter);
+	      refl != NULL;
+	      refl = next_refl(refl, iter) ) {
+
+		signed int h, k, l;
+
+		get_indices(refl, &h, &k, &l);
+		printf("%3i %3i %3i\n", h, k, l);
+
+		for ( i=0; i<num_items; i++ ) {
+			if ( (check[i].h == h)
+			  && (check[i].k == k)
+			  && (check[i].l == l)
+			  && (check[i].found == 0) ) {
+				check[i].found = 1;
+			}
+		}
+
+	}
+	for ( i=0; i<num_items; i++ ) {
+		if ( check[i].found == 0 ) {
+
+			Reflection *test;
+
+			fprintf(stderr, "Iteration didn't find %3i %3i %3i %i\n",
+			        check[i].h, check[i].k, check[i].l, i);
+			test = find_refl(list, check[i].h, check[i].k,
+			                 check[i].l);
+			if ( test == NULL ) {
+				fprintf(stderr, "Not in list\n");
+			} else {
+				fprintf(stderr, "But found in list.\n");
+			}
+			return 1;
+
+		}
+	}
+
 	/* Check that all the reflections can be found,
 	 * and delete the first few. */
 	for ( i=0; i<num_items; i++ ) {
@@ -187,7 +230,7 @@ int main(int argc, char *argv[])
 	printf("Running list test...");
 	fflush(stdout);
 
-	for ( i=0; i<100; i++ ) {
+	for ( i=0; i<1; i++ ) {
 		if ( test_lists(4096*random()/RAND_MAX) ) return 1;
 	}