Create function to update input profiles table

Transformed existing function to create a function that takes
in an outputs table containing the means at the end of the routine and
passes these to the inputs profiles table as the initial parameters.
The lightest mass is always unfixed (bounds are relaxed).
The h ratio is used depending on whether the sample contains hydrogen.

Author: GuiMacielPereira

src/mvesuvio/analysis_reduction.py

@@ -516,21 +516,23 @@ def _create_means_table(self):
             OutputWorkspace=self._workspace_being_fit.name() + "_means"
         )
         table.addColumn(type="str", name="label")
+        table.addColumn(type="float", name="mass")
         table.addColumn(type="float", name="mean_width")
         table.addColumn(type="float", name="std_width")
         table.addColumn(type="float", name="mean_intensity")
         table.addColumn(type="float", name="std_intensity")
 
         print("\nCreated Table with means and std:")
         print("\nMass    Mean \u00B1 Std Widths    Mean \u00B1 Std Intensities\n")
-        for label, mean_width, std_width, mean_intensity, std_intensity in zip(
+        for label, mass, mean_width, std_width, mean_intensity, std_intensity in zip(
             self._profiles_table.column("label"),
+            self._masses,
             self._mean_widths,
             self._std_widths,
             self._mean_intensity_ratios,
             self._std_intensity_ratios,
         ):
-            table.addRow([label, mean_width, std_width, mean_intensity, std_intensity])
+            table.addRow([label, mass, mean_width, std_width, mean_intensity, std_intensity])
             print(f"{label:5s}  {mean_width:10.5f} \u00B1 {std_width:7.5f}  \
                     {mean_intensity:10.5f} \u00B1 {std_intensity:7.5f}\n")
 

src/mvesuvio/analysis_routines.py

@@ -4,12 +4,11 @@
 from mantid.api import AlgorithmFactory, AlgorithmManager
 import numpy as np
 
-from mvesuvio.util.analysis_helpers import loadRawAndEmptyWsFromUserPath, cropAndMaskWorkspace
+from mvesuvio.util.analysis_helpers import loadRawAndEmptyWsFromUserPath, cropAndMaskWorkspace, NeutronComptonProfile
 from mvesuvio.analysis_reduction import AnalysisRoutine
-from mvesuvio.analysis_reduction import NeutronComptonProfile
 from tests.testhelpers.calibration.algorithms import create_algorithm
 
-def _create_analysis_object_from_current_interface(IC, running_tests=False):
+def _create_analysis_object_from_current_interface(IC, running_tests=False, overwrite_profiles_table=None):
     ws = loadRawAndEmptyWsFromUserPath(
         userWsRawPath=IC.userWsRawPath,
         userWsEmptyPath=IC.userWsEmptyPath,
@@ -27,17 +26,20 @@ def _create_analysis_object_from_current_interface(IC, running_tests=False):
         maskTOFRange=IC.maskTOFRange
     )
 
-    profiles = []
-    for mass, intensity, width, center, intensity_bound, width_bound, center_bound in zip(
-        IC.masses, IC.initPars[::3], IC.initPars[1::3], IC.initPars[2::3],
-        IC.bounds[::3], IC.bounds[1::3], IC.bounds[2::3]
-    ):
-        profiles.append(NeutronComptonProfile(
-            label=str(mass), mass=mass, intensity=intensity, width=width, center=center,
-            intensity_bounds=list(intensity_bound), width_bounds=list(width_bound), center_bounds=list(center_bound)
-        ))
-
-    profiles_table = create_profiles_table(cropedWs.name()+"_initial_parameters", profiles)
+    if overwrite_profiles_table:
+        profiles_table = overwrite_profiles_table
+    else:
+        profiles = []
+        for mass, intensity, width, center, intensity_bound, width_bound, center_bound in zip(
+            IC.masses, IC.initPars[::3], IC.initPars[1::3], IC.initPars[2::3],
+            IC.bounds[::3], IC.bounds[1::3], IC.bounds[2::3]
+        ):
+            profiles.append(NeutronComptonProfile(
+                label=str(mass), mass=mass, intensity=intensity, width=width, center=center,
+                intensity_bounds=list(intensity_bound), width_bounds=list(width_bound), center_bounds=list(center_bound)
+            ))
+
+        profiles_table = create_profiles_table(cropedWs.name()+"_initial_parameters", profiles)
 
     kwargs = {
         "InputWorkspace": cropedWs,
@@ -90,40 +92,6 @@ def create_profiles_table(name, profiles: list[NeutronComptonProfile]):
     return table
 
 
-def set_initial_profiles_from(self, source: 'AnalysisRoutine'):
-    
-    # Set intensities
-    for p in self._profiles.values():
-        if np.isclose(p.mass, 1, atol=0.1):    # Hydrogen present
-            p.intensity = source._h_ratio * source._get_lightest_profile().mean_intensity
-            continue
-        p.intensity = source.profiles[p.label].mean_intensity
-
-    # Normalise intensities
-    sum_intensities = sum([p.intensity for p in self._profiles.values()])
-    for p in self._profiles.values():
-        p.intensity /= sum_intensities
-        
-    # Set widths
-    for p in self._profiles.values():
-        try:
-            p.width = source.profiles[p.label].mean_width
-        except KeyError:
-            continue
-
-    # Fix all widths except lightest mass
-    for p in self._profiles.values():
-        if p == self._get_lightest_profile():
-            continue
-        p.width_bounds = [p.width, p.width]
-
-    return
-
-def _get_lightest_profile(self):
-    profiles = [p for p in self._profiles.values()]
-    masses = [p.mass for p in self._profiles.values()]
-    return profiles[np.argmin(masses)]
-
 def runIndependentIterativeProcedure(IC, clearWS=True, running_tests=False):
     """
     Runs the iterative fitting of NCP, cleaning any previously stored workspaces.
@@ -137,7 +105,6 @@ def runIndependentIterativeProcedure(IC, clearWS=True, running_tests=False):
 
     alg = _create_analysis_object_from_current_interface(IC, running_tests=running_tests)
     alg.execute()
-    means_table = alg.getPropertyValue("OutputMeansTable")
     return alg
 
 
@@ -212,12 +179,18 @@ def createTableWSHRatios():
 
 def runJoint(bckwdIC, fwdIC):
 
-    backRoutine = _create_analysis_object_from_current_interface(bckwdIC)
-    frontRoutine = _create_analysis_object_from_current_interface(fwdIC)
+    back_alg= _create_analysis_object_from_current_interface(bckwdIC)
+    front_alg= _create_analysis_object_from_current_interface(fwdIC)
 
-    backRoutine.execute()
-    frontRoutine.set_initial_profiles_from(backRoutine)
-    frontRoutine.execute()
+    back_alg.execute()
+    incoming_means_table = back_alg.getProperty("OutputMeansTable").value
+    receiving_profiles_table = front_alg.getProperty("InputProfiles").value
+
+    fixed_profiles_table = fix_profiles(incoming_means_table, receiving_profiles_table)
+
+    front_alg.setProperty("InputProfiles", fixed_profiles_table)
+
+    front_alg.execute()
     return
 
 

src/mvesuvio/util/analysis_helpers.py

@@ -6,6 +6,26 @@
 
 from mvesuvio.analysis_fitting import passDataIntoWS
 
+from dataclasses import dataclass
+
+@dataclass(frozen=False)
+class NeutronComptonProfile:
+    label: str
+    mass: float
+
+    intensity: float
+    width: float
+    center: float
+
+    intensity_bounds: list[float, float]
+    width_bounds: list[float, float]
+    center_bounds: list[float, float]
+
+    mean_intensity: float = None
+    mean_width: float = None
+    mean_center: float = None
+
+
 
 def loadRawAndEmptyWsFromUserPath(userWsRawPath, userWsEmptyPath, 
                                   tofBinning, name, scaleRaw, scaleEmpty, subEmptyFromRaw):
@@ -152,3 +172,62 @@ def createWS(dataX, dataY, dataE, wsName, parentWorkspace=None):
         ParentWorkspace=parentWorkspace
     )
     return ws
+
+
+def fix_profile_parameters(incoming_means_table, receiving_profiles_table, h_ratio):
+    means_dict = _convert_table_to_dict(incoming_means_table)
+    profiles_dict = _convert_table_to_dict(receiving_profiles_table)
+
+    # Set intensities
+    for p in profiles_dict.values():
+        if np.isclose(p['mass'], 1, atol=0.1):    # Hydrogen present
+            p['intensity'] = h_ratio * _get_lightest_profile(means_dict)['mean_intensity']
+            continue
+        p['intensity'] = means_dict[p['label']]['mean_intensity']
+
+    # Normalise intensities
+    sum_intensities = sum([p['intensity'] for p in profiles_dict.values()])
+    for p in profiles_dict.values():
+        p['intensity'] /= sum_intensities
+        
+    # Set widths
+    for p in profiles_dict.values():
+        try:
+            p['width'] = means_dict[p['label']]['mean_width']
+        except KeyError:
+            continue
+
+    # Fix all widths except lightest mass
+    for p in profiles_dict.values():
+        if p == _get_lightest_profile(profiles_dict):
+            continue
+        p['width_bounds'] = str([p['width'] , p['width']])
+
+    result_profiles_table = _convert_dict_to_table(profiles_dict)
+    return result_profiles_table
+
+
+def _convert_table_to_dict(table):
+    result_dict = {}
+    for i in range(table.rowCount()):
+        row_dict = table.row(i) 
+        result_dict[row_dict['label']] = row_dict
+    return result_dict
+
+
+def _convert_dict_to_table(m_dict):
+    table = CreateEmptyTableWorkspace()
+    for p in m_dict.values():
+        if table.columnCount() == 0:
+            for key, value in p.items():
+                value_type = 'str' if isinstance(value, str) else 'float'
+                table.addColumn(value_type, key)
+
+        table.addRow(p)
+    return table
+
+
+def _get_lightest_profile(p_dict):
+    profiles = [p for p in p_dict.values()]
+    masses = [p['mass'] for p in p_dict.values()]
+    return profiles[np.argmin(masses)]

tests/unit/analysis/test_analysis_functions.py

@@ -0,0 +1,109 @@
+import unittest
+import numpy as np
+import numpy.testing as nptest
+from mock import MagicMock
+from mvesuvio.util.analysis_helpers import extractWS, _convert_dict_to_table, fix_profile_parameters
+from mantid.simpleapi import CreateWorkspace, DeleteWorkspace
+
+
+class TestAnalysisHelpers(unittest.TestCase):
+    def setUp(self):
+        pass
+
+    def test_extract_ws(self):
+        data = [1, 2, 3]
+        ws = CreateWorkspace(DataX=data, DataY=data, DataE=data, NSpec=1, UnitX="some_unit")
+
+        dataX, dataY, dataE = extractWS(ws)
+        nptest.assert_array_equal([data], dataX)
+        nptest.assert_array_equal([data], dataY)
+        nptest.assert_array_equal([data], dataE)
+
+        DeleteWorkspace(ws)
+
+
+    def test_convert_dict_to_table(self):
+        d = {'H': {'label': 'H', 'mass': 1, 'intensity': 1}}
+        table = _convert_dict_to_table(d)
+        self.assertEqual(['label', 'mass', 'intensity'], table.getColumnNames())
+        self.assertEqual({'label': 'H', 'mass': 1, 'intensity': 1}, table.row(0))
+
+
+    def test_fix_profile_parameters_with_H(self):
+        means_table_mock = MagicMock()
+        means_table_mock.rowCount.return_value = 3 
+        means_table_mock.row.side_effect = [
+            {'label': '16.0', 'mass': 16.0, 'mean_width': 8.974, 'std_width': 1.401, 'mean_intensity': 0.176, 'std_intensity': 0.08722},
+            {'label': '27.0', 'mass': 27.0, 'mean_width': 15.397, 'std_width': 1.131, 'mean_intensity': 0.305, 'std_intensity': 0.04895},
+            {'label': '12.0', 'mass': 12.0, 'mean_width': 13.932, 'std_width': 0.314, 'mean_intensity': 0.517, 'std_intensity': 0.09531}
+        ]
+        profiles_table_mock = MagicMock()
+        profiles_table_mock.rowCount.return_value = 4
+        profiles_table_mock.row.side_effect = [
+            {'label': '1.0079', 'mass': 1.0078, 'intensity': 1.0, 'intensity_bounds': '[0, None]', 'width': 4.699, 'width_bounds': '[3, 6]', 'center': 0.0, 'center_bounds': '[-3, 1]'},
+            {'label': '16.0', 'mass': 16.0, 'intensity': 1.0, 'intensity_bounds': '[0, None]', 'width': 12, 'width_bounds': '[0, None]', 'center': 0.0, 'center_bounds': '[-3, 1]'},
+            {'label': '12.0', 'mass': 12.0, 'intensity': 1.0, 'intensity_bounds': '[0, None]', 'width': 8, 'width_bounds': '[0, None]', 'center': 0.0, 'center_bounds': '[-3, 1]'},
+            {'label': '27.0', 'mass': 27.0, 'intensity': 1.0, 'intensity_bounds': '[0, None]', 'width': 13, 'width_bounds': '[0, None]', 'center': 0.0, 'center_bounds': '[-3, 1]'}
+        ]
+
+        result_table = fix_profile_parameters(means_table_mock, profiles_table_mock, h_ratio=14.7)
+
+        # For some reason, reading floats from table introduces small variations in stored values
+        # TODO: Fix floating positions eg. 8.973999977111816 -> 8.974
+        self.assertEqual(
+            result_table.row(0),
+            {'label': '1.0079', 'mass': 1.0077999830245972, 'intensity': 0.8839251399040222, 'intensity_bounds': '[0, None]', 'width': 4.698999881744385, 'width_bounds': '[3, 6]', 'center': 0.0, 'center_bounds': '[-3, 1]'}
+        )
+        self.assertEqual(
+            result_table.row(1),
+            {'label': '16.0', 'mass': 16.0, 'intensity': 0.020470114424824715, 'intensity_bounds': '[0, None]', 'width': 8.973999977111816, 'width_bounds': '[8.974, 8.974]', 'center': 0.0, 'center_bounds': '[-3, 1]'}
+        )
+        self.assertEqual(
+            result_table.row(2),
+            {'label': '12.0', 'mass': 12.0, 'intensity': 0.06013096123933792, 'intensity_bounds': '[0, None]', 'width': 13.932000160217285, 'width_bounds': '[13.932, 13.932]', 'center': 0.0, 'center_bounds': '[-3, 1]'}
+        )
+        self.assertEqual(
+            result_table.row(3),
+            {'label': '27.0', 'mass': 27.0, 'intensity': 0.0354737788438797, 'intensity_bounds': '[0, None]', 'width': 15.397000312805176, 'width_bounds': '[15.397, 15.397]', 'center': 0.0, 'center_bounds': '[-3, 1]'}
+        )
+        
+
+    def test_fix_profile_parameters_without_H(self):
+        # TODO: Use a more physical example containing Deuterium
+        # Same profiles as before, except when H is not present,
+        # the first mass will be made free to vary and the intensities don't change
+
+        means_table_mock = MagicMock()
+        means_table_mock.rowCount.return_value = 3 
+        means_table_mock.row.side_effect = [
+            {'label': '16.0', 'mass': 16.0, 'mean_width': 8.974, 'std_width': 1.401, 'mean_intensity': 0.176, 'std_intensity': 0.08722},
+            {'label': '27.0', 'mass': 27.0, 'mean_width': 15.397, 'std_width': 1.131, 'mean_intensity': 0.305, 'std_intensity': 0.04895},
+            {'label': '12.0', 'mass': 12.0, 'mean_width': 13.932, 'std_width': 0.314, 'mean_intensity': 0.517, 'std_intensity': 0.09531}
+        ]
+        profiles_table_mock = MagicMock()
+        profiles_table_mock.rowCount.return_value = 3
+        profiles_table_mock.row.side_effect = [
+            {'label': '16.0', 'mass': 16.0, 'intensity': 1.0, 'intensity_bounds': '[0, None]', 'width': 12, 'width_bounds': '[0, None]', 'center': 0.0, 'center_bounds': '[-3, 1]'},
+            {'label': '12.0', 'mass': 12.0, 'intensity': 1.0, 'intensity_bounds': '[0, None]', 'width': 8, 'width_bounds': '[0, None]', 'center': 0.0, 'center_bounds': '[-3, 1]'},
+            {'label': '27.0', 'mass': 27.0, 'intensity': 1.0, 'intensity_bounds': '[0, None]', 'width': 13, 'width_bounds': '[0, None]', 'center': 0.0, 'center_bounds': '[-3, 1]'}
+        ]
+
+        result_table = fix_profile_parameters(means_table_mock, profiles_table_mock, h_ratio=14.7)
+
+        # For some reason, reading floats from table introduces small variations in stored values
+        # TODO: Fix floating positions eg. 8.973999977111816 -> 8.974
+        self.assertEqual(
+            result_table.row(0),
+            {'label': '16.0', 'mass': 16.0, 'intensity': 0.17635270953178406, 'intensity_bounds': '[0, None]', 'width': 8.973999977111816, 'width_bounds': '[8.974, 8.974]', 'center': 0.0, 'center_bounds': '[-3, 1]'}
+        )
+        self.assertEqual(
+            result_table.row(1),
+            {'label': '12.0', 'mass': 12.0, 'intensity': 0.5180360674858093, 'intensity_bounds': '[0, None]', 'width': 13.932000160217285, 'width_bounds': '[0, None]', 'center': 0.0, 'center_bounds': '[-3, 1]'}
+        )
+        self.assertEqual(
+            result_table.row(2),
+            {'label': '27.0', 'mass': 27.0, 'intensity': 0.3056112229824066, 'intensity_bounds': '[0, None]', 'width': 15.397000312805176, 'width_bounds': '[15.397, 15.397]', 'center': 0.0, 'center_bounds': '[-3, 1]'}
+        )
+
+if __name__ == "__main__":
+    unittest.main()