Implement "multi-indexing"

3 files changed, 196 insertions, 113 deletions

diff --git a/src/index.c b/src/index.c
index edfa33e5..4d2a158f 100644
--- a/src/index.c
+++ b/src/index.c
@@ -43,39 +43,67 @@ static IndexingPrivate *indexing_private(IndexingMethod indm)
 }
 
 
-IndexingPrivate *prepare_indexing(IndexingMethod indm, UnitCell *cell,
+IndexingPrivate **prepare_indexing(IndexingMethod *indm, UnitCell *cell,
                                   const char *filename, struct detector *det,
                                   double nominal_photon_energy)
 {
-	switch ( indm ) {
-	case INDEXING_NONE :
-		return indexing_private(indm);
-	case INDEXING_DIRAX :
-		return indexing_private(indm);
-	case INDEXING_MOSFLM :
-		return indexing_private(indm);
-	case INDEXING_TEMPLATE :
-		return generate_templates(cell, filename, det,
-		                          nominal_photon_energy);
+	int n;
+	int nm = 0;
+	IndexingPrivate **iprivs;
+
+	while ( indm[nm] != INDEXING_NONE ) nm++;
+	STATUS("Preparing %i indexing methods.\n", nm);
+	iprivs = malloc((nm+1) * sizeof(IndexingPrivate *));
+
+	for ( n=0; n<nm; n++ ) {
+
+		switch ( indm[n] ) {
+		case INDEXING_NONE :
+			ERROR("Tried to prepare INDEXING_NONE!\n");
+			break;
+		case INDEXING_DIRAX :
+			 iprivs[n] = indexing_private(indm[n]);
+			 break;
+		case INDEXING_MOSFLM :
+			 iprivs[n] = indexing_private(indm[n]);
+			 break;
+		case INDEXING_TEMPLATE :
+			iprivs[n] = generate_templates(cell, filename, det,
+				                  nominal_photon_energy);
+			break;
+		}
+
+		n++;
+
 	}
-	return 0;
+	iprivs[n] = NULL;
+
+	return iprivs;
 }
 
 
-void cleanup_indexing(IndexingPrivate *priv)
+void cleanup_indexing(IndexingPrivate **priv)
 {
-	switch ( priv->indm ) {
-	case INDEXING_NONE :
-		free(priv);
-		break;
-	case INDEXING_DIRAX :
-		free(priv);
-		break;
-	case INDEXING_MOSFLM :
-		free(priv);
-		break;
-	case INDEXING_TEMPLATE :
-		free_templates(priv);
+	int n = 0;
+
+	while ( priv[n] != NULL ) {
+
+		switch ( priv[n]->indm ) {
+		case INDEXING_NONE :
+			free(priv[n]);
+			break;
+		case INDEXING_DIRAX :
+			free(priv[n]);
+			break;
+		case INDEXING_MOSFLM :
+			free(priv[n]);
+			break;
+		case INDEXING_TEMPLATE :
+			free_templates(priv[n]);
+		}
+
+		n++;
+
 	}
 }
 
@@ -246,81 +274,93 @@ void map_all_peaks(struct image *image)
 }
 
 
-void index_pattern(struct image *image, UnitCell *cell, IndexingMethod indm,
-                   int cellr, int verbose, IndexingPrivate *ipriv)
+void index_pattern(struct image *image, UnitCell *cell, IndexingMethod *indm,
+                   int cellr, int verbose, IndexingPrivate **ipriv)
 {
 	int i;
+	int n = 0;
 
 	map_all_peaks(image);
-	write_drx(image);
 
-	image->ncells = 0;
+	while ( indm[n] != INDEXING_NONE ) {
 
-	/* Index (or not) as appropriate */
-	switch ( indm ) {
-	case INDEXING_NONE :
-		return;
-	case INDEXING_DIRAX :
-		run_dirax(image);
-		break;
-	case INDEXING_MOSFLM :
-		write_spt(image);
-		write_img(image); /* dummy image. not needed for each frame.*/
-		run_mosflm(image,cell);
-		break;
-	case INDEXING_TEMPLATE :
-		match_templates(image, ipriv);
-		break;
-	}
-
-	if ( image->ncells == 0 ) {
-		STATUS("No candidate cells found.\n");
-		return;
-	}
+		image->ncells = 0;
 
-	if ( (cellr == CELLR_NONE) || (indm == INDEXING_TEMPLATE) ) {
-		image->indexed_cell = image->candidate_cells[0];
-		if ( verbose ) {
-			STATUS("--------------------\n");
-			STATUS("The indexed cell (matching not performed):\n");
-			cell_print(image->indexed_cell);
-			STATUS("--------------------\n");
+		/* Index as appropriate */
+		switch ( indm[n] ) {
+		case INDEXING_NONE :
+			return;
+		case INDEXING_DIRAX :
+			STATUS("Running DirAx...\n");
+			write_drx(image);
+			run_dirax(image);
+			break;
+		case INDEXING_MOSFLM :
+			STATUS("Running MOSFLM...\n");
+			write_spt(image);
+			write_img(image); /* dummy image */
+			run_mosflm(image, cell);
+			break;
+		case INDEXING_TEMPLATE :
+			match_templates(image, ipriv[n]);
+			break;
 		}
-		return;
-	}
 
-	for ( i=0; i<image->ncells; i++ ) {
-
-		UnitCell *new_cell = NULL;
-
-		if ( verbose ) {
-			STATUS("--------------------\n");
-			STATUS("Candidate cell %i (before matching):\n", i);
-			cell_print(image->candidate_cells[i]);
-			STATUS("--------------------\n");
+		if ( image->ncells == 0 ) {
+			STATUS("No candidate cells found.\n");
+			return;
 		}
 
-		/* Match or reduce the cell as appropriate */
-		switch ( cellr ) {
-		case CELLR_NONE :
-			/* Never happens */
-			break;
-		case CELLR_REDUCE :
-			new_cell = match_cell(image->candidate_cells[i],
-			                      cell, verbose, 1);
-			break;
-		case CELLR_COMPARE :
-			new_cell = match_cell(image->candidate_cells[i],
-			                      cell, verbose, 0);
-			break;
+		if ( (cellr == CELLR_NONE) || (indm[n] == INDEXING_TEMPLATE) ) {
+			image->indexed_cell = image->candidate_cells[0];
+			if ( verbose ) {
+				STATUS("--------------------\n");
+				STATUS("The indexed cell (matching not"
+				       " performed):\n");
+				cell_print(image->indexed_cell);
+				STATUS("--------------------\n");
+			}
+			return;
 		}
 
-		image->indexed_cell = new_cell;
-		if ( new_cell != NULL ) {
-			STATUS("Matched on attempt %i.\n", i);
-			goto done;
+		for ( i=0; i<image->ncells; i++ ) {
+
+			UnitCell *new_cell = NULL;
+
+			if ( verbose ) {
+				STATUS("--------------------\n");
+				STATUS("Candidate cell %i (before matching):\n",
+				       i);
+				cell_print(image->candidate_cells[i]);
+				STATUS("--------------------\n");
+			}
+
+			/* Match or reduce the cell as appropriate */
+			switch ( cellr ) {
+			case CELLR_NONE :
+				/* Never happens */
+				break;
+			case CELLR_REDUCE :
+				new_cell = match_cell(image->candidate_cells[i],
+					              cell, verbose, 1);
+				break;
+			case CELLR_COMPARE :
+				new_cell = match_cell(image->candidate_cells[i],
+					              cell, verbose, 0);
+				break;
+			}
+
+			image->indexed_cell = new_cell;
+			if ( new_cell != NULL ) {
+				STATUS("Matched on attempt %i.\n", i);
+				goto done;
+			}
+
 		}
 
+		/* Move on to the next indexing method */
+		n++;
+
 	}
 
 done:
@@ -328,3 +368,38 @@ done:
 		cell_free(image->candidate_cells[i]);
 	}
 }
+
+
+IndexingMethod *build_indexer_list(const char *str, int *need_cell)
+{
+	int n, i;
+	char **methods;
+	IndexingMethod *list;
+	*need_cell = 0;
+
+	n = assplode(str, ",", &methods, ASSPLODE_NONE);
+	list = malloc((n+1)*sizeof(IndexingMethod));
+
+	for ( i=0; i<n; i++ ) {
+
+		if ( strcmp(methods[i], "dirax") == 0) {
+			list[i] = INDEXING_DIRAX;
+		} else if ( strcmp(methods[i], "mosflm") == 0) {
+			list[i] = INDEXING_MOSFLM;
+		} else if ( strcmp(methods[i], "template") == 0) {
+			list[i] = INDEXING_TEMPLATE;
+			*need_cell = 1;
+		} else {
+			ERROR("Unrecognised indexing method '%s'\n",
+			      methods[i]);
+			return NULL;
+		}
+
+		free(methods[i]);
+
+	}
+	free(methods);
+	list[i] = INDEXING_NONE;
+
+	return list;
+}
diff --git a/src/index.h b/src/index.h
index 29be9dcc..b58c9b53 100644
--- a/src/index.h
+++ b/src/index.h
@@ -42,7 +42,7 @@ enum {
 
 typedef struct _indexingprivate IndexingPrivate;
 
-extern IndexingPrivate *prepare_indexing(IndexingMethod indm, UnitCell *cell,
+extern IndexingPrivate **prepare_indexing(IndexingMethod *indm, UnitCell *cell,
                                          const char *filename,
                                          struct detector *det,
                                          double nominal_photon_energy);
@@ -50,9 +50,11 @@ extern IndexingPrivate *prepare_indexing(IndexingMethod indm, UnitCell *cell,
 extern void map_all_peaks(struct image *image);
 
 extern void index_pattern(struct image *image, UnitCell *cell,
-                          IndexingMethod indm, int cellr, int verbose,
-                          IndexingPrivate *priv);
+                          IndexingMethod *indm, int cellr, int verbose,
+                          IndexingPrivate **priv);
 
-extern void cleanup_indexing(IndexingPrivate *priv);
+extern void cleanup_indexing(IndexingPrivate **priv);
+
+extern IndexingMethod *build_indexer_list(const char *str, int *need_cell);
 
 #endif	/* INDEX_H */
diff --git a/src/indexamajig.c b/src/indexamajig.c
index 8d762fa4..a2b5e80d 100644
--- a/src/indexamajig.c
+++ b/src/indexamajig.c
@@ -68,8 +68,8 @@ struct static_index_args
 	float threshold;
 	float min_gradient;
 	struct detector *det;
-	IndexingMethod indm;
-	IndexingPrivate *ipriv;
+	IndexingMethod *indm;
+	IndexingPrivate **ipriv;
 	const double *intensities;
 	const unsigned char *flags;
 	const char *sym;  /* Symmetry of "intensities" and "flags" */
@@ -325,7 +325,7 @@ static void process_image(void *pp, int cookie)
 	int config_gpu = pargs->static_args.config_gpu;
 	int config_simulate = pargs->static_args.config_simulate;
 	int config_polar = pargs->static_args.config_polar;
-	IndexingMethod indm = pargs->static_args.indm;
+	IndexingMethod *indm = pargs->static_args.indm;
 	const double *intensities = pargs->static_args.intensities;
 	const unsigned char *flags = pargs->static_args.flags;
 	struct gpu_context *gctx = pargs->static_args.gctx;
@@ -543,7 +543,10 @@ int main(int argc, char *argv[])
 	float min_gradient = 100000.0;
 	struct detector *det;
 	char *geometry = NULL;
-	IndexingMethod indm;
+	IndexingMethod *indm;
+	IndexingPrivate **ipriv;
+	int indexer_needs_cell;
+	int reduction_needs_cell;
 	char *indm_str = NULL;
 	UnitCell *cell;
 	double *intensities = NULL;
@@ -561,7 +564,6 @@ int main(int argc, char *argv[])
 	pthread_mutex_t gpu_mutex = PTHREAD_MUTEX_INITIALIZER;
 	char prepare_line[1024];
 	char prepare_filename[1024];
-	IndexingPrivate *ipriv;
 	struct queue_args qargs;
 	struct beam_params *beam = NULL;
 	double nominal_photon_energy;
@@ -784,31 +786,30 @@ int main(int argc, char *argv[])
 		       " try to index anything.\n"
 		       "If that isn't what you wanted, re-run with"
 		       " --indexing=<method>.\n");
-		indm = INDEXING_NONE;
-	} else if ( strcmp(indm_str, "none") == 0 ) {
-		indm = INDEXING_NONE;
-	} else if ( strcmp(indm_str, "dirax") == 0) {
-		indm = INDEXING_DIRAX;
-	} else if ( strcmp(indm_str, "mosflm") == 0) {
-		indm = INDEXING_MOSFLM;
-	} else if ( strcmp(indm_str, "template") == 0) {
-		indm = INDEXING_TEMPLATE;
+		indm = NULL;
 	} else {
-		ERROR("Unrecognised indexing method '%s'\n", indm_str);
-		return 1;
+		indm = build_indexer_list(indm_str, &indexer_needs_cell);
+		if ( indm == NULL ) {
+			ERROR("Invalid indexer list '%s'\n", indm_str);
+			return 1;
+		}
+		free(indm_str);
 	}
-	free(indm_str);
 
+	reduction_needs_cell = 0;
 	if ( scellr == NULL ) {
 		STATUS("You didn't specify a cell reduction method, so I'm"
 		       " going to use 'reduce'.\n");
 		cellr = CELLR_REDUCE;
+		reduction_needs_cell = 1;
 	} else if ( strcmp(scellr, "none") == 0 ) {
 		cellr = CELLR_NONE;
 	} else if ( strcmp(scellr, "reduce") == 0) {
 		cellr = CELLR_REDUCE;
+		reduction_needs_cell = 1;
 	} else if ( strcmp(scellr, "compare") == 0) {
 		cellr = CELLR_COMPARE;
+		reduction_needs_cell = 1;
 	} else {
 		ERROR("Unrecognised cell reduction method '%s'\n", scellr);
 		return 1;
@@ -827,14 +828,15 @@ int main(int argc, char *argv[])
 	}
 	free(geometry);
 
-	if ( (cellr != CELLR_NONE) || (indm == INDEXING_TEMPLATE) ) {
+	if ( reduction_needs_cell || indexer_needs_cell ) {
 		cell = load_cell_from_pdb(pdb);
 		if ( cell == NULL ) {
 			ERROR("Couldn't read unit cell (from %s)\n", pdb);
 			return 1;
 		}
 	} else {
-		STATUS("No cell needed because --no-match was used.\n");
+		STATUS("No cell needed for these choices of indexing"
+		       " and reduction.\n");
 		cell = NULL;
 	}
 	free(pdb);
@@ -866,11 +868,15 @@ int main(int argc, char *argv[])
 	qargs.use_this_one_instead = prepare_line;
 
 	/* Prepare the indexer */
-	ipriv = prepare_indexing(indm, cell, prepare_filename, det,
-	                         nominal_photon_energy);
-	if ( ipriv == NULL ) {
-		ERROR("Failed to prepare indexing.\n");
-		return 1;
+	if ( indm != NULL ) {
+		ipriv = prepare_indexing(indm, cell, prepare_filename, det,
+		                         nominal_photon_energy);
+		if ( ipriv == NULL ) {
+			ERROR("Failed to prepare indexing.\n");
+			return 1;
+		}
+	} else {
+		ipriv = NULL;
 	}
 
 	gsl_set_error_handler_off();