From 35e1caf0d3aba000e60da977aa2bc6a37d7576a5 Mon Sep 17 00:00:00 2001
From: Richard Waite <richard.waite@stfc.ac.uk>
Date: Mon, 9 Oct 2023 18:56:07 +0100
Subject: [PATCH 01/24] Add method to return point group symbol of Laue class
 of point group

Point groups are defined here
https://github.com/mantidproject/mantid/blob/34301681ba037303d682493c4948f9a92f440d51/Framework/Geometry/src/Crystal/PointGroupFactory.cpp#L176-L229
---
 .../inc/MantidGeometry/Crystal/PointGroup.h   |  2 +
 Framework/Geometry/src/Crystal/PointGroup.cpp | 49 +++++++++++++++++++
 2 files changed, 51 insertions(+)

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..df9f8f45b52e 100644
--- a/Framework/Geometry/src/Crystal/PointGroup.cpp
+++ b/Framework/Geometry/src/Crystal/PointGroup.cpp
@@ -80,6 +80,55 @@ 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 {
+        return "-3m"; // not sure this works for -31m which is equivilent but different axis convention
+      }
+    }
+  case CrystalSystem::Cubic:
+    if (m_symbolHM == "23" || m_symbolHM == "m-3") {
+      return "m-3";
+    } else {
+      return "m-3m";
+    }
+  }
+}
+
 /**
  * Generates a set of hkls
  *

From 3a929fc8e55efffbd5ed4f0cc504cef7f7012f73 Mon Sep 17 00:00:00 2001
From: Richard Waite <richard.waite@stfc.ac.uk>
Date: Tue, 10 Oct 2023 18:11:30 +0100
Subject: [PATCH 02/24] Add getLauePointGroupSymbol to python export

---
 .../PythonInterface/mantid/geometry/src/Exports/PointGroup.cpp   | 1 +
 1 file changed, 1 insertion(+)

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))

From 52ce3cced431a2b5483d4d6b1a02459b5c594666 Mon Sep 17 00:00:00 2001
From: Richard Waite <richard.waite@stfc.ac.uk>
Date: Fri, 13 Oct 2023 15:49:24 +0100
Subject: [PATCH 03/24] Add algorithm to do Laue symmeterisation of
 MDHistoWorkspaces

---
 .../PythonInterface/plugins/CMakeLists.txt    |   1 +
 .../plugins/algorithms/SymmetriseMDHisto.py   | 117 ++++++++++++++++++
 2 files changed, 118 insertions(+)
 create mode 100644 Framework/PythonInterface/plugins/algorithms/SymmetriseMDHisto.py

diff --git a/Framework/PythonInterface/plugins/CMakeLists.txt b/Framework/PythonInterface/plugins/CMakeLists.txt
index abd30b0ef3db..81a5bf80bd9a 100644
--- a/Framework/PythonInterface/plugins/CMakeLists.txt
+++ b/Framework/PythonInterface/plugins/CMakeLists.txt
@@ -185,6 +185,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..1a9cfd345e32
--- /dev/null
+++ b/Framework/PythonInterface/plugins/algorithms/SymmetriseMDHisto.py
@@ -0,0 +1,117 @@
+# Mantid Repository : https://github.com/mantidproject/mantid
+#
+# Copyright &copy; 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
+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.",
+        )
+
+        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.isSubscribedSymbol(ptgr_sym):
+                issues["Spacegroup"] = "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
+
+        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
+
+        # 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())
+
+        # symmetrise data
+        signal = ws.getSignalArray().copy()
+        npix = np.zeros(signal.shape, dtype=int)
+        inonzero = abs(signal) > 1e-10
+        npix[inonzero] = 1
+
+        for sym_op in laue_ptgr.getSymmetryOperations():
+            transformed = sym_op.transformHKL([1, 2, 3])
+            ws_out = self.child_alg("PermuteMD", InputWorkspace=ws, Axes=[int(abs(iax)) - 1 for iax in transformed])
+            # could call TransformMD with Scaling=np.sign(transformed) but it's actually slower!
+            signs = np.sign(transformed)
+            signal_tmp = ws_out.getSignalArray().copy()[:: int(signs[0]), :: int(signs[1]), :: int(signs[2])]
+            npix_tmp = abs(signal_tmp) > 1e-10
+            signal += signal_tmp
+            npix += npix_tmp
+
+        # set symmetrised signal in ws
+        inonzero = abs(npix) > 0
+        signal[inonzero] = signal[inonzero] / npix[inonzero]
+
+        ws_out.setSignalArray(signal)
+
+        # 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
+
+
+AlgorithmFactory.subscribe(SymmetriseMDHisto)

From 39d1148e6f700a6e35c38e2306c645785149b578 Mon Sep 17 00:00:00 2001
From: Richard Waite <richard.waite@stfc.ac.uk>
Date: Thu, 19 Oct 2023 15:38:50 +0100
Subject: [PATCH 04/24] Add validation for ws dimensions

---
 .../plugins/algorithms/SymmetriseMDHisto.py          | 12 ++++++++++++
 1 file changed, 12 insertions(+)

diff --git a/Framework/PythonInterface/plugins/algorithms/SymmetriseMDHisto.py b/Framework/PythonInterface/plugins/algorithms/SymmetriseMDHisto.py
index 1a9cfd345e32..7158c8d796d2 100644
--- a/Framework/PythonInterface/plugins/algorithms/SymmetriseMDHisto.py
+++ b/Framework/PythonInterface/plugins/algorithms/SymmetriseMDHisto.py
@@ -63,6 +63,13 @@ def validateInputs(self):
         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
+        lo, hi, nbins = _get_dim_extents_and_nbins(ws, 0)
+        if not np.isclose(lo, -hi):
+            issues["InputWorkspace"] = "Workspace must have have dimensions centered on 0 (i.e. min = -max)."
+        for idim in range(1, ws.getNumDims()):
+            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."
 
         return issues
 
@@ -114,4 +121,9 @@ def child_alg(self, alg_name, **kwargs):
         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)

From 52ab0de4fbecdc57cff9688747adbf47ea6aeb60 Mon Sep 17 00:00:00 2001
From: Richard Waite <richard.waite@stfc.ac.uk>
Date: Thu, 26 Oct 2023 16:42:09 +0100
Subject: [PATCH 05/24] Fix typos in point group validation code

---
 .../PythonInterface/plugins/algorithms/SymmetriseMDHisto.py | 6 ++----
 1 file changed, 2 insertions(+), 4 deletions(-)

diff --git a/Framework/PythonInterface/plugins/algorithms/SymmetriseMDHisto.py b/Framework/PythonInterface/plugins/algorithms/SymmetriseMDHisto.py
index 7158c8d796d2..b89275663765 100644
--- a/Framework/PythonInterface/plugins/algorithms/SymmetriseMDHisto.py
+++ b/Framework/PythonInterface/plugins/algorithms/SymmetriseMDHisto.py
@@ -58,8 +58,8 @@ def validateInputs(self):
             if not SpaceGroupFactory.isSubscribedSymbol(spgr_sym):
                 issues["Spacegroup"] = "Not a valid spacegroup symbol."
         elif ptgr_sym and not spgr_sym:
-            if not PointGroupFactory.isSubscribedSymbol(ptgr_sym):
-                issues["Spacegroup"] = "Not a valid spacegroup symbol."
+            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
@@ -93,7 +93,6 @@ def PyExec(self):
         npix = np.zeros(signal.shape, dtype=int)
         inonzero = abs(signal) > 1e-10
         npix[inonzero] = 1
-
         for sym_op in laue_ptgr.getSymmetryOperations():
             transformed = sym_op.transformHKL([1, 2, 3])
             ws_out = self.child_alg("PermuteMD", InputWorkspace=ws, Axes=[int(abs(iax)) - 1 for iax in transformed])
@@ -107,7 +106,6 @@ def PyExec(self):
         # set symmetrised signal in ws
         inonzero = abs(npix) > 0
         signal[inonzero] = signal[inonzero] / npix[inonzero]
-
         ws_out.setSignalArray(signal)
 
         # assign output

From 9385fcea765477ee194b509c32875223eceae964 Mon Sep 17 00:00:00 2001
From: Richard Waite <richard.waite@stfc.ac.uk>
Date: Fri, 27 Oct 2023 15:09:54 +0100
Subject: [PATCH 06/24] Symmeterise errors and add option for weighted
 averaging

Also zero intialise arrays as symmetry operations include identity (was previously double counting the original data in the average)
---
 .../plugins/algorithms/SymmetriseMDHisto.py   | 49 ++++++++++++++-----
 1 file changed, 38 insertions(+), 11 deletions(-)

diff --git a/Framework/PythonInterface/plugins/algorithms/SymmetriseMDHisto.py b/Framework/PythonInterface/plugins/algorithms/SymmetriseMDHisto.py
index b89275663765..9605bba5bbb2 100644
--- a/Framework/PythonInterface/plugins/algorithms/SymmetriseMDHisto.py
+++ b/Framework/PythonInterface/plugins/algorithms/SymmetriseMDHisto.py
@@ -44,6 +44,13 @@ def PyInit(self):
             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)
@@ -78,6 +85,7 @@ def PyExec(self):
         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:
@@ -88,25 +96,44 @@ def PyExec(self):
             ptgr = PointGroupFactory.createPointGroup(ptgr_sym)
             laue_ptgr = PointGroupFactory.createPointGroup(ptgr.getLauePointGroupSymbol())
 
-        # symmetrise data
-        signal = ws.getSignalArray().copy()
-        npix = np.zeros(signal.shape, dtype=int)
-        inonzero = abs(signal) > 1e-10
-        npix[inonzero] = 1
+        # setup data array (note symmetry operations include identity so start sums from 0)
+        signal = np.zeros(3 * [ws.getDimension(0).getNBins()])
+        weights = np.zeros(signal.shape)
+        if not weighted_average:
+            error_sq = np.zeros(signal.shape)
+
         for sym_op in laue_ptgr.getSymmetryOperations():
             transformed = sym_op.transformHKL([1, 2, 3])
             ws_out = self.child_alg("PermuteMD", InputWorkspace=ws, Axes=[int(abs(iax)) - 1 for iax in transformed])
             # could call TransformMD with Scaling=np.sign(transformed) but it's actually slower!
             signs = np.sign(transformed)
             signal_tmp = ws_out.getSignalArray().copy()[:: int(signs[0]), :: int(signs[1]), :: int(signs[2])]
-            npix_tmp = abs(signal_tmp) > 1e-10
-            signal += signal_tmp
-            npix += npix_tmp
+            error_sq_tmp = ws_out.getErrorSquaredArray().copy()[:: int(signs[0]), :: int(signs[1]), :: int(signs[2])]
+            if weighted_average:
+                with np.errstate(divide="ignore", invalid="ignore"):
+                    weights_tmp = 1.0 / error_sq_tmp
+                weights_tmp[~np.isfinite(weights)] = 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
 
-        # set symmetrised signal in ws
-        inonzero = abs(npix) > 0
-        signal[inonzero] = signal[inonzero] / npix[inonzero]
+            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]
         ws_out.setSignalArray(signal)
+        ws_out.setErrorSquaredArray(error_sq)
 
         # assign output
         self.setProperty("OutputWorkspace", ws_out)

From f438f0a10d1fd7ab1ccca41c5c7baddbb65d90e7 Mon Sep 17 00:00:00 2001
From: Richard Waite <richard.waite@stfc.ac.uk>
Date: Fri, 27 Oct 2023 15:16:20 +0100
Subject: [PATCH 07/24] Add validation to ensure errors exist if
 WeightedAverage==True

Weighted average would produce all NaNs in signal
---
 .../PythonInterface/plugins/algorithms/SymmetriseMDHisto.py  | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/Framework/PythonInterface/plugins/algorithms/SymmetriseMDHisto.py b/Framework/PythonInterface/plugins/algorithms/SymmetriseMDHisto.py
index 9605bba5bbb2..c039d6032193 100644
--- a/Framework/PythonInterface/plugins/algorithms/SymmetriseMDHisto.py
+++ b/Framework/PythonInterface/plugins/algorithms/SymmetriseMDHisto.py
@@ -77,7 +77,10 @@ def validateInputs(self):
         for idim in range(1, ws.getNumDims()):
             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.getErrorArraySquared().any():
+                issues["WeightedAverage"] = "Cannot perform weighted average on data with no errors."
         return issues
 
     def PyExec(self):

From 35a547323fc7d048f814ed5318f228d80f0fd6b1 Mon Sep 17 00:00:00 2001
From: Richard Waite <richard.waite@stfc.ac.uk>
Date: Fri, 27 Oct 2023 15:32:10 +0100
Subject: [PATCH 08/24] Reset signal and error at origin to avoid duplicate
 counting in sum

Doesn't matter for unweighted average of signal, but does for errors
---
 .../plugins/algorithms/SymmetriseMDHisto.py              | 9 +++++++++
 1 file changed, 9 insertions(+)

diff --git a/Framework/PythonInterface/plugins/algorithms/SymmetriseMDHisto.py b/Framework/PythonInterface/plugins/algorithms/SymmetriseMDHisto.py
index c039d6032193..0505490de56c 100644
--- a/Framework/PythonInterface/plugins/algorithms/SymmetriseMDHisto.py
+++ b/Framework/PythonInterface/plugins/algorithms/SymmetriseMDHisto.py
@@ -135,6 +135,15 @@ def PyExec(self):
         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]
+
+        # if there is a bin at origin reset signal and error to avoid duplicate counting in sum
+        nbins = signal.shape[0]
+        if not nbins % 2:
+            icen = nbins // 2
+            signal[icen, icen, icen] = ws.getSignalArray()[icen, icen, icen].copy()
+            error_sq[icen, icen, icen] = ws.getErrorSquaredArray()[icen, icen, icen].copy()
+
+        # set data on workspace
         ws_out.setSignalArray(signal)
         ws_out.setErrorSquaredArray(error_sq)
 

From 440a8ad4b561ac911fe1788715ec54de0f9de471 Mon Sep 17 00:00:00 2001
From: Richard Waite <richard.waite@stfc.ac.uk>
Date: Fri, 27 Oct 2023 15:52:07 +0100
Subject: [PATCH 09/24] Reset error at origin to avoid duplicate counting in
 weighted avg sum

Note that both weighted and unweighted average (mean) the signal at the origin is unchanged from original value (a weighted average is equivilent to the mean if the error is the same for all contributing bins)
---
 .../PythonInterface/plugins/algorithms/SymmetriseMDHisto.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/Framework/PythonInterface/plugins/algorithms/SymmetriseMDHisto.py b/Framework/PythonInterface/plugins/algorithms/SymmetriseMDHisto.py
index 0505490de56c..cbdbeec78c84 100644
--- a/Framework/PythonInterface/plugins/algorithms/SymmetriseMDHisto.py
+++ b/Framework/PythonInterface/plugins/algorithms/SymmetriseMDHisto.py
@@ -136,11 +136,11 @@ def PyExec(self):
             # 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]
 
-        # if there is a bin at origin reset signal and error to avoid duplicate counting in sum
+        # if there is a bin at origin reset error to avoid duplicate counting in sum if weighted_average
+        # note in both weighted and unweighted average the signal at the origin is unchanged from original value
         nbins = signal.shape[0]
-        if not nbins % 2:
+        if weighted_average and nbins % 2:
             icen = nbins // 2
-            signal[icen, icen, icen] = ws.getSignalArray()[icen, icen, icen].copy()
             error_sq[icen, icen, icen] = ws.getErrorSquaredArray()[icen, icen, icen].copy()
 
         # set data on workspace

From 679017479f8beb53f04daddc24221029268b35c3 Mon Sep 17 00:00:00 2001
From: Richard Waite <richard.waite@stfc.ac.uk>
Date: Mon, 30 Oct 2023 13:25:30 +0000
Subject: [PATCH 10/24] Fix bug where wrong elements masked in weights of
 weighted avg

---
 .../PythonInterface/plugins/algorithms/SymmetriseMDHisto.py  | 5 ++---
 1 file changed, 2 insertions(+), 3 deletions(-)

diff --git a/Framework/PythonInterface/plugins/algorithms/SymmetriseMDHisto.py b/Framework/PythonInterface/plugins/algorithms/SymmetriseMDHisto.py
index cbdbeec78c84..174312be334d 100644
--- a/Framework/PythonInterface/plugins/algorithms/SymmetriseMDHisto.py
+++ b/Framework/PythonInterface/plugins/algorithms/SymmetriseMDHisto.py
@@ -79,7 +79,7 @@ def validateInputs(self):
                 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.getErrorArraySquared().any():
+            if not ws.getErrorSquaredArray().any():
                 issues["WeightedAverage"] = "Cannot perform weighted average on data with no errors."
         return issues
 
@@ -115,12 +115,11 @@ def PyExec(self):
             if weighted_average:
                 with np.errstate(divide="ignore", invalid="ignore"):
                     weights_tmp = 1.0 / error_sq_tmp
-                weights_tmp[~np.isfinite(weights)] = 0
+                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
 

From 3d0a1d4eeb0fd4048158008acac01814e3e4601c Mon Sep 17 00:00:00 2001
From: Richard Waite <richard.waite@stfc.ac.uk>
Date: Mon, 30 Oct 2023 16:00:22 +0000
Subject: [PATCH 11/24] Replace TransposeMD with np.transpose and use np.flip
 for speed

Reduced execution time from ~300s to ~30s
---
 .../plugins/algorithms/SymmetriseMDHisto.py          | 12 ++++++------
 1 file changed, 6 insertions(+), 6 deletions(-)

diff --git a/Framework/PythonInterface/plugins/algorithms/SymmetriseMDHisto.py b/Framework/PythonInterface/plugins/algorithms/SymmetriseMDHisto.py
index 174312be334d..55fceef1b94d 100644
--- a/Framework/PythonInterface/plugins/algorithms/SymmetriseMDHisto.py
+++ b/Framework/PythonInterface/plugins/algorithms/SymmetriseMDHisto.py
@@ -107,11 +107,10 @@ def PyExec(self):
 
         for sym_op in laue_ptgr.getSymmetryOperations():
             transformed = sym_op.transformHKL([1, 2, 3])
-            ws_out = self.child_alg("PermuteMD", InputWorkspace=ws, Axes=[int(abs(iax)) - 1 for iax in transformed])
-            # could call TransformMD with Scaling=np.sign(transformed) but it's actually slower!
-            signs = np.sign(transformed)
-            signal_tmp = ws_out.getSignalArray().copy()[:: int(signs[0]), :: int(signs[1]), :: int(signs[2])]
-            error_sq_tmp = ws_out.getErrorSquaredArray().copy()[:: int(signs[0]), :: int(signs[1]), :: int(signs[2])]
+            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 weighted_average:
                 with np.errstate(divide="ignore", invalid="ignore"):
                     weights_tmp = 1.0 / error_sq_tmp
@@ -119,7 +118,7 @@ def PyExec(self):
                 signal_tmp *= weights_tmp
             else:
                 error_sq += error_sq_tmp
-                weights_tmp = (abs(signal_tmp) > 1e-10).astype(int)  #  number of non-empty bins
+                weights_tmp = (abs(signal_tmp) > 1e-10).astype(int)  # number of non-empty bins
             signal += signal_tmp
             weights += weights_tmp
 
@@ -143,6 +142,7 @@ def PyExec(self):
             error_sq[icen, icen, icen] = ws.getErrorSquaredArray()[icen, icen, icen].copy()
 
         # set data on workspace
+        ws_out = self.child_alg("CloneWorkspace", InputWorkspace=ws)
         ws_out.setSignalArray(signal)
         ws_out.setErrorSquaredArray(error_sq)
 

From 70354bfa350b50171633d8c3253cb13410213814 Mon Sep 17 00:00:00 2001
From: Richard Waite <richard.waite@stfc.ac.uk>
Date: Mon, 30 Oct 2023 16:31:53 +0000
Subject: [PATCH 12/24] Zero bins invariant to symmetry operation (if not
 identity operation)

---
 .../plugins/algorithms/SymmetriseMDHisto.py    | 18 ++++++++++--------
 1 file changed, 10 insertions(+), 8 deletions(-)

diff --git a/Framework/PythonInterface/plugins/algorithms/SymmetriseMDHisto.py b/Framework/PythonInterface/plugins/algorithms/SymmetriseMDHisto.py
index 55fceef1b94d..b75ec0de0818 100644
--- a/Framework/PythonInterface/plugins/algorithms/SymmetriseMDHisto.py
+++ b/Framework/PythonInterface/plugins/algorithms/SymmetriseMDHisto.py
@@ -100,17 +100,26 @@ def PyExec(self):
             laue_ptgr = PointGroupFactory.createPointGroup(ptgr.getLauePointGroupSymbol())
 
         # setup data array (note symmetry operations include identity so start sums from 0)
-        signal = np.zeros(3 * [ws.getDimension(0).getNBins()])
+        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
         for sym_op in laue_ptgr.getSymmetryOperations():
             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
@@ -134,13 +143,6 @@ def PyExec(self):
             # 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]
 
-        # if there is a bin at origin reset error to avoid duplicate counting in sum if weighted_average
-        # note in both weighted and unweighted average the signal at the origin is unchanged from original value
-        nbins = signal.shape[0]
-        if weighted_average and nbins % 2:
-            icen = nbins // 2
-            error_sq[icen, icen, icen] = ws.getErrorSquaredArray()[icen, icen, icen].copy()
-
         # set data on workspace
         ws_out = self.child_alg("CloneWorkspace", InputWorkspace=ws)
         ws_out.setSignalArray(signal)

From f50194c1228e7dbedbc019d84b77d64584158cd5 Mon Sep 17 00:00:00 2001
From: Richard Waite <richard.waite@stfc.ac.uk>
Date: Mon, 30 Oct 2023 17:40:49 +0000
Subject: [PATCH 13/24] Validation to ensure 3D MDHisto workspaces only

---
 .../plugins/algorithms/SymmetriseMDHisto.py               | 8 ++++++--
 1 file changed, 6 insertions(+), 2 deletions(-)

diff --git a/Framework/PythonInterface/plugins/algorithms/SymmetriseMDHisto.py b/Framework/PythonInterface/plugins/algorithms/SymmetriseMDHisto.py
index b75ec0de0818..3ffe93ea2b32 100644
--- a/Framework/PythonInterface/plugins/algorithms/SymmetriseMDHisto.py
+++ b/Framework/PythonInterface/plugins/algorithms/SymmetriseMDHisto.py
@@ -71,10 +71,14 @@ def validateInputs(self):
             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 have dimensions centered on 0 (i.e. min = -max)."
-        for idim in range(1, ws.getNumDims()):
+            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)

From 3ba898692f2d284639ca03a50f7c77e2fd6fa887 Mon Sep 17 00:00:00 2001
From: Richard Waite <richard.waite@stfc.ac.uk>
Date: Mon, 30 Oct 2023 17:43:00 +0000
Subject: [PATCH 14/24] Add unit test for even nbin (i.e. no bins centers at 0)

---
 .../python/plugins/algorithms/CMakeLists.txt  |  1 +
 .../algorithms/SymmetriseMDHistoTest.py       | 75 +++++++++++++++++++
 2 files changed, 76 insertions(+)
 create mode 100644 Framework/PythonInterface/test/python/plugins/algorithms/SymmetriseMDHistoTest.py

diff --git a/Framework/PythonInterface/test/python/plugins/algorithms/CMakeLists.txt b/Framework/PythonInterface/test/python/plugins/algorithms/CMakeLists.txt
index da51be20860b..107fa1f8858b 100644
--- a/Framework/PythonInterface/test/python/plugins/algorithms/CMakeLists.txt
+++ b/Framework/PythonInterface/test/python/plugins/algorithms/CMakeLists.txt
@@ -127,6 +127,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..335f8262e2cb
--- /dev/null
+++ b/Framework/PythonInterface/test/python/plugins/algorithms/SymmetriseMDHistoTest.py
@@ -0,0 +1,75 @@
+# Mantid Repository : https://github.com/mantidproject/mantid
+#
+# Copyright &copy; 2022 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)
+        ws = self._make_md_ws(np.ones(8), "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(np.ones(8), "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, "ws4")
+
+        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="2,2,2",
+            Names="H,K,L",
+            Units="rlu,rlu,rlu",
+            OutputWorkspace=wsname,
+        )
+
+
+if __name__ == "__main__":
+    unittest.main()

From 6edc32fb2692b94fae86f0b88d0bfeca6e389731 Mon Sep 17 00:00:00 2001
From: Richard Waite <richard.waite@stfc.ac.uk>
Date: Tue, 31 Oct 2023 09:58:03 +0000
Subject: [PATCH 15/24] Add unit test for odd nbin (i.e. bins centers at 0)

---
 .../algorithms/SymmetriseMDHistoTest.py       | 54 +++++++++++++++++--
 1 file changed, 49 insertions(+), 5 deletions(-)

diff --git a/Framework/PythonInterface/test/python/plugins/algorithms/SymmetriseMDHistoTest.py b/Framework/PythonInterface/test/python/plugins/algorithms/SymmetriseMDHistoTest.py
index 335f8262e2cb..aac891d49d64 100644
--- a/Framework/PythonInterface/test/python/plugins/algorithms/SymmetriseMDHistoTest.py
+++ b/Framework/PythonInterface/test/python/plugins/algorithms/SymmetriseMDHistoTest.py
@@ -15,8 +15,8 @@ def tearDown(self):
         ADS.clear()
 
     def test_even_nbins_unweighted_avg_full_signal(self):
-        signal = np.ones(8)
-        ws = self._make_md_ws(np.ones(8), "ws1")
+        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")
 
@@ -25,7 +25,7 @@ def test_even_nbins_unweighted_avg_full_signal(self):
 
     def test_even_nbins_weighted_avg_full_signal(self):
         signal = np.ones(8)
-        ws = self._make_md_ws(np.ones(8), "ws2")
+        ws = self._make_md_ws(signal, "ws2")
 
         ws_sym = SymmetriseMDHisto(ws, PointGroup="-1", WeightedAverage=True, OutputWorkspace=ws.name() + "_sym")
 
@@ -47,7 +47,51 @@ def test_even_nbins_unweighted_avg_partial_signal(self):
     def test_even_nbins_weighted_avg_partial_signal(self):
         signal = np.zeros(8)
         signal[0] = 1.0
-        ws = self._make_md_ws(signal, "ws4")
+        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")
 
@@ -64,7 +108,7 @@ def _make_md_ws(signal, wsname):
             ErrorInput=signal,
             Dimensionality=3,
             Extents="-1,1,-1,1,-1,1",
-            NumberOfBins="2,2,2",
+            NumberOfBins=3 * [nbins],
             Names="H,K,L",
             Units="rlu,rlu,rlu",
             OutputWorkspace=wsname,

From c25fa3be3f27e2fee1a556636ea408df2d7bc604 Mon Sep 17 00:00:00 2001
From: Richard Waite <richard.waite@stfc.ac.uk>
Date: Tue, 31 Oct 2023 14:46:27 +0000
Subject: [PATCH 16/24] Add algorithm docs

---
 .../algorithms/SymmetriseMDHisto-v1.rst       | 45 +++++++++++++++++++
 1 file changed, 45 insertions(+)
 create mode 100644 docs/source/algorithms/SymmetriseMDHisto-v1.rst

diff --git a/docs/source/algorithms/SymmetriseMDHisto-v1.rst b/docs/source/algorithms/SymmetriseMDHisto-v1.rst
new file mode 100644
index 000000000000..a479aa37354d
--- /dev/null
+++ b/docs/source/algorithms/SymmetriseMDHisto-v1.rst
@@ -0,0 +1,45 @@
+.. algorithm::
+
+.. summary::
+
+.. relatedalgorithms::
+
+.. properties::
+
+Description
+-----------
+
+Function to symmeterise an :ref:`MDHistoWorkspace <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:**
+
+.. code-block:: python
+
+    from mantid.simpleapi import *
+
+    y = np.ones(8)
+    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)
+
+.. categories::
+
+.. sourcelink::
\ No newline at end of file

From 7da96b60ce0fd839ddf24f1d6c038a6aa42a6ce4 Mon Sep 17 00:00:00 2001
From: Richard Waite <richard.waite@stfc.ac.uk>
Date: Tue, 31 Oct 2023 15:27:21 +0000
Subject: [PATCH 17/24] Handle different axes conventions in hexagonal groups

---
 Framework/Geometry/src/Crystal/PointGroup.cpp | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/Framework/Geometry/src/Crystal/PointGroup.cpp b/Framework/Geometry/src/Crystal/PointGroup.cpp
index df9f8f45b52e..38c537975ebf 100644
--- a/Framework/Geometry/src/Crystal/PointGroup.cpp
+++ b/Framework/Geometry/src/Crystal/PointGroup.cpp
@@ -116,8 +116,10 @@ std::string PointGroup::getLauePointGroupSymbol() const {
       // Hexagonal
       if (m_symbolHM == "3" || m_symbolHM == "-3") {
         return "-3";
+      } else if (m_symbolHM.substr(1) == "1") {
+        return "-31m";
       } else {
-        return "-3m"; // not sure this works for -31m which is equivilent but different axis convention
+        return "-3m1";
       }
     }
   case CrystalSystem::Cubic:

From 8ccd407c6872ed2e9845454a531b3bf3ba445bf9 Mon Sep 17 00:00:00 2001
From: Richard Waite <richard.waite@stfc.ac.uk>
Date: Wed, 29 Nov 2023 09:53:27 +0000
Subject: [PATCH 18/24] Fix gcc warning

---
 Framework/Geometry/src/Crystal/PointGroup.cpp | 1 +
 1 file changed, 1 insertion(+)

diff --git a/Framework/Geometry/src/Crystal/PointGroup.cpp b/Framework/Geometry/src/Crystal/PointGroup.cpp
index 38c537975ebf..9ac2dfd20024 100644
--- a/Framework/Geometry/src/Crystal/PointGroup.cpp
+++ b/Framework/Geometry/src/Crystal/PointGroup.cpp
@@ -129,6 +129,7 @@ std::string PointGroup::getLauePointGroupSymbol() const {
       return "m-3m";
     }
   }
+  return "-1"; // never used but required for gcc warning
 }
 
 /**

From 94d1c51e3ea476ca2698c6a1294001e723f364a4 Mon Sep 17 00:00:00 2001
From: Richard Waite <richard.waite@stfc.ac.uk>
Date: Thu, 21 Dec 2023 11:13:32 +0000
Subject: [PATCH 19/24] Add default block and warning

---
 Framework/Geometry/src/Crystal/PointGroup.cpp | 9 ++++++++-
 1 file changed, 8 insertions(+), 1 deletion(-)

diff --git a/Framework/Geometry/src/Crystal/PointGroup.cpp b/Framework/Geometry/src/Crystal/PointGroup.cpp
index 9ac2dfd20024..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 <algorithm>
@@ -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.
@@ -128,8 +133,10 @@ std::string PointGroup::getLauePointGroupSymbol() const {
     } 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
   }
-  return "-1"; // never used but required for gcc warning
 }
 
 /**

From a67792557ccf2c5b1369f3ae8d2acb3598f9171b Mon Sep 17 00:00:00 2001
From: Richard Waite <richard.waite@stfc.ac.uk>
Date: Thu, 21 Dec 2023 11:22:06 +0000
Subject: [PATCH 20/24] Add doc test

---
 .../python/plugins/algorithms/SymmetriseMDHistoTest.py |  2 +-
 docs/source/algorithms/SymmetriseMDHisto-v1.rst        | 10 ++++++++--
 2 files changed, 9 insertions(+), 3 deletions(-)

diff --git a/Framework/PythonInterface/test/python/plugins/algorithms/SymmetriseMDHistoTest.py b/Framework/PythonInterface/test/python/plugins/algorithms/SymmetriseMDHistoTest.py
index aac891d49d64..08f592356855 100644
--- a/Framework/PythonInterface/test/python/plugins/algorithms/SymmetriseMDHistoTest.py
+++ b/Framework/PythonInterface/test/python/plugins/algorithms/SymmetriseMDHistoTest.py
@@ -1,6 +1,6 @@
 # Mantid Repository : https://github.com/mantidproject/mantid
 #
-# Copyright &copy; 2022 ISIS Rutherford Appleton Laboratory UKRI,
+# Copyright &copy; 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 +
diff --git a/docs/source/algorithms/SymmetriseMDHisto-v1.rst b/docs/source/algorithms/SymmetriseMDHisto-v1.rst
index a479aa37354d..96327409ae39 100644
--- a/docs/source/algorithms/SymmetriseMDHisto-v1.rst
+++ b/docs/source/algorithms/SymmetriseMDHisto-v1.rst
@@ -27,11 +27,11 @@ Usage
 
 **Example:**
 
-.. code-block:: python
+.. testcode:: SymmetriseMDHisto
 
     from mantid.simpleapi import *
 
-    y = np.ones(8)
+    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',
@@ -40,6 +40,12 @@ Usage
 
     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

From cd7276c309b15bdaccfde123229da955d4461e91 Mon Sep 17 00:00:00 2001
From: Richard Waite <richard.waite@stfc.ac.uk>
Date: Thu, 21 Dec 2023 11:41:07 +0000
Subject: [PATCH 21/24] Add progress reporting

---
 .../plugins/algorithms/SymmetriseMDHisto.py                | 7 +++++--
 1 file changed, 5 insertions(+), 2 deletions(-)

diff --git a/Framework/PythonInterface/plugins/algorithms/SymmetriseMDHisto.py b/Framework/PythonInterface/plugins/algorithms/SymmetriseMDHisto.py
index 3ffe93ea2b32..efac2da19aeb 100644
--- a/Framework/PythonInterface/plugins/algorithms/SymmetriseMDHisto.py
+++ b/Framework/PythonInterface/plugins/algorithms/SymmetriseMDHisto.py
@@ -4,7 +4,7 @@
 #   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
+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
@@ -112,7 +112,10 @@ def PyExec(self):
 
         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
-        for sym_op in laue_ptgr.getSymmetryOperations():
+        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)

From 41bd47b3f299339ad658d59df6e93b541819ab56 Mon Sep 17 00:00:00 2001
From: Richard Waite <richard.waite@stfc.ac.uk>
Date: Thu, 21 Dec 2023 11:47:57 +0000
Subject: [PATCH 22/24] Fix static analysis error

---
 docs/source/algorithms/SymmetriseMDHisto-v1.rst | 1 +
 1 file changed, 1 insertion(+)

diff --git a/docs/source/algorithms/SymmetriseMDHisto-v1.rst b/docs/source/algorithms/SymmetriseMDHisto-v1.rst
index 96327409ae39..d3087ed32129 100644
--- a/docs/source/algorithms/SymmetriseMDHisto-v1.rst
+++ b/docs/source/algorithms/SymmetriseMDHisto-v1.rst
@@ -44,6 +44,7 @@ Usage
 
 Output:
 .. testoutput:: SymmetriseMDHisto
+
     All HKL voxels have non-zero intensity? True
 
 .. categories::

From 2cc410c579a59fb3f8df1587495a739fcdc3a93d Mon Sep 17 00:00:00 2001
From: Richard Waite <richard.waite@stfc.ac.uk>
Date: Thu, 21 Dec 2023 11:53:58 +0000
Subject: [PATCH 23/24] Add release note

---
 .../v6.9.0/Diffraction/Single_Crystal/New_features/36342.rst     | 1 +
 1 file changed, 1 insertion(+)
 create mode 100644 docs/source/release/v6.9.0/Diffraction/Single_Crystal/New_features/36342.rst

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 <algm-SymmetriseMDHisto>` to symmetrise :ref:`MDHistoWorkspaces <MDHistoWorkspace>` by operations in the Laue class of the point group.
\ No newline at end of file

From bbf205037df705641f08e28155c0f00608387e09 Mon Sep 17 00:00:00 2001
From: Richard Waite <richard.waite@stfc.ac.uk>
Date: Thu, 21 Dec 2023 15:02:21 +0000
Subject: [PATCH 24/24] Fix doc test with new line

---
 docs/source/algorithms/SymmetriseMDHisto-v1.rst | 1 +
 1 file changed, 1 insertion(+)

diff --git a/docs/source/algorithms/SymmetriseMDHisto-v1.rst b/docs/source/algorithms/SymmetriseMDHisto-v1.rst
index d3087ed32129..c0d126dc0357 100644
--- a/docs/source/algorithms/SymmetriseMDHisto-v1.rst
+++ b/docs/source/algorithms/SymmetriseMDHisto-v1.rst
@@ -43,6 +43,7 @@ Usage
     print("All HKL voxels have non-zero intensity?", np.all(ws_sym.getSignalArray()))
 
 Output:
+
 .. testoutput:: SymmetriseMDHisto
 
     All HKL voxels have non-zero intensity? True