Merge branch 'tom/xgandalf'

5 files changed, 400 insertions, 0 deletions

diff --git a/libcrystfel/CMakeLists.txt b/libcrystfel/CMakeLists.txt
index f6f25ec9..25cf66bc 100644
--- a/libcrystfel/CMakeLists.txt
+++ b/libcrystfel/CMakeLists.txt
@@ -68,6 +68,7 @@ set(LIBCRYSTFEL_SOURCES
     src/felix.c
     src/peakfinder8.c
     src/taketwo.c
+    src/xgandalf.c
 )
 
 if (HAVE_FFTW)
@@ -106,6 +107,7 @@ set(LIBCRYSTFEL_HEADERS
     src/felix.h
     src/peakfinder8.h
     src/taketwo.h
+    src/xgandalf.h
 )
 
 add_library(${PROJECT_NAME} SHARED
diff --git a/libcrystfel/src/index.c b/libcrystfel/src/index.c
index ff84c629..4ff4405a 100644
--- a/libcrystfel/src/index.c
+++ b/libcrystfel/src/index.c
@@ -57,6 +57,7 @@
 #include "felix.h"
 #include "predict-refine.h"
 #include "taketwo.h"
+#include "xgandalf.h"
 
 
 struct _indexingprivate
@@ -66,6 +67,7 @@ struct _indexingprivate
 	float tolerance[4];
 
 	struct taketwo_options *ttopts;
+	struct xgandalf_options *xgandalf_opts;
 
 	int n_methods;
 	IndexingMethod *methods;
@@ -182,6 +184,10 @@ static char *base_indexer_str(IndexingMethod indm)
 		strcpy(str, "taketwo");
 		break;
 
+		case INDEXING_XGANDALF:
+		strcpy(str, "xgandalf");
+		break;
+
 		case INDEXING_SIMULATION :
 		strcpy(str, "simulation");
 		break;
@@ -219,6 +225,7 @@ static char *friendly_indexer_name(IndexingMethod m)
 
 
 static void *prepare_method(IndexingMethod *m, UnitCell *cell,
+                            struct xgandalf_options *xgandalf_opts,
                             struct felix_options *felix_opts)
 {
 	char *str;
@@ -259,6 +266,10 @@ static void *prepare_method(IndexingMethod *m, UnitCell *cell,
 		priv = taketwo_prepare(m, cell);
 		break;
 
+		case INDEXING_XGANDALF :
+		priv = xgandalf_prepare(m, cell, xgandalf_opts);
+		break;
+
 		default :
 		ERROR("Don't know how to prepare indexing method %i\n", *m);
 		break;
@@ -290,6 +301,7 @@ IndexingPrivate *setup_indexing(const char *method_list, UnitCell *cell,
                                 struct detector *det, float *ltl,
                                 IndexingFlags flags,
                                 struct taketwo_options *ttopts,
+                                struct xgandalf_options *xgandalf_opts,
                                 struct felix_options *felix_opts)
 {
 	int i, n;
@@ -386,6 +398,7 @@ IndexingPrivate *setup_indexing(const char *method_list, UnitCell *cell,
 		int j;
 
 		ipriv->engine_private[i] = prepare_method(&methods[i], cell,
+		                                          xgandalf_opts,
 		                                          felix_opts);
 
 		if ( ipriv->engine_private[i] == NULL ) return NULL;
@@ -411,6 +424,7 @@ IndexingPrivate *setup_indexing(const char *method_list, UnitCell *cell,
 	for ( i=0; i<4; i++ ) ipriv->tolerance[i] = ltl[i];
 
 	ipriv->ttopts = ttopts;
+	ipriv->xgandalf_opts = xgandalf_opts;
 
 	STATUS("List of indexing methods:\n");
 	for ( i=0; i<n; i++ ) {
@@ -467,6 +481,10 @@ void cleanup_indexing(IndexingPrivate *ipriv)
 			taketwo_cleanup(ipriv->engine_private[n]);
 			break;
 
+			case INDEXING_XGANDALF :
+			xgandalf_cleanup(ipriv->engine_private[n]);
+			break;
+
 			default :
 			ERROR("Don't know how to clean up indexing method %i\n",
 			      ipriv->methods[n]);
@@ -577,6 +595,10 @@ static int try_indexer(struct image *image, IndexingMethod indm,
 		r = taketwo_index(image, ipriv->ttopts, mpriv);
 		break;
 
+		case INDEXING_XGANDALF :
+		r = run_xgandalf(image, mpriv);
+		break;
+
 		default :
 		ERROR("Unrecognised indexing method: %i\n", indm);
 		return 0;
@@ -945,6 +967,11 @@ IndexingMethod get_indm_from_string_2(const char *str, int *err)
 			method = INDEXING_DEFAULTS_TAKETWO;
 			have_method = 1;
 
+		} else if ( strcmp(bits[i], "xgandalf") == 0) {
+			if ( have_method ) return warn_method(str);
+			method = INDEXING_DEFAULTS_XGANDALF;
+			have_method = 1;
+
 		} else if ( strcmp(bits[i], "none") == 0) {
 			if ( have_method ) return warn_method(str);
 			method = INDEXING_NONE;
@@ -1042,6 +1069,7 @@ char *detect_indexing_methods(UnitCell *cell)
 	do_probe(dirax_probe, cell, methods);
 	do_probe(asdf_probe, cell, methods);
 	do_probe(xds_probe, cell, methods);
+	do_probe(xgandalf_probe, cell, methods);
 	/* Don't automatically use TakeTwo or Felix (yet) */
 	//do_probe(taketwo_probe, cell, methods);
 	//do_probe(felix_probe, cell, methods);
diff --git a/libcrystfel/src/index.h b/libcrystfel/src/index.h
index 8ece9227..2099b4d7 100644
--- a/libcrystfel/src/index.h
+++ b/libcrystfel/src/index.h
@@ -56,6 +56,8 @@
 #define INDEXING_DEFAULTS_XDS (INDEXING_XDS | INDEXING_USE_LATTICE_TYPE \
                                      | INDEXING_USE_CELL_PARAMETERS)
 
+#define INDEXING_DEFAULTS_XGANDALF (INDEXING_XGANDALF | INDEXING_USE_CELL_PARAMETERS)
+
 /**
  * IndexingMethod:
  * @INDEXING_NONE: No indexing to be performed
@@ -67,6 +69,7 @@
  * @INDEXING_DEBUG: Results injector for debugging
  * @INDEXING_ASDF: Use in-built "asdf" indexer
  * @INDEXING_TAKETWO: Use in-built "taketwo" indexer
+ * @INDEXING_XGANDALF: Invoke XGANDALF
  * @INDEXING_ERROR: Special value for unrecognised indexing engine name
  * @INDEXING_USE_LATTICE_TYPE: Use lattice type and centering information to
  *   guide the indexing process.
@@ -90,6 +93,7 @@ typedef enum {
 	INDEXING_DEBUG = 7,
 	INDEXING_ASDF = 8,
 	INDEXING_TAKETWO = 9,
+	INDEXING_XGANDALF = 10,
 
 	INDEXING_ERROR = 255,  /* Unrecognised indexing engine */
 
@@ -136,6 +140,7 @@ extern IndexingMethod get_indm_from_string_2(const char *method, int *err);
 #include "cell.h"
 #include "image.h"
 #include "taketwo.h"
+#include "xgandalf.h"
 #include "felix.h"
 
 
@@ -143,6 +148,7 @@ extern IndexingPrivate *setup_indexing(const char *methods, UnitCell *cell,
                                        struct detector *det, float *ltl,
                                        IndexingFlags flags,
                                        struct taketwo_options *ttopts,
+                                       struct xgandalf_options *xgandalf_opts,
                                        struct felix_options *felix_opts);
 
 extern char *detect_indexing_methods(UnitCell *cell);
diff --git a/libcrystfel/src/xgandalf.c b/libcrystfel/src/xgandalf.c
new file mode 100644
index 00000000..9a2e54ec
--- /dev/null
+++ b/libcrystfel/src/xgandalf.c
@@ -0,0 +1,306 @@
+/*
+ * xgandalf.c
+ *
+ * Interface to XGANDALF indexer
+ *
+ * Copyright © 2017-2018 Deutsches Elektronen-Synchrotron DESY,
+ *                       a research centre of the Helmholtz Association.
+ *
+ * Authors:
+ *   2017-2018 Yaroslav Gevorkov <yaroslav.gevorkov@desy.de>
+ *
+ * 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/>.
+ *
+ */
+
+#include "xgandalf.h"
+
+#ifdef HAVE_XGANDALF
+#include <stdlib.h>
+
+#include "utils.h"
+#include "cell-utils.h"
+#include "peaks.h"
+
+#include "xgandalf/adaptions/crystfel/Lattice.h"
+#include "xgandalf/adaptions/crystfel/ExperimentSettings.h"
+#include "xgandalf/adaptions/crystfel/IndexerPlain.h"
+
+struct xgandalf_private_data {
+	IndexerPlain *indexer;
+	reciprocalPeaks_1_per_A_t reciprocalPeaks_1_per_A;
+
+	IndexingMethod indm;
+	UnitCell *cellTemplate;
+	Lattice_t sampleRealLatticeReduced_A;   //same as cellTemplate
+};
+
+#define FAKE_DETECTOR_DISTANCE (0.1)
+#define FAKE_DETECTOR_RADIUS (0.1)
+#define FAKE_BEAM_ENERGY (1)
+#define FAKE_DIVERGENCE_ANGLE_DEG (0.05)
+#define FAKE_NON_MONOCHROMATICITY (0.005)
+#define FAKE_REFLECTION_RADIUS (0.0001)
+
+#define MAX_ASSEMBLED_LATTICES_COUNT (10)
+
+static void reduceCell(UnitCell* cell);
+static void makeRightHanded(UnitCell* cell);
+
+
+int run_xgandalf(struct image *image, void *ipriv)
+{
+	struct xgandalf_private_data *xgandalf_private_data =
+	                                      (struct xgandalf_private_data*) ipriv;
+	reciprocalPeaks_1_per_A_t *reciprocalPeaks_1_per_A =
+	                          &(xgandalf_private_data->reciprocalPeaks_1_per_A);
+
+	int peakCountMax = image_feature_count(image->features);
+	reciprocalPeaks_1_per_A->peakCount = 0;
+	for (int i = 0; i < peakCountMax && i < MAX_PEAK_COUNT_FOR_INDEXER; i++) {
+		struct imagefeature *f;
+		f = image_get_feature(image->features, i);
+		if (f == NULL) {
+			continue;
+		}
+
+		reciprocalPeaks_1_per_A->coordinates_x[reciprocalPeaks_1_per_A->peakCount]
+		                                                        = f->rx * 1e-10;
+		reciprocalPeaks_1_per_A->coordinates_y[reciprocalPeaks_1_per_A->peakCount]
+		                                                        = f->ry * 1e-10;
+		reciprocalPeaks_1_per_A->coordinates_z[reciprocalPeaks_1_per_A->peakCount]
+		                                                        = f->rz * 1e-10;
+		reciprocalPeaks_1_per_A->peakCount++;
+	}
+
+	Lattice_t assembledLattices[MAX_ASSEMBLED_LATTICES_COUNT];
+	int assembledLatticesCount;
+	IndexerPlain_index(xgandalf_private_data->indexer,
+	                   assembledLattices,
+	                   &assembledLatticesCount,
+	                   MAX_ASSEMBLED_LATTICES_COUNT,
+	                   *reciprocalPeaks_1_per_A,
+	                   NULL);
+
+	if (assembledLatticesCount > 0) { //no multi-lattice at the moment
+		assembledLatticesCount = 1;
+	}
+
+	int goodLatticesCount = assembledLatticesCount;
+	for (int i = 0; i < assembledLatticesCount && i < 1; i++) {
+		reorderLattice(&(xgandalf_private_data->sampleRealLatticeReduced_A),
+		               &assembledLattices[i]);
+
+		UnitCell *uc;
+		uc = cell_new();
+
+		Lattice_t *l = &assembledLattices[i];
+
+		cell_set_cartesian(uc, l->ax * 1e-10, l->ay * 1e-10, l->az * 1e-10,
+		                       l->bx * 1e-10, l->by * 1e-10, l->bz * 1e-10,
+		                       l->cx * 1e-10, l->cy * 1e-10, l->cz * 1e-10);
+		makeRightHanded(uc);
+
+		if (validate_cell(uc)) {
+			STATUS("Problem with returned cell!\n");
+		}
+
+		Crystal *cr = crystal_new();
+		if (cr == NULL) {
+			ERROR("Failed to allocate crystal.\n");
+			return 0;
+		}
+		crystal_set_cell(cr, uc);
+		image_add_crystal(image, cr);
+
+	}
+
+	return goodLatticesCount;
+}
+
+
+void *xgandalf_prepare(IndexingMethod *indm, UnitCell *cell,
+                       struct xgandalf_options *xgandalf_opts)
+{
+	struct xgandalf_private_data *xgandalf_private_data =
+	                              malloc(sizeof(struct xgandalf_private_data));
+	allocReciprocalPeaks(&(xgandalf_private_data->reciprocalPeaks_1_per_A));
+	xgandalf_private_data->indm = *indm;
+	xgandalf_private_data->cellTemplate = NULL;
+
+	float tolerance = xgandalf_opts->tolerance;
+	samplingPitch_t samplingPitch = xgandalf_opts->sampling_pitch;
+	gradientDescentIterationsCount_t gradientDescentIterationsCount =
+	                                xgandalf_opts->grad_desc_iterations;
+
+	if (*indm & INDEXING_USE_CELL_PARAMETERS) {
+
+		xgandalf_private_data->cellTemplate = cell;
+
+		UnitCell* primitiveCell = uncenter_cell(cell, NULL);
+		reduceCell(primitiveCell);
+
+		double asx, asy, asz, bsx, bsy, bsz, csx, csy, csz;
+		int ret = cell_get_reciprocal(primitiveCell, &asx, &asy, &asz,
+		                                             &bsx, &bsy, &bsz,
+		                                             &csx, &csy, &csz);
+		if (ret != 0) {
+			ERROR("cell_get_reciprocal did not finish properly!");
+		}
+
+		Lattice_t sampleReciprocalLattice_1_per_A = {
+		        .ax = asx * 1e-10, .ay = asy * 1e-10, .az = asz * 1e-10,
+		        .bx = bsx * 1e-10, .by = bsy * 1e-10, .bz = bsz * 1e-10,
+		        .cx = csx * 1e-10, .cy = csy * 1e-10, .cz = csz * 1e-10 };
+
+		double ax, ay, az, bx, by, bz, cx, cy, cz;
+		ret = cell_get_cartesian(primitiveCell, &ax, &ay, &az,
+		                                        &bx, &by, &bz,
+		                                        &cx, &cy, &cz);
+		if (ret != 0) {
+			ERROR("cell_get_cartesian did not finish properly!");
+		}
+		Lattice_t sampleRealLatticeReduced_A = {
+		        .ax = ax * 1e10, .ay = ay * 1e10, .az = az * 1e10,
+		        .bx = bx * 1e10, .by = by * 1e10, .bz = bz * 1e10,
+		        .cx = cx * 1e10, .cy = cy * 1e10, .cz = cz * 1e10 };
+		xgandalf_private_data->sampleRealLatticeReduced_A = sampleRealLatticeReduced_A;
+
+		ExperimentSettings *experimentSettings =
+				ExperimentSettings_new(FAKE_BEAM_ENERGY,
+				                       FAKE_DETECTOR_DISTANCE,
+				                       FAKE_DETECTOR_RADIUS,
+				                       FAKE_DIVERGENCE_ANGLE_DEG,
+				                       FAKE_NON_MONOCHROMATICITY,
+				                       sampleReciprocalLattice_1_per_A,
+				                       tolerance,
+				                       FAKE_REFLECTION_RADIUS);
+
+		xgandalf_private_data->indexer = IndexerPlain_new(experimentSettings);
+		IndexerPlain_setSamplingPitch(xgandalf_private_data->indexer,
+		                              samplingPitch);
+		IndexerPlain_setGradientDescentIterationsCount(xgandalf_private_data->indexer,
+		                                               gradientDescentIterationsCount);
+
+		if (xgandalf_opts->no_deviation_from_provided_cell) {
+			IndexerPlain_setRefineWithExactLattice(xgandalf_private_data->indexer, 1);
+		}
+
+		ExperimentSettings_delete(experimentSettings);
+		cell_free(primitiveCell);
+
+	} else {
+
+		Lattice_t sampleRealLatticeReduced_A = { 0, 0, 0, 0, 0, 0, 0, 0, 0 };
+		xgandalf_private_data->sampleRealLatticeReduced_A = sampleRealLatticeReduced_A;
+
+		ExperimentSettings *experimentSettings =
+		   ExperimentSettings_new_nolatt(FAKE_BEAM_ENERGY,
+		                                 FAKE_DETECTOR_DISTANCE,
+		                                 FAKE_DETECTOR_RADIUS,
+		                                 FAKE_DIVERGENCE_ANGLE_DEG,
+		                                 FAKE_NON_MONOCHROMATICITY,
+		                                 xgandalf_opts->minLatticeVectorLength_A,
+		                                 xgandalf_opts->maxLatticeVectorLength_A,
+		                                 FAKE_REFLECTION_RADIUS);
+
+		xgandalf_private_data->indexer = IndexerPlain_new(experimentSettings);
+		IndexerPlain_setSamplingPitch(xgandalf_private_data->indexer,
+		                              samplingPitch);
+		IndexerPlain_setGradientDescentIterationsCount(xgandalf_private_data->indexer,
+		                                               gradientDescentIterationsCount);
+
+		ExperimentSettings_delete(experimentSettings);
+	}
+
+	/* Flags that XGANDALF knows about */
+	*indm &= INDEXING_METHOD_MASK | INDEXING_USE_CELL_PARAMETERS;
+
+	return xgandalf_private_data;
+}
+
+
+void xgandalf_cleanup(void *pp)
+{
+	struct xgandalf_private_data *xgandalf_private_data = pp;
+
+	freeReciprocalPeaks(xgandalf_private_data->reciprocalPeaks_1_per_A);
+	IndexerPlain_delete(xgandalf_private_data->indexer);
+}
+
+
+static void reduceCell(UnitCell *cell)
+{
+	double ax, ay, az, bx, by, bz, cx, cy, cz;
+	cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
+
+	Lattice_t l = { ax, ay, az, bx, by, bz, cx, cy, cz };
+
+	reduceLattice(&l);
+
+	cell_set_cartesian(cell, l.ax, l.ay, l.az,
+	                         l.bx, l.by, l.bz,
+	                         l.cx, l.cy, l.cz);
+
+	makeRightHanded(cell);
+}
+
+
+static void makeRightHanded(UnitCell *cell)
+{
+	double ax, ay, az, bx, by, bz, cx, cy, cz;
+	cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
+
+	if ( !right_handed(cell) ) {
+		cell_set_cartesian(cell, -ax, -ay, -az, -bx, -by, -bz, -cx, -cy, -cz);
+	}
+}
+
+
+const char *xgandalf_probe(UnitCell *cell)
+{
+	return "xgandalf";
+}
+
+#else
+
+int run_xgandalf(struct image *image, void *ipriv)
+{
+	ERROR("This copy of CrystFEL was compiled without XGANDALF support.\n");
+	return 0;
+}
+
+
+void *xgandalf_prepare(IndexingMethod *indm, UnitCell *cell,
+                       struct xgandalf_options *xgandalf_opts)
+{
+	ERROR("This copy of CrystFEL was compiled without XGANDALF support.\n");
+	ERROR("To use XGANDALF indexing, recompile with XGANDALF.\n");
+	return NULL;
+}
+
+
+void xgandalf_cleanup(void *pp)
+{
+}
+
+
+const char *xgandalf_probe(UnitCell *cell)
+{
+	return NULL;
+}
+
+#endif // HAVE_XGANDALF
diff --git a/libcrystfel/src/xgandalf.h b/libcrystfel/src/xgandalf.h
new file mode 100644
index 00000000..1aced417
--- /dev/null
+++ b/libcrystfel/src/xgandalf.h
@@ -0,0 +1,58 @@
+/*
+ * xgandalf.h
+ *
+ * Interface to XGANDALF indexer
+ *
+ * Copyright © 2017-2018 Deutsches Elektronen-Synchrotron DESY,
+ *                       a research centre of the Helmholtz Association.
+ *
+ * Authors:
+ *   2017-2018 Yaroslav Gevorkov <yaroslav.gevorkov@desy.de>
+ *
+ * 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/>.
+ *
+ */
+
+#ifndef LIBCRYSTFEL_SRC_XGANDALF_H
+#define LIBCRYSTFEL_SRC_XGANDALF_H
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <stddef.h>
+
+struct xgandalf_options {
+	unsigned int sampling_pitch;
+	unsigned int grad_desc_iterations;
+	float tolerance;
+	unsigned int no_deviation_from_provided_cell;
+	float minLatticeVectorLength_A;
+	float maxLatticeVectorLength_A;
+};
+
+#include "index.h"
+
+extern int run_xgandalf(struct image *image, void *ipriv);
+
+extern void *xgandalf_prepare(IndexingMethod *indm, UnitCell *cell,
+                              struct xgandalf_options *xgandalf_opts);
+
+extern void xgandalf_cleanup(void *pp);
+extern const char *xgandalf_probe(UnitCell *cell);
+
+
+#endif /* LIBCRYSTFEL_SRC_XGANDALF_H */