Merge pull request #37975 from robertapplin/37973-fix-calculation-of-…

…eisf-error

Fix calculation of EISF errors in BayesQuasi algorithm

Framework/PythonInterface/plugins/algorithms/WorkflowAlgorithms/BayesQuasi.py

@@ -8,6 +8,7 @@
 
 import os
 import numpy as np
+from typing import Tuple
 
 from mantid.api import (
     AlgorithmFactory,
@@ -24,6 +25,14 @@
 from IndirectCommon import *
 
 
+def _calculate_eisf(
+    height: np.ndarray, height_error: np.ndarray, amplitude: np.ndarray, amplitude_error: np.ndarray
+) -> Tuple[np.ndarray, np.ndarray]:
+    eisf = height / (height + amplitude)
+    eisf_error = (1 / (height + amplitude) ** 2) * np.sqrt((amplitude * height_error) ** 2 + (height * amplitude_error) ** 2)
+    return eisf, eisf_error
+
+
 class BayesQuasi(PythonAlgorithm):
     _program = None
     _samWS = None
@@ -654,39 +663,26 @@ def C2Fw(self, sname):
             height_data, height_error = np.asarray(height_data), np.asarray(height_error)
 
             # calculate EISF and EISF error
-            total = height_data + amplitude_data
-            EISF_data = height_data / total
-            total_error = height_error**2 + amplitude_error**2
-            EISF_error = EISF_data * np.sqrt((height_error**2 / height_data**2) + (total_error / total**2))
+            eisf_data, eisf_error = _calculate_eisf(height_data, height_error, amplitude_data, amplitude_error)
 
             # interlace amplitudes and widths of the peaks
-            y.append(np.asarray(height_data))
-            for amp, width, EISF in zip(amplitude_data, width_data, EISF_data):
-                y.append(amp)
-                y.append(width)
-                y.append(EISF)
+            y.extend(height_data)
+            y.extend(np.hstack((amplitude_data, width_data, eisf_data)).flatten("F"))
 
             # interlace amplitude and width errors of the peaks
-            e.append(np.asarray(height_error))
-            for amp, width, EISF in zip(amplitude_error, width_error, EISF_error):
-                e.append(amp)
-                e.append(width)
-                e.append(EISF)
+            e.extend(height_error)
+            e.extend(np.hstack((amplitude_error, width_error, eisf_error)).flatten("F"))
 
             # create x data and axis names for each function
             axis_names.append("f" + str(nl) + ".f0." + "Height")
-            x.append(x_data)
+            x.extend(x_data)
             for j in range(1, nl + 1):
                 axis_names.append("f" + str(nl) + ".f" + str(j) + ".Amplitude")
-                x.append(x_data)
+                x.extend(x_data)
                 axis_names.append("f" + str(nl) + ".f" + str(j) + ".FWHM")
-                x.append(x_data)
+                x.extend(x_data)
                 axis_names.append("f" + str(nl) + ".f" + str(j) + ".EISF")
-                x.append(x_data)
-
-        x = np.asarray(x).flatten()
-        y = np.asarray(y).flatten()
-        e = np.asarray(e).flatten()
+                x.extend(x_data)
 
         s_api.CreateWorkspace(
             OutputWorkspace=output_workspace,
@@ -763,12 +759,12 @@ def _read_ql_file(self, file_name, nl):
                 # SIGIK
                 line = next(line_pointer)
                 amp = AMAX * math.sqrt(math.fabs(line[0]) + 1.0e-20)
-                block_amplitude_e.append(amp)
+                block_amplitude_e.append(abs(amp))
 
                 # SIGFK
                 line = next(line_pointer)
                 FWHM = 2.0 * HWHM * math.sqrt(math.fabs(line[0]) + 1.0e-20)
-                block_FWHM_e.append(FWHM)
+                block_FWHM_e.append(abs(FWHM))
 
             # append data from block
             amp_data.append(block_amplitude)
@@ -778,7 +774,7 @@ def _read_ql_file(self, file_name, nl):
             # append error values from block
             amp_error.append(block_amplitude_e)
             FWHM_error.append(block_FWHM_e)
-            height_error.append(block_height_e)
+            height_error.append(abs(block_height_e))
 
         return q_data, (amp_data, FWHM_data, height_data), (amp_error, FWHM_error, height_error)
 

Framework/PythonInterface/test/python/plugins/algorithms/WorkflowAlgorithms/BayesQuasiTest.py

@@ -6,8 +6,9 @@
 # SPDX - License - Identifier: GPL - 3.0 +
 import unittest
 import numpy as np
-from mantid.simpleapi import *
+from mantid.simpleapi import BayesQuasi, CreateWorkspace, DeleteWorkspace, Load
 from mantid.api import MatrixWorkspace, WorkspaceGroup
+from plugins.algorithms.WorkflowAlgorithms.BayesQuasi import _calculate_eisf
 
 
 class BayesQuasiTest(unittest.TestCase):
@@ -114,6 +115,18 @@ def test_run_with_zero_at_the_end_of_data(self):
         self.assertTrue(isinstance(result, MatrixWorkspace))
         self.assertTrue(isinstance(prob, MatrixWorkspace))
 
+    def test_calculate_eisf_with_positive_height_and_amplitudes(self):
+        height = np.array([0.2, 0.4, 0.6, 0.8, 1.0])
+        height_error = np.full(5, 0.02)
+        amplitude = np.array([0.5, 0.4, 0.3, 0.2, 0.1])
+        amplitude_error = np.full(5, 0.01)
+        eisf, eisf_error = _calculate_eisf(height, height_error, amplitude, amplitude_error)
+
+        expected_eisf = [0.28571429, 0.5, 0.66666667, 0.8, 0.90909091]
+        expected_eisf_error = [0.02081232, 0.01397542, 0.01047566, 0.00894427, 0.00842813]
+        self.assertTrue(np.allclose(expected_eisf, eisf))
+        self.assertTrue(np.allclose(expected_eisf_error, eisf_error))
+
     # --------------------------------Validate results------------------------------------------------
 
     def _validate_QLr_shape(self, result, probability, group):

docs/source/release/v6.11.0/Inelastic/Bugfixes/37973.rst

@@ -0,0 +1 @@
+- Fixed an invalid calculation of the EISF errors calculated on the Quasi tab of the :ref:`Inelastic Bayes fitting <interface-inelastic-bayes-fitting>` interface.