Improve third derivative (#148)

Previous third derivative result was padded with 6 zeros on each side.
This improvement extends the range of the kinematic, y-space and resolution arrays by 6 indices on each side, so that when the third derivative is taken, the size of the resulting array is 12 indices less than the input (6 on each side).

So in essence this just means that I am extending the range before taking the derivative so that after the derivative the size matches the size before the extension. This makes sure the derivative is continuous and does not including zero padding, which contributes to a smoother ncp profile.

Added and updated unit tests. Updated benchmark for system tests.

src/mvesuvio/analysis_reduction.py

@@ -10,9 +10,10 @@
                             CloneWorkspace, DeleteWorkspace, VesuvioCalculateGammaBackground, \
                             VesuvioCalculateMS, Scale, RenameWorkspace, Minus, CreateSampleShape, \
                             VesuvioThickness, Integration, Divide, Multiply, DeleteWorkspaces, \
-                            CreateWorkspace, CreateSampleWorkspace
+                            CreateWorkspace
 
-from mvesuvio.util.analysis_helpers import numericalThirdDerivative, load_resolution, load_instrument_params
+from mvesuvio.util.analysis_helpers import numerical_third_derivative, load_resolution, load_instrument_params, \
+                                            extend_range_of_array
 
 np.set_printoptions(suppress=True, precision=4, linewidth=200)
 
@@ -207,8 +208,7 @@ def _update_workspace_data(self):
         self._set_kinematic_arrays(self._dataX)
         self._set_gaussian_resolution()
         self._set_lorentzian_resolution()
-        # Calculate y space after kinematics
-        self._set_y_space_arrays(self._dataX)
+        self._set_y_space_arrays()
 
         self._fit_parameters = np.zeros((len(self._dataY), 3 * self._profiles_table.rowCount() + 3))
         self._row_being_fit = 0 
@@ -331,6 +331,10 @@ def _fit_neutron_compton_profiles(self):
 
 
     def _set_kinematic_arrays(self, dataX):
+
+        # Extend range due to third derivative cutting edges
+        dataX = extend_range_of_array(dataX, 6)
+
         det, plick, angle, T0, L0, L1 = np.hsplit(self._instrument_params, 6)  # each is of len(dataX)
 
         # T0 is electronic delay due to instruments
@@ -349,9 +353,8 @@ def _set_kinematic_arrays(self, dataX):
         return
 
 
-    def _set_y_space_arrays(self, dataX):
-        # Prepare arrays to broadcast
-        dataX = dataX[np.newaxis, :, :]
+    def _set_y_space_arrays(self):
+
         delta_Q = self._deltaQ[np.newaxis, :, :]
         delta_E = self._deltaE[np.newaxis, :, :]
         masses = self._masses.reshape(-1, 1, 1)
@@ -543,27 +546,26 @@ def _neutron_compton_profiles(self, pars):
         Neutron Compton Profile distribution on TOF space for a single spectrum. 
         Calculated from kinematics, J(y) and resolution functions.
         """
-
         intensities = pars[::3].reshape(-1, 1)
         widths = pars[1::3].reshape(-1, 1)
         centers = pars[2::3].reshape(-1, 1)
         masses = self._masses.reshape(-1, 1)
 
-        E0 = self._E0[self._row_being_fit]
-        deltaQ = self._deltaQ[self._row_being_fit]
-
         gaussian_width = self._get_gaussian_resolution(centers)
         lorentzian_width = self._get_lorentzian_resolution(centers)
         total_gaussian_width = np.sqrt(widths**2 + gaussian_width**2)
 
         JOfY = scipy.special.voigt_profile(self._y_space_arrays[self._row_being_fit] - centers, total_gaussian_width, lorentzian_width)
 
-        FSE = (
-            -numericalThirdDerivative(self._y_space_arrays[self._row_being_fit], JOfY)
-            * widths**4
-            / deltaQ
-            * 0.72
-        )
+        # Third derivative cuts edges of array by 6 indices
+        JOfY_third_derivative = numerical_third_derivative(self._y_space_arrays[self._row_being_fit], JOfY)
+
+        deltaQ = self._deltaQ[self._row_being_fit, 6: -6]
+        E0 = self._E0[self._row_being_fit, 6: -6]
+        JOfY = JOfY[:, 6:-6]
+
+        FSE = - JOfY_third_derivative * widths**4 / deltaQ * 0.72
+
         NCP = intensities * (JOfY+FSE) * E0 * E0 ** (-0.92) * masses / deltaQ
         FSE = intensities * FSE * E0 * E0 ** (-0.92) * masses / deltaQ
         return NCP, FSE
@@ -572,6 +574,7 @@ def _neutron_compton_profiles(self, pars):
     def _get_gaussian_resolution(self, centers):
         proximity_to_y_centers = np.abs(self._y_space_arrays[self._row_being_fit] - centers)
         gauss_resolution = self._gaussian_resolution[self._row_being_fit]
+        assert proximity_to_y_centers.shape==gauss_resolution.shape
         return np.take_along_axis(gauss_resolution, proximity_to_y_centers.argmin(axis=1, keepdims=True), axis=1)
 
 
@@ -626,6 +629,7 @@ def _set_gaussian_resolution(self):
     def _get_lorentzian_resolution(self, centers):
         proximity_to_y_centers = np.abs(self._y_space_arrays[self._row_being_fit] - centers)
         lorentzian_resolution = self._lorentzian_resolution[self._row_being_fit]
+        assert proximity_to_y_centers.shape==lorentzian_resolution.shape
         return np.take_along_axis(lorentzian_resolution, proximity_to_y_centers.argmin(axis=1, keepdims=True), axis=1)
 
 

src/mvesuvio/util/analysis_helpers.py

@@ -237,7 +237,7 @@ def extractWS(ws):
     return ws.extractX(), ws.extractY(), ws.extractE()
 
 
-def numericalThirdDerivative(x, y):
+def numerical_third_derivative(x, y):
     k6 = (- y[:, 12:] + y[:, :-12]) * 1
     k5 = (+ y[:, 11:-1] - y[:, 1:-11]) * 24
     k4 = (- y[:, 10:-2] + y[:, 2:-10]) * 192
@@ -248,11 +248,7 @@ def numericalThirdDerivative(x, y):
     dev = k1 + k2 + k3 + k4 + k5 + k6
     dev /= np.power(x[:, 7:-5] - x[:, 6:-6], 3)
     dev /= 12**3
-
-    derivative = np.zeros_like(y)
-    # Padded with zeros left and right to return array with same shape
-    derivative[:, 6:-6] = dev
-    return derivative
+    return dev
 
 
 def load_resolution(instrument_params):
@@ -368,3 +364,9 @@ def calculate_h_ratio(means_table):
 
     return sorted_intensities[0] / sorted_intensities[1] 
 
+
+def extend_range_of_array(arr, n_columns):
+    arr = arr.copy()
+    left_extend = arr[:, :n_columns] + (arr[:, 0] - arr[:, n_columns]).reshape(-1, 1)
+    right_extend = arr[:, -n_columns:] + (arr[:, -1] - arr[:, -n_columns-1]).reshape(-1, 1)
+    return np.concatenate([left_extend, arr, right_extend], axis=-1)

tests/data/analysis/inputs/analysis_test.py

@@ -79,7 +79,7 @@ class ForwardAnalysisInputs(SampleParameters):
             [-3, 1],
         ]
     constraints = ()
-    noOfMSIterations = 3  # 4
+    noOfMSIterations = 1  # 4
     MSCorrectionFlag = True
     GammaCorrectionFlag = True
 

tests/system/analysis/test_analysis.py

@@ -37,7 +37,7 @@ def _run(cls):
     def _load_benchmark_results(cls):
         benchmarkPath = Path(__file__).absolute().parent.parent.parent / "data" / "analysis" / "benchmark"
         benchmarkResults = np.load(
-            str(benchmarkPath / "stored_spec_144-182_iter_3_GC_MS.npz")
+            str(benchmarkPath / "stored_spec_144-182_iter_1_GC_MS.npz")
         )
         cls._benchmarkResults = benchmarkResults
 

tests/unit/analysis/test_analysis_helpers.py

@@ -1,10 +1,11 @@
 import unittest
 import numpy as np
+import scipy
 import dill
 import numpy.testing as nptest
 from mock import MagicMock
 from mvesuvio.util.analysis_helpers import extractWS, _convert_dict_to_table,  \
-    fix_profile_parameters, calculate_h_ratio
+    fix_profile_parameters, calculate_h_ratio, extend_range_of_array, numerical_third_derivative
 from mantid.simpleapi import CreateWorkspace, DeleteWorkspace
 
 
@@ -126,5 +127,27 @@ def test_conversion_of_constraints(self):
         self.assertEqual(converted_constraints[1]['fun']([0, 0, 0, 2, 0, 0, 1]), 2-0.7234)
 
 
+    def test_extend_range_of_array_for_increasing_range(self):
+        x = np.arange(10)
+        x = np.vstack([x, 2*x])
+        x_extended = extend_range_of_array(x, 5)
+        np.testing.assert_array_equal(x_extended, np.vstack([np.arange(-5, 15, 1), np.arange(-10, 30, 2)]))
+
+
+    def test_extend_range_of_array_for_decreasing_range(self):
+        x = np.linspace(-5, 5, 21)
+        x = np.vstack([x, 2*x])
+        x_extended = extend_range_of_array(x, 5)
+        np.testing.assert_array_equal(x_extended, np.vstack([np.linspace(-7.5, 7.5, 31), np.linspace(-15, 15, 31)]))
+
+
+    def test_numerical_third_derivative(self):
+        x= np.linspace(-20, 20, 300)    # Workspaces are about 300 points of range
+        x = np.vstack([x, 2*x])
+        y = scipy.special.voigt_profile(x, 5, 5)
+        numerical_derivative = numerical_third_derivative(x, y)
+        expected_derivative = np.array([np.gradient(np.gradient(np.gradient(y_i, x_i), x_i), x_i)[6: -6] for y_i, x_i in zip(y, x) ])
+        np.testing.assert_allclose(numerical_derivative, expected_derivative, atol=1e-6)
+
 if __name__ == "__main__":
     unittest.main()