Implement smallcell indexing algorithm

5 files changed, 1132 insertions, 0 deletions

diff --git a/libcrystfel/meson.build b/libcrystfel/meson.build
index a8d43b4a..1fc89253 100644
--- a/libcrystfel/meson.build
+++ b/libcrystfel/meson.build
@@ -152,6 +152,7 @@ libcrystfel_sources = ['src/image.c',
                        'src/indexers/xgandalf.c',
                        'src/indexers/pinkindexer.c',
                        'src/indexers/fromfile.c',
+                       'src/indexers/smallcell.c',
                        symop_lex_ch,
                        symop_parse_ch]
 
diff --git a/libcrystfel/src/index.c b/libcrystfel/src/index.c
index d6c58c4a..c684cc1f 100644
--- a/libcrystfel/src/index.c
+++ b/libcrystfel/src/index.c
@@ -60,6 +60,7 @@
 #include "indexers/xgandalf.h"
 #include "indexers/pinkindexer.h"
 #include "indexers/fromfile.h"
+#include "indexers/smallcell.h"
 
 #include "uthash.h"
 
@@ -161,6 +162,10 @@ char *base_indexer_str(IndexingMethod indm)
 		strcpy(str, "file");
 		break;
 
+		case INDEXING_SMALLCELL :
+                strcpy(str, "smallcell");
+                break;
+
 		default :
 		strcpy(str, "(unknown)");
 		break;
@@ -197,6 +202,7 @@ static void *prepare_method(IndexingMethod *m, UnitCell *cell,
                             struct felix_options *felix_opts,
                             struct taketwo_options *taketwo_opts,
                             struct fromfile_options *fromfile_opts,
+			    struct smallcell_options *smallcell_opts,
                             struct asdf_options *asdf_opts)
 {
 	char *str;
@@ -229,6 +235,10 @@ static void *prepare_method(IndexingMethod *m, UnitCell *cell,
 		priv = fromfile_prepare(m, fromfile_opts);
 		break;
 
+		case INDEXING_SMALLCELL :
+                priv = smallcell_prepare(m, smallcell_opts, cell);
+                break;
+
 		case INDEXING_FELIX :
 		priv = felix_prepare(m, cell, felix_opts);
 		break;
@@ -328,6 +338,7 @@ IndexingPrivate *setup_indexing(const char *method_list,
                                 struct pinkindexer_options *pinkIndexer_opts,
                                 struct felix_options *felix_opts,
                                 struct fromfile_options *fromfile_opts,
+				struct smallcell_options *smallcell_opts,
                                 struct asdf_options *asdf_opts)
 {
 	IndexingPrivate *ipriv;
@@ -400,6 +411,7 @@ IndexingPrivate *setup_indexing(const char *method_list,
 		                                          felix_opts,
 		                                          ttopts,
 		                                          fromfile_opts,
+							  smallcell_opts,
 							  asdf_opts);
 
 		if ( ipriv->engine_private[i] == NULL ) {
@@ -489,6 +501,10 @@ void cleanup_indexing(IndexingPrivate *ipriv)
 			fromfile_cleanup(ipriv->engine_private[n]);
 			break;
 
+			case INDEXING_SMALLCELL :
+                        smallcell_cleanup(ipriv->engine_private[n]);
+                        break;
+
 			case INDEXING_TAKETWO :
 			taketwo_cleanup(ipriv->engine_private[n]);
 			break;
@@ -614,6 +630,13 @@ static int try_indexer(struct image *image, IndexingMethod indm,
 		profile_end("fromfile");
 		break;
 
+		case INDEXING_SMALLCELL :
+                set_last_task("indexing:smallcell");
+                profile_start("smallcell");
+                r = smallcell_index(image, mpriv);
+                profile_end("smallcell");
+                break;
+
 		case INDEXING_FELIX :
 		set_last_task("indexing:felix");
 		profile_start("felix");
@@ -1123,6 +1146,11 @@ IndexingMethod get_indm_from_string_2(const char *str, int *err)
 			method = INDEXING_FILE;
 			return method;
 
+		} else if ( strcmp(bits[i], "smallcell") == 0) {
+                        if ( have_method ) return warn_method(str);
+                        method = INDEXING_SMALLCELL;
+                        return method;
+
 		} else if ( strcmp(bits[i], "latt") == 0) {
 			method = set_lattice(method);
 
@@ -1225,6 +1253,7 @@ void default_method_options(struct taketwo_options **ttopts,
                             struct pinkindexer_options **pinkIndexer_opts,
                             struct felix_options **felix_opts,
                             struct fromfile_options **fromfile_opts,
+			    struct smallcell_options **smallcell_opts,
                             struct asdf_options **asdf_opts)
 {
 	taketwo_default_options(ttopts);
@@ -1232,5 +1261,6 @@ void default_method_options(struct taketwo_options **ttopts,
 	pinkIndexer_default_options(pinkIndexer_opts);
 	felix_default_options(felix_opts);
 	fromfile_default_options(fromfile_opts);
+	smallcell_default_options(smallcell_opts);
 	asdf_default_options(asdf_opts);
 }
diff --git a/libcrystfel/src/index.h b/libcrystfel/src/index.h
index e8167674..db6d77d0 100644
--- a/libcrystfel/src/index.h
+++ b/libcrystfel/src/index.h
@@ -54,6 +54,11 @@
                                    | INDEXING_USE_CELL_PARAMETERS \
                                    | INDEXING_USE_LATTICE_TYPE)
 
+#define INDEXING_DEFAULTS_SMALLCELL (INDEXING_SMALLCELL \
+                                   | INDEXING_USE_CELL_PARAMETERS \
+                                   | INDEXING_USE_LATTICE_TYPE)
+
+
 #define INDEXING_DEFAULTS_XDS (INDEXING_XDS | INDEXING_USE_LATTICE_TYPE \
                                      | INDEXING_USE_CELL_PARAMETERS)
 
@@ -78,6 +83,7 @@ typedef enum {
 	INDEXING_TAKETWO = 9,     /**< Use built-in TakeTwo algorithm */
 	INDEXING_XGANDALF = 10,   /**< Use XGANDALF (via optional library) */
 	INDEXING_PINKINDEXER = 11,/**< Use PinkIndexer (via optional library) */
+	INDEXING_SMALLCELL = 12,/**< Use Smallcell */
 
 	INDEXING_ERROR = 255,     /**< Special value for unrecognised indexing
 	                           *   engine */
@@ -172,6 +178,10 @@ struct fromfile_options
 	char *filename;
 };
 
+struct smallcell_options
+{
+       char *filename; 
+};
 
 struct felix_options
 {
@@ -217,6 +227,7 @@ extern struct argp pinkIndexer_argp;
 extern struct argp taketwo_argp;
 extern struct argp xgandalf_argp;
 extern struct argp fromfile_argp;
+extern struct argp smallcell_argp;
 extern struct argp asdf_argp;
 
 extern void default_method_options(struct taketwo_options **ttopts,
@@ -224,6 +235,7 @@ extern void default_method_options(struct taketwo_options **ttopts,
                                    struct pinkindexer_options **pinkIndexer_opts,
                                    struct felix_options **felix_opts,
                                    struct fromfile_options **fromfile_opts,
+				   struct smallcell_options **smallcell_opts,
                                    struct asdf_options **asdf_opts);
 
 extern IndexingPrivate *setup_indexing(const char *methods,
@@ -238,6 +250,7 @@ extern IndexingPrivate *setup_indexing(const char *methods,
                                        struct pinkindexer_options *pinkIndexer_opts,
                                        struct felix_options *felix_opts,
                                        struct fromfile_options *fromfile_opts,
+				       struct smallcell_options *smallcell_opts,
                                        struct asdf_options *asdf_opts);
 
 extern const IndexingMethod *indexing_methods(IndexingPrivate *p, int *n);
diff --git a/libcrystfel/src/indexers/smallcell.c b/libcrystfel/src/indexers/smallcell.c
new file mode 100644
index 00000000..ff29128f
--- /dev/null
+++ b/libcrystfel/src/indexers/smallcell.c
@@ -0,0 +1,1043 @@
+/*
+ * smallcell.c
+ *
+ * Perform indexing from solution file
+ *
+ * Copyright © 2020-2021 Max-Planck-Gesellschaft
+ *                       zur Förderung der Wissenschaften e.V.
+ * Copyright © 2021 Deutsches Elektronen-Synchrotron DESY,
+ *                  a research centre of the Helmholtz Association.
+ *
+ * Authors:
+ *   2020 Pascal Hogan-Lamarre <pascal.hogan.lamarre@mail.utoronto.ca>
+ *   2021 Thomas White <thomas.white@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 <libcrystfel-config.h>
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+#include <string.h>
+#include <assert.h>
+#include <fenv.h>
+#include <unistd.h>
+#include <argp.h>
+
+#include "image.h"
+#include "index.h"
+#include "cell-utils.h"
+#include "symmetry.h"
+#include "image.h"
+#include "utils.h"
+#include "peaks.h"
+#include "geometry.h"
+#include "detgeom.h"
+#include <gsl/gsl_matrix.h>
+#include <gsl/gsl_vector.h>
+#include <gsl/gsl_blas.h>
+#include <gsl/gsl_linalg.h>
+#include <gsl/gsl_multifit.h>
+#include <gsl/gsl_fit.h>
+#include "cell.h"
+#include "taketwo.h"
+
+/** \file smallcell.h */
+
+#define MAX_NEIGH (50)
+#define MAX_NODES (8072)
+#define DIFF_TOL (1e8)
+#define PIXEL_RING_TOL (6)
+struct sortmerefl
+{
+	signed int h;
+	signed int k;
+	signed int l;
+	double resolution;
+	int multi;
+	int ring_number;
+};
+
+
+static int cmpres(const void *av, const void *bv)
+{
+	const struct sortmerefl *a = av;
+	const struct sortmerefl *b = bv;
+	return a->resolution > b->resolution;
+}
+
+struct g_matrix
+{
+
+	double A;
+	double B;
+	double C;
+	double D;
+	double E;
+	double F;
+	double G;
+	double H;
+	double J;
+
+};
+
+struct smallcell_private
+{
+	struct sortmerefl *powderrings;
+	int num_rings;
+	SymOpList *sym;
+	struct g_matrix g9;
+	UnitCell *template;
+};
+
+
+void *smallcell_prepare(IndexingMethod * indm, struct smallcell_options *opts,
+                        UnitCell * cell)
+{
+	STATUS("\n");
+	STATUS("*******************************************************************\n");
+	STATUS("****                    Welcome to SmallCell                  ****\n");
+	STATUS("*******************************************************************\n");
+	STATUS("\n");
+
+
+	struct smallcell_private *dp;
+
+	char *add_unique_axis(const char *inp, char ua)
+	{
+		char *pg = cfmalloc(64);
+		if (pg == NULL)
+			return NULL;
+		snprintf(pg, 63, "%s_ua%c", inp, ua);
+		return pg;
+	}
+
+
+	char *get_chiral_holohedry(UnitCell * cell)
+	{
+		LatticeType lattice = cell_get_lattice_type(cell);
+		char *pg;
+		int add_ua = 1;
+
+		switch (lattice) {
+			case L_TRICLINIC:
+				pg = "1";
+				add_ua = 0;
+				break;
+
+			case L_MONOCLINIC:
+				pg = "2";
+				break;
+
+			case L_ORTHORHOMBIC:
+				pg = "222";
+				add_ua = 0;
+				break;
+
+			case L_TETRAGONAL:
+				pg = "422";
+				break;
+
+			case L_RHOMBOHEDRAL:
+				pg = "3_R";
+				add_ua = 0;
+				break;
+
+			case L_HEXAGONAL:
+				if (cell_get_centering(cell) == 'H') {
+					pg = "3_H";
+					add_ua = 0;
+				} else {
+					pg = "622";
+				}
+				break;
+
+			case L_CUBIC:
+				pg = "432";
+				add_ua = 0;
+				break;
+
+			default:
+				pg = "error";
+				break;
+		}
+
+		if (add_ua) {
+			return add_unique_axis(pg, cell_get_unique_axis(cell));
+		} else {
+			return cfstrdup(pg);
+		}
+	}
+
+
+	SymOpList *sym;
+
+	char *pg = get_chiral_holohedry(cell);
+	sym = get_pointgroup(pg);
+	cffree(pg);
+
+	double mres = 1 / (1e-10);      // 1/Angstrom
+
+	double ax, ay, az;
+	double bx, by, bz;
+	double cx, cy, cz;
+	int hmax, kmax, lmax;
+	signed int h, k, l;
+	RefList *reflist;
+	int i, n;
+	RefListIterator *iter;
+	Reflection *ring;
+	struct sortmerefl *sortus;
+
+	dp = cfmalloc(sizeof(struct smallcell_private));
+
+	cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
+	hmax = mres * modulus(ax, ay, az);
+	kmax = mres * modulus(bx, by, bz);
+	lmax = mres * modulus(cx, cy, cz);
+	reflist = reflist_new();
+	for (h = -hmax; h <= hmax; h++) {
+		for (k = -kmax; k <= kmax; k++) {
+			for (l = -lmax; l <= lmax; l++) {
+
+				signed int ha, ka, la;
+
+				if (forbidden_reflection(cell, h, k, l))
+					continue;
+				if (2.0 * resolution(cell, h, k, l) > mres)
+					continue;
+
+				if (sym != NULL) {
+
+					Reflection *refl;
+
+					get_asymm(sym, h, k, l, &ha, &ka, &la);
+					refl = find_refl(reflist, ha, ka, la);
+					if (refl == NULL) {
+						refl = add_refl(reflist, ha, ka,
+						                la);
+						set_redundancy(refl, 1);
+					} else {
+						set_redundancy(refl,
+						               get_redundancy
+						               (refl) + 1);
+					}
+
+				} else {
+					Reflection *refl;
+					refl = add_refl(reflist, h, k, l);
+					set_redundancy(refl, 1);
+				}
+
+			}
+		}
+	}
+
+	n = num_reflections(reflist);
+	sortus = cfmalloc(n * sizeof(struct sortmerefl));
+	i = 0;
+	for (ring = first_refl(reflist, &iter);
+	     ring != NULL; ring = next_refl(ring, iter)) {
+		signed int rh, rk, rl;
+
+		get_indices(ring, &rh, &rk, &rl);
+
+		sortus[i].ring_number = i;      // how many rings
+		sortus[i].h = rh;
+		sortus[i].k = rk;
+		sortus[i].l = rl;
+		sortus[i].resolution = 2.0 * resolution(cell, rh, rk, rl);      /* one over d */
+		sortus[i].multi = get_redundancy(ring);
+		i++;
+
+	}
+
+	qsort(sortus, n, sizeof(struct sortmerefl), cmpres);
+
+	STATUS("\nAll powder rings up to %f Ångstrøms.\n", 1e+10 / mres);
+	STATUS("Note that screw axis or glide plane absences are not "
+	       "omitted from this list.\n");
+	STATUS("\n  No.  d (Å)   1/d (m^-1)    h    k    l    multiplicity\n");
+	STATUS("------------------------------------------------------\n");
+	for (i = 0; i < n; i++) {
+		printf("%4i  %10.3f %e %4i %4i %4i    m = %i\n",
+		       i + 1, 1e10 / sortus[i].resolution, sortus[i].resolution,
+		       sortus[i].h, sortus[i].k, sortus[i].l, sortus[i].multi);
+	}
+
+
+	// get reciprocal unit cell elements in order to create G* matrix and store in private
+
+	double asx;
+	double bsx;
+	double csx;
+	double asy;
+	double bsy;
+	double csy;
+	double asz;
+	double bsz;
+	double csz;
+
+	cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
+
+	cell_get_reciprocal(cell,
+	                    &asx, &asy, &asz,
+	                    &bsx, &bsy, &bsz, &csx, &csy, &csz);
+
+	dp->g9.A = (asx * asx) + (asy * asy) + (asz * asz);
+	dp->g9.B = (asx * bsx) + (asy * bsy) + (asz * bsz);
+	dp->g9.C = (asx * csx) + (asy * csy) + (asz * csz);
+	dp->g9.D = (bsx * asx) + (bsy * asy) + (bsz * asz);
+	dp->g9.E = (bsx * bsx) + (bsy * bsy) + (bsz * bsz);
+	dp->g9.F = (bsx * csx) + (bsy * csy) + (bsz * csz);
+	dp->g9.G = (csx * asx) + (csy * asy) + (csz * asz);
+	dp->g9.H = (csx * bsx) + (csy * bsy) + (csz * bsz);
+	dp->g9.J = (csx * csx) + (csy * csy) + (csz * csz);
+
+
+	dp->sym = sym;
+	dp->powderrings = sortus;
+	dp->num_rings = n;
+
+	return dp;
+}
+
+/* PeakInfo structure for storing peak-related information for every match found*/
+typedef struct PeakInfo
+{
+	int peak_number;
+	double peak_res;
+	double x, y, z;
+	int h, k, l;
+	struct PeakInfo *neigh[MAX_NEIGH];
+	int n_neigh;
+	double weight_list[MAX_NEIGH];
+} PeakInfo;
+
+struct Nodelist
+{
+	int n_mem;
+	struct PeakInfo *mem[MAX_NODES];
+};
+
+struct Cliquelist
+{
+	int n;
+	struct Nodelist *list[MAX_NEIGH];
+};
+
+static struct Nodelist *CopyRlist(struct Nodelist *R)
+{
+	struct Nodelist *Rcopy = cfmalloc(sizeof(struct Nodelist));
+	Rcopy->n_mem = 0;
+	for (int i = 0; i < R->n_mem; i++) {
+		Rcopy->mem[i] = R->mem[i];
+		Rcopy->n_mem++;
+	}
+	return Rcopy;
+};
+
+//To make a list of neighbour nodes
+static struct Nodelist *neighbours(struct PeakInfo *in)
+{
+	int i;
+	struct Nodelist *list = cfmalloc(sizeof(struct Nodelist));
+	for (i = 0; i < in->n_neigh; i++) {
+		list->mem[i] = in->neigh[i];
+	}
+	list->n_mem = in->n_neigh;
+	return list;
+};
+
+//Function to see if a node is in a list
+static int isin(struct Nodelist *b, struct PeakInfo *test)
+{
+	int i;
+	for (i = 0; i < b->n_mem; i++) {
+		if (test == b->mem[i]) {
+			return 1;       //isin
+		}
+	}
+	return 0;               //is notin
+};
+
+static void add(struct Nodelist *c, struct PeakInfo *test)
+{
+
+	if (isin(c, test))
+		return;
+	assert(c->n_mem < MAX_NODES);
+	c->mem[c->n_mem] = test;
+	c->n_mem++;
+};
+
+//To create a new list of nodes from the union between two node lists (i.e. c = a U b)
+static struct Nodelist *Union(struct Nodelist *a, struct Nodelist *b)
+{
+	struct Nodelist *c = cfmalloc(sizeof(struct Nodelist));
+	c->n_mem = 0;
+	int i;
+	int j;
+	for (i = 0; i < a->n_mem; i++) {
+		add(c, a->mem[i]);
+	}
+	for (j = 0; j < b->n_mem; j++) {
+		add(c, b->mem[j]);
+	}
+	return c;
+};
+
+//To create a new list of nodes from the intersection of two lists (i.e. c = a intersection b)
+static struct Nodelist *intersection(struct Nodelist *a, struct Nodelist *b)
+{
+	struct Nodelist *c = cfmalloc(sizeof(struct Nodelist));
+	c->n_mem = 0;
+	int j;
+	for (j = 0; j < a->n_mem; j++) {
+		if (isin(b, a->mem[j]) == 1) {
+			add(c, a->mem[j]);
+		}
+	}
+	return c;
+};
+
+//To create a new list of nodes (c) from the exclusion of list b from list a (i.e. c = a\b)
+static struct Nodelist *exclusion(struct Nodelist *a, struct Nodelist *b)
+{
+	struct Nodelist *c = cfmalloc(sizeof(struct Nodelist));
+	c->n_mem = 0;
+	int j;
+	for (j = 0; j < a->n_mem; j++) {
+		if (isin(b, a->mem[j]) == 0) {
+			add(c, a->mem[j]);
+		}
+	}
+	return c;
+};
+
+//Function to exclude a single node
+static struct Nodelist *exclunode(struct Nodelist *a, struct PeakInfo *v)
+{
+	int j;
+	for (j = 0; j < a->n_mem; j++) {
+		if (a->mem[j] == v) {
+			int hole_index = j;
+			for (int i = hole_index; i < a->n_mem - 1; i++) {
+				a->mem[i] = a->mem[i + 1];
+			}
+		}
+	}
+	a->n_mem = a->n_mem - 1;
+	return a;
+};
+
+//Function to append a node to a list (creating a new list)
+static struct Nodelist *append(struct Nodelist *a, struct PeakInfo *v)
+{
+	struct Nodelist *c = cfmalloc(sizeof(struct Nodelist));
+	c->n_mem = 0;
+	int j;
+	for (j = 0; j < a->n_mem; j++) {
+		add(c, a->mem[j]);
+	}
+
+	add(c, v);
+	return c;
+};
+
+//Average weight function
+static double avg_weight(int num_neigh, double *weights)
+{
+
+	double avg;
+	double avg_add = 0.0;
+	for (int j = 0; j < num_neigh; j++) {
+		avg_add = avg_add + weights[j];
+	}
+
+	avg = avg_add / num_neigh;
+
+	return avg;
+};
+
+//Bron-Kerbosch algorithm, with pivoting
+void BK(struct Nodelist *R, struct Nodelist *P, struct Nodelist *X,
+        struct Cliquelist *Max_cliques)
+{
+
+	//If P & X are empty -> add R to the un-mapped clique array
+	if (P->n_mem == 0 && X->n_mem == 0) {
+		Max_cliques->list[Max_cliques->n] = CopyRlist(R);
+		Max_cliques->n++;
+
+
+		return;
+	}
+	//Find pivot u from set of nodes in P U X
+	struct Nodelist *piv_pool = Union(P, X);
+	struct PeakInfo *piv = NULL;
+	piv = piv_pool->mem[0];
+	int n_max = piv_pool->mem[0]->n_neigh;
+	double n_max_weight =
+	        avg_weight(piv_pool->mem[0]->n_neigh,
+	                   piv_pool->mem[0]->weight_list);
+	for (int i = 1; i < piv_pool->n_mem; i++) {
+		if (piv_pool->mem[i]->n_neigh > n_max) {
+			n_max = piv_pool->mem[i]->n_neigh;
+			piv = piv_pool->mem[i];
+			n_max_weight =
+			        avg_weight(piv_pool->mem[i]->n_neigh,
+			                   piv_pool->mem[i]->weight_list);
+		} else if (piv_pool->mem[i]->n_neigh == n_max) {
+			if (piv_pool->n_mem == 1) {
+				piv = piv_pool->mem[0];
+			} else if ((n_max > 0)
+			           &&
+			           (avg_weight
+			            (piv_pool->mem[i]->n_neigh,
+			             piv_pool->mem[i]->weight_list) <
+			            n_max_weight)) {
+				piv = piv_pool->mem[i];
+			}
+		}
+	}
+
+	if (piv == NULL) {
+		printf("Couldn't find pivot\n");
+		return;
+	}
+	//Get list of neighbours of the pivot
+	struct Nodelist *piv_neighbours = neighbours(piv);
+	//Remove pivot neighbours from P list
+	struct Nodelist *P_excl = exclusion(P, piv_neighbours);
+	cffree(piv_pool);
+	cffree(piv_neighbours);
+
+
+	for (int u = 0; u < P_excl->n_mem; u++) {
+		//For all members of the subset of P we are concered with (P_excl = P\N(pivot))
+		struct PeakInfo *v = P_excl->mem[u];
+
+		//create Nodelist for v using neighbour function
+		struct Nodelist *v_neighs = neighbours(v);
+
+
+		//add v to R
+		struct Nodelist *R_new = append(R, v);
+
+
+		//Find intersection of P_excl with N(v)
+		struct Nodelist *P_new = intersection(P, v_neighs);
+
+
+		//''                 '' X with N(v)    
+		struct Nodelist *X_new = intersection(X, v_neighs);
+		BK(R_new, P_new, X_new, Max_cliques);
+
+		cffree(R_new);
+		cffree(P_new);
+		cffree(X_new);
+		cffree(v_neighs);
+
+
+		//Redefine P and X as P\v and XUv
+
+		exclunode(P, v);
+
+		add(X, v);
+	}
+	cffree(P_excl);
+	return;
+};
+
+
+int smallcell_index(struct image *image, void *mpriv)
+{
+	struct smallcell_private *priv = (struct smallcell_private *) mpriv;
+	struct sortmerefl *powderrings = priv->powderrings;
+	int num_rings = priv->num_rings;
+
+
+	//Assigning image data
+	ImageFeatureList *peaks = image->features;
+	double lambda = image->lambda;
+	struct detgeom *det = image->detgeom;
+
+	//Calculating a tolerance value of PIXEL_RING_TOL  pixel lengths from a peak
+	double pixel_metres = det->panels[0].pixel_pitch;
+	double ten_pix_len = pixel_metres * PIXEL_RING_TOL;
+	double detector_len = det->panels[0].cnz;
+	double detector_len_m = pixel_metres * detector_len;
+	double tol = ten_pix_len / (lambda * detector_len_m);
+
+	int npk = image_feature_count(peaks);
+
+	// No peaks -> no resolution! 
+	if (npk <= 3) {
+		printf("not enough peaks\n");
+		return 0;
+	}
+
+	printf("Current Image %s, %s\n", image->filename, image->ev);
+
+	int peak_infos_size = 100;      // Arbitrary initial size allocation
+	PeakInfo *peak_infos = cfmalloc(peak_infos_size * sizeof(PeakInfo));
+
+	int num_peak_infos = 0;
+	int peaks_with_matches = 0;
+
+	//Loop through each peak, calculate d, then 1/d value (based on estimate_peak_resolution from peak.c), then use match_rings function to create structs for this peak with all possible h,k,l values
+	for (int i = 0; i < npk; i++) {
+		struct imagefeature *f = image_get_feature(peaks, i);
+		double r[3];
+		detgeom_transform_coords(&det->panels[f->pn], f->fs, f->ss,
+		                         lambda, 0.0, 0.0, r);
+		double x_res = r[0];
+		double y_res = (r[1]);
+		double z_res = (r[2]);
+		double rns = modulus(r[0], r[1], r[2]);
+		int init_num_peak_infos = num_peak_infos;
+
+		for (int j = 0; j < num_rings; j++) {
+			if (fabs(powderrings[j].resolution - rns) <= tol) {
+
+				signed int h = powderrings[j].h;
+				signed int k = powderrings[j].k;
+				signed int l = powderrings[j].l;
+
+				//Looking for symmetries and creating more PeakInfo structs with these symmetry indices
+				SymOpMask *m = new_symopmask(priv->sym);
+				special_position(priv->sym, m, h, k, l);
+				int n = num_equivs(priv->sym, m);
+				for (int y = 0; y < n; y++) {
+
+					if (num_peak_infos >= peak_infos_size) {        //Adding more structs if necessary
+						peak_infos_size *= 2;
+						peak_infos =
+						        cfrealloc(peak_infos,
+						                  peak_infos_size
+						                  *
+						                  sizeof
+						                  (PeakInfo));
+					}
+
+					signed int ha, ka, la;
+					get_equiv(priv->sym, m, y, h, k, l, &ha,
+					          &ka, &la);
+					// Add ha, ka, la to list of reflections 
+					peak_infos[num_peak_infos].peak_number =
+					        i;
+					peak_infos[num_peak_infos].peak_res =
+					        rns;
+					peak_infos[num_peak_infos].h = ha;
+					peak_infos[num_peak_infos].k = ka;
+					peak_infos[num_peak_infos].l = la;
+					peak_infos[num_peak_infos].x = x_res;
+					peak_infos[num_peak_infos].y = y_res;
+					peak_infos[num_peak_infos].z = z_res;
+					for (int w = 0; w < MAX_NEIGH; w++) {
+						peak_infos[num_peak_infos].
+						        neigh[w] = NULL;
+					}
+					peak_infos[num_peak_infos].n_neigh = 0; //Initialise
+					(num_peak_infos)++;
+				}
+				free_symopmask(m);
+			}
+
+		}
+
+		if (num_peak_infos != init_num_peak_infos) {
+			peaks_with_matches++;
+		}
+	}
+
+	printf("The number of matches in this image (including symmetric indices) is %d for %d/%d peaks\n", num_peak_infos, peaks_with_matches, npk);
+
+
+	// Now to connect the nodes using calculated and measured reciprocal distance
+
+
+	struct g_matrix g9 = priv->g9;
+	double dtol = DIFF_TOL;
+	int n_connected_nodes = 0;
+
+	//Loop through peak numbers
+	for (int j = 0; j < num_peak_infos; j++) {
+
+		int node_a_h = peak_infos[j].h;
+		int node_a_k = peak_infos[j].k;
+		int node_a_l = peak_infos[j].l;
+
+		//Loop through the rest of the peak infos
+		for (int y = j + 1; y < num_peak_infos - 1; y++) {
+
+			if (peak_infos[y].peak_number ==
+			    peak_infos[j].peak_number)
+				continue;
+
+			//Observed dist. d_1 = fabs(res.vec_j - res.vec_y)
+			double d_1 =
+			        modulus(peak_infos[j].x - peak_infos[y].x,
+			                peak_infos[j].y - peak_infos[y].y,
+			                peak_infos[j].z - peak_infos[y].z);
+
+			//Predicted d_2
+			int node_b_h = peak_infos[y].h;
+			int node_b_k = peak_infos[y].k;
+			int node_b_l = peak_infos[y].l;
+
+			//d_2 = sqrt(Transpose(MI_b - MI_a).G*.(MI_b - MI_a))
+			double d_2 =
+			        sqrt((node_b_h -
+			              node_a_h) * (g9.A * (node_b_h -
+			                                   node_a_h) +
+			                           g9.D * (node_b_k -
+			                                   node_a_k) +
+			                           g9.G * (node_b_l - node_a_l))
+			             + (node_b_k -
+			                node_a_k) * (g9.B * (node_b_h -
+			                                     node_a_h) +
+			                             g9.E * (node_b_k -
+			                                     node_a_k) +
+			                             g9.H * (node_b_l -
+			                                     node_a_l))
+			             + (node_b_l -
+			                node_a_l) * (g9.C * (node_b_h -
+			                                     node_a_h) +
+			                             g9.F * (node_b_k -
+			                                     node_a_k) +
+			                             g9.J * (node_b_l -
+			                                     node_a_l)));
+
+
+			double diff = fabs(d_2 - d_1);
+
+			//Now test difference
+			if (diff <= dtol) {
+
+				//Connect nodes
+
+				if (peak_infos[j].n_neigh <= MAX_NEIGH) {
+
+					peak_infos[j].neigh[peak_infos[j].n_neigh] = &peak_infos[y];    //Pointing to the info of the connected peak (adding node_b as neighbour due to its connection to node_a)
+					peak_infos[j].weight_list[peak_infos[j].n_neigh] = diff;        //add weight for later
+					peak_infos[j].n_neigh++;        //increasing number of connection/neighbours that node_a has
+
+				} else {
+
+					printf("The number of neighbours for this node has exceeded MAX_NEIGH, need to manage memory\n");
+				}
+
+				if (peak_infos[y].n_neigh <= MAX_NEIGH) {
+
+					peak_infos[y].neigh[peak_infos[y].n_neigh] = &peak_infos[j];    //Also pointing other way(i.e. adding node_a as a neghbour to node_b as we wont loop back through the upper part of the list when j = peak_no of node_b)
+					peak_infos[y].weight_list[peak_infos[y].n_neigh] = diff;        //add weight for later
+					peak_infos[y].n_neigh++;        //increasing number of connection/neighbours that node_b has
+				} else {
+					printf("The number of neighbours for this node has exceeded MAX_NEIGH, need to manage memory\n");
+				}
+
+			}
+
+		}
+
+		//Counting number of nodes with 1 or more connection (may be unnecessary)...
+		if (peak_infos[j].n_neigh != 0) {
+			n_connected_nodes++;
+		}
+
+	}
+
+	//Store the max.cliques found
+	struct Cliquelist *Max_cliques =
+	        cfmalloc(sizeof(struct Nodelist) * sizeof(struct Cliquelist));
+	Max_cliques->n = 0;
+	// R: array of nodes forming a clique (array of pointers to node infos)
+	// P: array of all prosepctive nodes that are connected to R which may be added to R. To begin, this is all nodes i.e all peak_infos
+	// X: exculsion set (same form as R but nodes that are NOT candidats for the max. clique, were originaly in P)
+	struct Nodelist *R = cfmalloc(sizeof(struct Nodelist));
+	struct Nodelist *X = cfmalloc(sizeof(struct Nodelist));
+	struct Nodelist *P = cfmalloc(sizeof(struct Nodelist));
+	R->n_mem = 0;
+	X->n_mem = 0;
+	P->n_mem = 0;
+	//To make P; create nodelist of all peak_infos
+	int i;
+	for (i = 0; i < num_peak_infos; i++) {
+		if (peak_infos[i].n_neigh != 0) {
+			add(P, &peak_infos[i]);
+		}
+	}
+
+	if (P->n_mem <= 2) {
+		printf("no peaks with neighbours\n");
+		return 0;
+	}
+	//Call BK using current peak info nodes for this image
+
+	printf("running BK\n");
+	BK(R, P, X, Max_cliques);
+	printf("done\n");
+	printf("The number of cliques found = %d.\n", Max_cliques->n);
+
+	//get the max. clique from list of Maximal cliques found
+	int Max_clique_len = Max_cliques->list[0]->n_mem;
+	struct Cliquelist *Max =
+	        cfmalloc(sizeof(struct Nodelist) * sizeof(struct Cliquelist));
+	Max->n = 0;
+	Max->list[0] = Max_cliques->list[0];
+	Max->n++;
+	for (int m = 1; m < Max_cliques->n; m++) {
+		if (Max_cliques->list[m]->n_mem > Max_clique_len) {
+			Max_clique_len = Max_cliques->list[m]->n_mem;
+			for (int t = 0; t < Max->n; t++) {
+				Max->list[t] = NULL;
+			}
+			Max->n = 0;
+			Max->list[0] = Max_cliques->list[m];
+			Max->n++;
+		} else if (Max_cliques->list[m]->n_mem == Max_clique_len) {
+			Max->list[Max->n] = Max_cliques->list[m];
+			Max->n++;
+		}
+	}
+	//If more than one max_clique with the same number of nodes is found, take only the right-handed solution
+	//This requires first getting the unit cell for each max_clique, and then using the right_handed function from cell-utils
+	for (int m = 0; m < Max->n; m++) {
+		if (Max->list[m]->n_mem < 5 && m == (Max->n) - 1)
+			return 0;
+		if (Max->list[m]->n_mem < 5)
+			continue;
+		gsl_matrix *h_mat =
+		        gsl_matrix_calloc(3 * (Max->list[m]->n_mem), 9);
+		gsl_vector *h_vec = gsl_vector_alloc(3 * (Max->list[m]->n_mem));
+
+		int count_node = 0;
+		int col_count = 0;
+		int have_a = 0, have_b = 0, have_c = 0;
+		for (int i = 0; i < 3 * (Max->list[m]->n_mem); i++) {
+			if (i > 0 && i % 3 == 0) {
+				count_node++;
+				col_count = 0;
+			}
+
+			if (Max->list[m]->mem[count_node]->h != 0)
+				have_a = 1;
+			if (Max->list[m]->mem[count_node]->k != 0)
+				have_b = 1;
+			if (Max->list[m]->mem[count_node]->l != 0)
+				have_c = 1;
+
+			gsl_matrix_set(h_mat, i, col_count,
+			               Max->list[m]->mem[count_node]->h);
+			gsl_matrix_set(h_mat, i, col_count + 3,
+			               Max->list[m]->mem[count_node]->k);
+			gsl_matrix_set(h_mat, i, col_count + 6,
+			               Max->list[m]->mem[count_node]->l);
+			col_count++;
+		}
+
+		int count_mem = 0;
+		int count_comp = 0;
+		for (int i = 0; i < 3 * (Max->list[m]->n_mem); i++) {
+
+			gsl_vector_set(h_vec, i,
+			               Max->list[m]->mem[count_mem]->x);
+			if (count_comp == 0) {
+				gsl_vector_set(h_vec, i,
+				               Max->list[m]->mem[count_mem]->x);
+				count_comp++;
+			} else if (count_comp == 1) {
+				gsl_vector_set(h_vec, i,
+				               Max->list[m]->mem[count_mem]->y);
+				count_comp++;
+			} else if (count_comp == 2) {
+				gsl_vector_set(h_vec, i,
+				               Max->list[m]->mem[count_mem]->z);
+				count_comp = 0;
+				count_mem++;
+			}
+		}
+		//Solve matrix-vector equation for unit-cell for this clique m
+		gsl_vector *cell_vecs = gsl_vector_alloc(9);
+		//cell_vec = [a*x a*y a*z b*x b*y b*z c*x c*y c*z] 
+		double chisq;
+		gsl_matrix *cov =
+		        gsl_matrix_alloc(3 * (Max->list[m]->n_mem), 9);
+		gsl_multifit_linear_workspace *work =
+		        gsl_multifit_linear_alloc(3 * (Max->list[m]->n_mem), 9);
+		if (gsl_multifit_linear
+		    (h_mat, h_vec, cell_vecs, cov, &chisq, work)) {
+			ERROR("Multifit failed\n");
+		}
+		//Use the following function to make a unit cell file then can use the checks directly to see if it's a viable solution
+
+		UnitCell *uc;
+
+		uc = cell_new();
+/*	UnitCell *cell = priv->template;
+	cell_set_lattice_type(uc, cell_get_lattice_type(cell));*/
+		cell_set_reciprocal(uc,
+		                    gsl_vector_get(cell_vecs, 0),
+		                    gsl_vector_get(cell_vecs, 1),
+		                    gsl_vector_get(cell_vecs, 2),
+		                    gsl_vector_get(cell_vecs, 3),
+		                    gsl_vector_get(cell_vecs, 4),
+		                    gsl_vector_get(cell_vecs, 5),
+		                    gsl_vector_get(cell_vecs, 6),
+		                    gsl_vector_get(cell_vecs, 7),
+		                    gsl_vector_get(cell_vecs, 8));
+
+
+		if (uc == NULL) {
+			printf("Unit Cell not created.. returned NULL\n");
+			continue;
+		}
+
+		//Free up matrix and vector memeories
+		gsl_multifit_linear_free(work);
+		gsl_vector_free(cell_vecs);
+		gsl_vector_free(h_vec);
+		gsl_matrix_free(cov);
+		gsl_matrix_free(h_mat);
+
+		if (!(have_a && have_b && have_c))
+			return 0;
+
+		printf("Unit Cell created, testing if valid solution..\n");
+		if (validate_cell(uc) == 0) {
+
+			Crystal *cr;
+			printf("unit cell valid!\n");
+			cell_print(uc);
+			cr = crystal_new();
+			if (cr == NULL) {
+				ERROR("Failed to allocate crystal.\n");
+				continue;
+			}
+			crystal_set_cell(cr, uc);
+			image_add_crystal(image, cr);
+			return 1;
+		}
+
+		cell_free(uc);
+
+	}
+
+
+	for (int t = 0; t < Max_cliques->n; t++) {
+		cffree(Max_cliques->list[t]);
+	}
+	for (int t = 0; t < Max->n; t++) {
+		cffree(Max->list[t]);
+	}
+
+	cffree(Max_cliques);
+	cffree(Max);
+	cffree(X);
+	cffree(P);
+	cffree(R);
+	cffree(peak_infos);
+	return 0;
+}
+
+
+void smallcell_cleanup(void *mpriv)
+{
+/*	struct smallcell_private *dp = mpriv;
+	struct smallcell_entry *item, *tmp;
+
+	HASH_ITER(hh, dp->sol_hash, item, tmp) {
+		int i;
+		HASH_DEL(dp->sol_hash, item);
+		for ( i=0; i<item->n_crystals; i++ ) {
+			Crystal *cr = item->crystals[i];
+			cell_free(crystal_get_cell(cr));
+			crystal_free(cr);
+		}
+	}
+
+	cffree(dp);*/
+}
+
+
+static void smallcell_show_help()
+{
+	printf("Parameters for 'smallcell' indexing:\n"
+	       "     --smallcell-input-file\n"
+	       "                           Filename of indexing solution file\n");
+}
+
+
+int smallcell_default_options(struct smallcell_options **opts_ptr)
+{
+	struct smallcell_options *opts;
+	opts = cfmalloc(sizeof(struct smallcell_options));
+	if (opts == NULL)
+		return ENOMEM;
+	opts->filename = NULL;
+	*opts_ptr = opts;
+	return 0;
+}
+
+
+static error_t smallcell_parse_arg(int key, char *arg, struct argp_state *state)
+{
+	struct smallcell_options **opts_ptr = state->input;
+	int r;
+
+	switch (key) {
+
+		case ARGP_KEY_INIT:
+			r = smallcell_default_options(opts_ptr);
+			if (r)
+				return r;
+			break;
+
+		case 1:
+			smallcell_show_help();
+			return EINVAL;
+
+		case 2:
+			(*opts_ptr)->filename = cfstrdup(arg);
+			break;
+
+		default:
+			return ARGP_ERR_UNKNOWN;
+
+	}
+
+	return 0;
+}
+
+
+static struct argp_option smallcell_options[] = {
+
+	{"help-smallcell", 1, NULL, OPTION_NO_USAGE,
+	 "Show options for 'from file' indexing", 99},
+
+	{"smallcell-input-file", 2, "filename", OPTION_HIDDEN, NULL},
+	{0}
+};
+
+
+struct argp smallcell_argp = { smallcell_options, smallcell_parse_arg,
+	NULL, NULL, NULL, NULL, NULL
+};
diff --git a/libcrystfel/src/indexers/smallcell.h b/libcrystfel/src/indexers/smallcell.h
new file mode 100644
index 00000000..726a998f
--- /dev/null
+++ b/libcrystfel/src/indexers/smallcell.h
@@ -0,0 +1,45 @@
+/*
+ * smallcell.h
+ *
+ * Perform indexing from solution file
+ *
+ * Copyright © 2020-2021 Max-Planck-Gesellschaft
+ *                       zur Förderung der Wissenschaften e.V.
+ * Copyright © 2021 Deutsches Elektronen-Synchrotron DESY,
+ *                  a research centre of the Helmholtz Association.
+ *
+ * Authors:
+ *   2020 Pascal Hogan-Lamarre <pascal.hogan.lamarre@mail.utoronto.ca>
+ *   2021 Thomas White <thomas.white@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 SMALLCELL_H
+#define SMALLCELL_H
+
+#include <argp.h>
+
+#include "image.h"
+
+extern int smallcell_default_options(struct smallcell_options **opts_ptr);
+extern void *smallcell_prepare(IndexingMethod *indm,
+                              struct smallcell_options *opts, UnitCell *cell);
+extern int smallcell_index(struct image *image, void *mpriv);
+extern void smallcell_cleanup(void *mpriv);
+
+#endif	/* SMALLCELL_H */