process_hkl: Use a different (array) indexing method to speed it up lots

src/likelihood.c

@@ -19,17 +19,20 @@
 
 void detwin_intensities(const double *model, double *new_pattern,
                         const unsigned int *model_counts,
-                        unsigned int *new_counts)
+                        ReflItemList *items)
 {
 	/* Placeholder... */
 }
 
 void scale_intensities(const double *model, double *new_pattern,
                        const unsigned int *model_counts,
-                       unsigned int *new_counts, double f0, int f0_valid)
+                       ReflItemList *items, double f0, int f0_valid)
 {
 	double s;
 	unsigned int i;
+	unsigned int *new_counts;
+
+	new_counts = items_to_counts(items);
 
 	s = stat_scale_intensity(model, model_counts, new_pattern, new_counts);
 	if ( f0_valid ) printf("%f %f\n", s, f0);
@@ -41,4 +44,6 @@ void scale_intensities(const double *model, double *new_pattern,
 	for ( i=0; i<LIST_SIZE; i++ ) {
 		new_counts[i] *= s;
 	}
+
+	free(new_counts);
 }

src/likelihood.h

@@ -17,13 +17,16 @@
 #define LIKELIHOOD_H
 
 
+#include "utils.h"
+
+
 extern void detwin_intensities(const double *model, double *new_pattern,
 	                       const unsigned int *model_counts,
-                               unsigned int *new_counts);
+                               ReflItemList *items);
 
 extern void scale_intensities(const double *model, double *new_pattern,
 	                      const unsigned int *model_counts,
-                              unsigned int *new_counts, double f0,
+                              ReflItemList *items, double f0,
                               int f0_valid);
 
 

src/process_hkl.c

@@ -188,17 +188,21 @@ static void process_reflections(double *ref, unsigned int *counts,
 
 static void merge_pattern(double *model, const double *new,
                           unsigned int *model_counts,
-                          const unsigned int *counts, int mo, int sum)
+                          ReflItemList *items, int mo, int sum)
 {
-	signed int h, k, l;
+	int i;
 
-	for ( l=-INDMAX; l<INDMAX; l++ ) {
-	for ( k=-INDMAX; k<INDMAX; k++ ) {
-	for ( h=-INDMAX; h<INDMAX; h++ ) {
+	for ( i=0; i<num_items(items); i++ ) {
 
 		double intensity;
+		signed int h, k, l;
+		struct refl_item *item;
 
-		if ( lookup_count(counts, h, k, l) == 0 ) continue;
+		item = get_item(items, i);
+
+		h = item->h;
+		k = item->k;
+		l = item->l;
 
 		intensity = lookup_intensity(new, h, k, l);
 
@@ -217,8 +221,6 @@ static void merge_pattern(double *model, const double *new,
 		}
 
 	}
-	}
-	}
 }
 
 
@@ -243,7 +245,6 @@ int main(int argc, char *argv[])
 	int config_scale = 0;
 	char *intfile = NULL;
 	double *new_pattern = NULL;
-	unsigned int *new_counts = NULL;
 	unsigned int n_total_patterns;
 	unsigned int *truecounts = NULL;
 	char *sym = NULL;
@@ -251,6 +252,7 @@ int main(int argc, char *argv[])
 	float f0;
 	int f0_valid;
 	int n_nof0 = 0;
+	ReflItemList *items;
 
 	/* Long options */
 	const struct option longopts[] = {
@@ -344,7 +346,7 @@ int main(int argc, char *argv[])
 	model_counts = new_list_count();
 	cell = load_cell_from_pdb(pdb);
 	new_pattern = new_list_intensity();
-	new_counts = new_list_count();
+	items = new_items();
 
 	if ( strcmp(filename, "-") == 0 ) {
 		fh = stdin;
@@ -393,7 +395,7 @@ int main(int argc, char *argv[])
 			/* Detwin before scaling */
 			if ( config_detwin ) {
 				detwin_intensities(model, new_pattern,
-				                   model_counts, new_counts);
+				                   model_counts, items);
 			}
 
 			/* Assume a default I0 if we don't have one by now */
@@ -405,13 +407,13 @@ int main(int argc, char *argv[])
 			/* Scale if requested */
 			if ( config_scale ) {
 				scale_intensities(model, new_pattern,
-				                  model_counts, new_counts, f0,
+				                  model_counts, items, f0,
 				                  f0_valid);
 			}
 
 			/* Start of second or later pattern */
 			merge_pattern(model, new_pattern, model_counts,
-			              new_counts, config_maxonly, config_sum);
+			              items, config_maxonly, config_sum);
 
 			if ( (trueref != NULL) && config_every
 			    && (n_patterns % config_every == 0) ) {
@@ -424,9 +426,8 @@ int main(int argc, char *argv[])
 
 			if ( n_patterns == config_stopafter ) break;
 
-			zero_list_count(new_counts);
-
 			n_patterns++;
+			clear_items(items);
 
 			progress_bar(n_patterns, n_total_patterns, "Merging");
 
@@ -449,12 +450,12 @@ int main(int argc, char *argv[])
 
 		if ( (h==0) && (k==0) && (l==0) ) continue;
 
-		if ( lookup_count(new_counts, h, k, l) != 0 ) {
+		if ( find_item(items, h, k, l) != 0 ) {
 			ERROR("More than one measurement for %i %i %i in"
 			      " pattern number %i\n", h, k, l, n_patterns);
 		}
 		set_intensity(new_pattern, h, k, l, intensity);
-		set_count(new_counts, h, k, l, 1);
+		add_item(items, h, k, l);
 
 	} while ( rval != NULL );
 
@@ -480,5 +481,7 @@ int main(int argc, char *argv[])
 	STATUS("There were %u patterns.\n", n_patterns);
 	STATUS("%i had no f0 valid value.\n", n_nof0);
 
+	delete_items(items);
+
 	return 0;
 }

src/utils.c

@@ -283,3 +283,93 @@ int assplode(const char *a, const char *delims, char ***pbits,
 	*pbits = bits;
 	return n;
 }
+
+
+struct _reflitemlist {
+	struct refl_item *items;
+	int n_items;
+	int max_items;
+};
+
+void clear_items(ReflItemList *items)
+{
+	items->n_items = 0;
+}
+
+static void alloc_items(ReflItemList *items)
+{
+	items->items = realloc(items->items,
+	                       items->max_items*sizeof(struct refl_item));
+}
+
+ReflItemList *new_items()
+{
+	ReflItemList *new;
+	new = malloc(sizeof(ReflItemList));
+	new->max_items = 1024;
+	new->n_items = 0;
+	new->items = NULL;
+	alloc_items(new);
+	return new;
+}
+
+void delete_items(ReflItemList *items)
+{
+	if ( items->items != NULL ) free(items->items);
+	free(items);
+}
+
+void add_item(ReflItemList *items,
+                     signed int h, signed int k, signed int l)
+{
+	if ( items->n_items == items->max_items ) {
+		items->max_items += 1024;
+		alloc_items(items);
+	}
+
+	items->items[items->n_items].h = h;
+	items->items[items->n_items].k = k;
+	items->items[items->n_items].l = l;
+	items->n_items++;
+}
+
+int find_item(ReflItemList *items,
+                     signed int h, signed int k, signed int l)
+{
+	int i;
+
+	for ( i=0; i<items->n_items; i++ ) {
+		if ( items->items[i].h != h ) continue;
+		if ( items->items[i].k != k ) continue;
+		if ( items->items[i].l != l ) continue;
+		return 1;
+	}
+	return 0;
+}
+
+struct refl_item *get_item(ReflItemList *items, int i)
+{
+	if ( i >= items->n_items ) return NULL;
+	return &items->items[i];
+}
+
+int num_items(ReflItemList *items)
+{
+	return items->n_items;
+}
+
+unsigned int *items_to_counts(ReflItemList *items)
+{
+	unsigned int *c;
+	int i;
+
+	c = new_list_count();
+
+	for ( i=0; i<num_items(items); i++ ) {
+		struct refl_item *r;
+		r = get_item(items, i);
+		set_count(c, r->h, r->k, r->l, 1);
+	}
+
+	return c;
+}

src/utils.h

@@ -176,6 +176,28 @@ static inline double angle_between(double x1, double y1, double z1,
 #include "list_tmp.h"
 
 
+/* ----------- Reflection lists indexed by sequence (not indices) ----------- */
+
+typedef struct _reflitemlist ReflItemList;  /* Opaque */
+
+struct refl_item {
+	signed int h;
+	signed int k;
+	signed int l;
+};
+
+extern void clear_items(ReflItemList *items);
+extern ReflItemList *new_items(void);
+extern void delete_items(ReflItemList *items);
+extern void add_item(ReflItemList *items,
+                     signed int h, signed int k, signed int l);
+extern int find_item(ReflItemList *items,
+                     signed int h, signed int k, signed int l);
+extern struct refl_item *get_item(ReflItemList *items, int i);
+extern int num_items(ReflItemList *items);
+extern unsigned int *items_to_counts(ReflItemList *items);
+
+
 /* ------------------------------ Message macros ---------------------------- */
 
 #define ERROR(...) fprintf(stderr, __VA_ARGS__)