2D Pawley refinement for POLDI

Framework/PythonInterface/plugins/algorithms/PoldiAutoCorrelation6.py

@@ -5,14 +5,17 @@
 #   Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
 # SPDX - License - Identifier: GPL - 3.0 +
 from mantid.api import AlgorithmFactory, MatrixWorkspaceProperty, PythonAlgorithm, Progress
-from mantid.kernel import Direction, FloatBoundedValidator, UnitConversion, DeltaEModeType, UnitParams, UnitParametersMap
+from mantid.kernel import Direction, FloatBoundedValidator, UnitParams, StringListValidator
 import numpy as np
 from plugins.algorithms.poldi_utils import (
     get_instrument_settings_from_log,
     get_max_tof_from_chopper,
     get_dspac_limits,
     get_final_dspac_array,
 )
+from joblib import Parallel, delayed
+from functools import reduce
+from operator import iadd
 
 
 class PoldiAutoCorrelation(PythonAlgorithm):
@@ -51,6 +54,14 @@ def PyInit(self):
             validator=positive_float_validator,
             doc="Maximum wavelength to consider.",
         )
+        self.declareProperty(
+            "InterpolationMethod",
+            defaultValue="Linear",
+            direction=Direction.Input,
+            validator=StringListValidator(["Linear", "Nearest"]),
+            doc="Interpolation used when adding intensity to a given bin in correlation function - "
+            "'Nearest' is quicker but potentially less accurate.",
+        )
 
     def validateInputs(self):
         issues = dict()
@@ -79,8 +90,9 @@ def PyExec(self):
         cycle_time, slit_offsets, t0_const, l1_chop = get_instrument_settings_from_log(ws)
         # get detector positions from IDF
         si = ws.spectrumInfo()
-        tths = np.array([si.twoTheta(ispec) for ispec in range(ws.getNumberHistograms())])
-        l2s = [si.l2(ispec) for ispec in range(ws.getNumberHistograms())]
+        nspec = ws.getNumberHistograms()
+        tths = np.array([si.twoTheta(ispec) for ispec in range(nspec)])
+        l2s = np.asarray([si.l2(ispec) for ispec in range(nspec)])
         l1 = si.l1()
         # determine npulses to include in calc
         time_max = get_max_tof_from_chopper(l1, l1_chop, l2s, tths, lambda_max) + slit_offsets[0]
@@ -90,35 +102,28 @@ def PyExec(self):
         # but this is a small effect as detector doesn't cover much two-theta range
         dspac_min, dspac_max = get_dspac_limits(tths.min(), tths.max(), lambda_min, lambda_max)
         bin_width = ws.readX(0)[1] - ws.readX(0)[0]
-        dspacs = get_final_dspac_array(bin_width, dspac_min, dspac_max, time_max)
+        dspacs = get_final_dspac_array(bin_width, dspac_min, dspac_max, time_max)[:, None]
         # perform auto-correlation (Eq. 7 in POLDI concept paper)
         ipulses = np.arange(npulses)[:, None]
-        offsets = (ipulses * cycle_time + slit_offsets).flatten()
-        inter_corr = np.zeros((len(dspacs), len(offsets)), dtype=float)  # npulses*nslits
-        # sum over all spectra
-        params = UnitParametersMap()
-        nspec = ws.getNumberHistograms()
+        offsets = (ipulses * cycle_time + slit_offsets + t0_const).flatten()
+        # get time-of-flight from chopper to detector for neutron corresponding to d=1Ang
+        path_length_ratio = (l2s + l1 - l1_chop) / (l2s + l1)
+        tof_d1Ang = np.asarray([si.diffractometerConstants(ispec)[UnitParams.difc] * path_length_ratio[ispec] for ispec in range(nspec)])
+        # loop over spectra and add to intermediate correlation
         progress = Progress(self, start=0.0, end=1.0, nreports=nspec)
-        for ispec in range(nspec):
-            params[UnitParams.l2] = l2s[ispec]
-            params[UnitParams.twoTheta] = tths[ispec]
-            # TOF from chopper to detector for wavelength corresponding to d=1Ang
-            # equivalent to DIFC * (Ltot - L1_source-chop) / Ltot
-            tof_d1Ang = UnitConversion.run("dSpacing", "TOF", 1.0, l1 - l1_chop, DeltaEModeType.Elastic, params)
-            arrival_time = tof_d1Ang * dspacs[:, None] + offsets + t0_const
-            # detector clock reset for each pulse
-            # equivalent to np.mod(arrival_time, cycle_time) but order of magnitude quicker
-            time_in_cycle_period = arrival_time - cycle_time * np.floor(arrival_time / cycle_time)
-            itimes = (time_in_cycle_period / bin_width) - 0.5  # correct for first time bin center at bin_width / 2
-            ibins = np.floor(itimes).astype(int)
-            ibins_plus = ibins + 1
-            ibins_plus[ibins_plus > ws.blocksize() - 1] = 0
-            y = ws.readY(ispec)
-            inter_corr += (ibins + 1 - itimes) * y[ibins] + (itimes - ibins) * y[ibins_plus]
-            progress.report()
+        if self.getProperty("InterpolationMethod").value == "Linear":
+            do_autocorr = _autocorr_spec_linear
+        else:
+            do_autocorr = _autocorr_spec_nearest
+        inter_corr = reduce(
+            iadd,
+            Parallel(n_jobs=-2, prefer="threads", return_as="generator_unordered")(
+                delayed(do_autocorr)(ws.readY(ispec), tof_d1Ang[ispec], dspacs, offsets, bin_width, progress) for ispec in range(nspec)
+            ),
+        )  # npulses*nslits
         # average of inverse intermediate correlation func (Eq.8 in POLDI concept paper)
         with np.errstate(divide="ignore", invalid="ignore"):
-            corr = 1 / np.sum(1 / inter_corr, axis=1)
+            corr = 1 / np.nansum(1 / inter_corr, axis=1)
         ws_corr = self.exec_child_alg("CreateWorkspace", DataX=dspacs, DataY=corr, UnitX="dSpacing", YUnitLabel="Intensity (a.u.)")
         ws_corr = self.exec_child_alg("ConvertUnits", InputWorkspace=ws_corr, Target="MomentumTransfer")
         self.setProperty("OutputWorkspace", ws_corr)
@@ -133,4 +138,23 @@ def exec_child_alg(self, alg_name: str, **kwargs):
         return out_props[0] if len(out_props) == 1 else out_props
 
 
+def _autocorr_spec_linear(y, tof_d1Ang, dspacs, offsets, bin_width, progress):
+    arrival_time = tof_d1Ang * dspacs + offsets
+    itimes = (arrival_time / bin_width) - 0.5  # correct for first time bin center at bin_width / 2
+    ibins = itimes.astype(int)  # quicker than np.floor
+    frac_bin = ibins + 1 - itimes
+    progress.report()
+    return frac_bin * np.take(y, ibins, mode="wrap") + (1 - frac_bin) * np.take(y, ibins + 1, mode="wrap")
+
+
+def _autocorr_spec_nearest(y, tof_d1Ang, dspacs, offsets, bin_width, progress):
+    arrival_time = tof_d1Ang * dspacs + offsets
+    # round to nearest bin - note this truncation is not equivalent to floor as
+    # implicitly subtracting -0.5 as first time bin center at bin_width / 2 then
+    # add 0.5 so that truncation results in rounding of input
+    ibins = (arrival_time / bin_width).astype(int)
+    progress.report()
+    return np.take(y, ibins, mode="wrap")
+
+
 AlgorithmFactory.subscribe(PoldiAutoCorrelation)

Framework/PythonInterface/plugins/algorithms/poldi_utils.py

@@ -9,7 +9,7 @@
 import numpy as np
 from mantid.kernel import UnitConversion, DeltaEModeType, UnitParams, UnitParametersMap
 from typing import Tuple, Sequence, Optional, TYPE_CHECKING
-
+from joblib import Parallel, delayed
 
 if TYPE_CHECKING:
     from mantid.dataobjects import Workspace2D
@@ -19,7 +19,7 @@ def exec_alg(alg_name: str, **kwargs):
     alg = AlgorithmManager.create(alg_name)
     alg.initialize()
     alg.setAlwaysStoreInADS(False)
-    alg.setLogging(True)
+    alg.setLogging(False)
     alg.setProperties(kwargs)
     alg.execute()
     out_props = tuple(alg.getProperty(prop).value for prop in alg.outputProperties())
@@ -46,6 +46,8 @@ def load_poldi(
     """
     ws = exec_alg("LoadEmptyInstrument", Filename=fpath_idf)
     dat = np.loadtxt(fpath_data)
+    if dat.shape[0] == 500:
+        dat = dat.T
     cycle_time = _calc_cycle_time_from_chopper_speed(chopper_speed)
     bin_width = cycle_time / dat.shape[-1]
     ws = exec_alg("Rebin", InputWorkspace=ws, Params=f"0,{bin_width},{cycle_time}")
@@ -73,6 +75,10 @@ def _calc_cycle_time_from_chopper_speed(chopper_speed):
     return 60.0 / (4.0 * chopper_speed) * 1.0e6  # mus
 
 
+def get_t0_parameters(chopper):
+    return chopper.getNumberParameter("t0")[0], chopper.getNumberParameter("t0_const")[0]
+
+
 def get_instrument_settings_from_log(ws: Workspace2D) -> Tuple[float, np.ndarray[float], float, float]:
     """
     Function to get instrument settings from logs stored on workspace
@@ -86,8 +92,7 @@ def get_instrument_settings_from_log(ws: Workspace2D) -> Tuple[float, np.ndarray
     source = inst.getSource()  # might not need these
     chopper = inst.getComponentByName("chopper")
     l1_chop = (chopper.getPos() - source.getPos()).norm()
-    t0 = chopper.getNumberParameter("t0")[0]
-    t0_const = chopper.getNumberParameter("t0_const")[0]
+    t0, t0_const = get_t0_parameters(chopper)
     chopper_speed = ws.run().getPropertyAsSingleValue("chopperspeed")  # rpm
     cycle_time = _calc_cycle_time_from_chopper_speed(chopper_speed)  # mus
     # get chopper offsets in time (stored as x position of child components of chopper)
@@ -124,6 +129,28 @@ def get_final_dspac_array(bin_width: float, dspac_min: float, dspac_max: float,
     return np.linspace(dspac_min, dspac_max, nbins_dspac)
 
 
+def get_dspac_array_from_ws(ws: Workspace2D, lambda_min: float = 1.1, lambda_max: float = 5.0) -> np.ndarray[float]:
+    """
+    Function to calculate d-spacing bins given workspace
+    :param ws_2d: MatrixWorkspace containing POLDI data (raw instrument)
+    :param lambda_min: maximum wavelength (Ang) to consider
+    :param lambda_max: maximum wavelength (Ang) to consider
+    :return np.ndarray: array of d-spacing bins
+    """
+    _, slit_offsets, _, l1_chop = get_instrument_settings_from_log(ws)
+    # get detector positions from IDF
+    si = ws.spectrumInfo()
+    nspec = ws.getNumberHistograms()
+    tths = np.array([si.twoTheta(ispec) for ispec in range(nspec)])
+    l2s = np.asarray([si.l2(ispec) for ispec in range(nspec)])
+    l1 = si.l1()
+    # determine npulses to include in calc
+    time_max = get_max_tof_from_chopper(l1, l1_chop, l2s, tths, lambda_max) + slit_offsets[0]
+    dspac_min, dspac_max = get_dspac_limits(tths.min(), tths.max(), lambda_min, lambda_max)
+    bin_width = ws.readX(0)[1] - ws.readX(0)[0]
+    return get_final_dspac_array(bin_width, dspac_min, dspac_max, time_max)
+
+
 def get_max_tof_from_chopper(l1: float, l1_chop: float, l2s: Sequence[float], tths: Sequence[float], lambda_max: float) -> float:
     """
     Function to calculate TOF of longest wavelength neutron to reach detector with the largest total path
@@ -166,22 +193,27 @@ def simulate_2d_data(
     # get instrument settings
     cycle_time, slit_offsets, t0_const, l1_chop = get_instrument_settings_from_log(ws_sim)
     si = ws_sim.spectrumInfo()
-    tths = np.array([si.twoTheta(ispec) for ispec in range(ws_sim.getNumberHistograms())])
-    l2s = [si.l2(ispec) for ispec in range(ws_sim.getNumberHistograms())]
+    nspec = ws_sim.getNumberHistograms()
+    tths = np.asarray([si.twoTheta(ispec) for ispec in range(nspec)])
+    l2s = np.asarray([si.l2(ispec) for ispec in range(nspec)])
     l1 = si.l1()
     # get npulses to include
     time_max = get_max_tof_from_chopper(l1, l1_chop, l2s, tths, lambda_max) + slit_offsets[0]
     npulses = int(time_max // cycle_time)
     # simulate detected spectra
     ipulses = np.arange(npulses)[:, None]
-    offsets = (ipulses * cycle_time - slit_offsets).flatten()  # note different sign to auto-corr!
-    tofs = ws_sim.readX(0)[:, None] + offsets - t0_const  # same for all spectra
-    params = UnitParametersMap()
-    for ispec in range(ws_sim.getNumberHistograms()):
-        params[UnitParams.l2] = l2s[ispec]
-        params[UnitParams.twoTheta] = tths[ispec]
-        tof_d1Ang = UnitConversion.run("dSpacing", "TOF", 1.0, l1 - l1_chop, DeltaEModeType.Elastic, params)
-        ds = tofs / tof_d1Ang
-        ws_sim.setY(ispec, np.interp(ds, ws_1d.readX(0), ws_1d.readY(0)).sum(axis=1))
+    offsets = (ipulses * cycle_time - slit_offsets - t0_const).flatten()  # note different sign to auto-corr!
+    tofs = ws_sim.readX(0)[:, None] + offsets  # same for all spectra
+    path_length_ratio = (l2s + l1 - l1_chop) / (l2s + l1)
+    tof_d1Ang = np.asarray([si.diffractometerConstants(ispec)[UnitParams.difc] * path_length_ratio[ispec] for ispec in range(nspec)])
+    out = Parallel(n_jobs=-2, prefer="threads", return_as="generator")(
+        delayed(_do_interp)(tofs / tof_d1Ang[ispec], ws_1d.readX(0), ws_1d.readY(0)) for ispec in range(nspec)
+    )
+    # set y values
+    [ws_sim.setY(ispec, yvec) for ispec, yvec in enumerate(out)]
     ADS.addOrReplace(output_workspace, ws_sim)
     return ws_sim
+
+
+def _do_interp(dtarget, d, intensity):
+    return np.interp(dtarget, d, intensity).sum(axis=1)

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

@@ -26,15 +26,11 @@ def tearDownClass(cls):
 
     def test_exec_default_wavelength_range(self):
         ws_corr = PoldiAutoCorrelation(InputWorkspace=self.ws, OutputWorkspace="ws_corr", Version=6)
+        self._assert_auto_corr_workspace(ws_corr)
 
-        # assert min/max Q
-        self.assertAlmostEqual(ws_corr.readX(0)[0], 1.5144, delta=1e-3)
-        self.assertAlmostEqual(ws_corr.readX(0)[-1], 9.0832, delta=1e-3)
-        # assert bin width/number bins
-        self.assertEqual(ws_corr.blocksize(), 2460)
-        # assert max y-value at Bragg peak Q
-        _, imax = ws_corr.findY(ws_corr.readY(0).max())
-        self.assertAlmostEqual(ws_corr.readX(0)[imax], 3.272, delta=1e-2)  # (220) peak @ d = 1.920 Ang
+    def test_exec_nearest_interpolation(self):
+        ws_corr = PoldiAutoCorrelation(InputWorkspace=self.ws, OutputWorkspace="ws_corr_nearest", InterpolationMethod="Nearest", Version=6)
+        self._assert_auto_corr_workspace(ws_corr)
 
     def test_exec_cropped_wavelength_range(self):
         ws_corr = PoldiAutoCorrelation(InputWorkspace=self.ws, OutputWorkspace="ws_corr", WavelengthMin=2, WavelengthMax=4, Version=6)
@@ -45,6 +41,16 @@ def test_exec_cropped_wavelength_range(self):
         # assert bin width/number bins
         self.assertEqual(ws_corr.blocksize(), 1470)
 
+    def _assert_auto_corr_workspace(self, ws_corr):
+        # assert min/max Q
+        self.assertAlmostEqual(ws_corr.readX(0)[0], 1.5144, delta=1e-3)
+        self.assertAlmostEqual(ws_corr.readX(0)[-1], 9.0832, delta=1e-3)
+        # assert bin width/number bins
+        self.assertEqual(ws_corr.blocksize(), 2460)
+        # assert max y-value at Bragg peak Q
+        _, imax = ws_corr.findY(ws_corr.readY(0).max())
+        self.assertAlmostEqual(ws_corr.readX(0)[imax], 3.272, delta=1e-2)  # (220) peak @ d = 1.920 Ang
+
 
 if __name__ == "__main__":
     unittest.main()

Testing/Data/UnitTest/SI.cif.md5

@@ -0,0 +1 @@
+5545bb67c49f25aca8efd53b5b15bbff

docs/source/algorithms/PoldiAutoCorrelation-v6.rst

@@ -34,6 +34,11 @@ d-spacing is the inverse of the sum of the inverse intensities at each possible
 A high correlation intensity implies a consistently high intensity was measured at all arrival times due to all
 possible chopper openings/pulses.
 
+In general the arrival time of a neutron for a given d-spacing will not coincide with a bin-center in the data -
+some interpolation is required. There are two interpolation methods which can be specified with the
+``InterpolationMethod`` parameter: ``Linear`` (default - as in [1]_) and ``Nearest``. The ``Nearest`` method is a
+roughly a factor of 2-3 quicker but is potentially less accurate.
+
 Note by convention the correlation spectrum is converted from d-spacing into momentum transfer.
 
 Further details of the POLDI instrument and the reduction can be found in [1]_.

docs/source/release/v6.14.0/Diffraction/Engineering/New_features/39776.rst

@@ -0,0 +1,2 @@
+- Performance improvements to :ref:`algm-PoldiAutoCorrelation-v6` and function to simulate POLDI 2D workspace.
+- New parameter ``InterpolationMethod`` in :ref:`algm-PoldiAutoCorrelation-v6`. The default value ``Linear`` preserves existing behaviour (linear interpolation), an additonal option method ``Nearest`` has been added for faster execution.

scripts/Engineering/CMakeLists.txt

@@ -1,6 +1,6 @@
 # Tests for Engineering scripts/common code
 
-set(TEST_PY_FILES test/test_EnggUtils.py)
+set(TEST_PY_FILES test/test_EnggUtils.py test/test_pawley_utils.py)
 
 check_tests_valid(${CMAKE_CURRENT_SOURCE_DIR} ${TEST_PY_FILES})