From 15556a0dd716d85b573284ab5f3244ddb120855a Mon Sep 17 00:00:00 2001 From: GuiMacielPereira Date: Fri, 13 Sep 2024 14:34:28 +0100 Subject: [PATCH] Refactor routine for estimating h ratio When h ratio is not provided and hydrogen is present in the sample, the routine for estimating the h ratio to the lowest mass is triggered. Had to update this routine with the changes from converting vesuvio analysis into a mantid algorithm. --- src/mvesuvio/analysis_reduction.py | 15 +----- src/mvesuvio/analysis_routines.py | 54 +++++++++----------- src/mvesuvio/run_routine.py | 18 ++++++- src/mvesuvio/util/analysis_helpers.py | 18 ++++++- tests/unit/analysis/test_analysis_helpers.py | 11 +++- 5 files changed, 68 insertions(+), 48 deletions(-) diff --git a/src/mvesuvio/analysis_reduction.py b/src/mvesuvio/analysis_reduction.py index 662b54d3..7e10109a 100644 --- a/src/mvesuvio/analysis_reduction.py +++ b/src/mvesuvio/analysis_reduction.py @@ -9,7 +9,7 @@ CloneWorkspace, DeleteWorkspace, VesuvioCalculateGammaBackground, \ VesuvioCalculateMS, Scale, RenameWorkspace, Minus, CreateSampleShape, \ VesuvioThickness, Integration, Divide, Multiply, DeleteWorkspaces, \ - CreateWorkspace + CreateWorkspace, CreateSampleWorkspace from mvesuvio.util.analysis_helpers import loadConstants, numericalThirdDerivative @@ -142,7 +142,7 @@ def PyExec(self): self.run() def _setup(self): - self._name = self.getProperty("InputWorkspace").value.name() + self._name = self.getPropertyValue("InputWorkspace") self._ip_file = self.getProperty("InstrumentParametersFile").value self._number_of_iterations = self.getProperty("NumberOfIterations").value self._mask_spectra = self.getProperty("InvalidDetectors").value @@ -191,17 +191,6 @@ def _setup(self): self._fit_profiles_workspaces = {} - def calculate_h_ratio(self): - - if not np.isclose(min(self._masses), 1, atol=0.1): # Hydrogen not present - return None - - # Hydrogen present - intensities = self._mean_intensity_ratios - sorted_intensities = intensities[np.argsort(self._masses)] - - return sorted_intensities[0] / sorted_intensities[1] - def _update_workspace_data(self): diff --git a/src/mvesuvio/analysis_routines.py b/src/mvesuvio/analysis_routines.py index 37adad87..6110f2e1 100644 --- a/src/mvesuvio/analysis_routines.py +++ b/src/mvesuvio/analysis_routines.py @@ -1,10 +1,12 @@ # from .analysis_reduction import iterativeFitForDataReduction from mantid.api import AnalysisDataService -from mantid.simpleapi import CreateEmptyTableWorkspace, mtd +from mantid.simpleapi import CreateEmptyTableWorkspace, mtd, RenameWorkspace from mantid.api import AlgorithmFactory, AlgorithmManager import numpy as np -from mvesuvio.util.analysis_helpers import fix_profile_parameters, loadRawAndEmptyWsFromUserPath, cropAndMaskWorkspace, NeutronComptonProfile +from mvesuvio.util.analysis_helpers import fix_profile_parameters, \ + loadRawAndEmptyWsFromUserPath, cropAndMaskWorkspace, \ + NeutronComptonProfile, calculate_h_ratio from mvesuvio.analysis_reduction import AnalysisRoutine from tests.testhelpers.calibration.algorithms import create_algorithm @@ -39,7 +41,7 @@ def _create_analysis_object_from_current_interface(IC, running_tests=False): profiles_table = create_profiles_table(cropedWs.name()+"_initial_parameters", profiles) kwargs = { - "InputWorkspace": cropedWs, + "InputWorkspace": cropedWs.name(), "InputProfiles": profiles_table.name(), "InstrumentParametersFile": str(IC.InstrParsPath), "HRatioToLowestMass": IC.HToMassIdxRatio, @@ -133,12 +135,14 @@ def run_joint_algs(back_alg, front_alg): assert incoming_means_table is not None, "Means table from backward routine not correctly accessed." assert h_ratio is not None, "H ratio from backward routine not correctly accesssed." - receiving_profiles_table = front_alg.getProperty("InputProfiles").value + receiving_profiles_table = mtd[front_alg.getPropertyValue("InputProfiles")] - fix_profile_parameters(incoming_means_table, receiving_profiles_table, h_ratio) + fixed_profiles_table = fix_profile_parameters(incoming_means_table, receiving_profiles_table, h_ratio) + # Update original profiles table + RenameWorkspace(fixed_profiles_table, receiving_profiles_table.name()) # Even if the name is the same, need to trigger update - front_alg.setProperty("InputProfiles", receiving_profiles_table.name()) + front_alg.setPropertyValue("InputProfiles", receiving_profiles_table.name()) front_alg.execute() return @@ -164,23 +168,26 @@ def runPreProcToEstHRatio(bckwdIC, fwdIC): table_h_ratios = createTableWSHRatios() - 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) + + front_alg.execute() + + means_table = mtd[front_alg.getPropertyValue("OutputMeansTable")] + current_ratio = calculate_h_ratio(means_table) - frontRoutine.execute() - current_ratio = frontRoutine.calculate_h_ratio() table_h_ratios.addRow([current_ratio]) previous_ratio = np.nan while not np.isclose(current_ratio, previous_ratio, rtol=0.01): - backRoutine._h_ratio = current_ratio - backRoutine.execute() - frontRoutine.set_initial_profiles_from(backRoutine) - frontRoutine.execute() + back_alg.setProperty("HRatioToLowestMass", current_ratio) + run_joint_algs(back_alg, front_alg) previous_ratio = current_ratio - current_ratio = frontRoutine.calculate_h_ratio() + + means_table = mtd[front_alg.getPropertyValue("OutputMeansTable")] + current_ratio = calculate_h_ratio(means_table) table_h_ratios.addRow([current_ratio]) @@ -192,22 +199,9 @@ def runPreProcToEstHRatio(bckwdIC, fwdIC): def createTableWSHRatios(): table = CreateEmptyTableWorkspace( - OutputWorkspace="H_Ratios_Estimates" + OutputWorkspace="hydrogen_intensity_ratios_estimates" ) - table.addColumn(type="float", name="H Ratio to lowest mass at each iteration") + table.addColumn(type="float", name="Hydrogen intensity ratio to lowest mass at each iteration") return table -def isHPresent(masses) -> bool: - Hmask = np.abs(masses - 1) / 1 < 0.1 # H mass whithin 10% of 1 au - - if np.any(Hmask): # H present - print("\nH mass detected.\n") - assert ( - len(Hmask) > 1 - ), "When H is only mass present, run independent forward procedure, not joint." - assert Hmask[0], "H mass needs to be the first mass in masses and initPars." - assert sum(Hmask) == 1, "More than one mass very close to H were detected." - return True - else: - return False diff --git a/src/mvesuvio/run_routine.py b/src/mvesuvio/run_routine.py index 1b6e7ca5..29545374 100644 --- a/src/mvesuvio/run_routine.py +++ b/src/mvesuvio/run_routine.py @@ -8,11 +8,10 @@ runIndependentIterativeProcedure, runJointBackAndForwardProcedure, runPreProcToEstHRatio, - createTableWSHRatios, - isHPresent, ) from mantid.api import mtd +import numpy as np def runRoutine( userCtr, @@ -115,3 +114,18 @@ def checkInputs(crtIC): if (crtIC.procedure != "JOINT") & (crtIC.fitInYSpace is not None): assert crtIC.procedure == crtIC.fitInYSpace + + +def isHPresent(masses) -> bool: + Hmask = np.abs(masses - 1) / 1 < 0.1 # H mass whithin 10% of 1 au + + if np.any(Hmask): # H present + print("\nH mass detected.\n") + assert ( + len(Hmask) > 1 + ), "When H is only mass present, run independent forward procedure, not joint." + assert Hmask[0], "H mass needs to be the first mass in masses and initPars." + assert sum(Hmask) == 1, "More than one mass very close to H were detected." + return True + else: + return False diff --git a/src/mvesuvio/util/analysis_helpers.py b/src/mvesuvio/util/analysis_helpers.py index 1efebfec..12476eaa 100644 --- a/src/mvesuvio/util/analysis_helpers.py +++ b/src/mvesuvio/util/analysis_helpers.py @@ -206,8 +206,7 @@ def fix_profile_parameters(incoming_means_table, receiving_profiles_table, h_rat p['width_bounds'] = str([p['width'] , p['width']]) result_profiles_table = _convert_dict_to_table(profiles_dict) - RenameWorkspace(result_profiles_table, receiving_profiles_table.name()) - return mtd[result_profiles_table.name()] + return result_profiles_table def _convert_table_to_dict(table): @@ -234,3 +233,18 @@ 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)] + + +def calculate_h_ratio(means_table): + + masses = means_table.column("mass") + intensities = np.array(means_table.column("mean_intensity")) + + if not np.isclose(min(masses), 1, atol=0.1): # Hydrogen not present + return None + + # Hydrogen present + sorted_intensities = intensities[np.argsort(masses)] + + return sorted_intensities[0] / sorted_intensities[1] + diff --git a/tests/unit/analysis/test_analysis_helpers.py b/tests/unit/analysis/test_analysis_helpers.py index d67ce4c8..29a1146d 100644 --- a/tests/unit/analysis/test_analysis_helpers.py +++ b/tests/unit/analysis/test_analysis_helpers.py @@ -2,7 +2,8 @@ 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 mvesuvio.util.analysis_helpers import extractWS, _convert_dict_to_table, \ + fix_profile_parameters, calculate_h_ratio from mantid.simpleapi import CreateWorkspace, DeleteWorkspace @@ -105,5 +106,13 @@ def test_fix_profile_parameters_without_H(self): {'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]'} ) + + def test_calculate_h_ratio(self): + means_table_mock = MagicMock() + means_table_mock.column.side_effect = lambda x: [16, 1, 12] if x is "mass" else [0.1, 0.85, 0.05] + h_ratio = calculate_h_ratio(means_table_mock) + self.assertEqual(h_ratio, 0.85 / 0.05) + + if __name__ == "__main__": unittest.main()