Merge pull request #21568 from mantidproject/21560-index-fractional-peaks

Add support for indexing fractional peaks

Framework/PythonInterface/plugins/algorithms/RefineSatellitePeaks.py → Framework/PythonInterface/plugins/algorithms/FindSatellitePeaks.py

@@ -2,18 +2,23 @@ from mantid.kernel import *
 from mantid.dataobjects import PeaksWorkspaceProperty
 from mantid.api import *
 from mantid.simpleapi import *
+import fractional_indexing as indexing
 
-import numpy as np
-from scipy.spatial import KDTree
-from scipy.cluster.vq import kmeans2
-import scipy.cluster.hierarchy as hcluster
 
-
-class RefineSatellitePeaks(DataProcessorAlgorithm):
+class FindSatellitePeaks(DataProcessorAlgorithm):
 
     def category(self):
         return 'Crystal\\Peaks'
 
+    def seeAlso(self):
+        return [ "IndexSatellitePeaks" ]
+
+    def name(self):
+        return "FindSatellitePeaks"
+
+    def summary(self):
+        return "Algorithm for finding satellite peaks in an MDWorkspace in the HKL frame."
+
     def PyInit(self):
         self.declareProperty(PeaksWorkspaceProperty(name="NuclearPeaks",
                                                     defaultValue="",
@@ -71,12 +76,12 @@ class RefineSatellitePeaks(DataProcessorAlgorithm):
         nuclear = self.getProperty("NuclearPeaks").value
         sats = self.getProperty("SatellitePeaks").value
 
-        nuclear_hkls = self.get_hkls(nuclear)
-        sats_hkls = self.get_hkls(sats)
+        nuclear_hkls = indexing.get_hkls(nuclear)
+        sats_hkls = indexing.get_hkls(sats)
 
-        qs = self.find_q_vectors(nuclear_hkls, sats_hkls)
-        clusters, k = self.cluster_qs(qs, threshold=cluster_threshold)
-        qs = self.average_clusters(qs, clusters)
+        qs = indexing.find_q_vectors(nuclear_hkls, sats_hkls)
+        clusters, k = indexing.cluster_qs(qs, threshold=cluster_threshold, k=k)
+        qs = indexing.average_clusters(qs, clusters)
         predicted_satellites = self.create_fractional_peaks_workspace(qs, nuclear)
 
         centroid_satellites = CentroidPeaksMD(InputWorkspace=md, PeaksWorkspace=predicted_satellites,
@@ -92,73 +97,6 @@ class RefineSatellitePeaks(DataProcessorAlgorithm):
         self.log().notice("Q vectors are: \n{}".format(qs))
         self.setProperty("OutputWorkspace", satellites_int_spherical)
 
-    def get_hkls(self, peaks_workspace):
-        """Return a 2D numpy array from a peaks workspace.
-
-        :param peaks_workpace: the PeaksWorkspace to extract HKL values from
-        :returns: np.ndarry -- 2D array of HKL values
-        """
-        return np.array([np.array([peak['h'], peak['k'], peak['l']]) for peak in peaks_workspace])
-
-    def cluster_qs(self, qs, k=None, threshold=1.5):
-        """Cluster q vectors into discrete groups.
-
-        Classifies each of the q vectors into a number of clusters. The number of clusters used is decided by the parameters passed:
-            * If the k parameter is supplied then the q vectors are grouped into k clusters using kmeans.
-            * If the threshold parameter is supplied then the q vectors a split into groups based on cophenetic distance.
-
-        :param qs: list of q vectors to cluster. Each element should be a numpy array of length three.
-        :param k: number of clusters to use (optional).
-        :param threshold: cophenetic distance cut off point for new clusters (optional)
-        :returns: tuple (clusters, k)
-            Where:
-                list -- clusters is a list of cluster indicies which each q belongs to
-                int -- k is the number of clusters used
-        """
-        if k is not None:
-            centroid, clusters = kmeans2(qs, k)
-        else:
-            clusters = hcluster.fclusterdata(qs, threshold, criterion="distance")
-        return clusters, len(set(clusters))
-
-    def average_clusters(self, qs, clusters):
-        """Find the centroid of the clusters.
-
-        For each q vector, group them by their designated cluster and then compute
-        the average of the group.
-
-        :param qs: list of q vectors. Each element should be a numpy array.
-        :param clusters: the indicies of the cluster that each q belongs to.
-        :returns: np.ndarry -- the list of centroids for each cluster.
-        """
-
-        averaged_qs = []
-        for cluster_index in set(clusters):
-            averaged_qs.append(np.mean(qs[clusters==cluster_index], axis=0))
-        return np.array(averaged_qs)
-
-    def find_q_vectors(self, nuclear_hkls, sats_hkls):
-        """Find the q vector between the nuclear HKLs and the satellite peaks
-
-        Given a list of HKL positions and a list of fractional HKL positions of
-        satellite peaks, find the difference between each satellite and its nearest
-        integer HKL.
-
-        :param nuclear_hkls: the positions of integer HKL peaks.
-        :param sats_hkl: the positions of fractional "satellite" HKL peaks.
-        :returns: np.ndarray -- array of q vectors.
-        """
-        peak_map = KDTree(nuclear_hkls)
-        qs = []
-        for sat in sats_hkls:
-            distance, index = peak_map.query(sat, k=1)
-            if distance > 2:
-                #peak to far away from satellite ignore
-                continue
-            nearest_peak = nuclear_hkls[index]
-            qs.append(sat - nearest_peak)
-        return np.array(qs)
-
     def create_fractional_peaks_workspace(self, qs, nuclear):
         """Generate a peaks workspace of possible satellite peaks from a list of q vectors.
 
@@ -179,4 +117,4 @@ class RefineSatellitePeaks(DataProcessorAlgorithm):
         return predicted_satellites
 
 
-AlgorithmFactory.subscribe(RefineSatellitePeaks)
+AlgorithmFactory.subscribe(FindSatellitePeaks)

Framework/PythonInterface/plugins/algorithms/IndexSatellitePeaks.py

+import numpy as np
+from scipy.spatial import KDTree
+import fractional_indexing as indexing
+
+from mantid.kernel import Direction
+from mantid.api import (IPeaksWorkspaceProperty,
+                        ITableWorkspaceProperty, PythonAlgorithm, AlgorithmFactory)
+import mantid.simpleapi as api
+
+
+class IndexSatellitePeaks(PythonAlgorithm):
+
+    def category(self):
+        return 'Crystal\\Peaks'
+
+    def seeAlso(self):
+        return [ "FindSatellitePeaks" ]
+
+    def name(self):
+        return "IndexSatellitePeaks"
+
+    def summary(self):
+        return "Algorithm for indexing satellite peaks in superspace"
+
+    def PyInit(self):
+        self.declareProperty(IPeaksWorkspaceProperty(name="NuclearPeaks",
+                                                     defaultValue="",
+                                                     direction=Direction.Input),
+                             doc="Main integer HKL peaks. Q vectors will be calculated relative to these peaks.")
+
+        self.declareProperty(IPeaksWorkspaceProperty(name="SatellitePeaks",
+                                                     defaultValue="",
+                                                     direction=Direction.Input),
+                             doc="Positions of satellite peaks with fractional \
+                             HKL coordinates")
+
+        self.declareProperty(ITableWorkspaceProperty(name="OutputWorkspace",
+                                                     defaultValue="",
+                                                     direction=Direction.Output),
+                             doc="The indexed satellite peaks. This will be a  \
+                             table workspace with miller indicies h, k, l, m1, \
+                             m2, ..., mn.")
+
+        self.declareProperty('Tolerance', 0.3, direction=Direction.Input,
+                             doc="Tolerance on the noise of the q vectors")
+
+        self.declareProperty('NumOfQs', 1, direction=Direction.Input,
+                             doc="Number of independant q vectors")
+
+        self.declareProperty('ClusterThreshold', 1.5, direction=Direction.Input,
+                             doc="Threshold for automaticallty deciding on the number of q vectors to use. If NumOfQs found is set then this \
+                             is property is ignored.")
+
+    def PyExec(self):
+        tolerance = self.getProperty("Tolerance").value
+        k = int(self.getProperty("NumOfQs").value)
+        nuclear = self.getProperty("NuclearPeaks").value
+        satellites = self.getProperty("SatellitePeaks").value
+        cluster_threshold = self.getProperty("ClusterThreshold").value
+        n_trunc_decimals = int(np.ceil(abs(np.log10(tolerance))))
+
+        if nuclear.getNumberPeaks() == 0:
+            raise RuntimeError("The NuclearPeaks parameter must have at least one peak")
+
+        if satellites.getNumberPeaks() == 0:
+            raise RuntimeError("The SatellitePeaks parameter must have at least one peak")
+
+        nuclear_hkls = indexing.get_hkls(nuclear)
+        sats_hkls = indexing.get_hkls(satellites)
+
+        qs = indexing.find_q_vectors(nuclear_hkls, sats_hkls)
+        self.log().notice("K value is {}".format(k))
+
+        k = None if k == -1 else k
+        clusters, k = indexing.cluster_qs(qs, k=k, threshold=cluster_threshold)
+
+        qs = indexing.average_clusters(qs, clusters)
+        qs = indexing.trunc_decimals(qs, n_trunc_decimals)
+        qs = indexing.sort_vectors_by_norm(qs)
+
+        self.log().notice("Q vectors are: \n{}".format(qs))
+
+        indices = indexing.index_q_vectors(qs, tolerance)
+        ndim = indices.shape[1] + 3
+
+        hkls = indexing.find_nearest_integer_peaks(nuclear_hkls, sats_hkls)
+
+        hklm = np.zeros((hkls.shape[0], ndim))
+        hklm[:, :3] = np.round(hkls)
+
+        raw_qs = hkls - sats_hkls
+        peak_map = KDTree(qs)
+        for i, q in enumerate(raw_qs):
+            distance, index = peak_map.query(q, k=1)
+            hklm[i, 3:] = indices[index]
+
+        indexed = self.create_indexed_workspace(satellites, ndim, hklm)
+        self.setProperty("OutputWorkspace", indexed)
+
+    def create_indexed_workspace(self, fractional_peaks, ndim, hklm):
+        """Create a TableWorkepace that contains indexed peak data.
+
+        This produces a TableWorkepace that looks like a PeaksWorkspace but
+        with the additional index columns included. In future releases support
+        for indexing should be added to the PeaksWorkspace data type itself.
+
+        :param fractional_peaks: the peaks workspace containing peaks with
+            fractional HKL values.
+        :param ndim: the number of additional indexing columns to add.
+        :param hklm: the new higher dimensional miller indicies to add.
+        :returns: a table workspace with the indexed peak data
+        """
+        # Create table with the number of columns we need
+        types = ['int', 'long64', 'double', 'double', 'double', 'double',  'double', 'double',
+                 'double', 'double', 'double', 'float', 'str', 'float', 'float', 'V3D', 'V3D']
+        name = self.getPropertyValue("OutputWorkspace")
+        indexed = api.CreateEmptyTableWorkspace(name)
+        names = fractional_peaks.getColumnNames()
+
+        # Insert the extra columns for the addtional indicies
+        for i in range(ndim - 3):
+            names.insert(5 + i, 'm{}'.format(i + 1))
+            types.insert(5 + i, 'double')
+
+        names = np.array(names)
+        types = np.array(types)
+
+        # Create columns in the table workspace
+        for name, column_type in zip(names, types):
+            indexed.addColumn(column_type, name)
+
+        # Copy all columns from original workspace, ignoring HKLs
+        column_data = []
+        idx = np.arange(0, names.size)
+        hkl_mask = (idx < 2) | (idx > 4 + (ndim - 3))
+        for name in names[hkl_mask]:
+            column_data.append(fractional_peaks.column(name))
+
+        # Insert the addtional HKL columns into the data
+        for i, col in enumerate(hklm.T.tolist()):
+            column_data.insert(i + 2, col)
+
+        # Insert the columns into the table workspace
+        for i in range(fractional_peaks.rowCount()):
+            row = [column_data[j][i] for j in range(indexed.columnCount())]
+            indexed.addRow(row)
+
+        return indexed
+
+
+AlgorithmFactory.subscribe(IndexSatellitePeaks)

Framework/PythonInterface/plugins/algorithms/fractional_indexing.py

+import itertools
+import numpy as np
+import scipy.cluster.hierarchy as hcluster
+from scipy.cluster.vq import kmeans2
+from scipy.spatial import KDTree
+
+_MAX_REFLECTIONS = 10
+
+
+def find_bases(qs, tolerance):
+    """Find bases for the higher dimensional indexing scheme
+
+    This will attempt to find a set of basis vectors for the
+    higher dimensional indexing scheme.
+
+    :param qs: list of q vectors in HKL space defining satellite peak offsets
+    :param tolerance: the tolerance how close each vector must be to be indexed
+                      using the generated scheme
+    :return: tuple of (number of dimensions, nx3 matrix of basis vectors)
+    """
+    qs = list(sort_vectors_by_norm(qs))
+    final_qs = [qs.pop(0)]
+    refs, hkls = generate_hkl_grid(np.array(final_qs), _MAX_REFLECTIONS)
+
+    for q in qs:
+        regenerate_grid = not is_indexed(q, refs, tolerance)
+
+        if regenerate_grid:
+            final_qs.append(q)
+            refs, hkls = generate_hkl_grid(np.array(final_qs), _MAX_REFLECTIONS)
+
+    return len(final_qs), np.vstack(final_qs)
+
+
+def generate_hkl_grid(bases, upper_bound):
+    """Generate a grid of reflections and hkl values
+
+    :param bases: the set of bases to generate the reflections and HKL grid.
+    :param upper_bound: the upper bound HKL index to generate
+    :return tuple of (nx3 matrix of reflections, nx3 matrix of HKL indices)
+    """
+    search_range = np.arange(-upper_bound, upper_bound)
+    ndim = len(bases)
+    hkls = np.array(list(itertools.product(*[search_range] * ndim)))
+    reflections = np.array([np.dot(bases.T, hkl) for hkl in hkls])
+    return reflections, hkls
+
+
+def find_nearest_hkl(qs, reflections, hkls, tolerance):
+    """Find the nearest HKL value that indexes each of the q vectors
+
+    :param qs: the q vectors for satellite peaks to index
+    :param reflections: the list of reflections to search
+    :param hkls: the list of hkl values mapping to each reflection
+    :param tolerance: the tolerance on the difference between a given q and
+                      the nearest reflection
+    :return: MxN matrix of integers indexing the q vectors. Unindexed vectors
+             will be returned as all zeros.
+    """
+    kdtree = KDTree(reflections)
+
+    final_indexing = []
+    for i, q in enumerate(qs):
+        result = kdtree.query_ball_point(q, r=tolerance)
+        if len(result) > 0:
+            smallest = sort_vectors_by_norm(hkls[result])
+            final_indexing.append(smallest[0])
+        else:
+            final_indexing.append(np.zeros_like(hkls.shape[1]))
+
+    return np.vstack(final_indexing)
+
+
+def is_indexed(q, reflections, tolerance):
+    """Check if a q vector is indexed by a set of reflections
+
+    :param q: the q vector to check
+    :param reflections: the list of reflections to find if this matches
+    :param tolerance: the tolerance on the difference between q and the nearest
+                      reflection.
+    :return: whether this q is indexed by this list of reflections
+    """
+    kdtree = KDTree(reflections)
+    result = kdtree.query_ball_point(q, r=tolerance)
+    return len(result) > 0
+
+
+def unique_rows(a):
+    """Return the unique rows of a numpy array
+
+    In numpy >= 1.13 np.unique(x, axis=0) can be user instead. But for versions of
+    numpy before 1.13 this function can be used to achieve the same result.
+
+    :param a: the array to find unique rows for.
+    :returns: the unique rows of a
+    """
+    a = np.ascontiguousarray(a)
+    unique_a = np.unique(a.view([('', a.dtype)] * a.shape[1]))
+    return unique_a.view(a.dtype).reshape((unique_a.shape[0], a.shape[1]))
+
+
+def index_q_vectors(qs, tolerance=.03):
+    """Find the n-dimensional index of a collection of q vectors
+
+    This function will automatically attempt infer the correct number of
+    additional dimensions required to index the given list of q vectors.
+
+    :param qs: the q vectors to find indicies for
+    :param tolerance: the tolerance on whether two equivilent or multiple qs
+        should be classifed as the same.
+    :return: ndarray of indicies for each q vector
+    """
+
+    ndim, bases = find_bases(qs, tolerance)
+    refs, hkls = generate_hkl_grid(bases, _MAX_REFLECTIONS)
+    return find_nearest_hkl(qs, refs, hkls, tolerance)
+
+
+def cluster_qs(qs, k=None, threshold=1.5):
+    """Cluster q vectors into discrete groups.
+
+    Classifies each of the q vectors into a number of clusters. The number of clusters used is decided by the parameters passed:
+        * If the k parameter is supplied then the q vectors are grouped into k clusters using kmeans.
+        * If the threshold parameter is supplied then the q vectors a split into groups based on cophenetic distance.
+
+    :param qs: list of q vectors to cluster. Each element should be a numpy array of length three.
+    :param k: number of clusters to use (optional).
+    :param threshold: cophenetic distance cut off point for new clusters (optional)
+    :returns: tuple (clusters, k)
+        Where:
+            list -- clusters is a list of cluster indicies which each q belongs to
+            int -- k is the number of clusters used
+    """
+    if k is not None:
+        centroids = kmeans_plus_plus(qs, k)
+        _, clusters = kmeans2(qs, centroids, minit='matrix')
+        if len(set(clusters)) != k:
+            raise ValueError("Could not group the satellite reflections "
+                             "into {} clusters. Please check that you have "
+                             "at least {} satellites.".format(k,k))
+    else:
+        clusters = hcluster.fclusterdata(qs, threshold, criterion="distance")
+    return clusters, len(set(clusters))
+
+
+def find_q_vectors(nuclear_hkls, sats_hkls):
+    """Find the q vector between the nuclear HKLs and the satellite peaks
+
+    Given a list of HKL positions and a list of fractional HKL positions of
+    satellite peaks, find the difference between each satellite and its nearest
+    integer HKL.
+
+    :param nuclear_hkls: the positions of integer HKL peaks.
+    :param sats_hkls: the positions of fractional "satellite" HKL peaks.
+    :returns: np.ndarray -- array of q vectors.
+    """
+    peak_map = KDTree(nuclear_hkls)
+    qs = []
+    for sat in sats_hkls:
+        distance, index = peak_map.query(sat, k=1)
+        if distance > 2:
+            # peak too far away from satellite ignore
+            continue
+        nearest_peak = nuclear_hkls[index]
+        qs.append(sat - nearest_peak)
+    return np.array(qs)
+
+
+def find_nearest_integer_peaks(nuclear_hkls, sat_hkls):
+    """Find the nearest integer peak for each fractional peak
+
+    This will perform a spatial search to find the intger peak which is nearest the fractional satellite peak.
+
+    :param nuclear_hkls: the hkl poistions of each of the nuclear peaks.
+    :param sat_hkls: the fractional hkl position of the satellite peaks.
+    :returns: np.ndarray -- 2D array of HKL integer values for each fractional peak
+    """
+    peak_map = KDTree(nuclear_hkls)
+    hkls = []
+    for sat in sat_hkls:
+        distance, index = peak_map.query(sat, k=1)
+        nearest_peak = nuclear_hkls[index]
+        hkls.append(nearest_peak)
+    return np.array(hkls)
+
+
+def average_clusters(qs, clusters):
+    """Compute the centroid of each cluster of vectors
+
+    :param qs: the q vectors to find the centroids of
+    :param clusters: the cluster index that each q belongs to
+    :return: an ndarray of cluster centroids
+    """
+    averaged_qs = []
+    for cluster_index in set(clusters):
+        averaged_qs.append(np.mean(qs[clusters == cluster_index], axis=0))
+    return np.array(averaged_qs)
+
+
+def get_hkls(peaks_workspace):
+    """Get the H, K, and L columns from a PeaksWorkspace.
+
+    :param peaks_workspace: the peaks workspace to extract HKL values from
+    :return: 2D numpy array of HKL values.
+    """
+    return np.array([np.array([peak['h'], peak['k'], peak['l']])
+                     for peak in peaks_workspace])
+
+
+def remove_noninteger(matrix):
+    """Remove any non integer values from a matrix
+
+    :param matrix: the matrix to remove non integer values from
+    :return: matrix with non integer elements set to zero
+    """
+    matrix[np.absolute(np.mod(matrix, 1)) > 1e-14] = 0
+    return matrix
+
+
+def trunc_decimals(vec, n_decimals=2):
+    """Truncate the elements of a vector below n decimals places
+
+    :param vec: the vector to truncate the elements of
+    :param n_decimals: the number of decimal places to truncate from
+    :return: the vector with truncated elements
+    """
+    decade = 10**n_decimals
+    return np.trunc(vec*decade)/decade
+
+
+def sort_vectors_by_norm(vecs):
+    """Sort a list of row vectors by the Euclidean norm
+    :param vecs: vectors to sort according to their norms
+    :returns: ndarray of sorted row vectors
+    """
+    idx = np.argsort(norm_along_axis(np.abs(vecs)))
+    vecs = vecs[idx]
+    return vecs
+
+
+def norm_along_axis(vecs, axis=-1):
+    """Compute the euclidean norm along the rows of a matrix
+
+    :param vecs: vectors to compute the Euclidean norm for
+    :param axis: the axis to compute the norm along
+    :return: a ndarray with the norms along the chosen axis
+    """
+    return np.sum(vecs**2, axis=axis)**(1./2)
+
+
+def kmeans_plus_plus(points, k):
+    """Generate centroids for kmeans using the kmeans++ algorithm
+
+    Select a good set of starting values for kmeans by weighting the choice of
+    the next centroid by the distance to the existing centroids.
+
+    This algorithm roughly works as follows:
+        1) The first centroid is selected uniformly from the data.
+        2) The squared euclidean distance between all points to all centroids is calculated.
+        3) Take the minimum computed distance for each point over all centroids.
+        4) Normalise the list of distances
+        5) Choose a new centroid weighted by the computed distances
+
+    :param points: the points to select centroids from
+    :param k: the desired number of centroids
+    :return: Kx3 matrix of centroids
+    """
+    if points.shape[0] < k:
+        raise RuntimeError("k is greater than the number of points! Please choose a smaller k value")
+
+    centroid_indices = []
+    centroids = []
+
+    # pick the first point uniformly at random from the data
+    indices = np.indices(points.shape)[0][:, 0]
+    centroid_index = np.random.choice(indices)
+    centroid = points[centroid_index]
+
+    centroid_indices.append(centroid_index)
+    centroids.append(centroid)
+
+    for i in range(k-1):
+        # choose all points that are not already centroids
+        mask = np.array([x for x in indices if x not in centroid_indices])
+        pts = points[mask]
+
+        # calculate probability distribution for being picked based on squared
+        # distance from each centroid
+        distance_squared = np.array([norm_along_axis(pts - c, axis=1)**2 for c in centroids])
+        distance_squared = np.min(distance_squared, axis=0)
+        distance_squared /= np.sum(distance_squared)
+
+        # choose a new random centroid weighted by how far it is from the other
+        # centroids
+        centroid_index = np.random.choice(mask, p=distance_squared)
+        centroid_index = indices[centroid_index]
+        centroid = points[centroid_index]
+
+        centroid_indices.append(centroid_index)
+        centroids.append(centroid)
+
+    return np.array(centroids)

Framework/PythonInterface/test/python/plugins/algorithms/CMakeLists.txt

@@ -45,9 +45,11 @@ set ( TEST_PY_FILES
   FindEPPTest.py
   FindReflectometryLinesTest.py
   FitGaussianTest.py
+  FractionalIndexingTest.py
   GetEiT0atSNSTest.py
   GetNegMuMuonicXRDTest.py
   IndirectTransmissionTest.py
+  IndexSatellitePeaksTest.py
   LoadAndMergeTest.py
   LoadDNSLegacyTest.py
   LoadEmptyVesuvioTest.py

Framework/PythonInterface/test/python/plugins/algorithms/FractionalIndexingTest.py

+import unittest
+import numpy as np
+import numpy.testing as npt
+import fractional_indexing as indexing
+
+from mantid.geometry import UnitCell
+from mantid.simpleapi import CreatePeaksWorkspace, CreateSimulationWorkspace
+
+
+class FractionIndexingTests(unittest.TestCase):
+
+    def setUp(self):
+        # Need to set the random seed because the scipy kmeans algorithm
+        # randomly initilizes the starting centroids. This can lead to a
+        # different but equivilent indexing.
+        np.random.seed(10)
+
+    def test_find_bases_with_1d_modulation(self):
+        qs = np.array([
+            [0, 0, .13],
+        ])
+
+        ndim, bases = indexing.find_bases(qs, tolerance=.02)
+
+        self.assertEqual(ndim, 1)
+        npt.assert_array_equal(bases[0], qs[0])
+
+    def test_find_bases_with_2d_modulation(self):
+        qs = np.array([
+            [0, 0, .1],
+            [0, .1, 0],
+        ])
+
+        ndim, bases = indexing.find_bases(qs, tolerance=.02)
+
+        self.assertEqual(ndim, 2)
+        npt.assert_array_equal(bases[0], qs[0])
+        npt.assert_array_equal(bases[1], qs[1])
+
+    def test_find_bases_with_3d_modulation(self):
+        qs = np.array([
+            [0, 0, .1],
+            [0, .1, 0],
+            [.1, 0, 0],
+        ])
+
+        ndim, bases = indexing.find_bases(qs, tolerance=.02)
+
+        self.assertEqual(ndim, 3)
+        npt.assert_array_equal(bases[0], qs[0])
+        npt.assert_array_equal(bases[1], qs[1])
+        npt.assert_array_equal(bases[2], qs[2])
+
+    def test_find_bases_with_4d_modulation(self):
+        qs = np.array([
+            [0, 0, .1],
+            [0, .1, 0],
+            [.1, 0, 0],
+            [.15, 0, 0],
+        ])
+
+        ndim, bases = indexing.find_bases(qs, tolerance=.02)
+
+        self.assertEqual(ndim, 4)
+        npt.assert_array_equal(bases[0], qs[0])
+        npt.assert_array_equal(bases[1], qs[1])
+        npt.assert_array_equal(bases[2], qs[2])
+
+    def test_find_bases_with_2d_modulation_with_linear_combination(self):
+        qs = np.array([
+            [0, .1, .1],
+            [0, 0, .1],
+            [0, .1, 0],
+            [0, 0, .15],
+            [0, .1, .2],
+        ])
+
+        ndim, bases = indexing.find_bases(qs, tolerance=.02)
+
+        self.assertEqual(ndim, 3)
+        npt.assert_array_equal(bases[0], np.array([0, 0, .1]))
+        npt.assert_array_equal(bases[1], np.array([0, .1, 0]))
+        npt.assert_array_equal(bases[2], np.array([0, 0, .15]))
+
+    def test_find_bases_with_non_orthogonal_cell(self):
+        cell = UnitCell(1, 1, 1, 90, 95, 103)
+        cart_cell = cell.getB()
+        qs = np.array([
+            [.5, 0, 0],
+            [0, .5, 0],
+            [0, 0, .5],
+            [0, 0, .25],
+        ])
+
+        qs = np.dot(qs, cart_cell)
+
+        ndim, bases = indexing.find_bases(qs, 1e-5)
+        self.assertEqual(ndim, 3, "Number of dimensions must be 3")
+
+        expected_bases = np.array([
+            [0., 0., 0.25],
+            [0., .5, 0.],
+            [0.51521732, 0.11589868, 0.04490415]
+        ])
+
+        npt.assert_almost_equal(bases, expected_bases, err_msg="Basis vectors do not match")
+
+    def test_find_indexing_with_non_orthogonal_cell(self):
+        cell = UnitCell(1, 1, 1, 90, 95, 103)
+        cart_cell = cell.getB()
+        qs = np.array([
+            [.5, 0, 0],
+            [0, .5, 0],
+            [0, 0, .5],
+            [0, 0, .25],
+        ])
+
+        qs = np.dot(qs, cart_cell)
+
+        indices = indexing.index_q_vectors(qs)
+
+        expected_indexing = np.array([
+            [0, 0, 1],
+            [0, 1, 0],
+            [2, 0, 0],
+            [1, 0, 0]
+        ])
+
+        npt.assert_equal(indices, expected_indexing, err_msg="Indexing does not match expected.")
+
+    def test_index_with_modulation(self):
+        qs = np.array([
+            [0, .1, .1],
+            [0, -.1, 0],
+            [0, 0, .1],
+            [0, 0, -.1],
+            [0, .1, .2],
+            [0, 0, .45],
+        ])
+
+        expected_indexing = np.array([
+            [-1,  1,  0],
+            [ 1,  0,  0],
+            [ 0,  1,  0],
+            [ 0, -1,  0],
+            [-1,  2,  0],
+            [ 0,  0,  1]
+        ])
+
+        actual_indexing = indexing.index_q_vectors(qs, tolerance=.03)
+        npt.assert_array_equal(actual_indexing, expected_indexing)
+
+    def test_norm_along_axis(self):
+        vecs = np.array([
+            [0, 0, 3],
+            [0, 2, 0],
+            [1, 0, 0],
+            [0, 5, 5],
+            [0, 4, 4],
+        ])
+
+        expected_output = np.array([3, 2, 1, 7.071068, 5.656854])
+        norms = indexing.norm_along_axis(vecs)
+        npt.assert_allclose(norms, expected_output)
+
+    def test_sort_vectors_by_norm(self):
+        vecs = np.array([
+            [0, 0, 3],
+            [0, 2, 0],
+            [1, 0, 0],
+            [0, 5, 5],
+            [0, 4, 4],
+        ])
+
+        expected_output = np.array([
+            [1, 0, 0],
+            [0, 2, 0],
+            [0, 0, 3],
+            [0, 4, 4],
+            [0, 5, 5],
+        ])
+
+        vecs_sorted = indexing.sort_vectors_by_norm(vecs)
+        npt.assert_allclose(vecs_sorted, expected_output)
+
+    def test_trunc_decimals(self):
+        reference = np.array([0, 0, 0.1])
+        test_input = reference + np.random.random(3) * 0.1
+        result = indexing.trunc_decimals(test_input, 1)
+        npt.assert_array_equal(result, reference)
+
+    def test_remove_noninteger(self):
+        test_input = np.array([1, 2, 3, 1.1, 3.4, -10, -1.8])
+        reference = np.array([1, 2, 3, 0, 0, -10, 0])
+        result = indexing.remove_noninteger(test_input)
+        npt.assert_array_equal(result, reference)
+
+    def test_get_hkls(self):
+        ws = CreateSimulationWorkspace("IRIS", BinParams="1,5,10")
+        peaks = CreatePeaksWorkspace(ws, 2)
+        reference = np.array([
+            [1, 1, 2],
+            [2, 1, 4],
+        ])
+
+        peak = peaks.getPeak(0)
+        peak.setHKL(1, 1, 2)
+        peak = peaks.getPeak(1)
+        peak.setHKL(2, 1, 4)
+
+        hkl = indexing.get_hkls(peaks)
+        npt.assert_array_equal(hkl, reference)
+
+    def test_cluster_qs_with_fixed_k(self):
+        qs = np.array([
+            [0, .1, .1],
+            [0, .1, .1],
+            [0, .0, .1],
+            [0, .0, .1],
+            [0, .1, .1],
+        ])
+
+        qs += np.random.random(qs.shape) * 0.01
+
+        k = 2
+        clusters, k = indexing.cluster_qs(qs, k)
+        self.assertEqual(k, 2)
+        npt.assert_array_equal(clusters, np.array([0, 0, 1, 1, 0]))
+
+    def test_cluster_qs_with_auto_k(self):
+        qs = np.array([
+            [0, .1, .1],
+            [0, .1, .1],
+            [0, .0, .1],
+            [0, .0, .1],
+            [0, .1, .1],
+        ])
+
+        qs += np.random.random(qs.shape) * 0.01
+
+        clusters, k = indexing.cluster_qs(qs, threshold=0.01)
+        self.assertEqual(k, 2)
+        npt.assert_array_equal(clusters, np.array([2, 2, 1, 1, 2]))
+
+if __name__ == "__main__":
+    unittest.main()

Framework/PythonInterface/test/python/plugins/algorithms/IndexSatellitePeaksTest.py

+from mantid.simpleapi import mtd, Load, IndexSatellitePeaks
+import unittest
+import numpy as np
+import numpy.testing as npt
+
+
+class IndexSatellitePeaksTest(unittest.TestCase):
+
+    def setUp(self):
+        # Need to set the random seed because the scipy kmeans algorithm
+        # randomly initilizes the starting centroids.
+        np.random.seed(100)
+        self._nuclear_peaks = Load(
+            "WISH_peak_hkl_small.nxs", OutputWorkspace="nuclear_peaks")
+        self._peaks = Load(
+            "refine_satellites_fixed_q_test.nxs", OutputWorkspace="peaks")
+
+    def tearDown(self):
+        mtd.clear()
+
+    def test_exec_with_cluster_threshold(self):
+        expected_values = np.array([-1, -1, 1, 1])
+        indexed_peaks = IndexSatellitePeaks(self._nuclear_peaks, self._peaks,
+                                            Tolerance=0.1,
+                                            ClusterThreshold=1.5, NumOfQs=-1)
+        index_values = np.array(indexed_peaks.column("m1"))
+
+        npt.assert_array_equal(index_values, expected_values)
+        self.assertRaises(RuntimeError, indexed_peaks.column, "m2")
+
+    def test_exec_with_number_of_qs(self):
+        expected_values = np.array([-1, -1, 1, 1])
+        indexed_peaks = IndexSatellitePeaks(self._nuclear_peaks, self._peaks,
+                                            Tolerance=0.1, NumOfQs=2)
+        index_values = np.array(indexed_peaks.column("m1"))
+
+        npt.assert_array_equal(index_values, expected_values)
+        self.assertRaises(RuntimeError, indexed_peaks.column, "m2")
+
+
+if __name__ == "__main__":
+    unittest.main()

Testing/SystemTests/scripts/CMakeLists.txt

@@ -17,6 +17,8 @@ list(APPEND DATA_DIRS ${SYSTEM_TEST_DATA_DIR}/LARMOR)
 list(APPEND DATA_DIRS ${SYSTEM_TEST_DATA_DIR}/INTER)
 # Reference results
 list(APPEND DATA_DIRS ${ExternalData_BINARY_ROOT}/Testing/SystemTests/tests/analysis/reference)
+# Unit tests
+list(APPEND DATA_DIRS ${ExternalData_BINARY_ROOT}/Testing/Data/UnitTest)
 # Doc tests
 list(APPEND DATA_DIRS ${ExternalData_BINARY_ROOT}/Testing/Data/DocTest)
 

Testing/SystemTests/tests/analysis/RefineSatellitePeaksTest.py → Testing/SystemTests/tests/analysis/FindSatellitePeaksTest.py

-from mantid.simpleapi import RefineSatellitePeaks, Load
+from mantid.simpleapi import FindSatellitePeaks, Load
 import unittest
 import stresstesting
 
 
 def load_files():
-    md_workspace = Load(Filename="WISH_md_small.nxs", OutputWorkspace='md_workspace')
+    md_workspace = Load(Filename="WISH_md_small.nxs",
+                        OutputWorkspace='md_workspace')
     main_peaks = Load(Filename="WISH_peak_hkl_small.nxs",
                       OutputWorkspace='main_peaks')
     satellite_peaks = Load(Filename="WISH_peak_hkl_frac_small.nxs",
@@ -12,7 +13,7 @@ def load_files():
     return md_workspace, main_peaks, satellite_peaks
 
 
-class RefineSatellitePeaksTestFixedNumQ(stresstesting.MantidStressTest):
+class FindSatellitePeaksTestFixedNumQ(stresstesting.MantidStressTest):
 
     def requiredFiles(self):
         return ["WISH_md_small.nxs", "WISH_peak_hkl_small.nxs", "WISH_peak_hkl_frac_small.nxs"]
@@ -29,14 +30,14 @@ class RefineSatellitePeaksTestFixedNumQ(stresstesting.MantidStressTest):
 
         k = 2
 
-        self._satellites_refined = RefineSatellitePeaks(NuclearPeaks=main_peaks, SatellitePeaks=satellite_peaks, MDWorkspace=md_workspace,
-                                                        NumOfQs=k, **fixed_params)
+        self._satellites_refined = FindSatellitePeaks(NuclearPeaks=main_peaks, SatellitePeaks=satellite_peaks, MDWorkspace=md_workspace,
+                                                      NumOfQs=k, **fixed_params)
 
     def validate(self):
-        return self._satellites_refined.name(),'refine_satellites_fixed_q_test.nxs'
+        return self._satellites_refined.name(), 'refine_satellites_fixed_q_test.nxs'
 
 
-class RefineSatellitePeaksTestAutoFindQ(stresstesting.MantidStressTest):
+class FindSatellitePeaksTestAutoFindQ(stresstesting.MantidStressTest):
 
     def requiredFiles(self):
         return ["WISH_md_small.nxs", "WISH_peak_hkl_small.nxs", "WISH_peak_hkl_frac_small.nxs"]
@@ -52,11 +53,11 @@ class RefineSatellitePeaksTestAutoFindQ(stresstesting.MantidStressTest):
         }
 
         threshold = 1.0
-        self._satellites_refined = RefineSatellitePeaks(NuclearPeaks=main_peaks, SatellitePeaks=satellite_peaks, MDWorkspace=md_workspace,
-                                                        ClusterThreshold=threshold, **fixed_params)
+        self._satellites_refined = FindSatellitePeaks(NuclearPeaks=main_peaks, SatellitePeaks=satellite_peaks, MDWorkspace=md_workspace,
+                                                      ClusterThreshold=threshold, **fixed_params)
 
     def validate(self):
-        return self._satellites_refined.name(),'refine_satellites_auto_q_test.nxs'
+        return self._satellites_refined.name(), 'refine_satellites_auto_q_test.nxs'
 
 
 if __name__ == "__main__":

docs/source/algorithms/RefineSatellitePeaks-v1.rst → docs/source/algorithms/FindSatellitePeaks-v1.rst

@@ -7,17 +7,17 @@
 Description
 -----------
 
-RefineSatellitePeaks can be used to refine the locations of "satellite" peaks
-occurring at fractional HKL locations in reciprocal space. RefineSatellitePeaks
+FindSatellitePeaks can be used to refine the locations of "satellite" peaks
+occurring at fractional HKL locations in reciprocal space. FindSatellitePeaks
 takes a :ref:`MDWorkspace <MDWorkspace>` of experimental data, a
 :ref:`PeaksWorkspace <PeaksWorkspace>` containing the locations of peaks with
 integer HKL, and another PeaksWorkspace containing a subset of peaks found at
 fractional HKL. 
 
-.. figure:: /images/RefineSatellitePeaks-satellites.png
+.. figure:: ../images/FindSatellitePeaks-satellites.png
    :align: center
    :width: 600px
-   :alt: RefineSatellitePeaks-satellites.png
+   :alt: FindSatellitePeaks-satellites.png
 
    Satellite peaks exist at fractional HKL coordinates. The vector `q` is the
    Euclidean distance from the nearest integer HKL peak to the satellite. A
@@ -33,10 +33,10 @@ this threshold.  The centroid of each cluster calculated for each group and is
 used as the offset to predict the location of fractional peaks everywhere in
 the :ref:`MDWorkspace <MDWorkspace>`.
 
-.. figure:: /images/RefineSatellitePeaks-clustering.png
+.. figure:: ../images/FindSatellitePeaks-clustering.png
    :align: center
    :width: 700px
-   :alt: RefineSatellitePeaks-clustering.png
+   :alt: FindSatellitePeaks-clustering.png
 
    Due to noise in the experimental data the `q` vector for every satellite may
    vary slightly. This algorithm calculates the `q` vectors for all peaks passed
@@ -54,10 +54,10 @@ Finally the found satellite peaks are integerated using the
 `BackgroundInnerRadius`, and `BackgroundOuterRadius`. Satellite peaks are
 discarded if there I/sigma value is less than the parameter `IOverSigma`.
 
-.. figure:: /images/RefineSatellitePeaks-centroids.png
+.. figure:: ../images/FindSatellitePeaks-centroids.png
    :align: center
    :width: 700px
-   :alt: RefineSatellitePeaks-centroid.png
+   :alt: FindSatellitePeaks-centroid.png
 
    Finally, using the predicted satellite peak positions (green) the local
    centroid is found and this is used as the true position of the experimental
@@ -102,15 +102,19 @@ Related Algorithms
 - :ref:`FilterPeaks <algm-FilterPeaks-v1>` is used to remove peaks with zero
   intensity or a I/sigma below the desired threshold.
 
+- :ref:`IndexSatellitePeaks <algm-IndexSatellitePeaks-v1>` can be used to
+  obtain a superspace indexing of the peaks workspace output from this
+  algorithm.
+
 Usage
 -----
 
 .. include:: ../usagedata-note.txt
 
 
-**Example - calling RefineSatellitePeaks:**
+**Example - calling FindSatellitePeaks:**
 
-.. testcode:: RefineSatellitePeaks
+.. testcode:: FindSatellitePeaks
 
     md_workspace = Load(Filename='WISH_md_small.nxs', OutputWorkspace='WISH_md_small')
     satellite_peaks = Load(Filename='WISH_peak_hkl_frac_small.nxs', OutputWorkspace='WISH_peak_hkl_frac_small')
@@ -125,11 +129,11 @@ Usage
     }
 
 
-    satellites_refined = RefineSatellitePeaks(NuclearPeaks=main_peaks, SatellitePeaks=satellite_peaks, MDWorkspace=md_workspace, **params)
+    satellites_refined = FindSatellitePeaks(NuclearPeaks=main_peaks, SatellitePeaks=satellite_peaks, MDWorkspace=md_workspace, **params)
 
     print (len(satellites_refined))
 
-.. testcleanup:: RefineSatellitePeaks
+.. testcleanup:: FindSatellitePeaks
 
    DeleteWorkspace('WISH_md_small')
    DeleteWorkspace('WISH_peak_hkl_small')
@@ -138,7 +142,7 @@ Usage
 
 Output:
 
-.. testoutput:: RefineSatellitePeaks
+.. testoutput:: FindSatellitePeaks
 
     4
 

docs/source/algorithms/IndexSatellitePeaks-v1.rst

+.. algorithm::
+
+.. summary::
+
+.. properties::
+
+Description
+-----------
+
+IndexSatellitePeaks is used to find the higher dimensional miller indices that
+index a collection of "satellite" fractional single crystal peaks. This
+algorithm takes a PeaksWorkspace, output from the algorithm
+FindSatellitePeaks, and returns a new table with integer indices for each
+peak.
+
+This algorithm works by first attempting to finding the minimum number of
+distinct modulation (`q`) vectors that are required to fully index the
+collection of peaks. The full list of q vectors is found using the same method
+described in the :ref:`FindSatellitePeaks <algm-FindSatellitePeaks-v1>`
+algorithm. The basis set is then chosen from this list. If there are multiple
+choices of vectors the algorithm will always choose the smallest one. This
+defines the number of additional dimensions required to index the crystal with
+integer indices. 
+
+Once a basis of `q` vectors has been chosen for indexing the system all integer
+multiples of the basis set are generated. Each `q` vector is then compared with
+the set of all possible integer reflections and is assigned to the closest one
+found within the radius of the `Tolerance` parameter (measured by Euclidean
+distance in HKL space). Reflections which cannot be indexed within the
+tolerance are set to (0,0,0).
+
+.. warning:: The current version of the algorithm returns a
+    :ref:`TableWorkspace <Table Workspaces>` and not a :ref:`PeaksWorkspace
+    <PeaksWorkspace>`. This means that the workspace cannot be overlaid on the
+    slice viewer or the instrument view. 
+
+.. seealso:: As well as being able to export the data to the nexus file format,
+    saving the data to the Jana format is supported via the :ref:`SaveReflections
+    <algm-SaveReflections-v1>` algorithm.
+
+For more information on superspace crystallography see:
+
+- Van Smaalen, Sander. "An elementary introduction to superspace 
+  crystallography." Zeitschrift für Kristallographie-Crystalline Materials 
+  219, no. 11 (2004): 681-691. 
+
+- Van Smaalen, Sander. "Incommensurate crystal structures." Crystallography 
+  Reviews 4, no. 2 (1995): 79-202. 
+
+Related Algorithms
+------------------
+
+- :ref:`FindSatellitePeaks <algm-FindSatellitePeaks-v1>` can be used to
+  obtain a workspace of satellite peak positions in fractional HKL that is
+  required for input to this algorithm.
+
+- :ref:`SaveReflections <algm-SaveReflections-v1>` algorithm can be used to save
+  the table workspace output from this algorithm to the Jana text file format.
+
+Usage
+-----
+
+.. include:: ../usagedata-note.txt
+
+
+**Example - calling IndexSatellitePeaks:**
+
+.. testcode:: IndexSatellitePeaks
+
+    import numpy as np
+    np.random.seed(1)
+
+    nuclear_peaks = Load('WISH_peak_hkl_small.nxs')
+    satellite_peaks = Load("refine_satellites_fixed_q_test.nxs")
+    indexed_peaks = IndexSatellitePeaks(nuclear_peaks, satellite_peaks, tolerance=0.1, NumOfQs=2)
+    index_values = np.array(indexed_peaks.column("m1"))
+
+    for peak in indexed_peaks:
+        print("H: {h:>5} K: {k:>5} L: {l:>5} M1: {m1:>5}".format(**peak))
+
+Output:
+
+.. testoutput:: IndexSatellitePeaks
+
+    H:   2.0 K:  -1.0 L:  -3.0 M1:  -1.0
+    H:   2.0 K:  -1.0 L:   0.0 M1:  -1.0
+    H:   2.0 K:  -1.0 L:  -3.0 M1:   1.0
+    H:   2.0 K:  -1.0 L:   0.0 M1:   1.0
+
+.. categories::
+
+.. sourcelink::
+  :cpp: None
+  :h: None

docs/source/release/v3.13.0/diffraction.rst

@@ -88,10 +88,12 @@ New
 
 - New algorithm :ref:`IntegratePeaksProfileFitting <algm-IntegratePeaksProfileFitting>` to integrate peaks using 3D profile fitting in reciprocal space.
 
-- New algorithm :ref:`RefineSatellitePeaks <algm-RefineSatellitePeaks>` to predict the location of fractional satellite peaks using a set of nuclear peaks and a set of seed satellite peaks.
+- New algorithm :ref:`FindSatellitePeaks <algm-FindSatellitePeaks>` to predict the location of fractional satellite peaks using a set of nuclear peaks and a set of seed satellite peaks.
 
 - New TOPAZ instrument geometry for 2018 run cycle
 
+- New algorithm :ref:`IndexSatellitePeaks <algm-IndexSatellitePeaks>` to index satellite peaks found using the :ref:`FindSatellitePeaks <algm-FindSatellitePeaks>` algorithm.
+
 Improvements
 ############