Merge pull request #459 from mantidproject/11418_poldi_create_peaks_from_cell_file

Algorithm to call PoldiCreatePeaksFromCell with parameters from a file

Code/Mantid/Framework/PythonInterface/plugins/algorithms/PoldiCreatePeaksFromFile.py

+# pylint: disable=no-init,invalid-name,too-few-public-methods
+from mantid.kernel import *
+from mantid.simpleapi import *
+from mantid.api import *
+from mantid.geometry import *
+
+from pyparsing import *
+
+import os
+
+
+class PoldiCompound(object):
+    """Small helper class to handle the results from PoldiCrystalFileParser."""
+    _name = ""
+    _spacegroup = ""
+    _atomString = ""
+    _cellDict = ""
+
+    def __init__(self, name, elements):
+        self._name = name
+
+        self.assign(elements)
+
+    def assign(self, elements):
+        for c in elements:
+            if c[0] == "atoms":
+                self._atomString = ';'.join(c[1:])
+            elif c[0] == "lattice":
+                cellNames = ['a', 'b', 'c', 'alpha', 'beta', 'gamma']
+                self._cellDict = dict(zip(cellNames, c[1:]))
+            elif c[0] == "spacegroup":
+                self._spacegroup = c[1]
+
+    def getAtomString(self):
+        return self._atomString
+
+    def getCellParameters(self):
+        return self._cellDict
+
+    def getSpaceGroup(self):
+        return self._spacegroup
+
+    def getName(self):
+        return self._name
+
+
+def raiseParseErrorException(message):
+    raise ParseException(message)
+
+
+class PoldiCrystalFileParser(object):
+    """Small parser for crystal structure files used at POLDI
+
+    This class encapsulates a small parser for crystal structure files that are used at
+    POLDI. The files contains information about the lattice, the space group and the basis (atoms
+    in the asymmetric unit).
+
+    The file format is defined as follows:
+
+        Compound_1 {
+            Lattice: [1 - 6 floats] => a, b, c, alpha, beta, gamma
+            Spacegroup: [valid space group symbol]
+            Atoms; {
+                Element x y z [occupancy [U_eq]]
+                Element x y z [occupancy [U_eq]]
+            }
+        }
+
+        Compound_2 {
+            ...
+        }
+
+    The parser returns a list of PoldiCompound objects with the compounds that were found
+    in the file. These are then processed by PoldiCreatePeaksFromFile to generate arguments
+    for calling PoldiCreatePeaksFromCell.
+    """
+    elementSymbol = Word(alphas, min=1, max=2).setFailAction(
+        lambda o, s, loc, token: raiseParseErrorException("Element symbol must be one or two characters."))
+    integerNumber = Word(nums)
+    decimalSeparator = Literal('.')
+    floatNumber = Combine(
+        integerNumber +
+        Optional(decimalSeparator + Optional(integerNumber))
+    )
+
+    whiteSpace = Suppress(White())
+
+    atomLine = Combine(
+        elementSymbol + whiteSpace +
+        delimitedList(floatNumber, delim=White()),
+        joinString=' '
+    )
+
+    keyValueSeparator = Suppress(Literal(":"))
+
+    groupOpener = Suppress(Literal('{'))
+    groupCloser = Suppress(Literal('}'))
+
+    atomsGroup = Group(CaselessLiteral("atoms") + keyValueSeparator +
+                       groupOpener + delimitedList(atomLine, delim=lineEnd) + groupCloser)
+
+    unitCell = Group(CaselessLiteral("lattice") + keyValueSeparator + delimitedList(
+        floatNumber, delim=White()))
+
+    spaceGroup = Group(CaselessLiteral("spacegroup") + keyValueSeparator + Word(
+        alphanums + "-" + ' '))
+
+    compoundContent = Each([atomsGroup, unitCell, spaceGroup]).setFailAction(
+        lambda o, s, loc, token: raiseParseErrorException(
+            "One of 'Lattice', 'SpaceGroup', 'Atoms' is missing or contains errors."))
+
+    compoundName = Word(alphanums + '_')
+
+    compound = Group(compoundName + Optional(whiteSpace) + \
+                     groupOpener + compoundContent + groupCloser)
+
+    comment = Suppress(Literal('#') + restOfLine)
+
+    compounds = Optional(comment) + OneOrMore(compound).ignore(comment) + stringEnd
+
+    def __call__(self, contentString):
+        parsedContent = None
+
+        if os.path.isfile(contentString):
+            parsedContent = self._parseFile(contentString)
+        else:
+            parsedContent = self._parseString(contentString)
+
+        return [PoldiCompound(x[0], x[1:]) for x in parsedContent]
+
+    def _parseFile(self, filename):
+        return self.compounds.parseFile(filename)
+
+    def _parseString(self, stringContent):
+        return self.compounds.parseString(stringContent)
+
+
+class PoldiCreatePeaksFromFile(PythonAlgorithm):
+    _parser = PoldiCrystalFileParser()
+
+    def category(self):
+        return "SINQ\\POLDI"
+
+    def name(self):
+        return "PoldiLoadCrystalData"
+
+    def summary(self):
+        return ("The algorithm reads a POLDI crystal structure file and creates a WorkspaceGroup that contains tables"
+                "with the expected reflections.")
+
+    def PyInit(self):
+        self.declareProperty(
+            FileProperty(name="InputFile",
+                         defaultValue="",
+                         action=FileAction.Load,
+                         extensions=["dat"]),
+            doc="A file with POLDI crystal data.")
+
+        self.declareProperty("LatticeSpacingMin", 0.5,
+                             direction=Direction.Input,
+                             doc="Lowest allowed lattice spacing.")
+
+        self.declareProperty("LatticeSpacingMax", 0.0,
+                             direction=Direction.Input,
+                             doc="Largest allowed lattice spacing.")
+
+        self.declareProperty(
+            WorkspaceProperty(name="OutputWorkspace",
+                              defaultValue="", direction=Direction.Output),
+            doc="WorkspaceGroup with reflection tables.")
+
+
+    def PyExec(self):
+        crystalFileName = self.getProperty("InputFile").value
+        try:
+            # Try parsing the supplied file using PoldiCrystalFileParser
+            compounds = self._parser(crystalFileName)
+
+            dMin = self.getProperty("LatticeSpacingMin").value
+            dMax = self.getProperty("LatticeSpacingMax").value
+
+            workspaces = []
+
+            # Go through found compounds and run "_createPeaksFromCell" for each of them
+            # If two compounds have the same name, a warning is written to the log.
+            for compound in compounds:
+                if compound.getName() in workspaces:
+                    self.log().warning("A compound with the name '" + compound.getName() + \
+                                       "' has already been created. Please check the file '" + crystalFileName + "'")
+                else:
+                    workspaces.append(self._createPeaksFromCell(compound, dMin, dMax))
+
+            self.setProperty("OutputWorkspace", GroupWorkspaces(workspaces))
+
+        # All parse errors are caught here and logged as errors
+        except ParseException as error:
+            errorString = "Could not parse input file '" + crystalFileName + "'.\n"
+            errorString += "The parser reported the following error:\n\t" + str(error)
+
+            self.log().error(errorString)
+
+
+    def _createPeaksFromCell(self, compound, dMin, dMax):
+        if not SpaceGroupFactory.isSubscribedSymbol(compound.getSpaceGroup()):
+            raise RuntimeError("SpaceGroup '" + compound.getSpaceGroup() + "' is not registered.")
+
+        PoldiCreatePeaksFromCell(SpaceGroup=compound.getSpaceGroup(),
+                                 Atoms=compound.getAtomString(),
+                                 LatticeSpacingMin=dMin,
+                                 LatticeSpacingMax=dMax,
+                                 OutputWorkspace=compound.getName(),
+                                 **compound.getCellParameters())
+
+        return compound.getName()
+
+
+AlgorithmFactory.subscribe(PoldiCreatePeaksFromFile)

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

@@ -59,6 +59,7 @@ set ( TEST_PY_FILES
   ExportExperimentLogTest.py
   PoldiMergeTest.py
   VesuvioResolutionTest.py
+  PoldiCreatePeaksFromFileTest.py
 )
 
 check_tests_valid ( ${CMAKE_CURRENT_SOURCE_DIR} ${TEST_PY_FILES} )

Code/Mantid/Framework/PythonInterface/test/python/plugins/algorithms/PoldiCreatePeaksFromFileTest.py

+# pylint: disable=no-init,invalid-name,too-many-public-methods
+import unittest
+from testhelpers import assertRaisesNothing
+from testhelpers.tempfile_wrapper import TemporaryFileHelper
+
+from mantid.kernel import *
+from mantid.api import *
+from mantid.simpleapi import *
+
+class PoldiCreatePeaksFromFileTest(unittest.TestCase):
+    testname = None
+
+    def __init__(self, *args):
+        unittest.TestCase.__init__(self, *args)
+
+    def test_Init(self):
+        assertRaisesNothing(self, AlgorithmManager.create, ("PoldiCreatePeaksFromFile"))
+
+    def test_FileOneCompoundOneAtom(self):
+        fileHelper = TemporaryFileHelper("""Silicon {
+                                                Lattice: 5.43 5.43 5.43 90.0 90.0 90.0
+                                                Spacegroup: F d -3 m
+                                                Atoms: {
+                                                    Si 0 0 0 1.0 0.05
+                                                }
+                                            }""")
+        ws = PoldiCreatePeaksFromFile(fileHelper.getName(), 0.7, 10.0)
+
+        # Check output GroupWorkspace
+        self.assertEquals(ws.getNumberOfEntries(), 1)
+        self.assertTrue(ws.contains("Silicon"))
+
+        # Check that the ouput is identical to what's expected
+        ws_expected = PoldiCreatePeaksFromCell("F d -3 m", "Si 0 0 0 1.0 0.05", a=5.43, LatticeSpacingMin=0.7)
+        si_ws = AnalysisDataService.retrieve("Silicon")
+        self._tablesAreEqual(si_ws, ws_expected)
+
+        # Clean up
+        self._cleanWorkspaces([ws, ws_expected])
+
+    def test_FileOneCompoundTwoAtoms(self):
+        # It's the same structure and the same reflections, just the structure factors are different
+        fileHelper = TemporaryFileHelper("""SiliconCarbon {
+                                                Lattice: 5.43 5.43 5.43 90.0 90.0 90.0
+                                                Spacegroup: F d -3 m
+                                                Atoms: {
+                                                    Si 0 0 0 0.9 0.05
+                                                    C 0 0 0 0.1 0.05
+                                                }
+                                                # Comment
+                                            }""")
+        ws = PoldiCreatePeaksFromFile(fileHelper.getName(), 0.7, 10.0)
+
+        self.assertEquals(ws.getNumberOfEntries(), 1)
+        self.assertTrue(ws.contains("SiliconCarbon"))
+
+        ws_expected = PoldiCreatePeaksFromCell("F d -3 m", "Si 0 0 0 0.9 0.05; C 0 0 0 0.1 0.05", a=5.43,
+                                               LatticeSpacingMin=0.7)
+        si_ws = AnalysisDataService.retrieve("SiliconCarbon")
+        self._tablesAreEqual(si_ws, ws_expected)
+
+        # Clean up
+        self._cleanWorkspaces([ws, ws_expected])
+
+    def test_FileTwoCompounds(self):
+        # It's the same structure and the same reflections, just the structure factors are different
+        fileHelper = TemporaryFileHelper("""SiliconCarbon {
+                                                Lattice: 5.43 5.43 5.43 90.0 90.0 90.0
+                                                Spacegroup: F d -3 m
+                                                Atoms: {
+                                                    Si 0 0 0 0.9 0.05
+                                                    C 0 0 0 0.1 0.05
+                                                }
+                                            }
+                                            Silicon {
+                                                Lattice: 5.43 5.43 5.43 90.0 90.0 90.0
+                                                Spacegroup: F d -3 m
+                                                Atoms: {
+                                                    Si 0 0 0 1.0 0.05
+                                                }
+                                            }""")
+        ws = PoldiCreatePeaksFromFile(fileHelper.getName(), 0.7, 10.0)
+
+        self.assertEquals(ws.getNumberOfEntries(), 2)
+        self.assertTrue(ws.contains("SiliconCarbon"))
+        self.assertTrue(ws.contains("Silicon"))
+
+        self._cleanWorkspaces([ws])
+
+    def test_FileFaultyLatticeStrings(self):
+        fhLatticeMissing = TemporaryFileHelper("""Silicon {
+                                                    Spacegroup: F d -3 m
+                                                    Atoms: {
+                                                        Si 0 0 0 1.0 0.05
+                                                    }
+                                                  }""")
+
+        fhNoLattice = TemporaryFileHelper("""Silicon {
+                                                Lattice:
+                                                Spacegroup: F d -3 m
+                                                Atoms: {
+                                                    Si 0 0 0 1.0 0.05
+                                                }
+                                             }""")
+
+        fhInvalidLattice = TemporaryFileHelper("""Silicon {
+                                                    Lattice: invalid
+                                                    Spacegroup: F d -3 m
+                                                    Atoms: {
+                                                        Si 0 0 0 1.0 0.05
+                                                    }
+                                                  }""")
+
+        self.assertRaises(RuntimeError, PoldiCreatePeaksFromFile, *(fhLatticeMissing.getName(), 0.7, 10.0, 'ws'))
+        self.assertRaises(RuntimeError, PoldiCreatePeaksFromFile, *(fhNoLattice.getName(), 0.7, 10.0, 'ws'))
+        self.assertRaises(RuntimeError, PoldiCreatePeaksFromFile, *(fhInvalidLattice.getName(), 0.7, 10.0, 'ws'))
+
+
+    def test_FileFaultySpaceGroupStrings(self):
+        fhSgMissing = TemporaryFileHelper("""Silicon {
+                                                Lattice: 5.43 5.43 5.43 90.0 90.0 90.0
+                                                Atoms: {
+                                                    Si 0 0 0 1.0 0.05
+                                                }
+                                             }""")
+
+        fhSgInvalid = TemporaryFileHelper("""Silicon {
+                                                    Lattice: 5.43 5.43 5.43 90.0 90.0 90.0
+                                                    Spacegroup: invalid
+                                                    Atoms: {
+                                                        Si 0 0 0 1.0 0.05
+                                                    }
+                                                  }""")
+
+        self.assertRaises(RuntimeError, PoldiCreatePeaksFromFile, *(fhSgMissing.getName(), 0.7, 10.0, 'ws'))
+        self.assertRaises(RuntimeError, PoldiCreatePeaksFromFile, *(fhSgInvalid.getName(), 0.7, 10.0, 'ws'))
+
+    def test_FileFaultyAtomStrings(self):
+        fhAtomsMissing = TemporaryFileHelper("""Silicon {
+                                                Lattice: 5.43 5.43 5.43 90.0 90.0 90.0
+                                                Spacegroup: F d -3 m
+                                             }""")
+
+        fhAtomsNoBraces = TemporaryFileHelper("""Silicon {
+                                                    Lattice: 5.43 5.43 5.43 90.0 90.0 90.0
+                                                    Spacegroup: invalid
+                                                    Atoms:
+                                                        Sis 0 0 0 1.0 0.05
+                                                  }""")
+        fhAtomsEmpty = TemporaryFileHelper("""Silicon {
+                                                    Lattice: 5.43 5.43 5.43 90.0 90.0 90.0
+                                                    Spacegroup: invalid
+                                                    Atoms: { }
+                                                  }""")
+
+        self.assertRaises(RuntimeError, PoldiCreatePeaksFromFile, *(fhAtomsMissing.getName(), 0.7, 10.0, 'ws'))
+        self.assertRaises(RuntimeError, PoldiCreatePeaksFromFile, *(fhAtomsNoBraces.getName(), 0.7, 10.0, 'ws'))
+        self.assertRaises(RuntimeError, PoldiCreatePeaksFromFile, *(fhAtomsEmpty.getName(), 0.7, 10.0, 'ws'))
+
+
+    def _tablesAreEqual(self, lhs, rhs):
+        self.assertEquals(lhs.rowCount(), rhs.rowCount(), msg="Row count of tables is different")
+
+        for r in range(lhs.rowCount()):
+            self.assertEquals(lhs.row(r), rhs.row(r), "Row " + str(r) + " of tables differ.")
+
+    def _cleanWorkspaces(self, wsList):
+        for ws in wsList:
+            DeleteWorkspace(ws)
+
+
+if __name__ == '__main__':
+    unittest.main()

Code/Mantid/Framework/PythonInterface/test/testhelpers/CMakeLists.txt

@@ -5,6 +5,7 @@
 set ( PY_FILES
   __init__.py
   algorithm_decorator.py
+  tempfile_wrapper.py
 )
 
 # Copy python files to output directory

Code/Mantid/Framework/PythonInterface/test/testhelpers/tempfile_wrapper.py

+from tempfile import NamedTemporaryFile
+import os
+
+
+class TemporaryFileHelper(object):
+    """Helper class for temporary files in unit tests
+
+    This class is a small helper for using temporary files for unit test. On instantiation, a temporary file will be
+    created (using NamedTemporaryFile from the tempfile module). If the string argument to the constructor is not empty,
+    its content will be written to that file. The getName()-method provides the name of the temporary file, which can
+    for example be passed to an algorithm that expects a FileProperty. On destruction of the TemporaryFileHelper object,
+    the temporary file is removed automatically using os.unlink().
+
+    Usage:
+      emptyFileHelper = TemporaryFileHelper()
+      fh = open(emptyFileHelper.getName(), 'r+')
+      fh.write("Something or other\n")
+      fh.close()
+
+      filledFileHelper = TemporaryFileHelper("Something or other\n")
+      other = open(filledFileHelper.getName(), 'r')
+      for line in other:
+        print line
+      other.close()
+
+      del emptyFileHelper
+      del filledFileHelper
+    """
+    tempFile = None
+
+    def __init__(self, fileContent=""):
+        self.tempFile = NamedTemporaryFile('r+', delete=False)
+
+        if fileContent:
+            self._setFileContent(fileContent)
+
+    def __del__(self):
+        os.unlink(self.tempFile.name)
+
+    def getName(self):
+        return self.tempFile.name
+
+    def _setFileContent(self, content):
+        fileHandle = open(self.getName(), 'r+')
+        fileHandle.write(content)
+        fileHandle.close()

Code/Mantid/Testing/Data/DocTest/PoldiCrystalFileExample.dat.md5

+0a93f7213e39cf02f7cb7ddb27f4d6f9

Code/Mantid/docs/source/algorithms/PoldiCreatePeaksFromFile-v1.rst

+
+.. algorithm::
+
+.. summary::
+
+.. alias::
+
+.. properties::
+
+Description
+-----------
+
+Some steps in the analysis of POLDI data require that detected peaks are indexed. This can be done by using
+:ref:`algm-PoldiIndexKnownCompounds`, which accepts a table with unindexed peaks and one or more workspaces with
+calculated peaks corresponding to the crystal structures that are expected in the sample. These can be calculated
+using the algorithm :ref:`algm-PoldiCreatePeaksFromCell`. Calling this algorithm over and over with the same
+parameters is not practical, but storing the tables is not practical either, since lattice parameters may change
+slightly from sample to sample.
+
+PoldiCreatePeaksFromFile reads a text file which contains one or more crystal structure definitions. Since the
+analysis requires information mainly about the lattice and the symmetry, the format is very simple. The following
+block shows how such a file would look when there are two compounds:
+
+.. code-block:: none
+
+    # The name may contain letters, numbers and _
+    Iron_FCC {
+        # Up to 6 values in the order a, b, c, alpha, beta, gamma.
+        # Lattice parameters are given in Angstrom.
+        Lattice: 3.65
+        Spacegroup: F m -3 m
+        Atoms: {
+            # Element x y z are mandatory. Optional occupancy and isotropic ADP (in Angstrom^2)
+            Fe 0.0 0.0 0.0
+        }
+    }
+
+    Iron_BCC {
+        Lattice: 2.88
+        Spacegroup: F m -3 m
+        Atoms: {
+            Fe 0.0 0.0 0.0
+        }
+    }
+
+Note that only the atoms in the asymmetric unit need to be specified, the space group is used to generate all
+equivalent atoms. This information is used to determine systematic absences, while the space group is also used by
+some POLDI algorithms to obtain the point group to get reflection multiplicities and more. Anything that follows the
+`#`-character is considered a comment and is ignored by the parser to allow documentation of the crystal structures
+if necessary.
+
+The algorithm will always produce a WorkspaceGroup which contains as many peak tables as compounds specified in the
+file.
+
+Usage
+-----
+
+.. include:: ../usagedata-note.txt
+
+The following usage example takes up the file showed above and passes it to the algorithm.
+
+.. testcode::
+
+    # Create two tables with expected peaks directly from a file
+    compounds = PoldiCreatePeaksFromFile('PoldiCrystalFileExample.dat', LatticeSpacingMin=0.7)
+
+    compound_count = compounds.getNumberOfEntries()
+    print 'Number of loaded compounds:', compound_count
+
+    for i in range(compound_count):
+        ws = compounds.getItem(i)
+        print 'Compound ' + str(i + 1) +':', ws.getName(), 'has', ws.rowCount(), 'reflections in the resolution range.'
+
+
+The script produces a WorkspaceGroup which contains a table with reflections for each compound in the file:
+
+.. testoutput::
+
+    Number of loaded compounds: 2
+    Compound 1: Iron_FCC has 11 reflections in the resolution range.
+    Compound 2: Iron_BCC has 8 reflections in the resolution range.
+
+.. testcleanup::
+
+    DeleteWorkspace('compounds')
+
+.. categories::