diff --git a/src/mvesuvio/analysis_reduction.py b/src/mvesuvio/analysis_reduction.py index 3e79743e..404b22d8 100644 --- a/src/mvesuvio/analysis_reduction.py +++ b/src/mvesuvio/analysis_reduction.py @@ -502,7 +502,8 @@ def _create_means_table(self): self._mean_intensity_ratios, self._std_intensity_ratios, ): - table.addRow([label, mass, mean_width, std_width, mean_intensity, std_intensity]) + # Explicit conversion to float required to match profiles table + table.addRow([label, float(mass), float(mean_width), float(std_width), float(mean_intensity), float(std_intensity)]) self.log().notice(f"{label:6s} {mean_width:10.5f} \u00B1 {std_width:7.5f}" + \ f"{mean_intensity:10.5f} \u00B1 {std_intensity:7.5f}\n") diff --git a/src/mvesuvio/main/run_routine.py b/src/mvesuvio/main/run_routine.py index 99be5f6d..8212bc0a 100644 --- a/src/mvesuvio/main/run_routine.py +++ b/src/mvesuvio/main/run_routine.py @@ -2,22 +2,21 @@ from mvesuvio.util import handle_config from mvesuvio.util.analysis_helpers import fix_profile_parameters, \ loadRawAndEmptyWsFromUserPath, cropAndMaskWorkspace, \ - NeutronComptonProfile, calculate_h_ratio, name_for_starting_ws + NeutronComptonProfile, calculate_h_ratio, name_for_starting_ws, \ + scattering_type, ws_history_matches_inputs, save_ws_from_load_vesuvio, \ + is_hydrogen_present, create_profiles_table, create_table_for_hydrogen_to_mass_ratios from mvesuvio.analysis_reduction import VesuvioAnalysisRoutine from mantid.api import mtd from mantid.api import AnalysisDataService -from mantid.simpleapi import CreateEmptyTableWorkspace, mtd, RenameWorkspace, \ - Load, LoadVesuvio, SaveNexus, DeleteWorkspace +from mantid.simpleapi import mtd, RenameWorkspace from mantid.api import AlgorithmFactory, AlgorithmManager -from mantid.kernel import logger import numpy as np from pathlib import Path import importlib import sys import dill # To convert constraints to string -import ntpath class Runner: @@ -38,9 +37,10 @@ def setup(self): self.yFitIC = ai.YSpaceFitInitialConditions self.userCtr = ai.UserScriptControls - checkInputs(self.userCtr) + for flag in [self.userCtr.procedure, self.userCtr.fitInYSpace]: + assert flag in ["BACKWARD", "FORWARD", "JOINT"], "Procedure not recognized" - # Names of workspaces to check if they exist to do fitting + # Names of workspaces to check if they exist to skip analysis self.ws_to_fit_y_space = [] self.classes_to_fit_y_space = [] for mode, i_cls in zip(["BACKWARD", "FORWARD"], [self.bckwdIC, self.fwdIC]): @@ -113,12 +113,12 @@ def runAnalysisRoutine(self): if (routine_type == "BACKWARD") | (routine_type== "JOINT"): - if isHPresent(self.fwdIC.masses) & (self.bckwdIC.HToMassIdxRatio==0): + if is_hydrogen_present(self.fwdIC.masses) & (self.bckwdIC.HToMassIdxRatio==0): self.run_estimate_h_ratio() return # TODO: make this automatic - assert isHPresent(self.fwdIC.masses) != ( + assert is_hydrogen_present(self.fwdIC.masses) != ( self.bckwdIC.HToMassIdxRatio==0 ), "No Hydrogen detected, HToMassIdxRatio has to be set to 0" @@ -143,7 +143,31 @@ def run_joint_analysis(self): AnalysisDataService.clear() back_alg = self._create_analysis_algorithm(self.wsBackIC, self.bckwdIC) front_alg = self._create_analysis_algorithm(self.wsFrontIC, self.fwdIC) - self.analysis_result = run_joint_algs(back_alg, front_alg) + self.run_joint_algs(back_alg, front_alg) + return + + + @classmethod + def run_joint_algs(cls, back_alg, front_alg): + + back_alg.execute() + + incoming_means_table = mtd[back_alg.getPropertyValue("OutputMeansTable")] + h_ratio = back_alg.getProperty("HRatioToLowestMass").value + + 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 = mtd[front_alg.getPropertyValue("InputProfiles")] + + 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.setPropertyValue("InputProfiles", receiving_profiles_table.name()) + + front_alg.execute() return @@ -165,7 +189,7 @@ def run_estimate_h_ratio(self): if not ((userInput == "y") or (userInput == "Y")): raise KeyboardInterrupt("Procedure interrupted.") - table_h_ratios = createTableWSHRatios() + table_h_ratios = create_table_for_hydrogen_to_mass_ratios() back_alg = self._create_analysis_algorithm(self.wsBackIC, self.bckwdIC) front_alg = self._create_analysis_algorithm(self.wsFrontIC, self.fwdIC) @@ -181,7 +205,7 @@ def run_estimate_h_ratio(self): while not np.isclose(current_ratio, previous_ratio, rtol=0.01): back_alg.setProperty("HRatioToLowestMass", current_ratio) - run_joint_algs(back_alg, front_alg) + self.run_joint_algs(back_alg, front_alg) previous_ratio = current_ratio @@ -217,16 +241,7 @@ def _create_analysis_algorithm(self, load_ai, ai): maskTOFRange=ai.maskTOFRange ) - profiles = [] - for mass, intensity, width, center, intensity_bound, width_bound, center_bound in zip( - ai.masses, ai.initPars[::3], ai.initPars[1::3], ai.initPars[2::3], - ai.bounds[::3], ai.bounds[1::3], ai.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 = self._make_neutron_profiles(ai) profiles_table = create_profiles_table(cropedWs.name()+"_initial_parameters", profiles) ipFilesPath = Path(handle_config.read_config_var("caching.ipfolder")) @@ -243,7 +258,7 @@ def _create_analysis_algorithm(self, load_ai, ai): "SampleHorizontalWidth": ai.horizontal_width, "SampleThickness": ai.thickness, "GammaCorrection": ai.GammaCorrectionFlag, - "ModeRunning": getRunningMode(load_ai), + "ModeRunning": scattering_type(load_ai), "TransmissionGuess": ai.transmission_guess, "MultipleScatteringOrder": int(ai.multiple_scattering_order), "NumberOfEvents": int(ai.number_of_events), @@ -264,24 +279,37 @@ def _create_analysis_algorithm(self, load_ai, ai): return alg + def _make_neutron_profiles(cls, ai): + profiles = [] + for mass, intensity, width, center, intensity_bound, width_bound, center_bound in zip( + ai.masses, ai.initPars[::3], ai.initPars[1::3], ai.initPars[2::3], + ai.bounds[::3], ai.bounds[1::3], ai.bounds[2::3] + ): + profiles.append(NeutronComptonProfile( + label=str(float(mass)), mass=float(mass), intensity=float(intensity), width=float(width), center=float(center), + intensity_bounds=list(intensity_bound), width_bounds=list(width_bound), center_bounds=list(center_bound) + )) + return profiles + + def _save_ws_if_not_on_path(self, load_ai): - runningMode = getRunningMode(load_ai).lower() + scatteringType = scattering_type(load_ai).lower() scriptName = handle_config.get_script_name() - rawWSName = scriptName + "_" + "raw" + "_" + runningMode + ".nxs" - emptyWSName = scriptName + "_" + "empty" + "_" + runningMode + ".nxs" + rawWSName = scriptName + "_" + "raw" + "_" + scatteringType + ".nxs" + emptyWSName = scriptName + "_" + "empty" + "_" + scatteringType + ".nxs" rawPath = self.input_ws_path / rawWSName emptyPath = self.input_ws_path / emptyWSName ipFilesPath = Path(handle_config.read_config_var("caching.ipfolder")) - if not wsHistoryMatchesInputs(load_ai.runs, load_ai.mode, load_ai.ipfile, rawPath): - saveWSFromLoadVesuvio(load_ai.runs, load_ai.mode, str(ipFilesPath/load_ai.ipfile), rawPath) + if not ws_history_matches_inputs(load_ai.runs, load_ai.mode, load_ai.ipfile, rawPath): + save_ws_from_load_vesuvio(load_ai.runs, load_ai.mode, str(ipFilesPath/load_ai.ipfile), rawPath) - if not wsHistoryMatchesInputs(load_ai.empty_runs, load_ai.mode, load_ai.ipfile, emptyPath): - saveWSFromLoadVesuvio(load_ai.empty_runs, load_ai.mode, str(ipFilesPath/load_ai.ipfile), emptyPath) + if not ws_history_matches_inputs(load_ai.empty_runs, load_ai.mode, load_ai.ipfile, emptyPath): + save_ws_from_load_vesuvio(load_ai.empty_runs, load_ai.mode, str(ipFilesPath/load_ai.ipfile), emptyPath) return rawPath, emptyPath @@ -312,155 +340,5 @@ def _set_output_paths(self, ai): return -def getRunningMode(wsIC): - if wsIC.__name__ == "LoadVesuvioBackParameters": - runningMode = "BACKWARD" - elif wsIC.__name__ == "LoadVesuvioFrontParameters": - runningMode = "FORWARD" - else: - raise ValueError( - f"Input class for loading workspace not valid: {wsIC.__name__}" - ) - return runningMode - - -def wsHistoryMatchesInputs(runs, mode, ipfile, localPath): - if not (localPath.is_file()): - logger.notice("Cached workspace not found") - return False - local_ws = Load(Filename=str(localPath)) - ws_history = local_ws.getHistory() - metadata = ws_history.getAlgorithmHistory(0) - - saved_runs = metadata.getPropertyValue("Filename") - if saved_runs != runs: - logger.notice( - f"Filename in saved workspace did not match: {saved_runs} and {runs}" - ) - return False - - saved_mode = metadata.getPropertyValue("Mode") - if saved_mode != mode: - logger.notice(f"Mode in saved workspace did not match: {saved_mode} and {mode}") - return False - - saved_ipfile_name = ntpath.basename(metadata.getPropertyValue("InstrumentParFile")) - if saved_ipfile_name != ipfile: - logger.notice( - f"IP files in saved workspace did not match: {saved_ipfile_name} and {ipfilename}" - ) - return False - - print("\nLocally saved workspace metadata matched with analysis inputs.\n") - DeleteWorkspace(local_ws) - return True - - -def saveWSFromLoadVesuvio(runs, mode, ipfile, localPath): - if "backward" in localPath.name: - spectra = "3-134" - elif "forward" in localPath.name: - spectra = "135-198" - else: - raise ValueError(f"Invalid name to save workspace: {localPath.name}") - - vesuvio_ws = LoadVesuvio( - Filename=runs, - SpectrumList=spectra, - Mode=mode, - InstrumentParFile=str(ipfile), - OutputWorkspace=localPath.name, - LoadLogFiles=False, - ) - - SaveNexus(vesuvio_ws, str(localPath.absolute())) - print(f"Workspace saved locally at: {localPath.absolute()}") - return - - -def checkInputs(crtIC): - try: - if ~crtIC.runRoutine: - return - except AttributeError: - if ~crtIC.runBootstrap: - return - - for flag in [crtIC.procedure, crtIC.fitInYSpace]: - assert ( - (flag == "BACKWARD") - | (flag == "FORWARD") - | (flag == "JOINT") - | (flag is None) - ), "Option not recognized." - - 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 - - -def create_profiles_table(name, profiles: list[NeutronComptonProfile]): - table = CreateEmptyTableWorkspace(OutputWorkspace=name) - table.addColumn(type="str", name="label") - table.addColumn(type="float", name="mass") - table.addColumn(type="float", name="intensity") - table.addColumn(type="str", name="intensity_bounds") - table.addColumn(type="float", name="width") - table.addColumn(type="str", name="width_bounds") - table.addColumn(type="float", name="center") - table.addColumn(type="str", name="center_bounds") - for p in profiles: - table.addRow([str(getattr(p, attr)) - if "bounds" in attr - else getattr(p, attr) - for attr in table.getColumnNames()]) - return table - - -def run_joint_algs(back_alg, front_alg): - - back_alg.execute() - - incoming_means_table = mtd[back_alg.getPropertyValue("OutputMeansTable")] - h_ratio = back_alg.getProperty("HRatioToLowestMass").value - - 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 = mtd[front_alg.getPropertyValue("InputProfiles")] - - 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.setPropertyValue("InputProfiles", receiving_profiles_table.name()) - - front_alg.execute() - return - - -def createTableWSHRatios(): - table = CreateEmptyTableWorkspace( - OutputWorkspace="hydrogen_intensity_ratios_estimates" - ) - table.addColumn(type="float", name="Hydrogen intensity ratio to lowest mass at each iteration") - return table - - if __name__ == "__main__": Runner().run() diff --git a/src/mvesuvio/util/analysis_helpers.py b/src/mvesuvio/util/analysis_helpers.py index b422d3ad..b808a7ca 100644 --- a/src/mvesuvio/util/analysis_helpers.py +++ b/src/mvesuvio/util/analysis_helpers.py @@ -1,7 +1,8 @@ from mantid.simpleapi import Load, Rebin, Scale, SumSpectra, Minus, CropWorkspace, \ - CloneWorkspace, MaskDetectors, CreateWorkspace, CreateEmptyTableWorkspace, \ - mtd, RenameWorkspace + MaskDetectors, CreateWorkspace, CreateEmptyTableWorkspace, \ + DeleteWorkspace, SaveNexus, LoadVesuvio, mtd +from mantid.kernel import logger import numpy as np import numbers @@ -9,6 +10,7 @@ from mvesuvio.util import handle_config from dataclasses import dataclass +import ntpath @dataclass(frozen=False) @@ -29,17 +31,124 @@ class NeutronComptonProfile: mean_center: float = None +def create_profiles_table(name, profiles: list[NeutronComptonProfile]): + table = CreateEmptyTableWorkspace(OutputWorkspace=name) + table.addColumn(type="str", name="label") + table.addColumn(type="float", name="mass") + table.addColumn(type="float", name="intensity") + table.addColumn(type="str", name="intensity_bounds") + table.addColumn(type="float", name="width") + table.addColumn(type="str", name="width_bounds") + table.addColumn(type="float", name="center") + table.addColumn(type="str", name="center_bounds") + for p in profiles: + table.addRow([str(getattr(p, attr)) + if "bounds" in attr + else getattr(p, attr) + for attr in table.getColumnNames()]) + return table + + +def create_table_for_hydrogen_to_mass_ratios(): + table = CreateEmptyTableWorkspace( + OutputWorkspace="hydrogen_intensity_ratios_estimates" + ) + table.addColumn(type="float", name="Hydrogen intensity ratio to lowest mass at each iteration") + return table + + +def is_hydrogen_present(masses) -> bool: + Hmask = np.abs(masses - 1) / 1 < 0.1 # H mass whithin 10% of 1 au + + if ~np.any(Hmask): # H not present + return False + + 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 + + +def ws_history_matches_inputs(runs, mode, ipfile, ws_path): + + if not (ws_path.is_file()): + logger.notice("Cached workspace not found") + return False + + ws = Load(Filename=str(ws_path)) + ws_history = ws.getHistory() + metadata = ws_history.getAlgorithmHistory(0) + + saved_runs = metadata.getPropertyValue("Filename") + if saved_runs != runs: + logger.notice( + f"Filename in saved workspace did not match: {saved_runs} and {runs}" + ) + return False + + saved_mode = metadata.getPropertyValue("Mode") + if saved_mode != mode: + logger.notice(f"Mode in saved workspace did not match: {saved_mode} and {mode}") + return False + + saved_ipfile_name = ntpath.basename(metadata.getPropertyValue("InstrumentParFile")) + if saved_ipfile_name != ipfile: + logger.notice( + f"IP files in saved workspace did not match: {saved_ipfile_name} and {ipfile}" + ) + return False + + print("\nLocally saved workspace metadata matched with analysis inputs.\n") + DeleteWorkspace(ws) + return True + + +def save_ws_from_load_vesuvio(runs, mode, ipfile, ws_path): + + if "backward" in ws_path.name: + spectra = "3-134" + elif "forward" in ws_path.name: + spectra = "135-198" + else: + raise ValueError(f"Invalid name to save workspace: {ws_path.name}") + + vesuvio_ws = LoadVesuvio( + Filename=runs, + SpectrumList=spectra, + Mode=mode, + InstrumentParFile=str(ipfile), + OutputWorkspace=ws_path.name, + LoadLogFiles=False, + ) + + SaveNexus(vesuvio_ws, str(ws_path.absolute())) + print(f"Workspace saved locally at: {ws_path.absolute()}") + return + + def name_for_starting_ws(load_ai): + name_suffix = scattering_type(load_ai, shorthand=True) + name = handle_config.get_script_name() + "_" + name_suffix + return name + + +def scattering_type(load_ai, shorthand=False): if load_ai.__name__ in ["LoadVesuvioBackParameters", "BackwardInitialConditions"]: - name_suffix = "bckwd" + scatteringType = "BACKWARD" + if shorthand: + scatteringType = "bckwd" elif load_ai.__name__ in ["LoadVesuvioFrontParameters", "ForwardInitialConditions"]: - name_suffix = "fwd" + scatteringType = "FORWARD" + if shorthand: + scatteringType = "fwd" else: raise ValueError( f"Input class for workspace not valid: {load_ai.__name__}" ) - name = handle_config.get_script_name() + "_" + name_suffix - return name + return scatteringType def loadRawAndEmptyWsFromUserPath(userWsRawPath, userWsEmptyPath,