diff --git a/Framework/Geometry/inc/MantidGeometry/Crystal/PointGroup.h b/Framework/Geometry/inc/MantidGeometry/Crystal/PointGroup.h index 1e186211dd0f..5d1bad2955ff 100644 --- a/Framework/Geometry/inc/MantidGeometry/Crystal/PointGroup.h +++ b/Framework/Geometry/inc/MantidGeometry/Crystal/PointGroup.h @@ -43,6 +43,8 @@ class MANTID_GEOMETRY_DLL PointGroup : public Group { CrystalSystem crystalSystem() const { return m_crystalSystem; } LatticeSystem latticeSystem() const { return m_latticeSystem; } + std::string getLauePointGroupSymbol() const; + /// Return true if the hkls are in same group bool isEquivalent(const Kernel::V3D &hkl, const Kernel::V3D &hkl2) const; diff --git a/Framework/Geometry/src/Crystal/PointGroup.cpp b/Framework/Geometry/src/Crystal/PointGroup.cpp index aeaa0090b4b8..4374abfeeb7b 100644 --- a/Framework/Geometry/src/Crystal/PointGroup.cpp +++ b/Framework/Geometry/src/Crystal/PointGroup.cpp @@ -8,6 +8,7 @@ #include "MantidGeometry/Crystal/PointGroupFactory.h" #include "MantidGeometry/Crystal/SymmetryElementFactory.h" #include "MantidGeometry/Crystal/SymmetryOperationFactory.h" +#include "MantidKernel/Logger.h" #include "MantidKernel/System.h" #include @@ -18,6 +19,10 @@ namespace Mantid::Geometry { using Kernel::IntMatrix; using Kernel::V3D; +namespace { +/// static logger object +Kernel::Logger g_log("PointGroup"); +} // namespace /** * Returns all equivalent reflections for the supplied hkl. @@ -80,6 +85,60 @@ bool PointGroup::isEquivalent(const Kernel::V3D &hkl, const Kernel::V3D &hkl2) c return (std::find(hklEquivalents.cbegin(), hklEquivalents.cend(), hkl2) != hklEquivalents.end()); } +std::string PointGroup::getLauePointGroupSymbol() const { + switch (m_crystalSystem) { + case CrystalSystem::Triclinic: + return "-1"; + case CrystalSystem::Monoclinic: + if (m_symbolHM.substr(0, 2) == "11") { + return "112/m"; // unique axis c + } else { + return "2/m"; // unique axis b + } + case CrystalSystem::Orthorhombic: + return "mmm"; + case CrystalSystem::Tetragonal: + if (m_symbolHM == "4" || m_symbolHM == "-4" || m_symbolHM == "4/m") { + return "4/m"; + } else { + return "4/mmm"; + } + case CrystalSystem::Hexagonal: + if (m_symbolHM == "6" || m_symbolHM == "-6" || m_symbolHM == "6/m") { + return "6/m"; + } else { + return "6/mmm"; + } + case CrystalSystem::Trigonal: + if (m_symbolHM.substr(m_symbolHM.size() - 1) == "r") { + // Rhombohedral + if (m_symbolHM == "3" || m_symbolHM == "-3") { + return "-3 r"; + } else { + return "-3m r"; + } + } else { + // Hexagonal + if (m_symbolHM == "3" || m_symbolHM == "-3") { + return "-3"; + } else if (m_symbolHM.substr(1) == "1") { + return "-31m"; + } else { + return "-3m1"; + } + } + case CrystalSystem::Cubic: + if (m_symbolHM == "23" || m_symbolHM == "m-3") { + return "m-3"; + } else { + return "m-3m"; + } + default: + g_log.warning() << "Invalid crystal system - returning group with lowest symmetry (inversion only).\n"; + return "-1"; // never used but required for gcc warning + } +} + /** * Generates a set of hkls * diff --git a/Framework/PythonInterface/mantid/geometry/src/Exports/PointGroup.cpp b/Framework/PythonInterface/mantid/geometry/src/Exports/PointGroup.cpp index 8ad1a9deb1c8..2286264a37e4 100644 --- a/Framework/PythonInterface/mantid/geometry/src/Exports/PointGroup.cpp +++ b/Framework/PythonInterface/mantid/geometry/src/Exports/PointGroup.cpp @@ -89,6 +89,7 @@ void export_PointGroup() { .def("getEquivalents", &getEquivalents, (arg("self"), arg("hkl")), "Returns an array with all symmetry equivalents of the supplied " "HKL.") + .def("getLauePointGroupSymbol", &PointGroup::getLauePointGroupSymbol, arg("self")) .def("getReflectionFamily", &getReflectionFamily, (arg("self"), arg("hkl")), "Returns the same HKL for all symmetry equivalents.") .def(str(self)) diff --git a/Framework/PythonInterface/plugins/CMakeLists.txt b/Framework/PythonInterface/plugins/CMakeLists.txt index 05883dbfd663..d2fdff595a9d 100644 --- a/Framework/PythonInterface/plugins/CMakeLists.txt +++ b/Framework/PythonInterface/plugins/CMakeLists.txt @@ -186,6 +186,7 @@ set(PYTHON_PLUGINS algorithms/CNCSSuggestTIB.py algorithms/HYSPECSuggestTIB.py algorithms/Symmetrise.py + algorithms/SymmetriseMDHisto.py algorithms/TOFTOFCropWorkspace.py algorithms/TOFTOFMergeRuns.py algorithms/TestWorkspaceGroupProperty.py diff --git a/Framework/PythonInterface/plugins/algorithms/SymmetriseMDHisto.py b/Framework/PythonInterface/plugins/algorithms/SymmetriseMDHisto.py new file mode 100644 index 000000000000..efac2da19aeb --- /dev/null +++ b/Framework/PythonInterface/plugins/algorithms/SymmetriseMDHisto.py @@ -0,0 +1,174 @@ +# Mantid Repository : https://github.com/mantidproject/mantid +# +# Copyright © 2023 ISIS Rutherford Appleton Laboratory UKRI, +# NScD Oak Ridge National Laboratory, European Spallation Source, +# Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS +# SPDX - License - Identifier: GPL - 3.0 + +from mantid.api import AlgorithmFactory, IMDHistoWorkspaceProperty, PythonAlgorithm, Progress +from mantid.kernel import Direction, EnabledWhenProperty, PropertyCriterion +from mantid.geometry import SpaceGroupFactory, PointGroupFactory +import numpy as np + + +class SymmetriseMDHisto(PythonAlgorithm): + def category(self): + return "Crystal\\DataHandling" + + def summary(self): + return "Symmetrise MDHistoWorkspace using symmetry operations of Laue class point group." + + def PyInit(self): + """Initilize the algorithms properties""" + + self.declareProperty( + IMDHistoWorkspaceProperty("InputWorkspace", "", Direction.Input), + doc="Input MDHistoWorkspace to symmetrise.", + ) + + self.declareProperty( + IMDHistoWorkspaceProperty("OutputWorkspace", "", Direction.Output), + doc="Symmetrised MDHistoWorkspace.", + ) + + self.declareProperty( + name="Pointgroup", + defaultValue="", + direction=Direction.Input, + doc="Point group Hermann–Mauguin symbol used to determine the point group of the Laue class.", + ) + + self.declareProperty( + name="Spacegroup", + defaultValue="", + direction=Direction.Input, + doc="Spacegroup Hermann–Mauguin symbol used to determine the point group of the Laue class.", + ) + + self.declareProperty( + name="WeightedAverage", + defaultValue=False, + direction=Direction.Input, + doc="Perform average weighted by the error squared.", + ) + + enable_spacegroup = EnabledWhenProperty("Pointgroup", PropertyCriterion.IsDefault) + self.setPropertySettings("Spacegroup", enable_spacegroup) + enable_pointgroup = EnabledWhenProperty("Spacegroup", PropertyCriterion.IsDefault) + self.setPropertySettings("Pointgroup", enable_pointgroup) + + def validateInputs(self): + issues = dict() + # check point group of laue class can be retrieved + spgr_sym = self.getProperty("Spacegroup").value + ptgr_sym = self.getProperty("Pointgroup").value + if spgr_sym and not ptgr_sym: + if not SpaceGroupFactory.isSubscribedSymbol(spgr_sym): + issues["Spacegroup"] = "Not a valid spacegroup symbol." + elif ptgr_sym and not spgr_sym: + if not PointGroupFactory.isSubscribed(ptgr_sym): + issues["Pointgroup"] = "Not a valid spacegroup symbol." + else: + issues["Spacegroup"] = "Please only provide one of Spacegroup or Pointgroup." + # check workspace has same extent and binning along all axes + ws = self.getProperty("InputWorkspace").value + # check nDims=3 + ndims = ws.getNumDims() + if ndims != 3: + issues["InputWorkspace"] = "Workspace must have 3 dimensions." + lo, hi, nbins = _get_dim_extents_and_nbins(ws, 0) + if not np.isclose(lo, -hi): + issues["InputWorkspace"] = "Workspace must have dimensions centered on 0 (i.e. min = -max)." + for idim in range(1, ndims): + if not all(np.isclose(_get_dim_extents_and_nbins(ws, idim), [lo, hi, nbins])): + issues["InputWorkspace"] = "Workspace must have same binning along all dimensions." + # check errors exist if WeightedAverage==True (would produce all NaNs in signal) + if self.getProperty("WeightedAverage").value: + if not ws.getErrorSquaredArray().any(): + issues["WeightedAverage"] = "Cannot perform weighted average on data with no errors." + return issues + + def PyExec(self): + """Execute the algorithm""" + ws = self.getProperty("InputWorkspace").value + spgr_sym = self.getProperty("Spacegroup").value + ptgr_sym = self.getProperty("Pointgroup").value + weighted_average = self.getProperty("WeightedAverage").value + + # get point group of Laue class + if spgr_sym: + spgr = SpaceGroupFactory.createSpaceGroup(spgr_sym) + ptgr = PointGroupFactory.createPointGroupFromSpaceGroup(spgr) + laue_ptgr = PointGroupFactory.createPointGroup(ptgr.getLauePointGroupSymbol()) + else: + ptgr = PointGroupFactory.createPointGroup(ptgr_sym) + laue_ptgr = PointGroupFactory.createPointGroup(ptgr.getLauePointGroupSymbol()) + + # setup data array (note symmetry operations include identity so start sums from 0) + nbins = ws.getDimension(0).getNBins() + signal = np.zeros(3 * [nbins]) + weights = np.zeros(signal.shape) + if not weighted_average: + error_sq = np.zeros(signal.shape) + + bin_index = np.arange(signal.size).reshape(signal.shape) + bin_at_zero = nbins % 2 # odd number of bins implies bin center at 0 along each axis + sym_ops = laue_ptgr.getSymmetryOperations() + prog_reporter = Progress(self, start=0.0, end=1.0, nreports=len(sym_ops)) + for sym_op in sym_ops: + prog_reporter.report("Symmetrising") + transformed = sym_op.transformHKL([1, 2, 3]) + axes = [int(abs(iax)) - 1 for iax in transformed] + iflip = np.flatnonzero(np.asarray(transformed) < 0) + signal_tmp = np.flip(np.transpose(ws.getSignalArray().copy(), axes=axes), axis=iflip) + error_sq_tmp = np.flip(np.transpose(ws.getErrorSquaredArray().copy(), axes=axes), axis=iflip) + if sym_op.getIdentifier() != "x,y,z" and bin_at_zero: + # zero bins invariant to symmetry operation (if not identity operation) + transformed_bin_index = np.flip(np.transpose(bin_index, axes=axes), axis=iflip) + i_invariant = transformed_bin_index == bin_index + signal_tmp[i_invariant] = 0 + error_sq_tmp[i_invariant] = 0 + if weighted_average: + with np.errstate(divide="ignore", invalid="ignore"): + weights_tmp = 1.0 / error_sq_tmp + weights_tmp[~np.isfinite(weights_tmp)] = 0 + signal_tmp *= weights_tmp + else: + error_sq += error_sq_tmp + weights_tmp = (abs(signal_tmp) > 1e-10).astype(int) # number of non-empty bins + signal += signal_tmp + weights += weights_tmp + + # normalise signal and errors by weights + inonzero = abs(weights) > 0 + signal[inonzero] = signal[inonzero] / weights[inonzero] + if weighted_average: + # in this case the weights = sum_i[1/error_i^2] + with np.errstate(divide="ignore", invalid="ignore"): + error_sq = 1.0 / weights + error_sq[~np.isfinite(error_sq)] = 0 + else: + # weights is the number of bins that contributed to each bin in the symmetrised sum - i.e. average is a mean + error_sq[inonzero] = error_sq[inonzero] / weights[inonzero] + + # set data on workspace + ws_out = self.child_alg("CloneWorkspace", InputWorkspace=ws) + ws_out.setSignalArray(signal) + ws_out.setErrorSquaredArray(error_sq) + + # assign output + self.setProperty("OutputWorkspace", ws_out) + + def child_alg(self, alg_name, **kwargs): + alg = self.createChildAlgorithm(alg_name, enableLogging=False, StoreInADS=False) + for prop, value in kwargs.items(): + alg.setProperty(prop, value) + alg.execute() + return alg.getProperty("OutputWorkspace").value + + +def _get_dim_extents_and_nbins(ws, idim): + dim = ws.getDimension(idim) + return dim.getMinimum(), dim.getMaximum(), dim.getNBins() + + +AlgorithmFactory.subscribe(SymmetriseMDHisto) diff --git a/Framework/PythonInterface/test/python/plugins/algorithms/CMakeLists.txt b/Framework/PythonInterface/test/python/plugins/algorithms/CMakeLists.txt index 559c114a293c..64e8a907f96a 100644 --- a/Framework/PythonInterface/test/python/plugins/algorithms/CMakeLists.txt +++ b/Framework/PythonInterface/test/python/plugins/algorithms/CMakeLists.txt @@ -128,6 +128,7 @@ set(TEST_PY_FILES CNCSSuggestTIBTest.py HYSPECSuggestTIBTest.py SymmetriseTest.py + SymmetriseMDHistoTest.py UpdatePeakParameterTableValueTest.py SANSSubtractTest.py TOFTOFCropWorkspaceTest.py diff --git a/Framework/PythonInterface/test/python/plugins/algorithms/SymmetriseMDHistoTest.py b/Framework/PythonInterface/test/python/plugins/algorithms/SymmetriseMDHistoTest.py new file mode 100644 index 000000000000..08f592356855 --- /dev/null +++ b/Framework/PythonInterface/test/python/plugins/algorithms/SymmetriseMDHistoTest.py @@ -0,0 +1,119 @@ +# Mantid Repository : https://github.com/mantidproject/mantid +# +# Copyright © 2023 ISIS Rutherford Appleton Laboratory UKRI, +# NScD Oak Ridge National Laboratory, European Spallation Source, +# Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS +# SPDX - License - Identifier: GPL - 3.0 + + +import unittest +from mantid.simpleapi import SymmetriseMDHisto, CreateMDHistoWorkspace, AnalysisDataService as ADS +import numpy as np + + +class SymmeteriseMDHistoTest(unittest.TestCase): + def tearDown(self): + ADS.clear() + + def test_even_nbins_unweighted_avg_full_signal(self): + signal = np.ones(8) # 2 x 2 x 2 + ws = self._make_md_ws(signal, "ws1") + + ws_sym = SymmetriseMDHisto(ws, PointGroup="-1", WeightedAverage=False, OutputWorkspace=ws.name() + "_sym") + + self.assertTrue(np.allclose(signal, ws_sym.getSignalArray().flatten())) + self.assertTrue(np.allclose(signal, ws_sym.getErrorSquaredArray().flatten())) + + def test_even_nbins_weighted_avg_full_signal(self): + signal = np.ones(8) + ws = self._make_md_ws(signal, "ws2") + + ws_sym = SymmetriseMDHisto(ws, PointGroup="-1", WeightedAverage=True, OutputWorkspace=ws.name() + "_sym") + + self.assertTrue(np.allclose(signal, ws_sym.getSignalArray().flatten())) + self.assertTrue(np.allclose(0.5 * signal, ws_sym.getErrorSquaredArray().flatten())) # 2 contributing bins + + def test_even_nbins_unweighted_avg_partial_signal(self): + signal = np.zeros(8) + signal[0] = 1.0 + ws = self._make_md_ws(signal, "ws4") + + ws_sym = SymmetriseMDHisto(ws, PointGroup="-1", WeightedAverage=False, OutputWorkspace=ws.name() + "_sym") + + expected_signal = signal.copy() + expected_signal[-1] = 1.0 + self.assertTrue(np.allclose(expected_signal, ws_sym.getSignalArray().flatten())) + self.assertTrue(np.allclose(expected_signal, ws_sym.getErrorSquaredArray().flatten())) # 1 contributing bin + + def test_even_nbins_weighted_avg_partial_signal(self): + signal = np.zeros(8) + signal[0] = 1.0 + ws = self._make_md_ws(signal, "ws5") + + ws_sym = SymmetriseMDHisto(ws, PointGroup="-1", WeightedAverage=True, OutputWorkspace=ws.name() + "_sym") + + expected_signal = signal.copy() + expected_signal[-1] = 1.0 + self.assertTrue(np.allclose(expected_signal, ws_sym.getSignalArray().flatten())) + self.assertTrue(np.allclose(expected_signal, ws_sym.getErrorSquaredArray().flatten())) # 1 contributing bin + + def test_odd_nbins_unweighted_avg_full_signal(self): + signal = np.ones(27) # 3 x 3 x 3 + ws = self._make_md_ws(signal, "ws6") + + ws_sym = SymmetriseMDHisto(ws, PointGroup="-1", WeightedAverage=False, OutputWorkspace=ws.name() + "_sym") + + self.assertTrue(np.allclose(signal, ws_sym.getSignalArray().flatten())) + self.assertTrue(np.allclose(signal, ws_sym.getErrorSquaredArray().flatten())) + + def test_odd_nbins_weighted_avg_full_signal(self): + signal = np.ones(27) + ws = self._make_md_ws(signal, "ws7") + + ws_sym = SymmetriseMDHisto(ws, PointGroup="-1", WeightedAverage=True, OutputWorkspace=ws.name() + "_sym") + + self.assertTrue(np.allclose(signal, ws_sym.getSignalArray().flatten())) + expected_error_sq = 0.5 * signal + expected_error_sq[len(signal) // 2] = 1 # inversion center at origin not double counted + self.assertTrue(np.allclose(expected_error_sq, ws_sym.getErrorSquaredArray().flatten())) + + def test_odd_nbins_unweighted_avg_partial_signal(self): + signal = np.zeros(27) + signal[0] = 1.0 + ws = self._make_md_ws(signal, "ws8") + + ws_sym = SymmetriseMDHisto(ws, PointGroup="-1", WeightedAverage=False, OutputWorkspace=ws.name() + "_sym") + + expected_signal = signal.copy() + expected_signal[-1] = 1.0 + self.assertTrue(np.allclose(expected_signal, ws_sym.getSignalArray().flatten())) + self.assertTrue(np.allclose(expected_signal, ws_sym.getErrorSquaredArray().flatten())) # 1 contributing bin + + def test_odd_nbins_weighted_avg_partial_signal(self): + signal = np.zeros(27) + signal[0] = 1.0 + ws = self._make_md_ws(signal, "ws9") + + ws_sym = SymmetriseMDHisto(ws, PointGroup="-1", WeightedAverage=True, OutputWorkspace=ws.name() + "_sym") + + expected_signal = signal.copy() + expected_signal[-1] = 1.0 + self.assertTrue(np.allclose(expected_signal, ws_sym.getSignalArray().flatten())) + self.assertTrue(np.allclose(expected_signal, ws_sym.getErrorSquaredArray().flatten())) # 1 contributing bin + + @staticmethod + def _make_md_ws(signal, wsname): + nbins = int(len(signal) ** (1 / 3)) + return CreateMDHistoWorkspace( + SignalInput=signal, + ErrorInput=signal, + Dimensionality=3, + Extents="-1,1,-1,1,-1,1", + NumberOfBins=3 * [nbins], + Names="H,K,L", + Units="rlu,rlu,rlu", + OutputWorkspace=wsname, + ) + + +if __name__ == "__main__": + unittest.main() diff --git a/docs/source/algorithms/SymmetriseMDHisto-v1.rst b/docs/source/algorithms/SymmetriseMDHisto-v1.rst new file mode 100644 index 000000000000..c0d126dc0357 --- /dev/null +++ b/docs/source/algorithms/SymmetriseMDHisto-v1.rst @@ -0,0 +1,53 @@ +.. algorithm:: + +.. summary:: + +.. relatedalgorithms:: + +.. properties:: + +Description +----------- + +Function to symmeterise an :ref:`MDHistoWorkspace ` with 3 dimensions (HKL) using the +symmetry operations of the point group of the Laue class (inferred from the ``PointGroup`` or ``Spacegroup`` provided). + +Two types of average can be performed: + + 1. ``WeightedAverage=False`` (simple mean of the signal values and errors are added in quadrature then divided by the number of contributing bins) + 2. ``WeightedAverage=True`` (weighted average performed using inverse square of errors as weights). + +If (2) ``WeightedAverage=True`` then the signal is given by :math:`s = \Sigma_i w_i s_i` where +:math:`s_i` is the signal in symmetrically equivalent bin :math:`i` and the weights are the inverse of the +error squared in each bin, :math:`w_i = 1/e_i^2`. The final error is given by :math:`1/e^2 = \Sigma_i w_i`. + + +Usage +----- + +**Example:** + +.. testcode:: SymmetriseMDHisto + + from mantid.simpleapi import * + + y = np.r_[4*[1], 4*[0]] + nbins = int(len(y) ** (1 / 3)) + ws = CreateMDHistoWorkspace(SignalInput=y, ErrorInput=y, + Dimensionality=3, Extents='-1,1,-1,1,-1,1', + NumberOfBins=3 * [nbins], Names='H,K,L', + Units='rlu,rlu,rlu') + + ws_sym = SymmetriseMDHisto(ws, PointGroup="-1", WeightedAverage=False) + + print("All HKL voxels have non-zero intensity?", np.all(ws_sym.getSignalArray())) + +Output: + +.. testoutput:: SymmetriseMDHisto + + All HKL voxels have non-zero intensity? True + +.. categories:: + +.. sourcelink:: \ No newline at end of file diff --git a/docs/source/release/v6.9.0/Diffraction/Single_Crystal/New_features/36342.rst b/docs/source/release/v6.9.0/Diffraction/Single_Crystal/New_features/36342.rst new file mode 100644 index 000000000000..2fe6a0a23a83 --- /dev/null +++ b/docs/source/release/v6.9.0/Diffraction/Single_Crystal/New_features/36342.rst @@ -0,0 +1 @@ +- New algorithm :ref:`SymmetriseMDHisto ` to symmetrise :ref:`MDHistoWorkspaces ` by operations in the Laue class of the point group. \ No newline at end of file