diff --git a/src/mvesuvio/analysis_reduction.py b/src/mvesuvio/analysis_reduction.py index 930b3284..1aa1ca59 100644 --- a/src/mvesuvio/analysis_reduction.py +++ b/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) diff --git a/src/mvesuvio/util/analysis_helpers.py b/src/mvesuvio/util/analysis_helpers.py index 93b74709..7b4cb07d 100644 --- a/src/mvesuvio/util/analysis_helpers.py +++ b/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) diff --git a/tests/data/analysis/benchmark/stored_spec_144-182_iter_1_GC_MS.npz b/tests/data/analysis/benchmark/stored_spec_144-182_iter_1_GC_MS.npz new file mode 100644 index 00000000..13f2d344 Binary files /dev/null and b/tests/data/analysis/benchmark/stored_spec_144-182_iter_1_GC_MS.npz differ diff --git a/tests/data/analysis/benchmark/stored_spec_144-182_iter_3_GC_MS.npz b/tests/data/analysis/benchmark/stored_spec_144-182_iter_3_GC_MS.npz deleted file mode 100644 index 433fda56..00000000 Binary files a/tests/data/analysis/benchmark/stored_spec_144-182_iter_3_GC_MS.npz and /dev/null differ diff --git a/tests/data/analysis/inputs/analysis_test.py b/tests/data/analysis/inputs/analysis_test.py index dbf660ec..2ae8121c 100644 --- a/tests/data/analysis/inputs/analysis_test.py +++ b/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 diff --git a/tests/data/analysis/inputs/analysis_test/output_files/spec_144-182_iter_1_GC_MS.npz b/tests/data/analysis/inputs/analysis_test/output_files/spec_144-182_iter_1_GC_MS.npz new file mode 100644 index 00000000..13f2d344 Binary files /dev/null and b/tests/data/analysis/inputs/analysis_test/output_files/spec_144-182_iter_1_GC_MS.npz differ diff --git a/tests/data/analysis/inputs/analysis_test/output_files/spec_144-182_iter_3_GC_MS.npz b/tests/data/analysis/inputs/analysis_test/output_files/spec_144-182_iter_3_GC_MS.npz deleted file mode 100644 index 2206a5f0..00000000 Binary files a/tests/data/analysis/inputs/analysis_test/output_files/spec_144-182_iter_3_GC_MS.npz and /dev/null differ diff --git a/tests/system/analysis/test_analysis.py b/tests/system/analysis/test_analysis.py index 64a7eb1e..fe31b224 100644 --- a/tests/system/analysis/test_analysis.py +++ b/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 diff --git a/tests/unit/analysis/test_analysis_helpers.py b/tests/unit/analysis/test_analysis_helpers.py index d6b20e5b..ab6f3211 100644 --- a/tests/unit/analysis/test_analysis_helpers.py +++ b/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() diff --git a/tests/unit/analysis/test_analysis_reduction.py b/tests/unit/analysis/test_analysis_reduction.py index f93c8fd1..43129e8b 100644 --- a/tests/unit/analysis/test_analysis_reduction.py +++ b/tests/unit/analysis/test_analysis_reduction.py @@ -25,27 +25,35 @@ def test_calculate_kinematics(self): [5, 5, 133.69, -0.2, 11.005, 0.5696], ]) dataX = np.array([ - [110.5, 111.5, 112.5, 113.5, 114.5], - [110.5, 111.5, 112.5, 113.5, 114.5], - [110.5, 111.5, 112.5, 113.5, 114.5], + [110.5, 111.5, 112.5, 113.5, 114.5, 115.5, 116.5], + [110.5, 111.5, 112.5, 113.5, 114.5, 115.5, 116.5], + [110.5, 111.5, 112.5, 113.5, 114.5, 115.5, 116.5], ]) alg._set_kinematic_arrays(dataX) - np.testing.assert_allclose(alg._v0, np.array([[0.12095, 0.11964, 0.11835, 0.11709, 0.11586], - [0.11988, 0.11858, 0.11732, 0.11608, 0.11487], - [0.11948, 0.1182 , 0.11694, 0.11571, 0.11451], - ]), atol=1e-4) - np.testing.assert_allclose(alg._E0, np.array([[76475.65823, 74821.94729, 73221.30191, 71671.47572, 70170.33998], - [75122.06378, 73511.83023, 71952.81999, 70442.88326, 68979.98182], - [74627.68443, 73033.23536, 71489.34475, 69993.89741, 68544.88764], - ]), atol=1e-4) - np.testing.assert_allclose(alg._deltaE, np.array([[71569.65823, 69915.94729, 68315.30191, 66765.47572, 65264.33998], - [70216.06378, 68605.83023, 67046.81999, 65536.88326, 64073.98182], - [69721.68443, 68127.23536, 66583.34475, 65087.89741, 63638.88764], - ]), atol=1e-4) - np.testing.assert_allclose(alg._deltaQ, np.array([[227.01905, 224.95887, 222.94348, 220.97148, 219.04148], - [226.21278, 224.18766, 222.20618, 220.26696, 218.36867], - [226.07138, 224.05877, 222.08939, 220.16185, 218.27485], - ]), atol=1e-4) + + # Kinematic arrays are calculated from exteded range for dataX + + np.testing.assert_allclose(alg._v0, np.array( + [[0.12949, 0.127984, 0.126513, 0.125075, 0.12367, 0.122295, 0.120951, 0.119636, 0.11835, 0.117091, 0.115858, 0.114651, 0.113469, 0.112311, 0.111176, 0.110064, 0.108974, 0.107906, 0.106858], + [0.128259, 0.126781, 0.125337, 0.123926, 0.122546, 0.121196, 0.119876, 0.118584, 0.11732, 0.116083, 0.114871, 0.113684, 0.112522, 0.111383, 0.110267, 0.109173, 0.108101, 0.107049, 0.106018], + [0.127807, 0.126339, 0.124905, 0.123504, 0.122133, 0.120792, 0.119481, 0.118198, 0.116942, 0.115712, 0.114508, 0.113329, 0.112174, 0.111042, 0.109933, 0.108845, 0.107779, 0.106734, 0.105709]] + ), atol=1e-6) + np.testing.assert_allclose(alg._E0, np.array( + [[87655.042194, 85628.1004, 83670.660984, 81779.582305, 79951.898249, 78184.806586, 76475.658228, 74821.947293, 73221.30191, 71671.475719, 70170.339977, 68715.876247, 67306.169611, 65939.402362, 64613.848144, 63327.866492, 62079.897763, 60868.458397, 59692.136511], + [85995.604273, 84025.62244, 82122.565654, 80283.436416, 78505.403187, 76785.789479, 75122.063781, 73511.830231, 71952.819995, 70442.883259, 68979.981824, 67562.182219, 66187.649309, 64854.640355, 63561.499488, 62306.652564, 61088.602368, 59905.924152, 58757.261461], + [85390.300476, 83440.962859, 81557.622084, 79737.332238, 77977.309967, 76274.923823, 74627.684428, 73033.235365, 71489.344748, 69993.897408, 68544.887641, 67140.412484, 65778.665459, 64457.930766, 63176.577872, 61933.056479, 60725.891829, 59553.680333, 58415.085488]] + ), atol=1e-6) + np.testing.assert_allclose(alg._deltaQ, np.array( + [[240.409925, 238.046957, 235.738868, 233.483784, 231.279914, 229.12555, 227.019055, 224.958865, 222.943485, 220.97148, 219.041479, 217.152166, 215.30228, 213.490611, 211.715999, 209.97733, 208.273532, 206.603579, 204.966482], + [239.365494, 237.045827, 234.779507, 232.564729, 230.399768, 228.282976, 226.212776, 224.18766, 222.206182, 220.26696, 218.368667, 216.510034, 214.689841, 212.906918, 211.160143, 209.448438, 207.770767, 206.126134, 204.513584], + [239.138969, 236.834796, 234.583415, 232.383047, 230.23199, 228.128619, 226.071376, 224.058774, 222.089387, 220.16185, 218.274853, 216.427142, 214.617513, 212.844811, 211.107927, 209.405797, 207.737397, 206.101744, 204.497891]] + ), atol=1e-6) + np.testing.assert_allclose(alg._deltaE, np.array( + [[82749.042194, 80722.1004, 78764.660984, 76873.582305, 75045.898249, 73278.806586, 71569.658228, 69915.947293, 68315.30191, 66765.475719, 65264.339977, 63809.876247, 62400.169611, 61033.402362, 59707.848144, 58421.866492, 57173.897763, 55962.458397, 54786.136511], + [81089.604273, 79119.62244, 77216.565654, 75377.436416, 73599.403187, 71879.789479, 70216.063781, 68605.830231, 67046.819995, 65536.883259, 64073.981824, 62656.182219, 61281.649309, 59948.640355, 58655.499488, 57400.652564, 56182.602368, 54999.924152, 53851.261461], + [80484.300476, 78534.962859, 76651.622084, 74831.332238, 73071.309967, 71368.923823, 69721.684428, 68127.235365, 66583.344748, 65087.897408, 63638.887641, 62234.412484, 60872.665459, 59551.930766, 58270.577872, 57027.056479, 55819.891829, 54647.680333, 53509.085488]] + ), atol=1e-6) + def test_set_y_space_arrays(self): alg = VesuvioAnalysisRoutine() @@ -63,7 +71,7 @@ def test_set_y_space_arrays(self): [70216.06378, 68605.83023, 67046.81999, 65536.88326, 64073.98182], [69721.68443, 68127.23536, 66583.34475, 65087.89741, 63638.88764], ]) - alg._set_y_space_arrays(dataX) + alg._set_y_space_arrays() y_spaces_expected = np.array( [[[-38.0885,-38.12637,-38.16414,-38.20185,-38.23948], [791.54274,779.75738,768.21943,756.92089,745.85436], @@ -99,16 +107,16 @@ def test_neutron_compton_profile_fit_function(self): alg._set_kinematic_arrays(alg._dataX) alg._set_gaussian_resolution() alg._set_lorentzian_resolution() - alg._set_y_space_arrays(alg._dataX) + alg._set_y_space_arrays() example_fit_parameters = np.array([7.1, 5.05, 0.02, 0.22, 12.71, 1.0, 0.0, 8.76, -1.1, 0.3, 13.897, 0.64]) alg._row_being_fit = 0 NCP, FSE = alg._neutron_compton_profiles(example_fit_parameters) - expected_NCP = np.array([ - [0.00004, 0.000059, 0.00009, 0.000138, 0.000216, 0.000336, 0.000648, 0.000944, 0.001349, 0.001891, 0.002596, 0.003491, 0.004595, 0.005921, 0.007464, 0.009194, 0.011051, 0.012936, 0.014711, 0.016212, 0.017269, 0.017734, 0.017518, 0.016611, 0.01509, 0.013112, 0.010879, 0.008601, 0.006465, 0.004605, 0.003092, 0.001939, 0.001119, 0.000574, 0.00024, 0.000055, -0.000033, -0.000064, -0.000064, -0.00005, -0.000033, -0.000017, -0.000004, 0.000005, 0.000012, 0.000016, 0.000018, 0.000027, 0.000025, 0.000023, 0.000022, 0.000021, 0.00002], - [0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.000001, 0.000001, 0.000001, 0.000001, 0.000001, 0.000001, 0.000001, 0.000002, 0.000002, 0.000003, 0.000004, 0.000005, 0.000006, 0.000008, 0.000011, 0.000043, 0.000149, 0.00034, 0.000926, 0.002246, 0.003149, 0.00241, 0.000926, 0.00007, 0.000067, 0.00003, 0.000019, 0.000014, 0.000011, 0.000009], - [0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.], - [ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.000001, 0.000001, 0.000001, 0.000001, 0.000001, 0.000001, 0.000001, 0.000002, 0.000002, 0.000002, 0.000003, 0.000004, 0.000005, 0.000006, 0.000008, 0.000011, 0.000017, 0.000021, 0.000052, 0.000152, 0.000455, 0.003186, 0.006654, 0.003274, 0.000527, 0.000117, 0.000058, 0.000037, 0.000026, 0.000019, 0.000015] - ]) + expected_NCP = np.array( + [[0.000055, 0.000083, 0.000126, 0.000192, 0.000292, 0.000438, 0.000648, 0.000944, 0.001349, 0.001891, 0.002596, 0.003491, 0.004595, 0.005921, 0.007464, 0.009194, 0.011051, 0.012936, 0.014711, 0.016212, 0.017269, 0.017734, 0.017518, 0.016611, 0.01509, 0.013112, 0.010879, 0.008601, 0.006465, 0.004605, 0.003092, 0.001939, 0.001119, 0.000574, 0.00024, 0.000055, -0.000033, -0.000064, -0.000064, -0.00005, -0.000033, -0.000017, -0.000004, 0.000005, 0.000012, 0.000016, 0.000018, 0.000019, 0.00002, 0.00002, 0.000019, 0.000019, 0.000018], + [0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.000001, 0.000001, 0.000001, 0.000001, 0.000001, 0.000001, 0.000001, 0.000002, 0.000002, 0.000003, 0.000004, 0.000005, 0.000006, 0.000008, 0.000011, 0.000043, 0.000149, 0.00034, 0.000926, 0.002246, 0.003149, 0.00241, 0.000926, 0.00007, -0.000056, 0.000005, 0.00002, 0.000014, 0.00001, 0.000008], + [0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., ], + [0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.000001, 0.000001, 0.000001, 0.000001, 0.000001, 0.000001, 0.000001, 0.000002, 0.000002, 0.000002, 0.000003, 0.000004, 0.000005, 0.000006, 0.000008, 0.000011, 0.000017, 0.000021, 0.000052, 0.000152, 0.000455, 0.003186, 0.006654, 0.003274, 0.000527, 0.00001, 0.000045, 0.000045, 0.000025, 0.000019, 0.000015]] + ) np.testing.assert_allclose(NCP, expected_NCP, atol=1e-6) @@ -123,7 +131,7 @@ def test_error_function(self): alg._set_kinematic_arrays(alg._dataX) alg._set_gaussian_resolution() alg._set_lorentzian_resolution() - alg._set_y_space_arrays(alg._dataX) + alg._set_y_space_arrays() example_fit_parameters = np.array([7.1, 5.05, 0.02, 0.22, 12.71, 1.0, 0.0, 8.76, -1.1, 0.3, 13.897, 0.64]) alg._row_being_fit = 0 @@ -139,8 +147,8 @@ def test_error_function(self): alg._dataY[:, 5:15] = 0 chi2_with_errors_with_tof_masked = alg._error_function(example_fit_parameters) - self.assertEqual(chi2_without_errors, chi2_with_errors * 0.0015**2) - self.assertEqual(chi2_without_errors, 1.762715850011478e-05) + self.assertAlmostEqual(chi2_without_errors, chi2_with_errors * 0.0015**2, places=15) + self.assertEqual(chi2_without_errors, 1.7627158500114776e-05) self.assertEqual(chi2_with_errors, 7.834292666717679) self.assertEqual(chi2_with_errors_with_tof_masked, 5.5864483595799195) @@ -163,7 +171,7 @@ def test_fit_neutron_compton_profiles_to_row(self): alg._set_kinematic_arrays(alg._dataX) alg._set_gaussian_resolution() alg._set_lorentzian_resolution() - alg._set_y_space_arrays(alg._dataX) + alg._set_y_space_arrays() synthetic_ncp = np.zeros_like(alg._dataX) example_fit_parameters = np.array([ @@ -241,7 +249,7 @@ def test_fit_neutron_compton_profiles_to_row_with_masked_tof(self): alg._set_kinematic_arrays(alg._dataX) alg._set_gaussian_resolution() alg._set_lorentzian_resolution() - alg._set_y_space_arrays(alg._dataX) + alg._set_y_space_arrays() # Create mock for storing ncp total result ncp_array_masked = np.zeros_like(alg._dataY) @@ -254,28 +262,28 @@ def test_fit_neutron_compton_profiles_to_row_with_masked_tof(self): fit_parameters_masked = alg._fit_parameters.copy() # Now cut range so that zeros are not part of dataY - # (Still need to leave a padding with 6 zeros due to numerical derivative in ncp) - alg._dataY = alg._dataY[:, cut_off_idx - 6:].reshape(1, -1) - alg._dataX = alg._dataX[:, cut_off_idx - 6:].reshape(1, -1) - alg._dataE = alg._dataE[:, cut_off_idx - 6:].reshape(1, -1) + alg._dataY = alg._dataY[:, cut_off_idx:].reshape(1, -1) + alg._dataX = alg._dataX[:, cut_off_idx:].reshape(1, -1) + alg._dataE = alg._dataE[:, cut_off_idx:].reshape(1, -1) # Run methods that depend on dataX alg._set_kinematic_arrays(alg._dataX) alg._set_gaussian_resolution() alg._set_lorentzian_resolution() - alg._set_y_space_arrays(alg._dataX) + alg._set_y_space_arrays() - ncp_array_cut = ncp_array_masked[:, cut_off_idx - 6:].reshape(1, -1) + ncp_array_cut = ncp_array_masked[:, cut_off_idx:].reshape(1, -1) alg._fit_profiles_workspaces = {"total": MagicMock(dataY=lambda arg: ncp_array_cut), "1": MagicMock()} alg._fit_neutron_compton_profiles_to_row() fit_parameters_cut = alg._fit_parameters.copy() np.testing.assert_allclose(fit_parameters_cut[:, 1:-2], fit_parameters_masked[:, 1:-2], atol=1e-6) - np.testing.assert_allclose(ncp_array_masked[:, cut_off_idx-6:], ncp_array_cut, atol=1e-6) + np.testing.assert_allclose(ncp_array_masked[:, cut_off_idx:], ncp_array_cut, atol=1e-6) def test_set_gaussian_resolution(self): + alg = VesuvioAnalysisRoutine() alg._instrument_params = np.array([ [144, 144, 54.1686, -0.41, 11.005, 0.720184], @@ -285,16 +293,16 @@ def test_set_gaussian_resolution(self): alg._masses = np.array([1, 12]) alg._dataX = np.array([ - [227., 233., 239.], - [227., 233., 239.], + [110.5, 111.5, 112.5, 113.5, 114.5, 115.5, 116.5], + [110.5, 111.5, 112.5, 113.5, 114.5, 115.5, 116.5], ]) alg._set_kinematic_arrays(alg._dataX) alg._set_gaussian_resolution() expected_gaussian_resolution = np.array( - [[[1.066538, 1.046284, 1.028226], - [7.070845, 6.854225, 6.662274]], - [[1.076188, 1.055946, 1.037921], - [7.219413, 7.003598, 6.812526]]] + [[[3.709857, 3.62982, 3.552662, 3.478251, 3.406461, 3.337172, 3.270274, 3.205661, 3.143232, 3.082894, 3.024556, 2.968136, 2.913551, 2.860728, 2.809592, 2.760078, 2.712118, 2.665653, 2.620622], + [31.128982, 30.438141, 29.771605, 29.128243, 28.50699, 27.906841, 27.326849, 26.766118, 26.223805, 25.699109, 25.191275, 24.699587, 24.223369, 23.761977, 23.314804, 22.881272, 22.460833, 22.052966, 21.657176]], + [[3.716323, 3.636666, 3.559864, 3.485785, 3.414305, 3.345307, 3.278679, 3.214318, 3.152124, 3.092005, 3.033872, 2.977641, 2.923233, 2.870573, 2.81959, 2.770216, 2.722386, 2.67604, 2.631121], + [31.204563, 30.517381, 29.854298, 29.214196, 28.596021, 27.998779, 27.421535, 26.863402, 26.323544, 25.801171, 25.295536, 24.805929, 24.33168, 23.872154, 23.426748, 22.99489, 22.576037, 22.169674, 21.77531, ]]] ) np.testing.assert_allclose(expected_gaussian_resolution, alg._gaussian_resolution, atol=1e-6) @@ -309,16 +317,16 @@ def test_set_lorentzian_resolution(self): alg._masses = np.array([1, 12]) alg._dataX = np.array([ - [227., 233., 239.], - [227., 233., 239.], + [110.5, 111.5, 112.5, 113.5, 114.5, 115.5, 116.5], + [110.5, 111.5, 112.5, 113.5, 114.5, 115.5, 116.5], ]) alg._set_kinematic_arrays(alg._dataX) alg._set_lorentzian_resolution() expected_lorentzian_resolution = np.array( - [[[0.119899, 0.119167, 0.118724], - [1.346864, 1.356234, 1.366692]], - [[0.124002, 0.123212, 0.122712], - [1.376647, 1.387298, 1.399038]]] + [[[0.313762, 0.308099, 0.302633, 0.297353, 0.292252, 0.287321, 0.282552, 0.277939, 0.273474, 0.269152, 0.264966, 0.26091, 0.256979, 0.253168, 0.249472, 0.245887, 0.242408, 0.239031, 0.235751], + [2.410374, 2.368446, 2.32816, 2.289439, 2.252211, 2.21641, 2.18197, 2.148833, 2.11694, 2.086238, 2.056676, 2.028205, 2.000781, 1.974359, 1.948898, 1.92436, 1.900708, 1.877906, 1.855921]], + [[0.317283, 0.311673, 0.306255, 0.301023, 0.295966, 0.291078, 0.28635, 0.281775, 0.277348, 0.273061, 0.268909, 0.264886, 0.260986, 0.257205, 0.253538, 0.24998, 0.246526, 0.243174, 0.239918], + [2.41149, 2.370058, 2.330248, 2.291987, 2.255202, 2.219827, 2.1858, 2.15306, 2.121553, 2.091224, 2.062023, 2.033904, 2.00682, 1.980729, 1.955591, 1.931366, 1.908019, 1.885515, 1.863821]]] ) np.testing.assert_allclose(expected_lorentzian_resolution, alg._lorentzian_resolution, atol=1e-6)