From 0d2c608b6ebf71bcd363812fc50420ecd11568cc Mon Sep 17 00:00:00 2001 From: GuiMacielPereira Date: Tue, 10 Dec 2024 11:45:42 +0000 Subject: [PATCH 1/7] Rename variables on input script --- src/mvesuvio/analysis_fitting.py | 13 ++-- src/mvesuvio/config/analysis_inputs.py | 62 +++++++++---------- src/mvesuvio/main/run_routine.py | 57 ++++++++--------- src/mvesuvio/util/analysis_helpers.py | 4 +- tests/data/analysis/inputs/analysis_test.py | 60 +++++++++--------- .../data/analysis/inputs/yspace_gauss_test.py | 4 -- tests/data/analysis/inputs/yspace_gc_test.py | 4 -- 7 files changed, 97 insertions(+), 107 deletions(-) diff --git a/src/mvesuvio/analysis_fitting.py b/src/mvesuvio/analysis_fitting.py index 2b6c13b0..0b811d44 100644 --- a/src/mvesuvio/analysis_fitting.py +++ b/src/mvesuvio/analysis_fitting.py @@ -54,7 +54,7 @@ def find_ws_name_fse_first_mass(ic): # Some dirty implementation to be deleted in the future ws_names_fse = [] ws_masses = [] - prefix = ic.name+'_'+str(ic.noOfMSIterations) + prefix = ic.name+'_'+str(ic.number_of_iterations_for_corrections) for ws_name in mtd.getObjectNames(): if ws_name.startswith(prefix) and ws_name.endswith('fse'): name_ending = ws_name.replace(prefix, "") @@ -103,8 +103,8 @@ def subtract_profiles_except_lightest(ic, ws): if len(ic.masses) == 1: return - ws_name_lightest_mass = ic.name + '_' + str(ic.noOfMSIterations) + '_' + str(min(ic.masses)) + '_ncp' - ws_name_profiles = ic.name + '_' + str(ic.noOfMSIterations) + '_total_ncp' + ws_name_lightest_mass = ic.name + '_' + str(ic.number_of_iterations_for_corrections) + '_' + str(min(ic.masses)) + '_ncp' + ws_name_profiles = ic.name + '_' + str(ic.number_of_iterations_for_corrections) + '_total_ncp' wsNcpExceptFirst = Minus(mtd[ws_name_profiles], mtd[ws_name_lightest_mass], OutputWorkspace=ws_name_profiles + '_except_lightest') @@ -126,7 +126,7 @@ def switchFirstTwoAxis(A): def ySpaceReduction(wsTOF, mass0, yFitIC, ic): """Seperate procedures depending on masking specified.""" - ws_name_lightest_mass = ic.name + '_' + str(ic.noOfMSIterations) + '_' + str(min(ic.masses)) + '_ncp' + ws_name_lightest_mass = ic.name + '_' + str(ic.number_of_iterations_for_corrections) + '_' + str(min(ic.masses)) + '_ncp' ncp = mtd[ws_name_lightest_mass].extractY() rebinPars = yFitIC.rebinParametersForYSpaceFit @@ -1503,9 +1503,10 @@ def extractData(ws, wsRes, ic): def loadInstrParsFileIntoArray(ic): ipFilesPath = Path(handle_config.read_config_var("caching.ipfolder")) - data = np.loadtxt(str(ipFilesPath / ic.ipfile), dtype=str)[1:].astype(float) + data = np.loadtxt(str(ipFilesPath / ic.instrument_parameters_file), dtype=str)[1:].astype(float) spectra = data[:, 0] - select_rows = np.where((spectra >= ic.firstSpec) & (spectra <= ic.lastSpec)) + firstSpec, lastSpec = [int(d) for d in ic.detectors.split('-')] + select_rows = np.where((spectra >= firstSpec) & (spectra <= lastSpec)) instrPars = data[select_rows] return instrPars diff --git a/src/mvesuvio/config/analysis_inputs.py b/src/mvesuvio/config/analysis_inputs.py index af287c78..06bacb16 100644 --- a/src/mvesuvio/config/analysis_inputs.py +++ b/src/mvesuvio/config/analysis_inputs.py @@ -11,39 +11,38 @@ class SampleParameters: @dataclass class BackwardAnalysisInputs(SampleParameters): - run_this_scattering_type = False + run_this_scattering_type = True fit_in_y_space = False runs = "43066-43076" empty_runs = "41876-41923" mode = "DoubleDifference" - ipfile = "ip2019.par" - firstSpec = 3 # 3 - lastSpec = 134 # 134 - maskedSpecAllNo = [18, 34, 42, 43, 59, 60, 62, 118, 119, 133] - tofBinning = "275.,1.,420" # Binning of ToF spectra - maskTOFRange = None - subEmptyFromRaw = True # Flag to control wether empty ws gets subtracted from raw - scaleEmpty = 1 # None or scaling factor - scaleRaw = 1 + instrument_parameters_file = "ip2019.par" + detectors = "3-134" + mask_detectors = [18, 34, 42, 43, 59, 60, 62, 118, 119, 133] + time_of_flight_binning = "275.,1.,420" + mask_time_of_flight_range = None + subtract_empty_workspace_from_raw = True + scale_empty_workspace = 1 # None or scaling factor + scale_raw_workspace = 1 masses = [12, 16, 27] # Intensities, NCP widths, NCP centers - initPars = [1, 12, 0.0, 1, 12, 0.0, 1, 12.5, 0.0] - bounds = [ + initial_fitting_parameters = [1, 12, 0.0, 1, 12, 0.0, 1, 12.5, 0.0] + fitting_bounds = [ [0, None], [8, 16], [-3, 1], [0, None], [8, 16], [-3, 1], [0, None], [11, 14], [-3, 1], ] constraints = () - noOfMSIterations = 0 # 4 - MSCorrectionFlag = True - HToMassIdxRatio = 19.0620008206 # Set to zero when H is not present + number_of_iterations_for_corrections = 0 # 4 + do_multiple_scattering_correction = True + intensity_ratio_of_hydrogen_to_lowest_mass = 19.0620008206 # Set to zero to disable transmission_guess = 0.8537 # Experimental value from VesuvioTransmission multiple_scattering_order = 2 - number_of_events = 1.0e5 - GammaCorrectionFlag = False + multiple_scattering_number_of_events = 1.0e5 + do_gamma_correction = False @dataclass @@ -54,20 +53,19 @@ class ForwardAnalysisInputs(SampleParameters): runs = "43066-43076" empty_runs = "43868-43911" mode = "SingleDifference" - ipfile = "ip2018_3.par" - firstSpec = 144 # 144 - lastSpec = 182 # 182 - maskedSpecAllNo = [173, 174, 179] - tofBinning = "110,1,430" # Binning of ToF spectra - maskTOFRange = None - subEmptyFromRaw = False # Flag to control wether empty ws gets subtracted from raw - scaleEmpty = 1 # None or scaling factor - scaleRaw = 1 + instrument_parameters_file = "ip2018_3.par" + detectors = '144-182' + mask_detectors = [173, 174, 179] + time_of_flight_binning = "110,1,430" + mask_time_of_flight_range = None + subtract_empty_workspace_from_raw = False + scale_empty_workspace = 1 # None or scaling factor + scale_raw_workspace = 1 masses = 1.0079, 12, 16, 27 # Intensities, NCP widths, NCP centers - initPars = [1, 4.7, 0, 1, 12.71, 0.0, 1, 8.76, 0.0, 1, 13.897, 0.0] - bounds = [ + initial_fitting_parameters = [1, 4.7, 0, 1, 12.71, 0.0, 1, 8.76, 0.0, 1, 13.897, 0.0] + fitting_bounds = [ [0, None], [3, 6], [-3, 1], [0, None], [12.71, 12.71], [-3, 1], [0, None], [8.76, 8.76], [-3, 1], @@ -75,12 +73,12 @@ class ForwardAnalysisInputs(SampleParameters): ] constraints = () - noOfMSIterations = 0 # 4 - MSCorrectionFlag = True + number_of_iterations_for_corrections = 0 # 4 + do_multiple_scattering_correction = True transmission_guess = 0.8537 # Experimental value from VesuvioTransmission multiple_scattering_order = 2 - number_of_events = 1.0e5 - GammaCorrectionFlag = True + multiple_scattering_number_of_events = 1.0e5 + do_gamma_correction = True @dataclass diff --git a/src/mvesuvio/main/run_routine.py b/src/mvesuvio/main/run_routine.py index 43e335c3..eeb58584 100644 --- a/src/mvesuvio/main/run_routine.py +++ b/src/mvesuvio/main/run_routine.py @@ -39,7 +39,7 @@ def setup(self): self.classes_to_fit_y_space = [] for ai_cls in [self.bckwd_ai, self.fwd_ai]: if ai_cls.fit_in_y_space: - self.ws_to_fit_y_space.append(name_for_starting_ws(ai_cls) + '_' + str(ai_cls.noOfMSIterations)) + self.ws_to_fit_y_space.append(name_for_starting_ws(ai_cls) + '_' + str(ai_cls.number_of_iterations_for_corrections)) self.classes_to_fit_y_space.append(ai_cls) self.analysis_result = None @@ -103,14 +103,14 @@ def runAnalysisRoutine(self): if self.bckwd_ai.run_this_scattering_type: - if is_hydrogen_present(self.fwd_ai.masses) & (self.bckwd_ai.HToMassIdxRatio==0): + if is_hydrogen_present(self.fwd_ai.masses) & (self.bckwd_ai.intensity_ratio_of_hydrogen_to_lowest_mass==0): self.run_estimate_h_ratio() return # TODO: make this automatic assert is_hydrogen_present(self.fwd_ai.masses) != ( - self.bckwd_ai.HToMassIdxRatio==0 - ), "No Hydrogen detected, HToMassIdxRatio has to be set to 0" + self.bckwd_ai.intensity_ratio_of_hydrogen_to_lowest_mass==0 + ), "No Hydrogen detected, intensity_ratio_of_hydrogen_to_lowest_mass has to be set to 0" if self.bckwd_ai.run_this_scattering_type and self.fwd_ai.run_this_scattering_type: self.run_joint_analysis() @@ -220,37 +220,38 @@ def _create_analysis_algorithm(self, ai): ws = loadRawAndEmptyWsFromUserPath( userWsRawPath=raw_path, userWsEmptyPath=empty_path, - tofBinning=ai.tofBinning, + tofBinning=ai.time_of_flight_binning, name=name_for_starting_ws(ai), - scaleRaw=ai.scaleRaw, - scaleEmpty=ai.scaleEmpty, - subEmptyFromRaw=ai.subEmptyFromRaw + scaleRaw=ai.scale_raw_workspace, + scaleEmpty=ai.scale_empty_workspace, + subEmptyFromRaw=ai.subtract_empty_workspace_from_raw ) + first_detector, last_detector = [int(s) for s in ai.detectors.split('-')] cropedWs = cropAndMaskWorkspace( ws, - firstSpec=ai.firstSpec, - lastSpec=ai.lastSpec, - maskedDetectors=ai.maskedSpecAllNo, - maskTOFRange=ai.maskTOFRange + firstSpec=first_detector, + lastSpec=last_detector, + maskedDetectors=ai.mask_detectors, + maskTOFRange=ai.mask_time_of_flight_range ) profiles_table = create_profiles_table(cropedWs.name()+"_initial_parameters", ai) ipFilesPath = Path(handle_config.read_config_var("caching.ipfolder")) kwargs = { "InputWorkspace": cropedWs.name(), "InputProfiles": profiles_table.name(), - "InstrumentParametersFile": str(ipFilesPath / ai.ipfile), - "HRatioToLowestMass": ai.HToMassIdxRatio if hasattr(ai, 'HRatioToLowestMass') else 0, - "NumberOfIterations": int(ai.noOfMSIterations), - "InvalidDetectors": [int(det) for det in ai.maskedSpecAllNo], - "MultipleScatteringCorrection": ai.MSCorrectionFlag, + "InstrumentParametersFile": str(ipFilesPath / ai.instrument_parameters_file), + "HRatioToLowestMass": ai.intensity_ratio_of_hydrogen_to_lowest_mass if hasattr(ai, 'intensity_ratio_of_hydrogen_to_lowest_mass') else 0, + "NumberOfIterations": int(ai.number_of_iterations_for_corrections), + "InvalidDetectors": [int(det) for det in ai.mask_detectors], + "MultipleScatteringCorrection": ai.do_multiple_scattering_correction, "SampleVerticalWidth": ai.vertical_width, "SampleHorizontalWidth": ai.horizontal_width, "SampleThickness": ai.thickness, - "GammaCorrection": ai.GammaCorrectionFlag, + "GammaCorrection": ai.do_gamma_correction, "ModeRunning": scattering_type(ai), "TransmissionGuess": ai.transmission_guess, "MultipleScatteringOrder": int(ai.multiple_scattering_order), - "NumberOfEvents": int(ai.number_of_events), + "NumberOfEvents": int(ai.multiple_scattering_number_of_events), "Constraints": str(dill.dumps(ai.constraints)), "ResultsPath": str(self.results_save_path), "FiguresPath": str(self.fig_save_path), @@ -271,8 +272,8 @@ def _create_analysis_algorithm(self, ai): 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] + ai.masses, ai.initial_fitting_parameters[::3], ai.initial_fitting_parameters[1::3], ai.initial_fitting_parameters[2::3], + ai.fitting_bounds[::3], ai.fitting_bounds[1::3], ai.fitting_bounds[2::3] ): profiles.append(NeutronComptonProfile( label=str(float(mass)), mass=float(mass), intensity=float(intensity), width=float(width), center=float(center), @@ -294,11 +295,11 @@ def _save_ws_if_not_on_path(self, load_ai): ipFilesPath = Path(handle_config.read_config_var("caching.ipfolder")) - 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 ws_history_matches_inputs(load_ai.runs, load_ai.mode, load_ai.instrument_parameters_file, rawPath): + save_ws_from_load_vesuvio(load_ai.runs, load_ai.mode, str(ipFilesPath/load_ai.instrument_parameters_file), rawPath) - 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) + if not ws_history_matches_inputs(load_ai.empty_runs, load_ai.mode, load_ai.instrument_parameters_file, emptyPath): + save_ws_from_load_vesuvio(load_ai.empty_runs, load_ai.mode, str(ipFilesPath/load_ai.instrument_parameters_file), emptyPath) return rawPath, emptyPath @@ -309,13 +310,13 @@ def _set_output_paths(self, ai): # Build Filename based on ic corr = "" - if ai.GammaCorrectionFlag & (ai.noOfMSIterations > 0): + if ai.do_gamma_correction & (ai.number_of_iterations_for_corrections > 0): corr += "_GC" - if ai.MSCorrectionFlag & (ai.noOfMSIterations > 0): + if ai.do_multiple_scattering_correction & (ai.number_of_iterations_for_corrections > 0): corr += "_MS" fileName = ( - f"spec_{ai.firstSpec}-{ai.lastSpec}_iter_{ai.noOfMSIterations}{corr}" + ".npz" + f"spec_{ai.detectors.strip()}_iter_{ai.number_of_iterations_for_corrections}{corr}" + ".npz" ) fileNameYSpace = fileName + "_ySpaceFit" + ".npz" diff --git a/src/mvesuvio/util/analysis_helpers.py b/src/mvesuvio/util/analysis_helpers.py index 7b4cb07d..300a5670 100644 --- a/src/mvesuvio/util/analysis_helpers.py +++ b/src/mvesuvio/util/analysis_helpers.py @@ -23,8 +23,8 @@ def create_profiles_table(name, ai): table.addColumn(type="float", name="center") table.addColumn(type="str", name="center_bounds") 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] + ai.masses, ai.initial_fitting_parameters[::3], ai.initial_fitting_parameters[1::3], ai.initial_fitting_parameters[2::3], + ai.fitting_bounds[::3], ai.fitting_bounds[1::3], ai.fitting_bounds[2::3] ): table.addRow( [str(float(mass)), float(mass), float(intensity), str(list(intensity_bound)), diff --git a/tests/data/analysis/inputs/analysis_test.py b/tests/data/analysis/inputs/analysis_test.py index 46975a20..180be197 100644 --- a/tests/data/analysis/inputs/analysis_test.py +++ b/tests/data/analysis/inputs/analysis_test.py @@ -1,30 +1,29 @@ class SampleParameters: transmission_guess = 0.8537 # Experimental value from VesuvioTransmission - multiple_scattering_order, number_of_events = 2, 1.0e5 + multiple_scattering_order = 2 + multiple_scattering_number_of_events = 1.0e5 vertical_width, horizontal_width, thickness = 0.1, 0.1, 0.001 # Expressed in meters class BackwardAnalysisInputs(SampleParameters): run_this_scattering_type = False fit_in_y_space = False - ipfile = "ip2019.par" + instrument_parameters_file = "ip2019.par" runs = "43066-43076" empty_runs = "41876-41923" - spectra = "3-134" + detectors = "3-134" mode = "DoubleDifference" - tofBinning = "275.,1.,420" # Binning of ToF spectra - maskTOFRange = None - maskedSpecAllNo = [18, 34, 42, 43, 59, 60, 62, 118, 119, 133] - firstSpec = 3 # 3 - lastSpec = 134 # 134 - subEmptyFromRaw = True - scaleEmpty = 1 - scaleRaw = 1 + time_of_flight_binning = "275.,1.,420" # Binning of ToF spectra + mask_time_of_flight_range = None + mask_detectors = [18, 34, 42, 43, 59, 60, 62, 118, 119, 133] + subtract_empty_workspace_from_raw = True + scale_empty_workspace = 1 + scale_raw_workspace = 1 masses = [12, 16, 27] - initPars = [1, 12, 0.0, 1, 12, 0.0, 1, 12.5, 0.0] - bounds =[ + initial_fitting_parameters = [1, 12, 0.0, 1, 12, 0.0, 1, 12.5, 0.0] + fitting_bounds =[ [0, None], [8, 16], [-3, 1], @@ -37,34 +36,33 @@ class BackwardAnalysisInputs(SampleParameters): ] constraints = () - noOfMSIterations = 3 # 4 - MSCorrectionFlag = True - HToMassIdxRatio = 19.0620008206 # Set to zero or None when H is not present - GammaCorrectionFlag = False + number_of_iterations_for_corrections = 3 # 4 + do_multiple_scattering_correction = True + intensity_ratio_of_hydrogen_to_lowest_mass = 19.0620008206 # Set to zero or None when H is not present + do_gamma_correction = False class ForwardAnalysisInputs(SampleParameters): run_this_scattering_type = True fit_in_y_space = False - ipfile = "ip2018_3.par" + instrument_parameters_file = "ip2018_3.par" runs = "43066-43076" empty_runs = "43868-43911" spectra = "144-182" mode = "SingleDifference" - tofBinning = "110,1,430" # Binning of ToF spectra - maskTOFRange = None - maskedSpecAllNo = [173, 174, 179] - firstSpec = 144 # 144 - lastSpec = 182 # 182 - subEmptyFromRaw = False - scaleEmpty = 1 - scaleRaw = 1 + time_of_flight_binning = "110,1,430" # Binning of ToF spectra + mask_time_of_flight_range = None + mask_detectors = [173, 174, 179] + detectors = '144-182' + subtract_empty_workspace_from_raw= False + scale_empty_workspace= 1 + scale_raw_workspace= 1 masses = [1.0079, 12, 16, 27] - initPars = [1, 4.7, 0, 1, 12.71, 0.0, 1, 8.76, 0.0, 1, 13.897, 0.0] - bounds =[ + initial_fitting_parameters= [1, 4.7, 0, 1, 12.71, 0.0, 1, 8.76, 0.0, 1, 13.897, 0.0] + fitting_bounds =[ [0, None], [3, 6], [-3, 1], @@ -79,9 +77,9 @@ class ForwardAnalysisInputs(SampleParameters): [-3, 1], ] constraints = () - noOfMSIterations = 1 # 4 - MSCorrectionFlag = True - GammaCorrectionFlag = True + number_of_iterations_for_corrections= 1 # 4 + do_multiple_scattering_correction= True + do_gamma_correction= True class YSpaceFitInputs: diff --git a/tests/data/analysis/inputs/yspace_gauss_test.py b/tests/data/analysis/inputs/yspace_gauss_test.py index efe7c7ae..72cb1cca 100644 --- a/tests/data/analysis/inputs/yspace_gauss_test.py +++ b/tests/data/analysis/inputs/yspace_gauss_test.py @@ -3,10 +3,6 @@ ForwardAnalysisInputs, YSpaceFitInputs, ) -ForwardAnalysisInputs.noOfMSIterations = 1 -ForwardAnalysisInputs.firstSpec = 164 -ForwardAnalysisInputs.lastSpec = 175 -ForwardAnalysisInputs.maskedSpecAllNo = [173, 174] ForwardAnalysisInputs.fit_in_y_space = True BackwardAnalysisInputs.fit_in_y_space = False ForwardAnalysisInputs.run_this_scattering_type = True diff --git a/tests/data/analysis/inputs/yspace_gc_test.py b/tests/data/analysis/inputs/yspace_gc_test.py index cbd11323..d5fa66cd 100644 --- a/tests/data/analysis/inputs/yspace_gc_test.py +++ b/tests/data/analysis/inputs/yspace_gc_test.py @@ -4,10 +4,6 @@ YSpaceFitInputs, ) -ForwardAnalysisInputs.noOfMSIterations = 1 -ForwardAnalysisInputs.firstSpec = 164 -ForwardAnalysisInputs.lastSpec = 175 -ForwardAnalysisInputs.maskedSpecAllNo = [173, 174] ForwardAnalysisInputs.fit_in_y_space = True BackwardAnalysisInputs.fit_in_y_space = False ForwardAnalysisInputs.run_this_scattering_type = True From 91f77cfad2ed1c46d37053a1f6272ed5c8324554 Mon Sep 17 00:00:00 2001 From: GuiMacielPereira Date: Tue, 10 Dec 2024 13:35:37 +0000 Subject: [PATCH 2/7] Change range to mask in tof to use '-' Previously the range in tof to be masked was being specified by ','. It makes more sense to use '-' to describe the range. --- src/mvesuvio/util/analysis_helpers.py | 6 ++-- tests/unit/analysis/test_analysis_helpers.py | 31 +++++++++++++++++++- 2 files changed, 33 insertions(+), 4 deletions(-) diff --git a/src/mvesuvio/util/analysis_helpers.py b/src/mvesuvio/util/analysis_helpers.py index 300a5670..b2453bce 100644 --- a/src/mvesuvio/util/analysis_helpers.py +++ b/src/mvesuvio/util/analysis_helpers.py @@ -203,13 +203,13 @@ def cropAndMaskWorkspace(ws, firstSpec, lastSpec, maskedDetectors, maskTOFRange) OutputWorkspace=newWsName, ) - maskBinsWithZeros(wsCrop, maskTOFRange) # Used to mask resonance peaks + mask_time_of_flight_bins_with_zeros(wsCrop, maskTOFRange) # Used to mask resonance peaks MaskDetectors(Workspace=wsCrop, SpectraList=maskedDetectors) return wsCrop -def maskBinsWithZeros(ws, maskTOFRange): +def mask_time_of_flight_bins_with_zeros(ws, maskTOFRange): """ Masks a given TOF range on ws with zeros on dataY. Leaves errors dataE unchanged, as they are used by later treatments. @@ -220,7 +220,7 @@ def maskBinsWithZeros(ws, maskTOFRange): return dataX, dataY, dataE = extractWS(ws) - start, end = [int(s) for s in maskTOFRange.split(",")] + start, end = [float(s) for s in maskTOFRange.split("-")] assert ( start <= end ), "Start value for masking needs to be smaller or equal than end." diff --git a/tests/unit/analysis/test_analysis_helpers.py b/tests/unit/analysis/test_analysis_helpers.py index ab6f3211..84c7280d 100644 --- a/tests/unit/analysis/test_analysis_helpers.py +++ b/tests/unit/analysis/test_analysis_helpers.py @@ -5,7 +5,8 @@ import numpy.testing as nptest from mock import MagicMock from mvesuvio.util.analysis_helpers import extractWS, _convert_dict_to_table, \ - fix_profile_parameters, calculate_h_ratio, extend_range_of_array, numerical_third_derivative + fix_profile_parameters, calculate_h_ratio, extend_range_of_array, numerical_third_derivative, \ + mask_time_of_flight_bins_with_zeros from mantid.simpleapi import CreateWorkspace, DeleteWorkspace @@ -149,5 +150,33 @@ def test_numerical_third_derivative(self): expected_derivative = np.array([np.gradient(np.gradient(np.gradient(y_i, x_i), x_i), x_i)[6: -6] for y_i, x_i in zip(y, x) ]) np.testing.assert_allclose(numerical_derivative, expected_derivative, atol=1e-6) + + def test_mask_time_of_flight_bins_with_zeros(self): + data_x = np.arange(10).reshape(1, -1) * np.ones((3, 1)) + data_y = np.ones((3, 10)) + data_e = np.ones((3, 10)) + workspace_mock = MagicMock() + workspace_mock.extractX.return_value = data_x + workspace_mock.extractY.return_value = data_y + workspace_mock.extractE.return_value = data_e + + actual_data_x = np.zeros((3, 10)) + actual_data_y = np.zeros((3, 10)) + actual_data_e = np.zeros((3, 10)) + + workspace_mock.dataY.side_effect = lambda i: actual_data_y[i] + workspace_mock.dataX.side_effect = lambda i: actual_data_x[i] + workspace_mock.dataE.side_effect = lambda i: actual_data_e[i] + + workspace_mock.getNumberHistograms.return_value = 3 + mask_time_of_flight_bins_with_zeros(workspace_mock, '4.5-7.3') + + np.testing.assert_allclose(actual_data_x, data_x) + np.testing.assert_allclose(actual_data_e, data_e) + expected_data_y = np.ones((3, 10)) + expected_data_y[(data_x >= 4.5) & (data_x <= 7.3)] = 0 + np.testing.assert_allclose(actual_data_y, expected_data_y) + + if __name__ == "__main__": unittest.main() From c3c62e9115f6531f57e92544a7c625bd90466fba Mon Sep 17 00:00:00 2001 From: GuiMacielPereira Date: Tue, 10 Dec 2024 17:05:19 +0000 Subject: [PATCH 3/7] Rename variables for inputs of y space --- src/mvesuvio/analysis_fitting.py | 100 +++++++++--------- src/mvesuvio/config/analysis_inputs.py | 32 ++++-- tests/data/analysis/inputs/analysis_test.py | 18 ++-- .../data/analysis/inputs/yspace_gauss_test.py | 2 +- tests/data/analysis/inputs/yspace_gc_test.py | 4 +- 5 files changed, 84 insertions(+), 72 deletions(-) diff --git a/src/mvesuvio/analysis_fitting.py b/src/mvesuvio/analysis_fitting.py index 0b811d44..3002800f 100644 --- a/src/mvesuvio/analysis_fitting.py +++ b/src/mvesuvio/analysis_fitting.py @@ -27,13 +27,13 @@ def fitInYSpaceProcedure(yFitIC, IC, wsTOF): wsTOFMass0 = subtract_profiles_except_lightest(IC, wsTOF) - if yFitIC.subtractFSE: + if yFitIC.subtract_calculated_fse_from_data: wsFSEMass0 = find_ws_name_fse_first_mass(IC) wsTOFMass0 = Minus(wsTOFMass0, wsFSEMass0, OutputWorkspace=wsTOFMass0.name() + "_fse") wsJoY, wsJoYAvg = ySpaceReduction(wsTOFMass0, IC.masses[0], yFitIC, IC) - if yFitIC.symmetrisationFlag: + if yFitIC.do_symmetrisation: wsJoYAvg = symmetrizeWs(wsJoYAvg) fitProfileMinuit(yFitIC, wsJoYAvg, wsResSum) @@ -44,7 +44,7 @@ def fitInYSpaceProcedure(yFitIC, IC, wsTOF): yfitResults = ResultsYFitObject(IC, yFitIC, wsTOF.name(), wsJoYAvg.name()) yfitResults.save() - if yFitIC.globalFit: + if yFitIC.do_global_fit: runGlobalFit(wsJoY, wsRes, IC, yFitIC) return yfitResults @@ -76,7 +76,7 @@ def calculateMantidResolutionFirstMass(IC, yFitIC, ws): ) Rebin( InputWorkspace="tmp", - Params=yFitIC.rebinParametersForYSpaceFit, + Params=yFitIC.range_for_rebinning_in_y_space, OutputWorkspace="tmp", ) @@ -129,10 +129,10 @@ def ySpaceReduction(wsTOF, mass0, yFitIC, ic): ws_name_lightest_mass = ic.name + '_' + str(ic.number_of_iterations_for_corrections) + '_' + str(min(ic.masses)) + '_ncp' ncp = mtd[ws_name_lightest_mass].extractY() - rebinPars = yFitIC.rebinParametersForYSpaceFit + rebinPars = yFitIC.range_for_rebinning_in_y_space if np.any(np.all(wsTOF.extractY() == 0, axis=0)): # Masked columns present - if yFitIC.maskTypeProcedure == "NAN": + if yFitIC.mask_zeros_with == "nan": # Build special workspace to store accumulated points wsJoY = convertToYSpace(wsTOF, mass0) xp = buildXRangeFromRebinPars(yFitIC) @@ -152,14 +152,14 @@ def ySpaceReduction(wsTOF, mass0, yFitIC, ic): wsJoYAvg = weightedAvgXBins(wsJoYN, xp) return wsJoYN, wsJoYAvg - elif yFitIC.maskTypeProcedure == "NCP": + elif yFitIC.mask_zeros_with == "ncp": wsTOF = replaceZerosWithNCP(wsTOF, ncp) else: raise ValueError( """ Masked TOF bins were found but no valid procedure in y-space fit was selected. - Options: 'NAN', 'NCP' + Options: 'nan', 'ncp' """ ) @@ -206,7 +206,7 @@ def replaceZerosWithNCP(ws, ncp): def buildXRangeFromRebinPars(yFitIC): # Range used in case mask is set to NAN first, step, last = [ - float(s) for s in yFitIC.rebinParametersForYSpaceFit.split(",") + float(s) for s in yFitIC.range_for_rebinning_in_y_space.split(",") ] xp = np.arange(first, last, step) + step / 2 # Correction to match Mantid range return xp @@ -519,7 +519,7 @@ def fitProfileMinuit(yFitIC, wsYSpaceSym, wsRes): resX, resY, resE = extractFirstSpectra(wsRes) assert np.all(dataX == resX), "Resolution should operate on the same range as DataX" - model, defaultPars, sharedPars = selectModelAndPars(yFitIC.fitModel) + model, defaultPars, sharedPars = selectModelAndPars(yFitIC.fitting_model) xDelta, resDense = oddPointsRes(resX, resY) @@ -544,11 +544,11 @@ def convolvedModel(x, y0, *pars): m = Minuit(costFun, **defaultPars) m.limits["A"] = (0, None) - if yFitIC.fitModel == "DOUBLE_WELL": + if yFitIC.fitting_model== "doublewell": m.limits["d"] = (0, None) m.limits["R"] = (0, None) - if yFitIC.fitModel == "SINGLE_GAUSSIAN": + if yFitIC.fitting_model== "gauss": m.simplex() m.migrad() @@ -616,7 +616,7 @@ def selectModelAndPars(modelFlag): The defaultPars should be in the same order as the signature of the function """ - if modelFlag == "SINGLE_GAUSSIAN": + if modelFlag == "gauss": def model(x, A, x0, sigma): return ( @@ -629,7 +629,7 @@ def model(x, A, x0, sigma): defaultPars = {"A": 1, "x0": 0, "sigma": 5} sharedPars = ["sigma"] # Used only in Global fit - elif modelFlag == "GC_C4_C6": + elif modelFlag == "gcc4c6": def model(x, A, x0, sigma1, c4, c6): return ( @@ -659,7 +659,7 @@ def model(x, A, x0, sigma1, c4, c6): defaultPars = {"A": 1, "x0": 0, "sigma1": 6, "c4": 0, "c6": 0} sharedPars = ["sigma1", "c4", "c6"] # Used only in Global fit - elif modelFlag == "GC_C4": + elif modelFlag == "gcc4": def model(x, A, x0, sigma1, c4): return ( @@ -681,7 +681,7 @@ def model(x, A, x0, sigma1, c4): defaultPars = {"A": 1, "x0": 0, "sigma1": 6, "c4": 0} sharedPars = ["sigma1", "c4"] # Used only in Global fit - elif modelFlag == "GC_C6": + elif modelFlag == "gcc6": def model(x, A, x0, sigma1, c6): return ( @@ -704,7 +704,7 @@ def model(x, A, x0, sigma1, c6): defaultPars = {"A": 1, "x0": 0, "sigma1": 6, "c6": 0} sharedPars = ["sigma1", "c6"] # Used only in Global fit - elif modelFlag == "DOUBLE_WELL": + elif modelFlag == "doublewell": def model(x, A, d, R, sig1, sig2): # h = 2.04 @@ -744,7 +744,7 @@ def model(x, A, d, R, sig1, sig2): } # TODO: Starting parameters and bounds? sharedPars = ["d", "R", "sig1", "sig2"] # Only varying parameter is amplitude A - elif modelFlag == "ANSIO_GAUSSIAN": + elif modelFlag == "ansiogauss": # Ansiotropic case def model(x, A, sig1, sig2): # h = 2.04 @@ -765,7 +765,7 @@ def model(x, A, sig1, sig2): defaultPars = {"A": 1, "sig1": 3, "sig2": 5} sharedPars = ["sig1", "sig2"] - elif modelFlag == "Gaussian3D": + elif modelFlag == "gauss3d": def model(x, A, sig_x, sig_y, sig_z): y = x[:, np.newaxis, np.newaxis] @@ -795,7 +795,7 @@ def model(x, A, sig_x, sig_y, sig_z): else: raise ValueError( - "Fitting Model not recognized, available options: 'SINGLE_GAUSSIAN', 'GC_C4_C6', 'GC_C4'" + "Fitting Model not recognized, available options: 'gauss', 'gcc4c6', 'gcc4', 'gcc6', 'ansiogauss' gauss3d'" ) print("\nShared Parameters: ", [key for key in sharedPars]) @@ -915,7 +915,7 @@ def runMinos(mObj, yFitIC, constrFunc, wsName): errors = list(mObj.errors) # If minos is set not to run, ouput columns with zeros on minos errors - if not (yFitIC.runMinos): + if not (yFitIC.run_minos): minosAutoErr = list(np.zeros((len(parameters), 2))) minosManErr = list(np.zeros((len(parameters), 2))) return parameters, values, errors, minosAutoErr, minosManErr @@ -927,7 +927,7 @@ def runMinos(mObj, yFitIC, constrFunc, wsName): bestFitErrs[p] = e if ( - yFitIC.fitModel == "SINGLE_GAUSSIAN" + yFitIC.fitting_model== "gauss" ): # Case with no positivity constraint, can use automatic minos() mObj.minos() me = mObj.merrors @@ -938,13 +938,13 @@ def runMinos(mObj, yFitIC, constrFunc, wsName): minosAutoErr.append([me[p].lower, me[p].upper]) minosManErr = list(np.zeros(np.array(minosAutoErr).shape)) - if yFitIC.showPlots: + if yFitIC.show_plots: plotAutoMinos(mObj, wsName) else: # Case with positivity constraint on function, use manual implementation # Changes values of minuit obj m, do not use m below this point merrors, fig = runAndPlotManualMinos( - mObj, constrFunc, bestFitVals, bestFitErrs, yFitIC.showPlots + mObj, constrFunc, bestFitVals, bestFitErrs, yFitIC.show_plots ) # Same as above, but the other way around @@ -953,7 +953,7 @@ def runMinos(mObj, yFitIC, constrFunc, wsName): minosManErr.append(merrors[p]) minosAutoErr = list(np.zeros(np.array(minosManErr).shape)) - if yFitIC.showPlots: + if yFitIC.show_plots: fig.canvas.manager.set_window_title(wsName + "_Manual_Implementation_MINOS") fig.show() @@ -990,7 +990,6 @@ def runAndPlotManualMinos(minuitObj, constrFunc, bestFitVals, bestFitErrs, showP ) merrors[p] = np.array([lerr, uerr]) - # if showPlots: # Hide plots not in use: for ax in axs.flat: if not ax.lines: # If empty list @@ -999,7 +998,6 @@ def runAndPlotManualMinos(minuitObj, constrFunc, bestFitVals, bestFitErrs, showP # ALl axes share same legend, so set figure legend to first axis handle, label = axs[0, 0].get_legend_handles_labels() fig.legend(handle, label, loc="lower right") - # fig.show() return merrors, fig @@ -1226,13 +1224,13 @@ def oddPointsRes(x, res): def fitProfileMantidFit(yFitIC, wsYSpaceSym, wsRes): print("\nFitting on the sum of spectra in the West domain ...\n") for minimizer in ["Levenberg-Marquardt", "Simplex"]: - if yFitIC.fitModel == "SINGLE_GAUSSIAN": + if yFitIC.fitting_model== "gauss": function = f"""composite=Convolution,FixResolution=true,NumDeriv=true; name=Resolution,Workspace={wsRes.name()},WorkspaceIndex=0; name=UserFunction,Formula=y0 + A*exp( -(x-x0)^2/2/sigma^2)/(2*3.1415*sigma^2)^0.5, y0=0,A=1,x0=0,sigma=5, ties=()""" - elif yFitIC.fitModel == "GC_C4_C6": + elif yFitIC.fitting_model== "gcc4c6": function = f""" composite=Convolution,FixResolution=true,NumDeriv=true; name=Resolution,Workspace={wsRes.name()},WorkspaceIndex=0,X=(),Y=(); @@ -1241,7 +1239,7 @@ def fitProfileMantidFit(yFitIC, wsYSpaceSym, wsRes): (64*((x-x0)/sqrt(2)/sigma1)^6 - 480*((x-x0)/sqrt(2)/sigma1)^4 + 720*((x-x0)/sqrt(2)/sigma1)^2 - 120)), y0=0, A=1,x0=0,sigma1=4.0,c4=0.0,c6=0.0,ties=(),constraints=(07s} = {v:>8.4f} \u00B1 {e:<8.4f}") print("\n") - if yFitIC.showPlots: + if yFitIC.show_plots: plotGlobalFit(dataX, dataY, dataE, m, totCost, wsYSpace.name()) return np.array(m.values), np.array( @@ -1533,7 +1533,7 @@ def groupDetectors(ipData, yFitIC): checkNGroupsValid(yFitIC, ipData) - print(f"\nNumber of gropus: {yFitIC.nGlobalFitGroups}") + print(f"\nNumber of gropus: {yFitIC.number_of_global_fit_groups}") L1 = ipData[:, -1].copy() theta = ipData[:, 2].copy() @@ -1547,7 +1547,7 @@ def groupDetectors(ipData, yFitIC): points = np.vstack((L1, theta)).T assert points.shape == (len(L1), 2), "Wrong shape." # Initial centers of groups - startingIdxs = np.linspace(0, len(points) - 1, yFitIC.nGlobalFitGroups).astype(int) + startingIdxs = np.linspace(0, len(points) - 1, yFitIC.number_of_global_fit_groups).astype(int) centers = points[ startingIdxs, : ] # Centers of cluster groups, NOT fitting parameter @@ -1558,7 +1558,7 @@ def groupDetectors(ipData, yFitIC): clusters = kMeansClustering(points, centers) idxList = formIdxList(clusters) - if yFitIC.showPlots: + if yFitIC.show_plots: fig, ax = plt.subplots(tight_layout=True, subplot_kw={"projection": "mantid"}) fig.canvas.manager.set_window_title("Grouping of detectors") plotFinalGroups(ax, ipData, idxList) @@ -1571,17 +1571,17 @@ def checkNGroupsValid(yFitIC, ipData): nSpectra = len(ipData) # Number of spectra in the workspace - if yFitIC.nGlobalFitGroups == "ALL": - yFitIC.nGlobalFitGroups = nSpectra + if yFitIC.number_of_global_fit_groups == "all": + yFitIC.number_of_global_fit_groups = nSpectra else: assert isinstance( - yFitIC.nGlobalFitGroups, int + yFitIC.number_of_global_fit_groups, int ), "Number of global groups needs to be an integer." assert ( - yFitIC.nGlobalFitGroups <= nSpectra + yFitIC.number_of_global_fit_groups <= nSpectra ), "Number of global groups needs to be less or equal to the no of unmasked spectra." assert ( - yFitIC.nGlobalFitGroups > 0 + yFitIC.number_of_global_fit_groups > 0 ), "NUmber of global groups needs to be bigger than zero" return @@ -1721,7 +1721,7 @@ def avgWeightDetGroups(dataX, dataY, dataE, dataRes, idxList, yFitIC): np.all(dataY == 0, axis=1) ), f"Input data should not include masked spectra at: {np.argwhere(np.all(dataY==0, axis=1))}" - if yFitIC.maskTypeProcedure == "NAN": + if yFitIC.mask_zeros_with == "nan": return avgGroupsWithBins(dataX, dataY, dataE, dataRes, idxList, yFitIC) # Use Default for unmasked or NCP masked diff --git a/src/mvesuvio/config/analysis_inputs.py b/src/mvesuvio/config/analysis_inputs.py index 06bacb16..0bc40390 100644 --- a/src/mvesuvio/config/analysis_inputs.py +++ b/src/mvesuvio/config/analysis_inputs.py @@ -11,7 +11,7 @@ class SampleParameters: @dataclass class BackwardAnalysisInputs(SampleParameters): - run_this_scattering_type = True + run_this_scattering_type = False fit_in_y_space = False runs = "43066-43076" @@ -83,15 +83,27 @@ class ForwardAnalysisInputs(SampleParameters): @dataclass class YSpaceFitInputs: - showPlots = True - symmetrisationFlag = False - subtractFSE = True - rebinParametersForYSpaceFit = "-25, 0.5, 25" # Needs to be symetric - fitModel = "SINGLE_GAUSSIAN" # Options: 'SINGLE_GAUSSIAN', 'GC_C4', 'GC_C6', 'GC_C4_C6', 'DOUBLE_WELL', 'ANSIO_GAUSSIAN', 'Gaussian3D' - runMinos = True - globalFit = True # Performs global fit with Minuit by default - nGlobalFitGroups = 4 # Number or string "ALL" - maskTypeProcedure = "NAN" # Options: 'NCP', 'NAN', None + show_plots = True + do_symmetrisation = False + subtract_calculated_fse_from_data = True + range_for_rebinning_in_y_space = "-25, 0.5, 25" # Needs to be symetric + # Fitting model options + # 'gauss': Single Gaussian + # 'gcc4': Gram-Charlier with C4 parameter + # 'gcc6': Gram-Charlier with C6 parameter + # 'doublewell': Double Well function + # 'gauss3D': 3-Dimensional Gaussian + fitting_model = "gauss" + run_minos = True + do_global_fit = True # Performs global fit with Minuit by default + # Number of groups of detectors to perform global (simultaneous) fit on + # Either an integer less than the number of detectors + # or option 'all', which does not form groups and fits all spectra simultaneously and individualy + number_of_global_fit_groups = 4 + # Type of masking + # 'nan': Zeros in workspace being fit are ignored + # 'ncp': Zeros in workspace being fit are replaced by the fitted neutron compton profile + mask_zeros_with = "nan" ######################## diff --git a/tests/data/analysis/inputs/analysis_test.py b/tests/data/analysis/inputs/analysis_test.py index 180be197..8cb9a915 100644 --- a/tests/data/analysis/inputs/analysis_test.py +++ b/tests/data/analysis/inputs/analysis_test.py @@ -83,15 +83,15 @@ class ForwardAnalysisInputs(SampleParameters): class YSpaceFitInputs: - showPlots = False - symmetrisationFlag = True - subtractFSE = False - rebinParametersForYSpaceFit = "-20, 0.5, 20" # Needs to be symetric - fitModel = "SINGLE_GAUSSIAN" - runMinos = True - globalFit = None - nGlobalFitGroups = 4 - maskTypeProcedure = None + show_plots = False + do_symmetrisation = True + subtract_calculated_fse_from_data = False + range_for_rebinning_in_y_space = "-20, 0.5, 20" # Needs to be symetric + fitting_model = "SINGLE_GAUSSIAN" + run_minos = True + do_global_fit = None + number_of_global_fit_groups = 4 + mask_zeros_with = None if (__name__ == "__main__") or (__name__ == "mantidqt.widgets.codeeditor.execution"): diff --git a/tests/data/analysis/inputs/yspace_gauss_test.py b/tests/data/analysis/inputs/yspace_gauss_test.py index 72cb1cca..2ec5c952 100644 --- a/tests/data/analysis/inputs/yspace_gauss_test.py +++ b/tests/data/analysis/inputs/yspace_gauss_test.py @@ -7,7 +7,7 @@ BackwardAnalysisInputs.fit_in_y_space = False ForwardAnalysisInputs.run_this_scattering_type = True BackwardAnalysisInputs.run_this_scattering_type = False -YSpaceFitInputs.fitModel = "SINGLE_GAUSSIAN" +YSpaceFitInputs.fitting_model = "gauss" if (__name__ == "__main__") or (__name__ == "mantidqt.widgets.codeeditor.execution"): import mvesuvio diff --git a/tests/data/analysis/inputs/yspace_gc_test.py b/tests/data/analysis/inputs/yspace_gc_test.py index d5fa66cd..d2b2a7d4 100644 --- a/tests/data/analysis/inputs/yspace_gc_test.py +++ b/tests/data/analysis/inputs/yspace_gc_test.py @@ -8,8 +8,8 @@ BackwardAnalysisInputs.fit_in_y_space = False ForwardAnalysisInputs.run_this_scattering_type = True BackwardAnalysisInputs.run_this_scattering_type = False -YSpaceFitInputs.fitModel = "GC_C4_C6" -YSpaceFitInputs.symmetrisationFlag = False +YSpaceFitInputs.fitting_model = "gcc4c6" +YSpaceFitInputs.do_symmetrisation = False if (__name__ == "__main__") or (__name__ == "mantidqt.widgets.codeeditor.execution"): From 88f8a01f5c295a1331193d00700941bb0a00ee74 Mon Sep 17 00:00:00 2001 From: GuiMacielPereira Date: Tue, 10 Dec 2024 18:03:21 +0000 Subject: [PATCH 4/7] Prevent new fit model from overriding old ws Previously the workspaces for the y-space fit were being overwritten when the fitting model changed. This change changes the name of the workspace which makes sure the workspaces do not get overwritten. --- src/mvesuvio/analysis_fitting.py | 26 ++++++++++---------------- 1 file changed, 10 insertions(+), 16 deletions(-) diff --git a/src/mvesuvio/analysis_fitting.py b/src/mvesuvio/analysis_fitting.py index 3002800f..8f052fc1 100644 --- a/src/mvesuvio/analysis_fitting.py +++ b/src/mvesuvio/analysis_fitting.py @@ -580,9 +580,10 @@ def constrFunc(*pars): # Constrain physical model before convolution dataYSigma *= chi2 # Weight the confidence band Residuals = dataY - dataYFit + wsYSpaceSym = CloneWorkspace(wsYSpaceSym, OutputWorkspace=wsYSpaceSym.name()+'_minuit_'+yFitIC.fitting_model) # Create workspace to store best fit curve and errors on the fit wsMinFit = createFitResultsWorkspace( - wsYSpaceSym, dataX, dataY, dataE, dataYFit, dataYSigma, Residuals + wsYSpaceSym, dataX, dataY, dataE, dataYFit, dataYSigma, Residuals ) saveMinuitPlot(yFitIC, wsMinFit, m) @@ -851,23 +852,16 @@ def ndata(self): def createFitResultsWorkspace( - wsYSpaceSym, dataX, dataY, dataE, dataYFit, dataYSigma, Residuals + wsYSpaceSym, dataX, dataY, dataE, dataYFit, dataYSigma, Residuals ): """Creates workspace similar to the ones created by Mantid Fit.""" - ws_fit_profile = CreateWorkspace( - DataX=dataX, - DataY=dataY, - DataE=dataE, - NSpec=1, - OutputWorkspace=wsYSpaceSym.name() + "_minuit", - ) ws_fit_complete = CreateWorkspace( DataX=np.concatenate((dataX, dataX, dataX)), DataY=np.concatenate((dataY, dataYFit, Residuals)), DataE=np.concatenate((dataE, dataYSigma, np.zeros(len(dataE)))), NSpec=3, - OutputWorkspace=wsYSpaceSym.name() + "_minuit_Workspace", + OutputWorkspace=wsYSpaceSym.name() + "_Workspace", ) return ws_fit_complete @@ -896,7 +890,7 @@ def saveMinuitPlot(yFitIC, wsMinuitFit, mObj): def createCorrelationTableWorkspace(wsYSpaceSym, parameters, corrMatrix): tableWS = CreateEmptyTableWorkspace( - OutputWorkspace=wsYSpaceSym.name() + "_minuit_NormalizedCovarianceMatrix" + OutputWorkspace=wsYSpaceSym.name() + "_NormalizedCovarianceMatrix" ) tableWS.setTitle("Minuit Fit") tableWS.addColumn(type="str", name="Name") @@ -1178,7 +1172,7 @@ def createFitParametersTableWorkspace( ): # Create Parameters workspace tableWS = CreateEmptyTableWorkspace( - OutputWorkspace=wsYSpaceSym.name() + "_minuit_Parameters" + OutputWorkspace=wsYSpaceSym.name() + "_Parameters" ) tableWS.setTitle("Minuit Fit") tableWS.addColumn(type="str", name="Name") @@ -1267,7 +1261,7 @@ def fitProfileMantidFit(yFitIC, wsYSpaceSym, wsRes): ) suffix = 'lm' if minimizer=="Levenberg-Marquardt" else minimizer.lower() - outputName = wsYSpaceSym.name() + "_" + suffix + outputName = wsYSpaceSym.name() + "_" + suffix + "_" + yFitIC.fitting_model CloneWorkspace(InputWorkspace=wsYSpaceSym, OutputWorkspace=outputName) Fit( @@ -1337,19 +1331,19 @@ def __init__(self, ic, yFitIC, wsFinalName, wsYSpaceAvgName): poptList = [] perrList = [] try: - wsFitMinuit = mtd[wsJoYAvg.name() + "_minuit_Parameters"] + wsFitMinuit = mtd[wsJoYAvg.name() + "_minuit_" + yFitIC.fitting_model + "_Parameters"] poptList.append(wsFitMinuit.column("Value")) perrList.append(wsFitMinuit.column("Error")) except: pass try: - wsFitLM = mtd[wsJoYAvg.name() + "_lm_Parameters"] + wsFitLM = mtd[wsJoYAvg.name() + "_lm_" + yFitIC.fitting_model + "_Parameters"] poptList.append(wsFitLM.column("Value")) perrList.append(wsFitLM.column("Error")) except: pass try: - wsFitSimplex = mtd[wsJoYAvg.name() + "_simplex_Parameters"] + wsFitSimplex = mtd[wsJoYAvg.name() + "_simplex_" + yFitIC.fitting_model + "_Parameters"] poptList.append(wsFitSimplex.column("Value")) perrList.append(wsFitSimplex.column("Error")) except: From 6cfcba84e827350a4f0d07a63d9e5454122afdf4 Mon Sep 17 00:00:00 2001 From: GuiMacielPereira Date: Wed, 11 Dec 2024 11:41:42 +0000 Subject: [PATCH 5/7] Save global fit and minos plots Save some plots by default that were not being saved --- src/mvesuvio/analysis_fitting.py | 26 ++++++++++++++++---------- 1 file changed, 16 insertions(+), 10 deletions(-) diff --git a/src/mvesuvio/analysis_fitting.py b/src/mvesuvio/analysis_fitting.py index 8f052fc1..68c44b3b 100644 --- a/src/mvesuvio/analysis_fitting.py +++ b/src/mvesuvio/analysis_fitting.py @@ -876,6 +876,7 @@ def saveMinuitPlot(yFitIC, wsMinuitFit, mObj): fig, ax = plt.subplots(subplot_kw={"projection": "mantid"}) ax.errorbar(wsMinuitFit, "k.", wkspIndex=0, label="Weighted Avg") ax.errorbar(wsMinuitFit, "r-", wkspIndex=1, label=leg) + ax.plot(wsMinuitFit, "b.", wkspIndex=2, label="Residuals") ax.set_xlabel("YSpace") ax.set_ylabel("Counts") ax.set_title("Minuit Fit") @@ -933,12 +934,12 @@ def runMinos(mObj, yFitIC, constrFunc, wsName): minosManErr = list(np.zeros(np.array(minosAutoErr).shape)) if yFitIC.show_plots: - plotAutoMinos(mObj, wsName) + plotAutoMinos(mObj, wsName, yFitIC) else: # Case with positivity constraint on function, use manual implementation # Changes values of minuit obj m, do not use m below this point merrors, fig = runAndPlotManualMinos( - mObj, constrFunc, bestFitVals, bestFitErrs, yFitIC.show_plots + mObj, constrFunc, bestFitVals, bestFitErrs, yFitIC.show_plots, yFitIC, wsName ) # Same as above, but the other way around @@ -948,13 +949,12 @@ def runMinos(mObj, yFitIC, constrFunc, wsName): minosAutoErr = list(np.zeros(np.array(minosManErr).shape)) if yFitIC.show_plots: - fig.canvas.manager.set_window_title(wsName + "_Manual_Implementation_MINOS") fig.show() return parameters, values, errors, minosAutoErr, minosManErr -def runAndPlotManualMinos(minuitObj, constrFunc, bestFitVals, bestFitErrs, showPlots): +def runAndPlotManualMinos(minuitObj, constrFunc, bestFitVals, bestFitErrs, showPlots, yFitIC, wsName): """ Runs brute implementation of minos algorithm and plots the profile for each parameter along the way. @@ -975,7 +975,7 @@ def runAndPlotManualMinos(minuitObj, constrFunc, bestFitVals, bestFitErrs, showP figsize=figsize, subplot_kw={"projection": "mantid"}, ) # subplot_kw={'projection':'mantid'} - # fig.canvas.manager.set_window_title("Plot of Manual Implementation MINOS") + fig.canvas.manager.set_window_title(wsName + "_minos") merrors = {} for p, ax in zip(minuitObj.parameters, axs.flat): @@ -992,6 +992,8 @@ def runAndPlotManualMinos(minuitObj, constrFunc, bestFitVals, bestFitErrs, showP # ALl axes share same legend, so set figure legend to first axis handle, label = axs[0, 0].get_legend_handles_labels() fig.legend(handle, label, loc="lower right") + savePath = yFitIC.figSavePath / fig.canvas.manager.get_window_title() + plt.savefig(savePath, bbox_inches="tight") return merrors, fig @@ -1104,7 +1106,7 @@ def errsFromMinosCurve(varSpace, varVal, fValsScipy, fValsMin, dChi2=1): return lerr, uerr -def plotAutoMinos(minuitObj, wsName): +def plotAutoMinos(minuitObj, wsName, yFitIC): # Set format of subplots height = 2 width = int(np.ceil(len(minuitObj.parameters) / 2)) @@ -1117,7 +1119,7 @@ def plotAutoMinos(minuitObj, wsName): figsize=figsize, subplot_kw={"projection": "mantid"}, ) - fig.canvas.manager.set_window_title(wsName + "_Plot_Automatic_MINOS") + fig.canvas.manager.set_window_title(wsName + "_autominos") for p, ax in zip(minuitObj.parameters, axs.flat): loc, fvals, status = minuitObj.mnprofile(p, bound=2) @@ -1137,6 +1139,8 @@ def plotAutoMinos(minuitObj, wsName): # ALl axes share same legend, so set figure legend to first axis handle, label = axs[0, 0].get_legend_handles_labels() fig.legend(handle, label, loc="lower right") + savePath = yFitIC.figSavePath / fig.canvas.manager.get_window_title() + plt.savefig(savePath, bbox_inches="tight") fig.show() @@ -1476,7 +1480,7 @@ def constr(*pars): print("\n") if yFitIC.show_plots: - plotGlobalFit(dataX, dataY, dataE, m, totCost, wsYSpace.name()) + plotGlobalFit(dataX, dataY, dataE, m, totCost, wsYSpace.name(), yFitIC) return np.array(m.values), np.array( m.errors @@ -1847,7 +1851,7 @@ def minuitInitialParameters(defaultPars, sharedPars, nSpec): return initPars -def plotGlobalFit(dataX, dataY, dataE, mObj, totCost, wsName): +def plotGlobalFit(dataX, dataY, dataE, mObj, totCost, wsName, yFitIC): if len(dataY) > 10: print("\nToo many axes to show in figure, skipping the plot ...\n") return @@ -1860,7 +1864,7 @@ def plotGlobalFit(dataX, dataY, dataE, mObj, totCost, wsName): tight_layout=True, subplot_kw={"projection": "mantid"}, ) - fig.canvas.manager.set_window_title(wsName + "_Plot_of_Global_Fit") + fig.canvas.manager.set_window_title(wsName + "_fitglobal") # Data used in Global Fit for i, (x, y, yerr, ax) in enumerate(zip(dataX, dataY, dataE, axs.flat)): @@ -1882,5 +1886,7 @@ def plotGlobalFit(dataX, dataY, dataE, mObj, totCost, wsName): ax.fill_between(x, yfit, label="\n".join(leg), alpha=0.4) ax.legend() + savePath = yFitIC.figSavePath / fig.canvas.manager.get_window_title() + plt.savefig(savePath, bbox_inches="tight") fig.show() return From 6e733e97433c8a9c07d15d761510b21f78a10bcd Mon Sep 17 00:00:00 2001 From: GuiMacielPereira Date: Wed, 11 Dec 2024 14:59:27 +0000 Subject: [PATCH 6/7] Save fit workspaces as ascii Added function to save workspaces containing the result in y space fit in Ascii format, and made minor esthetic improvements to plots that are also saved. --- src/mvesuvio/analysis_fitting.py | 48 +++++++++++++----- src/mvesuvio/main/run_routine.py | 2 +- .../analysis_test_empty_backward.nxs | Bin .../analysis_test_empty_forward.nxs | Bin .../analysis_test_raw_backward.nxs | Bin .../analysis_test_raw_forward.nxs | Bin .../yspace_gauss_test_empty_backward.nxs | Bin .../yspace_gauss_test_empty_forward.nxs | Bin .../yspace_gauss_test_raw_backward.nxs | Bin .../yspace_gauss_test_raw_forward.nxs | Bin .../yspace_gc_test_empty_backward.nxs | Bin .../yspace_gc_test_empty_forward.nxs | Bin .../yspace_gc_test_raw_backward.nxs | Bin .../yspace_gc_test_raw_forward.nxs | Bin 14 files changed, 36 insertions(+), 14 deletions(-) rename tests/data/analysis/inputs/analysis_test/{input_ws => input_workspaces}/analysis_test_empty_backward.nxs (100%) rename tests/data/analysis/inputs/analysis_test/{input_ws => input_workspaces}/analysis_test_empty_forward.nxs (100%) rename tests/data/analysis/inputs/analysis_test/{input_ws => input_workspaces}/analysis_test_raw_backward.nxs (100%) rename tests/data/analysis/inputs/analysis_test/{input_ws => input_workspaces}/analysis_test_raw_forward.nxs (100%) rename tests/data/analysis/inputs/yspace_gauss_test/{input_ws => input_workspaces}/yspace_gauss_test_empty_backward.nxs (100%) rename tests/data/analysis/inputs/yspace_gauss_test/{input_ws => input_workspaces}/yspace_gauss_test_empty_forward.nxs (100%) rename tests/data/analysis/inputs/yspace_gauss_test/{input_ws => input_workspaces}/yspace_gauss_test_raw_backward.nxs (100%) rename tests/data/analysis/inputs/yspace_gauss_test/{input_ws => input_workspaces}/yspace_gauss_test_raw_forward.nxs (100%) rename tests/data/analysis/inputs/yspace_gc_test/{input_ws => input_workspaces}/yspace_gc_test_empty_backward.nxs (100%) rename tests/data/analysis/inputs/yspace_gc_test/{input_ws => input_workspaces}/yspace_gc_test_empty_forward.nxs (100%) rename tests/data/analysis/inputs/yspace_gc_test/{input_ws => input_workspaces}/yspace_gc_test_raw_backward.nxs (100%) rename tests/data/analysis/inputs/yspace_gc_test/{input_ws => input_workspaces}/yspace_gc_test_raw_forward.nxs (100%) diff --git a/src/mvesuvio/analysis_fitting.py b/src/mvesuvio/analysis_fitting.py index 68c44b3b..725e8dd2 100644 --- a/src/mvesuvio/analysis_fitting.py +++ b/src/mvesuvio/analysis_fitting.py @@ -47,6 +47,7 @@ def fitInYSpaceProcedure(yFitIC, IC, wsTOF): if yFitIC.do_global_fit: runGlobalFit(wsJoY, wsRes, IC, yFitIC) + save_workspaces(yFitIC) return yfitResults @@ -880,8 +881,11 @@ def saveMinuitPlot(yFitIC, wsMinuitFit, mObj): ax.set_xlabel("YSpace") ax.set_ylabel("Counts") ax.set_title("Minuit Fit") + ax.grid() ax.legend() + + fileName = wsMinuitFit.name() + ".pdf" savePath = yFitIC.figSavePath / fileName plt.savefig(savePath, bbox_inches="tight") @@ -1468,23 +1472,16 @@ def constr(*pars): # Explicitly calculate errors m.hesse() - chi2 = m.fval / ( - np.sum(dataE != 0) - m.nfit - ) # Number of non zero points (considered in the fit) minus no of parameters - print(f"Value of Chi2/ndof: {chi2:.2f}") - print(f"Migrad Minimum valid: {m.valid}") + # Number of non zero points (considered in the fit) minus no of parameters + chi2 = m.fval / (np.sum(dataE != 0) - m.nfit) - print("\nResults of Global Fit:\n") - for p, v, e in zip(m.parameters, m.values, m.errors): - print(f"{p:>7s} = {v:>8.4f} \u00B1 {e:<8.4f}") - print("\n") + create_table_for_global_fit_parameters(wsYSpace.name(), m, chi2) if yFitIC.show_plots: plotGlobalFit(dataX, dataY, dataE, m, totCost, wsYSpace.name(), yFitIC) - return np.array(m.values), np.array( - m.errors - ) # Pass into array to store values in variable + # Pass into array to store values in variable + return np.array(m.values), np.array(m.errors) def extractData(ws, wsRes, ic): @@ -1850,6 +1847,24 @@ def minuitInitialParameters(defaultPars, sharedPars, nSpec): initPars[up + str(i)] = defaultPars[up] return initPars +def create_table_for_global_fit_parameters(wsName, m, chi2): + t = CreateEmptyTableWorkspace( + OutputWorkspace=wsName + "_gloabalfit_Parameters" + ) + t.setTitle("Global Fit Parameters") + t.addColumn(type="str", name="Name") + t.addColumn(type="float", name="Value") + t.addColumn(type="float", name="Error") + + print(f"Value of Chi2/ndof: {chi2:.2f}") + print(f"Migrad Minimum valid: {m.valid}") + print("\nResults of Global Fit:\n") + for p, v, e in zip(m.parameters, m.values, m.errors): + print(f"{p:>7s} = {v:>8.4f} \u00B1 {e:<8.4f}") + t.addRow([p, v, e]) + + t.addRow(["Cost function", chi2, 0]) + def plotGlobalFit(dataX, dataY, dataE, mObj, totCost, wsName, yFitIC): if len(dataY) > 10: @@ -1884,9 +1899,16 @@ def plotGlobalFit(dataX, dataY, dataE, mObj, totCost, wsName, yFitIC): for p, v, e in zip(signature, values, errors): leg.append(f"${p} = {v:.3f} \\pm {e:.3f}$") - ax.fill_between(x, yfit, label="\n".join(leg), alpha=0.4) + ax.plot(x, yfit, "r-", label="\n".join(leg)) + ax.plot(x, y-yfit, "b.", label=f"Residuals") + ax.grid() ax.legend() savePath = yFitIC.figSavePath / fig.canvas.manager.get_window_title() plt.savefig(savePath, bbox_inches="tight") fig.show() return + +def save_workspaces(yFitIC): + for ws_name in mtd.getObjectNames(): + if ws_name.endswith('Parameters') or ws_name.endswith('Workspace'): + SaveAscii(ws_name, str(yFitIC.figSavePath.parent / "output_files" / ws_name)) diff --git a/src/mvesuvio/main/run_routine.py b/src/mvesuvio/main/run_routine.py index eeb58584..354a9c0c 100644 --- a/src/mvesuvio/main/run_routine.py +++ b/src/mvesuvio/main/run_routine.py @@ -49,7 +49,7 @@ def setup(self): inputs_script_path = Path(handle_config.read_config_var("caching.inputs")) script_name = handle_config.get_script_name() self.experiment_path = inputs_script_path.parent / script_name - self.input_ws_path = self.experiment_path / "input_ws" + self.input_ws_path = self.experiment_path / "input_workspaces" self.input_ws_path.mkdir(parents=True, exist_ok=True) # TODO: Output paths should probably not be set like this diff --git a/tests/data/analysis/inputs/analysis_test/input_ws/analysis_test_empty_backward.nxs b/tests/data/analysis/inputs/analysis_test/input_workspaces/analysis_test_empty_backward.nxs similarity index 100% rename from tests/data/analysis/inputs/analysis_test/input_ws/analysis_test_empty_backward.nxs rename to tests/data/analysis/inputs/analysis_test/input_workspaces/analysis_test_empty_backward.nxs diff --git a/tests/data/analysis/inputs/analysis_test/input_ws/analysis_test_empty_forward.nxs b/tests/data/analysis/inputs/analysis_test/input_workspaces/analysis_test_empty_forward.nxs similarity index 100% rename from tests/data/analysis/inputs/analysis_test/input_ws/analysis_test_empty_forward.nxs rename to tests/data/analysis/inputs/analysis_test/input_workspaces/analysis_test_empty_forward.nxs diff --git a/tests/data/analysis/inputs/analysis_test/input_ws/analysis_test_raw_backward.nxs b/tests/data/analysis/inputs/analysis_test/input_workspaces/analysis_test_raw_backward.nxs similarity index 100% rename from tests/data/analysis/inputs/analysis_test/input_ws/analysis_test_raw_backward.nxs rename to tests/data/analysis/inputs/analysis_test/input_workspaces/analysis_test_raw_backward.nxs diff --git a/tests/data/analysis/inputs/analysis_test/input_ws/analysis_test_raw_forward.nxs b/tests/data/analysis/inputs/analysis_test/input_workspaces/analysis_test_raw_forward.nxs similarity index 100% rename from tests/data/analysis/inputs/analysis_test/input_ws/analysis_test_raw_forward.nxs rename to tests/data/analysis/inputs/analysis_test/input_workspaces/analysis_test_raw_forward.nxs diff --git a/tests/data/analysis/inputs/yspace_gauss_test/input_ws/yspace_gauss_test_empty_backward.nxs b/tests/data/analysis/inputs/yspace_gauss_test/input_workspaces/yspace_gauss_test_empty_backward.nxs similarity index 100% rename from tests/data/analysis/inputs/yspace_gauss_test/input_ws/yspace_gauss_test_empty_backward.nxs rename to tests/data/analysis/inputs/yspace_gauss_test/input_workspaces/yspace_gauss_test_empty_backward.nxs diff --git a/tests/data/analysis/inputs/yspace_gauss_test/input_ws/yspace_gauss_test_empty_forward.nxs b/tests/data/analysis/inputs/yspace_gauss_test/input_workspaces/yspace_gauss_test_empty_forward.nxs similarity index 100% rename from tests/data/analysis/inputs/yspace_gauss_test/input_ws/yspace_gauss_test_empty_forward.nxs rename to tests/data/analysis/inputs/yspace_gauss_test/input_workspaces/yspace_gauss_test_empty_forward.nxs diff --git a/tests/data/analysis/inputs/yspace_gauss_test/input_ws/yspace_gauss_test_raw_backward.nxs b/tests/data/analysis/inputs/yspace_gauss_test/input_workspaces/yspace_gauss_test_raw_backward.nxs similarity index 100% rename from tests/data/analysis/inputs/yspace_gauss_test/input_ws/yspace_gauss_test_raw_backward.nxs rename to tests/data/analysis/inputs/yspace_gauss_test/input_workspaces/yspace_gauss_test_raw_backward.nxs diff --git a/tests/data/analysis/inputs/yspace_gauss_test/input_ws/yspace_gauss_test_raw_forward.nxs b/tests/data/analysis/inputs/yspace_gauss_test/input_workspaces/yspace_gauss_test_raw_forward.nxs similarity index 100% rename from tests/data/analysis/inputs/yspace_gauss_test/input_ws/yspace_gauss_test_raw_forward.nxs rename to tests/data/analysis/inputs/yspace_gauss_test/input_workspaces/yspace_gauss_test_raw_forward.nxs diff --git a/tests/data/analysis/inputs/yspace_gc_test/input_ws/yspace_gc_test_empty_backward.nxs b/tests/data/analysis/inputs/yspace_gc_test/input_workspaces/yspace_gc_test_empty_backward.nxs similarity index 100% rename from tests/data/analysis/inputs/yspace_gc_test/input_ws/yspace_gc_test_empty_backward.nxs rename to tests/data/analysis/inputs/yspace_gc_test/input_workspaces/yspace_gc_test_empty_backward.nxs diff --git a/tests/data/analysis/inputs/yspace_gc_test/input_ws/yspace_gc_test_empty_forward.nxs b/tests/data/analysis/inputs/yspace_gc_test/input_workspaces/yspace_gc_test_empty_forward.nxs similarity index 100% rename from tests/data/analysis/inputs/yspace_gc_test/input_ws/yspace_gc_test_empty_forward.nxs rename to tests/data/analysis/inputs/yspace_gc_test/input_workspaces/yspace_gc_test_empty_forward.nxs diff --git a/tests/data/analysis/inputs/yspace_gc_test/input_ws/yspace_gc_test_raw_backward.nxs b/tests/data/analysis/inputs/yspace_gc_test/input_workspaces/yspace_gc_test_raw_backward.nxs similarity index 100% rename from tests/data/analysis/inputs/yspace_gc_test/input_ws/yspace_gc_test_raw_backward.nxs rename to tests/data/analysis/inputs/yspace_gc_test/input_workspaces/yspace_gc_test_raw_backward.nxs diff --git a/tests/data/analysis/inputs/yspace_gc_test/input_ws/yspace_gc_test_raw_forward.nxs b/tests/data/analysis/inputs/yspace_gc_test/input_workspaces/yspace_gc_test_raw_forward.nxs similarity index 100% rename from tests/data/analysis/inputs/yspace_gc_test/input_ws/yspace_gc_test_raw_forward.nxs rename to tests/data/analysis/inputs/yspace_gc_test/input_workspaces/yspace_gc_test_raw_forward.nxs From 67e12554e6ed89963c863aed34ff4808514d4086 Mon Sep 17 00:00:00 2001 From: GuiMacielPereira Date: Fri, 13 Dec 2024 17:49:43 +0000 Subject: [PATCH 7/7] Fix bug with rerunning y-space fit Re-running y-space fit with subtraction of fse turned on was picking up the wrong workspace that also ended in 'fse'. This correction guarantees the right workspace is chosen. --- src/mvesuvio/analysis_fitting.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/mvesuvio/analysis_fitting.py b/src/mvesuvio/analysis_fitting.py index 725e8dd2..d866764c 100644 --- a/src/mvesuvio/analysis_fitting.py +++ b/src/mvesuvio/analysis_fitting.py @@ -59,9 +59,9 @@ def find_ws_name_fse_first_mass(ic): for ws_name in mtd.getObjectNames(): if ws_name.startswith(prefix) and ws_name.endswith('fse'): name_ending = ws_name.replace(prefix, "") - match = re.search(r'\d+\.?\d*', name_ending) + match = re.search(r'_(\d+(?:\.\d+)?)_', name_ending) if match: # If float found in ws name ending - ws_masses.append(float(match.group())) + ws_masses.append(float(match.group().replace('_', ''))) ws_names_fse.append(ws_name) return ws_names_fse[np.argmin(ws_masses)]