diff --git a/src/mvesuvio/analysis_reduction.py b/src/mvesuvio/analysis_reduction.py index a18966f7..97dd19e6 100644 --- a/src/mvesuvio/analysis_reduction.py +++ b/src/mvesuvio/analysis_reduction.py @@ -516,6 +516,7 @@ 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") @@ -523,14 +524,15 @@ def _create_means_table(self): 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") diff --git a/src/mvesuvio/analysis_routines.py b/src/mvesuvio/analysis_routines.py index 27e53195..8bc603b3 100644 --- a/src/mvesuvio/analysis_routines.py +++ b/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 diff --git a/src/mvesuvio/util/analysis_helpers.py b/src/mvesuvio/util/analysis_helpers.py index 5125bc36..0e29c786 100644 --- a/src/mvesuvio/util/analysis_helpers.py +++ b/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)] diff --git a/tests/unit/analysis/test_analysis_functions.py b/tests/unit/analysis/test_analysis_functions.py deleted file mode 100644 index 15ef2324..00000000 --- a/tests/unit/analysis/test_analysis_functions.py +++ /dev/null @@ -1,26 +0,0 @@ -import unittest -import numpy as np -import numpy.testing as nptest -from mock import MagicMock -from mvesuvio.util.analysis_helpers import extractWS -from mantid.simpleapi import CreateWorkspace, DeleteWorkspace - - -class TestAnalysisFunctions(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) - - -if __name__ == "__main__": - unittest.main() diff --git a/tests/unit/analysis/test_analysis_helpers.py b/tests/unit/analysis/test_analysis_helpers.py new file mode 100644 index 00000000..d67ce4c8 --- /dev/null +++ b/tests/unit/analysis/test_analysis_helpers.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()