diff --git a/EVSVesuvio/__init__.py b/EVSVesuvio/__init__.py index 61aa99df..10699baf 100644 --- a/EVSVesuvio/__init__.py +++ b/EVSVesuvio/__init__.py @@ -7,5 +7,5 @@ """ version_info = (0, 0, 0) -__version__ = '.'.join(map(str, version_info)) -__project_url__ = 'https://github.com/mantidproject/vesuvio' +__version__ = ".".join(map(str, version_info)) +__project_url__ = "https://github.com/mantidproject/vesuvio" diff --git a/EVSVesuvio/analysis_runner.py b/EVSVesuvio/analysis_runner.py index a80e1f56..51c79a93 100644 --- a/EVSVesuvio/analysis_runner.py +++ b/EVSVesuvio/analysis_runner.py @@ -7,10 +7,14 @@ def run(yes_to_all=False): - scriptName = handle_config.read_config_var('caching.experiment') - experimentsPath = Path(handle_config.read_config_var('caching.location')) / "experiments" / scriptName # Path to the repository + scriptName = handle_config.read_config_var("caching.experiment") + experimentsPath = ( + Path(handle_config.read_config_var("caching.location")) + / "experiments" + / scriptName + ) # Path to the repository inputs_path = experimentsPath / "analysis_inputs.py" - ipFilesPath = Path(handle_config.read_config_var('caching.ipfolder')) + ipFilesPath = Path(handle_config.read_config_var("caching.ipfolder")) ai = import_from_path(inputs_path, "analysis_inputs") start_time = time.time() @@ -23,10 +27,20 @@ def run(yes_to_all=False): bootIC = ai.BootstrapInitialConditions userCtr = ai.UserScriptControls - runScript(userCtr, scriptName, wsBackIC, wsFrontIC, bckwdIC, fwdIC, yFitIC, bootIC, yes_to_all) + runScript( + userCtr, + scriptName, + wsBackIC, + wsFrontIC, + bckwdIC, + fwdIC, + yFitIC, + bootIC, + yes_to_all, + ) end_time = time.time() - print("\nRunning time: ", end_time-start_time, " seconds") + print("\nRunning time: ", end_time - start_time, " seconds") def import_from_path(path, name): @@ -37,5 +51,5 @@ def import_from_path(path, name): return module -if __name__ == '__main__': +if __name__ == "__main__": run() diff --git a/EVSVesuvio/config/analysis_inputs.py b/EVSVesuvio/config/analysis_inputs.py index c1a61e7f..8690a8ca 100644 --- a/EVSVesuvio/config/analysis_inputs.py +++ b/EVSVesuvio/config/analysis_inputs.py @@ -3,39 +3,44 @@ class LoadVesuvioBackParameters: def __init__(self, ipFilesPath): - self.ipfile=ipFilesPath / "ip2019.par" + self.ipfile = ipFilesPath / "ip2019.par" - runs="43066-43076" # 77K # The numbers of the runs to be analysed - empty_runs="41876-41923" # 77K # The numbers of the empty runs to be subtracted - mode='DoubleDifference' + runs = "43066-43076" # 77K # The numbers of the runs to be analysed + empty_runs = ( + "41876-41923" # 77K # The numbers of the empty runs to be subtracted + ) + mode = "DoubleDifference" - subEmptyFromRaw = True # Flag to control wether empty ws gets subtracted from raw - scaleEmpty = 1 # None or scaling factor + subEmptyFromRaw = True # Flag to control wether empty ws gets subtracted from raw + scaleEmpty = 1 # None or scaling factor scaleRaw = 1 class LoadVesuvioFrontParameters: def __init__(self, ipFilesPath): - self.ipfile=ipFilesPath / "ip2018_3.par" + self.ipfile = ipFilesPath / "ip2018_3.par" - runs='43066-43076' # 100K # The numbers of the runs to be analysed - empty_runs='43868-43911' # 100K # The numbers of the empty runs to be subtracted - mode='SingleDifference' + runs = "43066-43076" # 100K # The numbers of the runs to be analysed + empty_runs = ( + "43868-43911" # 100K # The numbers of the empty runs to be subtracted + ) + mode = "SingleDifference" - subEmptyFromRaw = False # Flag to control wether empty ws gets subtracted from raw - scaleEmpty = 1 # None or scaling factor + subEmptyFromRaw = False # Flag to control wether empty ws gets subtracted from raw + scaleEmpty = 1 # None or scaling factor scaleRaw = 1 class GeneralInitialConditions: """Used to define initial conditions shared by both Back and Forward scattering""" + # Sample slab parameters vertical_width, horizontal_width, thickness = 0.1, 0.1, 0.001 # Expressed in meters class BackwardInitialConditions(GeneralInitialConditions): def __init__(self, ipFilesPath): - self.InstrParsPath=ipFilesPath / "ip2018_3.par" + self.InstrParsPath = ipFilesPath / "ip2018_3.par" HToMassIdxRatio = 19.0620008206 # Set to zero or None when H is not present massIdx = 0 @@ -45,21 +50,25 @@ def __init__(self, ipFilesPath): # noOfMasses = len(masses) # Intensities, NCP widths, NCP centers - initPars = np.array([ - 1, 12, 0., - 1, 12, 0., - 1, 12.5, 0. - ]) - bounds = np.array([ - [0, np.nan], [8, 16], [-3, 1], - [0, np.nan], [8, 16], [-3, 1], - [0, np.nan], [11, 14], [-3, 1] - ]) + initPars = np.array([1, 12, 0.0, 1, 12, 0.0, 1, 12.5, 0.0]) + bounds = np.array( + [ + [0, np.nan], + [8, 16], + [-3, 1], + [0, np.nan], + [8, 16], + [-3, 1], + [0, np.nan], + [11, 14], + [-3, 1], + ] + ) constraints = () - noOfMSIterations = 3 #4 - firstSpec = 3 #3 - lastSpec = 134 #134 + noOfMSIterations = 3 # 4 + firstSpec = 3 # 3 + lastSpec = 134 # 134 maskedSpecAllNo = np.array([18, 34, 42, 43, 59, 60, 62, 118, 119, 133]) @@ -68,12 +77,12 @@ def __init__(self, ipFilesPath): GammaCorrectionFlag = False # # Parameters of workspaces in input_ws - tofBinning='275.,1.,420' # Binning of ToF spectra + tofBinning = "275.,1.,420" # Binning of ToF spectra maskTOFRange = None - transmission_guess = 0.8537 # Experimental value from VesuvioTransmission + transmission_guess = 0.8537 # Experimental value from VesuvioTransmission multiple_scattering_order = 2 - number_of_events = 1.e5 + number_of_events = 1.0e5 # Original data uses histogram data instead of point data runHistData = True @@ -82,28 +91,33 @@ def __init__(self, ipFilesPath): class ForwardInitialConditions(GeneralInitialConditions): def __init__(self, ipFilesPath): - self.InstrParsPath=ipFilesPath / "ip2018_3.par" + self.InstrParsPath = ipFilesPath / "ip2018_3.par" masses = np.array([1.0079, 12, 16, 27]) # Intensities, NCP widths, NCP centers - initPars = np.array([ - 1, 4.7, 0, - 1, 12.71, 0., - 1, 8.76, 0., - 1, 13.897, 0. - ]) - bounds = np.array([ - [0, np.nan], [3, 6], [-3, 1], - [0, np.nan], [12.71, 12.71], [-3, 1], - [0, np.nan], [8.76, 8.76], [-3, 1], - [0, np.nan], [13.897, 13.897], [-3, 1] - ]) + initPars = np.array([1, 4.7, 0, 1, 12.71, 0.0, 1, 8.76, 0.0, 1, 13.897, 0.0]) + bounds = np.array( + [ + [0, np.nan], + [3, 6], + [-3, 1], + [0, np.nan], + [12.71, 12.71], + [-3, 1], + [0, np.nan], + [8.76, 8.76], + [-3, 1], + [0, np.nan], + [13.897, 13.897], + [-3, 1], + ] + ) constraints = () - noOfMSIterations = 0 #4 - firstSpec = 144 #144 - lastSpec = 182 #182 + noOfMSIterations = 0 # 4 + firstSpec = 144 # 144 + lastSpec = 182 # 182 # Boolean Flags to control script MSCorrectionFlag = True @@ -111,12 +125,12 @@ def __init__(self, ipFilesPath): maskedSpecAllNo = np.array([173, 174, 179]) - tofBinning="110,1,430" # Binning of ToF spectra + tofBinning = "110,1,430" # Binning of ToF spectra maskTOFRange = None - transmission_guess = 0.8537 # Experimental value from VesuvioTransmission + transmission_guess = 0.8537 # Experimental value from VesuvioTransmission multiple_scattering_order = 2 - number_of_events = 1.e5 + number_of_events = 1.0e5 # Original data uses histogram data instead of point data runHistData = True @@ -126,19 +140,19 @@ def __init__(self, ipFilesPath): class YSpaceFitInitialConditions: showPlots = True symmetrisationFlag = True - rebinParametersForYSpaceFit = "-25, 0.5, 25" # Needs to be symetric - fitModel = "Gaussian3D" # Options: 'SINGLE_GAUSSIAN', 'GC_C4', 'GC_C6', 'GC_C4_C6', 'DOUBLE_WELL', 'ANSIO_GAUSSIAN', 'Gaussian3D' + rebinParametersForYSpaceFit = "-25, 0.5, 25" # Needs to be symetric + fitModel = "Gaussian3D" # 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 + globalFit = True # Performs global fit with Minuit by default + nGlobalFitGroups = 4 # Number or string "ALL" + maskTypeProcedure = "NAN" # Options: 'NCP', 'NAN', None class UserScriptControls: runRoutine = True # Choose main procedure to run - procedure = "FORWARD" # Options: None, "BACKWARD", "FORWARD", "JOINT" + procedure = "FORWARD" # Options: None, "BACKWARD", "FORWARD", "JOINT" # Choose on which ws to perform the fit in y space fitInYSpace = "FORWARD" # Options: None, "BACKWARD", "FORWARD", "JOINT" diff --git a/EVSVesuvio/scripts/__init__.py b/EVSVesuvio/scripts/__init__.py index 5e5164d3..18c7ddc7 100644 --- a/EVSVesuvio/scripts/__init__.py +++ b/EVSVesuvio/scripts/__init__.py @@ -18,15 +18,33 @@ def main(): def __set_up_parser(): - parser = argparse.ArgumentParser(description="Package to analyse Vesuvio instrument data") - subparsers = parser.add_subparsers(dest='command', required=True) + parser = argparse.ArgumentParser( + description="Package to analyse Vesuvio instrument data" + ) + subparsers = parser.add_subparsers(dest="command", required=True) config_parser = subparsers.add_parser("config", help="set mvesuvio configuration") - config_parser.add_argument("--set-cache", "-c", help="set the cache directory", default="", type=str) - config_parser.add_argument("--set-ipfolder", "-p", help="set the intrument parameters directory", default="", type=str) - config_parser.add_argument("--set-experiment", "-e", help="set the current experiment", default="", type=str) + config_parser.add_argument( + "--set-cache", "-c", help="set the cache directory", default="", type=str + ) + config_parser.add_argument( + "--set-ipfolder", + "-p", + help="set the intrument parameters directory", + default="", + type=str, + ) + config_parser.add_argument( + "--set-experiment", + "-e", + help="set the current experiment", + default="", + type=str, + ) run_parser = subparsers.add_parser("run", help="run mvesuvio analysis") - run_parser.add_argument("--yes", "-y", help="Say yes to all input prompts", action='store_true') + run_parser.add_argument( + "--yes", "-y", help="Say yes to all input prompts", action="store_true" + ) return parser @@ -36,28 +54,51 @@ def __setup_config(args): ipfolder_dir = handle_config.VESUVIO_IPFOLDER_PATH if handle_config.config_set(): - cache_dir = handle_config.read_config_var('caching.location') if not args or not args.set_cache else args.set_cache - experiment = handle_config.read_config_var('caching.experiment') if not args or not args.set_experiment else args.set_experiment - ipfolder_dir = handle_config.read_config_var('caching.ipfolder') if not args or not args.set_ipfolder else args.set_ipfolder + cache_dir = ( + handle_config.read_config_var("caching.location") + if not args or not args.set_cache + else args.set_cache + ) + experiment = ( + handle_config.read_config_var("caching.experiment") + if not args or not args.set_experiment + else args.set_experiment + ) + ipfolder_dir = ( + handle_config.read_config_var("caching.ipfolder") + if not args or not args.set_ipfolder + else args.set_ipfolder + ) else: cache_dir = config_dir if not args or not args.set_cache else args.set_cache - experiment = "default" if not args or not args.set_experiment else args.set_experiment - ipfolder_dir = ipfolder_dir if not args or not args.set_ipfolder else args.set_ipfolder + experiment = ( + "default" if not args or not args.set_experiment else args.set_experiment + ) + ipfolder_dir = ( + ipfolder_dir if not args or not args.set_ipfolder else args.set_ipfolder + ) handle_config.setup_default_ipfile_dir() - handle_config.set_config_vars({'caching.location': cache_dir, - 'caching.experiment': experiment, - 'caching.ipfolder': ipfolder_dir}) + handle_config.set_config_vars( + { + "caching.location": cache_dir, + "caching.experiment": experiment, + "caching.ipfolder": ipfolder_dir, + } + ) handle_config.setup_expr_dir(cache_dir, experiment) handle_config.check_dir_exists("IP folder", ipfolder_dir) def __run_analysis(yes_to_all): - environ['MANTIDPROPERTIES'] = path.join(handle_config.VESUVIO_CONFIG_PATH, "Mantid.user.properties") + environ["MANTIDPROPERTIES"] = path.join( + handle_config.VESUVIO_CONFIG_PATH, "Mantid.user.properties" + ) from EVSVesuvio import analysis_runner + analysis_runner.run(yes_to_all) -if __name__ == '__main__': +if __name__ == "__main__": main() diff --git a/EVSVesuvio/scripts/handle_config.py b/EVSVesuvio/scripts/handle_config.py index 89c1e82f..c3cb394c 100644 --- a/EVSVesuvio/scripts/handle_config.py +++ b/EVSVesuvio/scripts/handle_config.py @@ -1,7 +1,7 @@ import os from shutil import copyfile, copytree, ignore_patterns -VESUVIO_CONFIG_PATH = os.path.join(os.path.expanduser("~"), '.mvesuvio') +VESUVIO_CONFIG_PATH = os.path.join(os.path.expanduser("~"), ".mvesuvio") VESUVIO_CONFIG_FILE = "vesuvio.user.properties" VESUVIO_INPUTS_FILE = "analysis_inputs.py" VESUVIO_PACKAGE_PATH = os.path.dirname(os.path.dirname(os.path.abspath(__file__))) @@ -12,11 +12,13 @@ def __read_config(config_file_path, throw_on_not_found=True): lines = "" try: - with open(config_file_path, 'r') as file: + with open(config_file_path, "r") as file: lines = file.readlines() except IOError: if throw_on_not_found: - raise RuntimeError(f"Could not read from vesuvio config file: {config_file_path}") + raise RuntimeError( + f"Could not read from vesuvio config file: {config_file_path}" + ) return lines @@ -29,50 +31,63 @@ def set_config_vars(var_dict): match = False for var in var_dict: if line.startswith(var): - new_line = f'{var}={var_dict[var]}' - updated_lines.append(f'{new_line}\n') + new_line = f"{var}={var_dict[var]}" + updated_lines.append(f"{new_line}\n") match = True - print(f'Setting: {new_line}') + print(f"Setting: {new_line}") break if not match: updated_lines.append(line) - with open(file_path, 'w') as file: + with open(file_path, "w") as file: file.writelines(updated_lines) def read_config_var(var, throw_on_not_found=True): - file_path = f'{VESUVIO_CONFIG_PATH}/{VESUVIO_CONFIG_FILE}' + file_path = f"{VESUVIO_CONFIG_PATH}/{VESUVIO_CONFIG_FILE}" lines = __read_config(file_path, throw_on_not_found) result = "" for line in lines: if line.startswith(var): - result = line.split("=", 2)[1].strip('\n') + result = line.split("=", 2)[1].strip("\n") break if not result and throw_on_not_found: - raise ValueError(f'{var} was not found in the vesuvio config') + raise ValueError(f"{var} was not found in the vesuvio config") return result def setup_config_dir(config_dir): - success = __mk_dir('config', config_dir) + success = __mk_dir("config", config_dir) if success: - copyfile(os.path.join(VESUVIO_PACKAGE_PATH, "config", VESUVIO_CONFIG_FILE), os.path.join(config_dir, VESUVIO_CONFIG_FILE)) - copyfile(os.path.join(VESUVIO_PACKAGE_PATH, "config", MANTID_CONFIG_FILE), os.path.join(config_dir, MANTID_CONFIG_FILE)) + copyfile( + os.path.join(VESUVIO_PACKAGE_PATH, "config", VESUVIO_CONFIG_FILE), + os.path.join(config_dir, VESUVIO_CONFIG_FILE), + ) + copyfile( + os.path.join(VESUVIO_PACKAGE_PATH, "config", MANTID_CONFIG_FILE), + os.path.join(config_dir, MANTID_CONFIG_FILE), + ) def setup_expr_dir(cache_dir, experiment): expr_path = os.path.join(cache_dir, "experiments", experiment) - success = __mk_dir('experiment', expr_path) + success = __mk_dir("experiment", expr_path) if success: - copyfile(os.path.join(VESUVIO_PACKAGE_PATH, "config", VESUVIO_INPUTS_FILE), input_file_path(cache_dir, experiment)) + copyfile( + os.path.join(VESUVIO_PACKAGE_PATH, "config", VESUVIO_INPUTS_FILE), + input_file_path(cache_dir, experiment), + ) def setup_default_ipfile_dir(): if not os.path.isdir(VESUVIO_IPFOLDER_PATH): - copytree(os.path.join(VESUVIO_PACKAGE_PATH, "config", "ip_files"), VESUVIO_IPFOLDER_PATH, ignore=ignore_patterns('__*')) + copytree( + os.path.join(VESUVIO_PACKAGE_PATH, "config", "ip_files"), + VESUVIO_IPFOLDER_PATH, + ignore=ignore_patterns("__*"), + ) def __mk_dir(type, path): @@ -82,12 +97,12 @@ def __mk_dir(type, path): os.makedirs(path) success = True except: - print(f'Unable to make {type} directory at location: {path}') + print(f"Unable to make {type} directory at location: {path}") return success def config_set(): - if(read_config_var('caching.location', False)): + if read_config_var("caching.location", False): return True else: return False diff --git a/EVSVesuvio/system_tests/test_analysis.py b/EVSVesuvio/system_tests/test_analysis.py index c49a6881..7e0dac11 100644 --- a/EVSVesuvio/system_tests/test_analysis.py +++ b/EVSVesuvio/system_tests/test_analysis.py @@ -3,7 +3,14 @@ import numpy as np import numpy.testing as nptest from pathlib import Path -from EVSVesuvio.system_tests.tests_IC import scriptName, wsBackIC, wsFrontIC, bckwdIC, fwdIC, yFitIC +from EVSVesuvio.system_tests.tests_IC import ( + scriptName, + wsBackIC, + wsFrontIC, + bckwdIC, + fwdIC, + yFitIC, +) class BootstrapInitialConditions: # Not used, but still need to pass as arg @@ -37,7 +44,9 @@ def _run(cls): bootIC = BootstrapInitialConditions userCtr = UserScriptControls - scattRes, yfitRes = runScript(userCtr, scriptName, wsBackIC, wsFrontIC, bckwdIC, fwdIC, yFitIC, bootIC) + scattRes, yfitRes = runScript( + userCtr, scriptName, wsBackIC, wsFrontIC, bckwdIC, fwdIC, yFitIC, bootIC + ) wsFinal, forwardScatteringResults = scattRes @@ -57,9 +66,13 @@ def _load_benchmark_results(cls): def displayMask(mask, rtol, string): noDiff = np.sum(mask) maskSize = mask.size - print("\nNo of different "+string+f", rtol={rtol}:\n", - noDiff, " out of ", maskSize, - f"ie {100*noDiff/maskSize:.1f} %") + print( + "\nNo of different " + string + f", rtol={rtol}:\n", + noDiff, + " out of ", + maskSize, + f"ie {100*noDiff/maskSize:.1f} %", + ) class TestFitParameters(unittest.TestCase): @@ -181,5 +194,5 @@ def test_FinalWS(self): nptest.assert_array_equal(self.optws, self.oriws) -if __name__ == '__main__': +if __name__ == "__main__": unittest.main() diff --git a/EVSVesuvio/system_tests/test_bootstrap.py b/EVSVesuvio/system_tests/test_bootstrap.py index 4c93a09f..be83a439 100644 --- a/EVSVesuvio/system_tests/test_bootstrap.py +++ b/EVSVesuvio/system_tests/test_bootstrap.py @@ -3,7 +3,14 @@ import numpy as np import numpy.testing as nptest from pathlib import Path -from EVSVesuvio.system_tests.tests_IC import scriptName, wsBackIC, wsFrontIC, bckwdIC, fwdIC, yFitIC +from EVSVesuvio.system_tests.tests_IC import ( + scriptName, + wsBackIC, + wsFrontIC, + bckwdIC, + fwdIC, + yFitIC, +) class BootstrapInitialConditions: @@ -37,7 +44,9 @@ def _run_analysis(cls): yFitIC.fitModel = "SINGLE_GAUSSIAN" yFitIC.symmetrisationFlag = True - bootRes, noneRes = runScript(userCtr, scriptName, wsBackIC, wsFrontIC, bckwdIC, fwdIC, yFitIC, bootIC) + bootRes, noneRes = runScript( + userCtr, scriptName, wsBackIC, wsFrontIC, bckwdIC, fwdIC, yFitIC, bootIC + ) # TODO: Figure out why doing the two tests simultaneously fails the testing # Probably because running bootstrap alters the initial conditions of forward scattering @@ -58,9 +67,15 @@ def _run_analysis(cls): @classmethod def _load_benchmark_results(cls): testPath = Path(__file__).absolute().parent - cls._oriJointBack = np.load(str(testPath / "stored_boot_back.npz"))["boot_samples"] - cls._oriJointFront = np.load(str(testPath / "stored_boot_front.npz"))["boot_samples"] - cls._oriJointYFit = np.load(str(testPath / "stored_boot_yfit.npz"))["boot_samples"] + cls._oriJointBack = np.load(str(testPath / "stored_boot_back.npz"))[ + "boot_samples" + ] + cls._oriJointFront = np.load(str(testPath / "stored_boot_front.npz"))[ + "boot_samples" + ] + cls._oriJointYFit = np.load(str(testPath / "stored_boot_yfit.npz"))[ + "boot_samples" + ] @classmethod def setUpClass(cls): @@ -109,5 +124,5 @@ def testYFit(self): # nptest.assert_array_almost_equal(self._bootSingleYFitSamples, self._oriYFit) # -if __name__ == '__main__': +if __name__ == "__main__": unittest.main() diff --git a/EVSVesuvio/system_tests/test_jackknife.py b/EVSVesuvio/system_tests/test_jackknife.py index 5c432bde..d8f11bf8 100644 --- a/EVSVesuvio/system_tests/test_jackknife.py +++ b/EVSVesuvio/system_tests/test_jackknife.py @@ -3,7 +3,14 @@ import numpy as np import numpy.testing as nptest from pathlib import Path -from EVSVesuvio.system_tests.tests_IC import scriptName, wsBackIC, wsFrontIC, bckwdIC, fwdIC, yFitIC +from EVSVesuvio.system_tests.tests_IC import ( + scriptName, + wsBackIC, + wsFrontIC, + bckwdIC, + fwdIC, + yFitIC, +) class BootstrapInitialConditions: @@ -13,7 +20,7 @@ class BootstrapInitialConditions: fitInYSpace = None bootstrapType = "JACKKNIFE" - nSamples = 3 # Overwritten by running Jackknife + nSamples = 3 # Overwritten by running Jackknife skipMSIterations = False runningTest = True userConfirmation = False @@ -34,7 +41,9 @@ def _run_analysis(cls): bootIC = BootstrapInitialConditions userCtr = UserScriptControls - bootRes, noneRes = runScript(userCtr, scriptName, wsBackIC, wsFrontIC, bckwdIC, fwdIC, yFitIC, bootIC) + bootRes, noneRes = runScript( + userCtr, scriptName, wsBackIC, wsFrontIC, bckwdIC, fwdIC, yFitIC, bootIC + ) cls._jackBackSamples = bootRes["bckwdScat"].bootSamples cls._jackFrontSamples = bootRes["fwdScat"].bootSamples @@ -50,8 +59,12 @@ def _run_analysis(cls): @classmethod def _load_benchmark_results(cls): testPath = Path(__file__).absolute().parent - cls._oriJointBack = np.load(str(testPath / "stored_joint_jack_back.npz"))["boot_samples"] - cls._oriJointFront = np.load(str(testPath / "stored_joint_jack_front.npz"))["boot_samples"] + cls._oriJointBack = np.load(str(testPath / "stored_joint_jack_back.npz"))[ + "boot_samples" + ] + cls._oriJointFront = np.load(str(testPath / "stored_joint_jack_front.npz"))[ + "boot_samples" + ] @classmethod def setUpClass(cls): @@ -65,5 +78,5 @@ def testFront(self): nptest.assert_array_almost_equal(self._jackFrontSamples, self._oriJointFront) -if __name__ == '__main__': +if __name__ == "__main__": unittest.main() diff --git a/EVSVesuvio/system_tests/test_yspace_fit.py b/EVSVesuvio/system_tests/test_yspace_fit.py index fe931203..b1430e47 100644 --- a/EVSVesuvio/system_tests/test_yspace_fit.py +++ b/EVSVesuvio/system_tests/test_yspace_fit.py @@ -5,12 +5,19 @@ import numpy as np import unittest import numpy.testing as nptest -from EVSVesuvio.system_tests.tests_IC import scriptName, wsBackIC, wsFrontIC, bckwdIC, fwdIC, yFitIC +from EVSVesuvio.system_tests.tests_IC import ( + scriptName, + wsBackIC, + wsFrontIC, + bckwdIC, + fwdIC, + yFitIC, +) np.set_printoptions(suppress=True, precision=8, linewidth=150) -class BootstrapInitialConditions: # Not used, but still need to pass as arg +class BootstrapInitialConditions: # Not used, but still need to pass as arg runBootstrap = False @@ -47,10 +54,16 @@ def get_current_result(cls): @classmethod def _load_workspaces(cls): AnalysisDataService.clear() - wsFinal = Load(str(cls._test_path / "wsFinal.nxs"), OutputWorkspace=scriptName + "_FORWARD_1") + wsFinal = Load( + str(cls._test_path / "wsFinal.nxs"), + OutputWorkspace=scriptName + "_FORWARD_1", + ) for i in range(len(fwdIC.masses)): fileName = "wsFinal_ncp_" + str(i) + ".nxs" - Load(str(cls._test_path / fileName), OutputWorkspace=wsFinal.name() + "_TOF_Fitted_Profile_" + str(i)) + Load( + str(cls._test_path / fileName), + OutputWorkspace=wsFinal.name() + "_TOF_Fitted_Profile_" + str(i), + ) @classmethod def _run(cls): @@ -58,7 +71,9 @@ def _run(cls): userCtr = UserScriptControls yFitIC.fitModel = "SINGLE_GAUSSIAN" - scattRes, yfitRes = runScript(userCtr, scriptName, wsBackIC, wsFrontIC, bckwdIC, fwdIC, yFitIC, bootIC) + scattRes, yfitRes = runScript( + userCtr, scriptName, wsBackIC, wsFrontIC, bckwdIC, fwdIC, yFitIC, bootIC + ) cls._currentResults = yfitRes @classmethod @@ -159,7 +174,7 @@ def setUp(self): def test_perr(self): # print("\norierr:\n", self.oriperr, "\nopterr:\n", self.optperr) - nptest.assert_array_equal( self.oriperr, self.optperr) + nptest.assert_array_equal(self.oriperr, self.optperr) if __name__ == "__main__": diff --git a/EVSVesuvio/system_tests/test_yspace_fit_GC.py b/EVSVesuvio/system_tests/test_yspace_fit_GC.py index be29ad4b..4b448be3 100644 --- a/EVSVesuvio/system_tests/test_yspace_fit_GC.py +++ b/EVSVesuvio/system_tests/test_yspace_fit_GC.py @@ -5,12 +5,19 @@ import numpy as np import unittest import numpy.testing as nptest -from EVSVesuvio.system_tests.tests_IC import scriptName, wsBackIC, wsFrontIC, bckwdIC, fwdIC, yFitIC +from EVSVesuvio.system_tests.tests_IC import ( + scriptName, + wsBackIC, + wsFrontIC, + bckwdIC, + fwdIC, + yFitIC, +) np.set_printoptions(suppress=True, precision=8, linewidth=150) -class BootstrapInitialConditions: # Not used, but still need to pass as arg +class BootstrapInitialConditions: # Not used, but still need to pass as arg runBootstrap = False @@ -47,10 +54,16 @@ def get_current_result(cls): @classmethod def _load_workspaces(cls): AnalysisDataService.clear() - wsFinal = Load(str(cls._test_path / "wsFinal.nxs"), OutputWorkspace=scriptName + "_FORWARD_1") + wsFinal = Load( + str(cls._test_path / "wsFinal.nxs"), + OutputWorkspace=scriptName + "_FORWARD_1", + ) for i in range(len(fwdIC.masses)): fileName = "wsFinal_ncp_" + str(i) + ".nxs" - Load(str(cls._test_path / fileName), OutputWorkspace=wsFinal.name() + "_TOF_Fitted_Profile_" + str(i)) + Load( + str(cls._test_path / fileName), + OutputWorkspace=wsFinal.name() + "_TOF_Fitted_Profile_" + str(i), + ) @classmethod def _run(cls): @@ -58,12 +71,16 @@ def _run(cls): userCtr = UserScriptControls yFitIC.fitModel = "GC_C4_C6" - scattRes, yfitRes = runScript(userCtr, scriptName, wsBackIC, wsFrontIC, bckwdIC, fwdIC, yFitIC, bootIC) + scattRes, yfitRes = runScript( + userCtr, scriptName, wsBackIC, wsFrontIC, bckwdIC, fwdIC, yFitIC, bootIC + ) cls._currentResults = yfitRes @classmethod def _load_benchmark_results(cls): - cls._benchmarkResults = np.load(str(cls._test_path / "stored_yspace_fit_GC.npz")) + cls._benchmarkResults = np.load( + str(cls._test_path / "stored_yspace_fit_GC.npz") + ) class TestSymSumYSpace(unittest.TestCase): @@ -159,7 +176,7 @@ def setUp(self): def test_perr(self): # print("\norierr:\n", self.oriperr, "\nopterr:\n", self.optperr) - nptest.assert_array_equal( self.oriperr, self.optperr) + nptest.assert_array_equal(self.oriperr, self.optperr) if __name__ == "__main__": diff --git a/EVSVesuvio/system_tests/tests_IC.py b/EVSVesuvio/system_tests/tests_IC.py index 9d97742a..70d4fbbb 100644 --- a/EVSVesuvio/system_tests/tests_IC.py +++ b/EVSVesuvio/system_tests/tests_IC.py @@ -2,30 +2,34 @@ from pathlib import Path from EVSVesuvio.scripts import handle_config -ipFilesPath = Path(handle_config.read_config_var('caching.ipfolder')) +ipFilesPath = Path(handle_config.read_config_var("caching.ipfolder")) ipFilePath = ipFilesPath / "ip2018_3.par" class LoadVesuvioBackParameters: - runs="43066-43076" # 77K # The numbers of the runs to be analysed - empty_runs="41876-41923" # 77K # The numbers of the empty runs to be subtracted - spectra='3-134' # Spectra to be analysed - mode='DoubleDifference' - ipfile=str(ipFilesPath / "ip2019.par") - - subEmptyFromRaw = True # Flag to control wether empty ws gets subtracted from raw + runs = "43066-43076" # 77K # The numbers of the runs to be analysed + empty_runs = ( + "41876-41923" # 77K # The numbers of the empty runs to be subtracted + ) + spectra = "3-134" # Spectra to be analysed + mode = "DoubleDifference" + ipfile = str(ipFilesPath / "ip2019.par") + + subEmptyFromRaw = True # Flag to control wether empty ws gets subtracted from raw scaleEmpty = 1 scaleRaw = 1 class LoadVesuvioFrontParameters: - runs='43066-43076' # 100K # The numbers of the runs to be analysed - empty_runs='43868-43911' # 100K # The numbers of the empty runs to be subtracted - spectra='144-182' # Spectra to be analysed - mode='SingleDifference' - ipfile=str(ipFilesPath / "ip2018_3.par") - - subEmptyFromRaw = False # Flag to control wether empty ws gets subtracted from raw + runs = "43066-43076" # 100K # The numbers of the runs to be analysed + empty_runs = ( + "43868-43911" # 100K # The numbers of the empty runs to be subtracted + ) + spectra = "144-182" # Spectra to be analysed + mode = "SingleDifference" + ipfile = str(ipFilesPath / "ip2018_3.par") + + subEmptyFromRaw = False # Flag to control wether empty ws gets subtracted from raw scaleEmpty = 1 scaleRaw = 1 @@ -33,8 +37,8 @@ class LoadVesuvioFrontParameters: class GeneralInitialConditions: """Used to define initial conditions shared by both Back and Forward scattering""" - transmission_guess = 0.8537 # Experimental value from VesuvioTransmission - multiple_scattering_order, number_of_events = 2, 1.e5 + transmission_guess = 0.8537 # Experimental value from VesuvioTransmission + multiple_scattering_order, number_of_events = 2, 1.0e5 # Sample slab parameters vertical_width, horizontal_width, thickness = 0.1, 0.1, 0.001 # Expressed in meters @@ -43,27 +47,43 @@ class BackwardInitialConditions(GeneralInitialConditions): InstrParsPath = ipFilesPath / "ip2018_3.par" HToMassIdxRatio = 19.0620008206 # Set to zero or None when H is not present - massIdx = 0 # Idx of mass to take the ratio with + massIdx = 0 # Idx of mass to take the ratio with # Masses, instrument parameters and initial fitting parameters masses = np.array([12, 16, 27]) # noOfMasses = len(masses) - initPars = np.array([ - # Intensities, NCP widths, NCP centers - 1, 12, 0., - 1, 12, 0., - 1, 12.5, 0. - ]) - bounds = np.array([ - [0, np.nan], [8, 16], [-3, 1], - [0, np.nan], [8, 16], [-3, 1], - [0, np.nan], [11, 14], [-3, 1] - ]) + initPars = np.array( + [ + # Intensities, NCP widths, NCP centers + 1, + 12, + 0.0, + 1, + 12, + 0.0, + 1, + 12.5, + 0.0, + ] + ) + bounds = np.array( + [ + [0, np.nan], + [8, 16], + [-3, 1], + [0, np.nan], + [8, 16], + [-3, 1], + [0, np.nan], + [11, 14], + [-3, 1], + ] + ) constraints = () - noOfMSIterations = 0 #4 - firstSpec = 3 #3 - lastSpec = 13 #134 + noOfMSIterations = 0 # 4 + firstSpec = 3 # 3 + lastSpec = 13 # 134 maskedSpecAllNo = np.array([18, 34, 42, 43, 59, 60, 62, 118, 119, 133]) @@ -72,7 +92,7 @@ class BackwardInitialConditions(GeneralInitialConditions): GammaCorrectionFlag = False # # Parameters of workspaces in input_ws - tofBinning='275.,1.,420' # Binning of ToF spectra + tofBinning = "275.,1.,420" # Binning of ToF spectra maskTOFRange = None # Only when running tests, to match original script @@ -81,30 +101,49 @@ class BackwardInitialConditions(GeneralInitialConditions): class ForwardInitialConditions(GeneralInitialConditions): - InstrParsPath = ipFilePath masses = np.array([1.0079, 12, 16, 27]) # noOfMasses = len(masses) - initPars = np.array([ - # Intensities, NCP widths, NCP centers - 1, 4.7, 0, - 1, 12.71, 0., - 1, 8.76, 0., - 1, 13.897, 0. - ]) - bounds = np.array([ - [0, np.nan], [3, 6], [-3, 1], - [0, np.nan], [12.71, 12.71], [-3, 1], - [0, np.nan], [8.76, 8.76], [-3, 1], - [0, np.nan], [13.897, 13.897], [-3, 1] - ]) + initPars = np.array( + [ + # Intensities, NCP widths, NCP centers + 1, + 4.7, + 0, + 1, + 12.71, + 0.0, + 1, + 8.76, + 0.0, + 1, + 13.897, + 0.0, + ] + ) + bounds = np.array( + [ + [0, np.nan], + [3, 6], + [-3, 1], + [0, np.nan], + [12.71, 12.71], + [-3, 1], + [0, np.nan], + [8.76, 8.76], + [-3, 1], + [0, np.nan], + [13.897, 13.897], + [-3, 1], + ] + ) constraints = () - noOfMSIterations = 1 #4 - firstSpec = 164 #144 - lastSpec = 175 #182 + noOfMSIterations = 1 # 4 + firstSpec = 164 # 144 + lastSpec = 175 # 182 # Boolean Flags to control script MSCorrectionFlag = True @@ -112,7 +151,7 @@ class ForwardInitialConditions(GeneralInitialConditions): maskedSpecAllNo = np.array([173, 174, 179]) - tofBinning="110,1.,430" # Binning of ToF spectra + tofBinning = "110,1.,430" # Binning of ToF spectra maskTOFRange = None # Only when running tests @@ -123,18 +162,19 @@ class ForwardInitialConditions(GeneralInitialConditions): class YSpaceFitInitialConditions: showPlots = False symmetrisationFlag = True - rebinParametersForYSpaceFit = "-20, 0.5, 20" # Needs to be symetric + rebinParametersForYSpaceFit = "-20, 0.5, 20" # Needs to be symetric fitModel = "SINGLE_GAUSSIAN" runMinos = True globalFit = None nGlobalFitGroups = 4 maskTypeProcedure = None + # userControls and bootIC defined in corresponding test scripts wsFrontIC = LoadVesuvioFrontParameters -wsBackIC = LoadVesuvioBackParameters # THIS WAS SET TO FRONT +wsBackIC = LoadVesuvioBackParameters # THIS WAS SET TO FRONT fwdIC = ForwardInitialConditions bckwdIC = BackwardInitialConditions diff --git a/EVSVesuvio/vesuvio_analysis/ICHelpers.py b/EVSVesuvio/vesuvio_analysis/ICHelpers.py index aa6d61a8..12f12005 100644 --- a/EVSVesuvio/vesuvio_analysis/ICHelpers.py +++ b/EVSVesuvio/vesuvio_analysis/ICHelpers.py @@ -2,15 +2,20 @@ from pathlib import Path from EVSVesuvio.scripts import handle_config -experimentsPath = Path(handle_config.read_config_var('caching.location')) / "experiments" +experimentsPath = ( + Path(handle_config.read_config_var("caching.location")) / "experiments" +) def completeICFromInputs(IC, scriptName, wsIC): """Assigns new methods to the initial conditions class from the inputs of that class""" - assert IC.lastSpec > IC.firstSpec, "Last spectrum needs to be bigger than first spectrum" - assert ((IC.lastSpec<135) & (IC.firstSpec<135)) | ((IC.lastSpec>=135) & (IC.firstSpec>=135) - ), "First and last spec need to be both in Back or Front scattering." + assert ( + IC.lastSpec > IC.firstSpec + ), "Last spectrum needs to be bigger than first spectrum" + assert ((IC.lastSpec < 135) & (IC.firstSpec < 135)) | ( + (IC.lastSpec >= 135) & (IC.firstSpec >= 135) + ), "First and last spec need to be both in Back or Front scattering." if IC.lastSpec <= 134: IC.modeRunning = "BACKWARD" @@ -19,12 +24,14 @@ def completeICFromInputs(IC, scriptName, wsIC): else: raise ValueError("Invalid first and last spectra input.") - IC.name = scriptName+"_"+IC.modeRunning+"_" + IC.name = scriptName + "_" + IC.modeRunning + "_" IC.masses = IC.masses.astype(float) IC.noOfMasses = len(IC.masses) - IC.maskedSpecNo = IC.maskedSpecAllNo[(IC.maskedSpecAllNo>=IC.firstSpec) & (IC.maskedSpecAllNo<=IC.lastSpec)] + IC.maskedSpecNo = IC.maskedSpecAllNo[ + (IC.maskedSpecAllNo >= IC.firstSpec) & (IC.maskedSpecAllNo <= IC.lastSpec) + ] IC.maskedDetectorIdx = IC.maskedSpecNo - IC.firstSpec # Extract some attributes from wsIC @@ -35,7 +42,7 @@ def completeICFromInputs(IC, scriptName, wsIC): # When attribute InstrParsPath is not present, set it equal to path from wsIC try: - IC.InstrParsPath # If present, leave it unaltered + IC.InstrParsPath # If present, leave it unaltered except AttributeError: IC.InstrParsPath = wsIC.ipfile @@ -57,7 +64,7 @@ def completeICFromInputs(IC, scriptName, wsIC): IC.runningPreliminary = False # Set directories for figures - figSavePath = experimentsPath / scriptName /"figures" + figSavePath = experimentsPath / scriptName / "figures" figSavePath.mkdir(exist_ok=True) IC.figSavePath = figSavePath @@ -103,7 +110,9 @@ def getRunningMode(wsIC): elif wsIC.__class__.__name__ == "LoadVesuvioFrontParameters": runningMode = "forward" else: - raise ValueError(f"Input class for loading workspace not valid: {wsIC.__class__.__name__}") + raise ValueError( + f"Input class for loading workspace not valid: {wsIC.__class__.__name__}" + ) return runningMode @@ -113,13 +122,15 @@ def setOutputDirsForSample(IC, sampleName): # Build Filename based on ic corr = "" - if IC.GammaCorrectionFlag & (IC.noOfMSIterations>0): - corr+="_GC" - if IC.MSCorrectionFlag & (IC.noOfMSIterations>0): - corr+="_MS" + if IC.GammaCorrectionFlag & (IC.noOfMSIterations > 0): + corr += "_GC" + if IC.MSCorrectionFlag & (IC.noOfMSIterations > 0): + corr += "_MS" - fileName = f"spec_{IC.firstSpec}-{IC.lastSpec}_iter_{IC.noOfMSIterations}{corr}"+".npz" - fileNameYSpace = fileName + "_ySpaceFit"+".npz" + fileName = ( + f"spec_{IC.firstSpec}-{IC.lastSpec}_iter_{IC.noOfMSIterations}{corr}" + ".npz" + ) + fileNameYSpace = fileName + "_ySpaceFit" + ".npz" IC.resultsSavePath = outputPath / fileName IC.ySpaceFitSavePath = outputPath / fileNameYSpace @@ -127,14 +138,13 @@ def setOutputDirsForSample(IC, sampleName): def wsHistoryMatchesInputs(runs, mode, ipfile, localPath): - - if not(localPath.is_file()): + if not (localPath.is_file()): 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: print(f"Filename in saved workspace did not match: {saved_runs} and {runs}") @@ -147,18 +157,19 @@ def wsHistoryMatchesInputs(runs, mode, ipfile, localPath): saved_ipfile_name = Path(metadata.getPropertyValue("InstrumentParFile")).name if saved_ipfile_name != ipfile.name: - print(f"IP files in saved workspace did not match: {saved_ipfile_name} and {ipfile.name}") + print( + f"IP files in saved workspace did not match: {saved_ipfile_name} and {ipfile.name}" + ) return False return True def saveWSFromLoadVesuvio(runs, mode, ipfile, localPath): - if "backward" in localPath.name: - spectra = '3-134' + spectra = "3-134" elif "forward" in localPath.name: - spectra = '135-198' + spectra = "135-198" else: raise ValueError(f"Invalid name to save workspace: {localPath.name}") @@ -168,8 +179,8 @@ def saveWSFromLoadVesuvio(runs, mode, ipfile, localPath): Mode=mode, InstrumentParFile=str(ipfile), OutputWorkspace=localPath.name, - LoadLogFiles=False - ) + LoadLogFiles=False, + ) SaveNexus(vesuvio_ws, str(localPath)) print(f"Workspace saved locally at: {localPath}") @@ -177,10 +188,10 @@ def saveWSFromLoadVesuvio(runs, mode, ipfile, localPath): def completeBootIC(bootIC, bckwdIC, fwdIC, yFitIC): - if not(bootIC.runBootstrap): + if not (bootIC.runBootstrap): return - try: # Assume it is not running a test if atribute is not found + try: # Assume it is not running a test if atribute is not found bootIC.runningTest except AttributeError: bootIC.runningTest = False @@ -193,13 +204,13 @@ def setBootstrapDirs(bckwdIC, fwdIC, bootIC, yFitIC): """Form bootstrap output data paths""" # Select script name and experiments path - sampleName = bckwdIC.scriptName # Name of sample currently running + sampleName = bckwdIC.scriptName # Name of sample currently running # Used to store running times required to estimate Bootstrap total run time. bootIC.runTimesPath = experimentsPath / sampleName / "running_times.txt" # Make bootstrap and jackknife data directories - if bootIC.bootstrapType=="JACKKNIFE": + if bootIC.bootstrapType == "JACKKNIFE": bootPath = experimentsPath / sampleName / "jackknife_data" else: bootPath = experimentsPath / sampleName / "bootstrap_data" @@ -214,39 +225,43 @@ def setBootstrapDirs(bckwdIC, fwdIC, bootIC, yFitIC): # Create text file for logs logFilePath = dataPath / "data_files_log.txt" - if not(logFilePath.is_file()): + if not (logFilePath.is_file()): with open(logFilePath, "w") as txtFile: txtFile.write(header_string()) - for IC in [bckwdIC, fwdIC]: # Make save paths for .npz files + for IC in [bckwdIC, fwdIC]: # Make save paths for .npz files bootName, bootNameYFit = genBootFilesName(IC, bootIC) - IC.bootSavePath = dataPath / bootName # works because modeRunning has same strings as procedure + IC.bootSavePath = ( + dataPath / bootName + ) # works because modeRunning has same strings as procedure IC.bootYFitSavePath = dataPath / bootNameYFit IC.logFilePath = logFilePath IC.bootSavePathLog = logString(bootName, IC, yFitIC, bootIC, isYFit=False) - IC.bootYFitSavePathLog = logString(bootNameYFit, IC, yFitIC, bootIC, isYFit=True) + IC.bootYFitSavePathLog = logString( + bootNameYFit, IC, yFitIC, bootIC, isYFit=True + ) return -def genBootFilesName (IC, bootIC): +def genBootFilesName(IC, bootIC): """Generates save file name for either BACKWARD or FORWARD class""" nSamples = bootIC.nSamples - if bootIC.bootstrapType=="JACKKNIFE": + if bootIC.bootstrapType == "JACKKNIFE": nSamples = 3 if bootIC.runningTest else noOfHistsFromTOFBinning(IC) # Build Filename based on ic corr = "" - if IC.MSCorrectionFlag & (IC.noOfMSIterations>0): - corr+="_MS" - if IC.GammaCorrectionFlag & (IC.noOfMSIterations>0): - corr+="_GC" + if IC.MSCorrectionFlag & (IC.noOfMSIterations > 0): + corr += "_MS" + if IC.GammaCorrectionFlag & (IC.noOfMSIterations > 0): + corr += "_GC" fileName = f"spec_{IC.firstSpec}-{IC.lastSpec}_iter_{IC.noOfMSIterations}{corr}" - bootName = fileName + f"_nsampl_{nSamples}"+".npz" - bootNameYFit = fileName + "_ySpaceFit" + f"_nsampl_{nSamples}"+".npz" + bootName = fileName + f"_nsampl_{nSamples}" + ".npz" + bootNameYFit = fileName + "_ySpaceFit" + f"_nsampl_{nSamples}" + ".npz" return bootName, bootNameYFit @@ -260,24 +275,40 @@ def header_string(): def logString(bootDataName, IC, yFitIC, bootIC, isYFit): if isYFit: - log = (bootDataName+" : "+bootIC.bootstrapType - + " | "+str(bootIC.fitInYSpace) - + " | "+str(yFitIC.symmetrisationFlag) - + " | "+yFitIC.rebinParametersForYSpaceFit - + " | "+yFitIC.fitModel - + " | "+str(yFitIC.maskTypeProcedure)) + log = ( + bootDataName + + " : " + + bootIC.bootstrapType + + " | " + + str(bootIC.fitInYSpace) + + " | " + + str(yFitIC.symmetrisationFlag) + + " | " + + yFitIC.rebinParametersForYSpaceFit + + " | " + + yFitIC.fitModel + + " | " + + str(yFitIC.maskTypeProcedure) + ) else: - log = (bootDataName+" : "+bootIC.bootstrapType - + " | "+str(bootIC.procedure) - + " | "+IC.tofBinning - + " | "+str(IC.maskTOFRange)) + log = ( + bootDataName + + " : " + + bootIC.bootstrapType + + " | " + + str(bootIC.procedure) + + " | " + + IC.tofBinning + + " | " + + str(IC.maskTOFRange) + ) return log def noOfHistsFromTOFBinning(IC): # Convert first to float and then to int because of decimal points start, spacing, end = [int(float(s)) for s in IC.tofBinning.split(",")] - return int((end-start)/spacing) - 1 # To account for last column being ignored + return int((end - start) / spacing) - 1 # To account for last column being ignored def buildFinalWSName(scriptName: str, procedure: str, IC): @@ -289,7 +320,7 @@ def buildFinalWSName(scriptName: str, procedure: str, IC): def completeYFitIC(yFitIC, sampleName): # Set directories for figures - figSavePath = experimentsPath / sampleName / "figures" + figSavePath = experimentsPath / sampleName / "figures" figSavePath.mkdir(exist_ok=True) yFitIC.figSavePath = figSavePath return diff --git a/EVSVesuvio/vesuvio_analysis/analysis_functions.py b/EVSVesuvio/vesuvio_analysis/analysis_functions.py index 53ca9d6e..b56c5000 100644 --- a/EVSVesuvio/vesuvio_analysis/analysis_functions.py +++ b/EVSVesuvio/vesuvio_analysis/analysis_functions.py @@ -12,22 +12,32 @@ def iterativeFitForDataReduction(ic): createTableInitialParameters(ic) - initialWs = loadRawAndEmptyWsFromUserPath(ic) # Do this before alternative bootstrap to extract name() + initialWs = loadRawAndEmptyWsFromUserPath( + ic + ) # Do this before alternative bootstrap to extract name() if ic.runningSampleWS: - initialWs = RenameWorkspace(InputWorkspace=ic.sampleWS, OutputWorkspace=initialWs.name()) + initialWs = RenameWorkspace( + InputWorkspace=ic.sampleWS, OutputWorkspace=initialWs.name() + ) cropedWs = cropAndMaskWorkspace(ic, initialWs) - wsToBeFitted = CloneWorkspace(InputWorkspace=cropedWs, OutputWorkspace=cropedWs.name()+"0") + wsToBeFitted = CloneWorkspace( + InputWorkspace=cropedWs, OutputWorkspace=cropedWs.name() + "0" + ) for iteration in range(ic.noOfMSIterations + 1): # Workspace from previous iteration - wsToBeFitted = mtd[ic.name+str(iteration)] + wsToBeFitted = mtd[ic.name + str(iteration)] ncpTotal = fitNcpToWorkspace(ic, wsToBeFitted) - mWidths, stdWidths, mIntRatios, stdIntRatios = extractMeans(wsToBeFitted.name(), ic) - createMeansAndStdTableWS(wsToBeFitted.name(), ic, mWidths, stdWidths, mIntRatios, stdIntRatios) + mWidths, stdWidths, mIntRatios, stdIntRatios = extractMeans( + wsToBeFitted.name(), ic + ) + createMeansAndStdTableWS( + wsToBeFitted.name(), ic, mWidths, stdWidths, mIntRatios, stdIntRatios + ) # When last iteration, skip MS and GC if iteration == ic.noOfMSIterations: @@ -42,16 +52,22 @@ def iterativeFitForDataReduction(ic): if ic.GammaCorrectionFlag: wsGC = createWorkspacesForGammaCorrection(ic, mWidths, mIntRatios, wsNCPM) - Minus(LHSWorkspace="tmpNameWs", RHSWorkspace=wsGC, OutputWorkspace="tmpNameWs") + Minus( + LHSWorkspace="tmpNameWs", RHSWorkspace=wsGC, OutputWorkspace="tmpNameWs" + ) if ic.MSCorrectionFlag: wsMS = createWorkspacesForMSCorrection(ic, mWidths, mIntRatios, wsNCPM) - Minus(LHSWorkspace="tmpNameWs", RHSWorkspace=wsMS, OutputWorkspace="tmpNameWs") + Minus( + LHSWorkspace="tmpNameWs", RHSWorkspace=wsMS, OutputWorkspace="tmpNameWs" + ) - remaskValues(ic.name, "tmpNameWS") # Masks cols in the same place as in ic.name - RenameWorkspace(InputWorkspace="tmpNameWs", OutputWorkspace=ic.name+str(iteration+1)) + remaskValues(ic.name, "tmpNameWS") # Masks cols in the same place as in ic.name + RenameWorkspace( + InputWorkspace="tmpNameWs", OutputWorkspace=ic.name + str(iteration + 1) + ) - wsFinal = mtd[ic.name+str(ic.noOfMSIterations)] + wsFinal = mtd[ic.name + str(ic.noOfMSIterations)] fittingResults = resultsObject(ic) fittingResults.save() return wsFinal, fittingResults @@ -64,12 +80,12 @@ def remaskValues(wsName, wsToMaskName): """ ws = mtd[wsName] dataX, dataY, dataE = extractWS(ws) - mask = np.all(dataY==0, axis=0) + mask = np.all(dataY == 0, axis=0) wsM = mtd[wsToMaskName] dataXM, dataYM, dataEM = extractWS(wsM) dataYM[:, mask] = 0 - if np.all(dataE==0): + if np.all(dataE == 0): dataEM = np.zeros(dataEM.shape) passDataIntoWS(dataXM, dataYM, dataEM, wsM) @@ -81,18 +97,27 @@ def createTableInitialParameters(ic): if ic.modeRunning == "BACKWARD": print(f"\nH ratio to mass with idx={ic.massIdx}: {ic.HToMassIdxRatio}\n") - meansTableWS = CreateEmptyTableWorkspace(OutputWorkspace=ic.name+"_Initial_Parameters") - meansTableWS.addColumn(type='float', name="Mass") - meansTableWS.addColumn(type='float', name="Initial Widths") - meansTableWS.addColumn(type='str', name="Bounds Widths") - meansTableWS.addColumn(type='float', name="Initial Intensities") - meansTableWS.addColumn(type='str', name="Bounds Intensities") - meansTableWS.addColumn(type='float', name="Initial Centers") - meansTableWS.addColumn(type='str', name="Bounds Centers") + meansTableWS = CreateEmptyTableWorkspace( + OutputWorkspace=ic.name + "_Initial_Parameters" + ) + meansTableWS.addColumn(type="float", name="Mass") + meansTableWS.addColumn(type="float", name="Initial Widths") + meansTableWS.addColumn(type="str", name="Bounds Widths") + meansTableWS.addColumn(type="float", name="Initial Intensities") + meansTableWS.addColumn(type="str", name="Bounds Intensities") + meansTableWS.addColumn(type="float", name="Initial Centers") + meansTableWS.addColumn(type="str", name="Bounds Centers") print("\nCreated Table with Initial Parameters:") - for m, iw, bw, ii, bi, inc, bc in zip(ic.masses.astype( - float), ic.initPars[1::3], ic.bounds[1::3], ic.initPars[0::3], ic.bounds[0::3], ic.initPars[2::3], ic.bounds[2::3]): + for m, iw, bw, ii, bi, inc, bc in zip( + ic.masses.astype(float), + ic.initPars[1::3], + ic.bounds[1::3], + ic.initPars[0::3], + ic.bounds[0::3], + ic.initPars[2::3], + ic.bounds[2::3], + ): meansTableWS.addRow([m, iw, str(bw), ii, str(bi), inc, str(bc)]) print("\nMass: ", m) print(f"{'Initial Intensity:':>20s} {ii:<8.3f} Bounds: {bi}") @@ -102,31 +127,55 @@ def createTableInitialParameters(ic): def loadRawAndEmptyWsFromUserPath(ic): + print("\nLoading local workspaces ...\n") + Load(Filename=str(ic.userWsRawPath), OutputWorkspace=ic.name + "raw") + Rebin( + InputWorkspace=ic.name + "raw", + Params=ic.tofBinning, + OutputWorkspace=ic.name + "raw", + ) - print('\nLoading local workspaces ...\n') - Load(Filename=str(ic.userWsRawPath), OutputWorkspace=ic.name+"raw") - Rebin(InputWorkspace=ic.name+'raw', Params=ic.tofBinning, - OutputWorkspace=ic.name+'raw') - - assert (isinstance(ic.scaleRaw, float)) | (isinstance(ic.scaleRaw, int)), "Scaling factor of raw ws needs to be float or int." - Scale(InputWorkspace=ic.name+'raw', OutputWorkspace=ic.name+'raw', Factor=str(ic.scaleRaw)) + assert (isinstance(ic.scaleRaw, float)) | ( + isinstance(ic.scaleRaw, int) + ), "Scaling factor of raw ws needs to be float or int." + Scale( + InputWorkspace=ic.name + "raw", + OutputWorkspace=ic.name + "raw", + Factor=str(ic.scaleRaw), + ) - SumSpectra(InputWorkspace=ic.name+'raw', OutputWorkspace=ic.name+'raw'+'_sum') - wsToBeFitted = CloneWorkspace(InputWorkspace=ic.name+'raw', OutputWorkspace=ic.name+"uncroped_unmasked") + SumSpectra(InputWorkspace=ic.name + "raw", OutputWorkspace=ic.name + "raw" + "_sum") + wsToBeFitted = CloneWorkspace( + InputWorkspace=ic.name + "raw", OutputWorkspace=ic.name + "uncroped_unmasked" + ) # if ic.mode=="DoubleDifference": if ic.subEmptyFromRaw: - Load(Filename=str(ic.userWsEmptyPath), OutputWorkspace=ic.name+"empty") - Rebin(InputWorkspace=ic.name+'empty', Params=ic.tofBinning, - OutputWorkspace=ic.name+'empty') + Load(Filename=str(ic.userWsEmptyPath), OutputWorkspace=ic.name + "empty") + Rebin( + InputWorkspace=ic.name + "empty", + Params=ic.tofBinning, + OutputWorkspace=ic.name + "empty", + ) - assert (isinstance(ic.scaleEmpty, float)) | (isinstance(ic.scaleEmpty, int)), "Scaling factor of empty ws needs to be float or int" - Scale(InputWorkspace=ic.name+'empty', OutputWorkspace=ic.name+'empty', Factor=str(ic.scaleEmpty)) + assert (isinstance(ic.scaleEmpty, float)) | ( + isinstance(ic.scaleEmpty, int) + ), "Scaling factor of empty ws needs to be float or int" + Scale( + InputWorkspace=ic.name + "empty", + OutputWorkspace=ic.name + "empty", + Factor=str(ic.scaleEmpty), + ) - SumSpectra(InputWorkspace=ic.name+'empty', OutputWorkspace=ic.name+'empty'+'_sum') + SumSpectra( + InputWorkspace=ic.name + "empty", OutputWorkspace=ic.name + "empty" + "_sum" + ) - wsToBeFitted = Minus(LHSWorkspace=ic.name+'raw', RHSWorkspace=ic.name+'empty', - OutputWorkspace=ic.name+"uncroped_unmasked") + wsToBeFitted = Minus( + LHSWorkspace=ic.name + "raw", + RHSWorkspace=ic.name + "empty", + OutputWorkspace=ic.name + "uncroped_unmasked", + ) return wsToBeFitted @@ -134,7 +183,9 @@ def cropAndMaskWorkspace(ic, ws): """Returns cloned and cropped workspace with modified name""" # Read initial Spectrum number wsFirstSpec = ws.getSpectrumNumbers()[0] - assert ic.firstSpec >= wsFirstSpec, "Can't crop workspace, firstSpec < first spectrum in workspace." + assert ( + ic.firstSpec >= wsFirstSpec + ), "Can't crop workspace, firstSpec < first spectrum in workspace." initialIdx = ic.firstSpec - wsFirstSpec lastIdx = ic.lastSpec - wsFirstSpec @@ -142,11 +193,12 @@ def cropAndMaskWorkspace(ic, ws): newWsName = ws.name().split("uncroped")[0] # Retrieve original name wsCrop = CropWorkspace( InputWorkspace=ws, - StartWorkspaceIndex=initialIdx, EndWorkspaceIndex=lastIdx, - OutputWorkspace=newWsName - ) + StartWorkspaceIndex=initialIdx, + EndWorkspaceIndex=lastIdx, + OutputWorkspace=newWsName, + ) - maskBinsWithZeros(wsCrop, ic) # Used to mask resonance peaks + maskBinsWithZeros(wsCrop, ic) # Used to mask resonance peaks MaskDetectors(Workspace=wsCrop, WorkspaceIndexList=ic.maskedDetectorIdx) return wsCrop @@ -159,13 +211,15 @@ def maskBinsWithZeros(ws, IC): Used to mask resonance peaks. """ - if IC.maskTOFRange is None: # Masked TOF bins not found, skip + if IC.maskTOFRange is None: # Masked TOF bins not found, skip return dataX, dataY, dataE = extractWS(ws) start, end = [int(s) for s in IC.maskTOFRange.split(",")] - assert start <= end, "Start value for masking needs to be smaller or equal than end." - mask = (dataX >= start) & (dataX <= end) # TOF region to mask + assert ( + start <= end + ), "Start value for masking needs to be smaller or equal than end." + mask = (dataX >= start) & (dataX <= end) # TOF region to mask dataY[mask] = 0 @@ -181,20 +235,31 @@ def fitNcpToWorkspace(IC, ws): """ dataX, dataY, dataE = extractWS(ws) - if IC.runHistData: # Converts point data from workspaces to histogram data + if IC.runHistData: # Converts point data from workspaces to histogram data dataY, dataX, dataE = histToPointData(dataY, dataX, dataE) - resolutionPars, instrPars, kinematicArrays, ySpacesForEachMass = prepareFitArgs(IC, dataX) + resolutionPars, instrPars, kinematicArrays, ySpacesForEachMass = prepareFitArgs( + IC, dataX + ) print("\nFitting NCP:\n") - arrFitPars = fitNcpToArray(IC, dataY, dataE, resolutionPars, instrPars, kinematicArrays, ySpacesForEachMass) + arrFitPars = fitNcpToArray( + IC, dataY, dataE, resolutionPars, instrPars, kinematicArrays, ySpacesForEachMass + ) createTableWSForFitPars(ws.name(), IC.noOfMasses, arrFitPars) arrBestFitPars = arrFitPars[:, 1:-2] - ncpForEachMass, ncpTotal = calculateNcpArr(IC, arrBestFitPars, resolutionPars, instrPars, kinematicArrays, ySpacesForEachMass) + ncpForEachMass, ncpTotal = calculateNcpArr( + IC, + arrBestFitPars, + resolutionPars, + instrPars, + kinematicArrays, + ySpacesForEachMass, + ) ncpSumWSs = createNcpWorkspaces(ncpForEachMass, ncpTotal, ws, IC) - wsDataSum = SumSpectra(InputWorkspace=ws, OutputWorkspace=ws.name()+"_Sum") + wsDataSum = SumSpectra(InputWorkspace=ws, OutputWorkspace=ws.name() + "_Sum") plotSumNCPFits(wsDataSum, *ncpSumWSs, IC) return ncpTotal @@ -213,11 +278,13 @@ def histToPointData(dataY, dataX, dataE): """ histWidths = dataX[:, 1:] - dataX[:, :-1] - assert np.min(histWidths) == np.max(histWidths), "Histogram widhts need to be the same length" + assert np.min(histWidths) == np.max( + histWidths + ), "Histogram widhts need to be the same length" dataYp = dataY[:, :-1] dataEp = dataE[:, :-1] - dataXp = dataX[:, :-1] + histWidths[0, 0]/2 + dataXp = dataX[:, :-1] + histWidths[0, 0] / 2 return dataYp, dataXp, dataEp @@ -227,7 +294,9 @@ def prepareFitArgs(ic, dataX): v0, E0, delta_E, delta_Q = calculateKinematicsArrays(dataX, instrPars) kinematicArrays = np.array([v0, E0, delta_E, delta_Q]) - ySpacesForEachMass = convertDataXToYSpacesForEachMass(dataX, ic.masses, delta_Q, delta_E) + ySpacesForEachMass = convertDataXToYSpacesForEachMass( + dataX, ic.masses, delta_Q, delta_E + ) kinematicArrays = reshapeArrayPerSpectrum(kinematicArrays) ySpacesForEachMass = reshapeArrayPerSpectrum(ySpacesForEachMass) @@ -247,17 +316,17 @@ def loadInstrParsFileIntoArray(InstrParsPath, firstSpec, lastSpec): def loadResolutionPars(instrPars): """Resolution of parameters to propagate into TOF resolution - Output: matrix with each parameter in each column""" + Output: matrix with each parameter in each column""" spectrums = instrPars[:, 0] L = len(spectrums) # For spec no below 135, back scattering detectors, mode is double difference # For spec no 135 or above, front scattering detectors, mode is single difference - dE1 = np.where(spectrums < 135, 88.7, 73) #meV, STD - dE1_lorz = np.where(spectrums < 135, 40.3, 24) #meV, HFHM - dTOF = np.repeat(0.37, L) #us - dTheta = np.repeat(0.016, L) #rad - dL0 = np.repeat(0.021, L) #meters - dL1 = np.repeat(0.023, L) #meters + dE1 = np.where(spectrums < 135, 88.7, 73) # meV, STD + dE1_lorz = np.where(spectrums < 135, 40.3, 24) # meV, HFHM + dTOF = np.repeat(0.37, L) # us + dTheta = np.repeat(0.016, L) # rad + dL0 = np.repeat(0.021, L) # meters + dL1 = np.repeat(0.023, L) # meters resolutionPars = np.vstack((dE1, dTOF, dTheta, dL0, dL1, dE1_lorz)).transpose() return resolutionPars @@ -267,15 +336,22 @@ def calculateKinematicsArrays(dataX, instrPars): """Kinematics quantities calculated from TOF data""" mN, Ef, en_to_vel, vf, hbar = loadConstants() - det, plick, angle, T0, L0, L1 = np.hsplit(instrPars, 6) #each is of len(dataX) - t_us = dataX - T0 #T0 is electronic delay due to instruments - v0 = vf * L0 / ( vf * t_us - L1 ) - E0 = np.square( v0 / en_to_vel ) #en_to_vel is a factor used to easily change velocity to energy and vice-versa + det, plick, angle, T0, L0, L1 = np.hsplit(instrPars, 6) # each is of len(dataX) + t_us = dataX - T0 # T0 is electronic delay due to instruments + v0 = vf * L0 / (vf * t_us - L1) + E0 = np.square( + v0 / en_to_vel + ) # en_to_vel is a factor used to easily change velocity to energy and vice-versa delta_E = E0 - Ef - delta_Q2 = 2. * mN / hbar**2 * ( E0 + Ef - 2. * np.sqrt(E0*Ef) * np.cos(angle/180.*np.pi) ) - delta_Q = np.sqrt( delta_Q2 ) - return v0, E0, delta_E, delta_Q #shape(no of spectrums, no of bins) + delta_Q2 = ( + 2.0 + * mN + / hbar**2 + * (E0 + Ef - 2.0 * np.sqrt(E0 * Ef) * np.cos(angle / 180.0 * np.pi)) + ) + delta_Q = np.sqrt(delta_Q2) + return v0, E0, delta_E, delta_Q # shape(no of spectrums, no of bins) def reshapeArrayPerSpectrum(A): @@ -297,17 +373,20 @@ def convertDataXToYSpacesForEachMass(dataX, masses, delta_Q, delta_E): mN, Ef, en_to_vel, vf, hbar = loadConstants() masses = masses.reshape(masses.size, 1, 1) - energyRecoil = np.square( hbar * delta_Q ) / 2. / masses - ySpacesForEachMass = masses / hbar**2 /delta_Q * (delta_E - energyRecoil) #y-scaling + energyRecoil = np.square(hbar * delta_Q) / 2.0 / masses + ySpacesForEachMass = ( + masses / hbar**2 / delta_Q * (delta_E - energyRecoil) + ) # y-scaling return ySpacesForEachMass -def fitNcpToArray(ic, dataY, dataE, resolutionPars, instrPars, kinematicArrays, ySpacesForEachMass): +def fitNcpToArray( + ic, dataY, dataE, resolutionPars, instrPars, kinematicArrays, ySpacesForEachMass +): """Takes dataY as a 2D array and returns the 2D array best fit parameters.""" - arrFitPars = np.zeros((len(dataY), len(ic.initPars)+3)) + arrFitPars = np.zeros((len(dataY), len(ic.initPars) + 3)) for i in range(len(dataY)): - specFitPars = fitNcpToSingleSpec( dataY[i], dataE[i], @@ -315,51 +394,58 @@ def fitNcpToArray(ic, dataY, dataE, resolutionPars, instrPars, kinematicArrays, resolutionPars[i], instrPars[i], kinematicArrays[i], - ic - ) + ic, + ) arrFitPars[i] = specFitPars - if np.all(specFitPars==0): + if np.all(specFitPars == 0): print("Skipped spectra.") else: - with np.printoptions(suppress=True, precision=4, linewidth=200, threshold=sys.maxsize): + with np.printoptions( + suppress=True, precision=4, linewidth=200, threshold=sys.maxsize + ): print(specFitPars) - assert ~np.all(arrFitPars==0), "Either Fits are all zero or assignment of fitting not working" + assert ~np.all( + arrFitPars == 0 + ), "Either Fits are all zero or assignment of fitting not working" return arrFitPars def createTableWSForFitPars(wsName, noOfMasses, arrFitPars): - tableWS = CreateEmptyTableWorkspace(OutputWorkspace=wsName+"_Best_Fit_NCP_Parameters") + tableWS = CreateEmptyTableWorkspace( + OutputWorkspace=wsName + "_Best_Fit_NCP_Parameters" + ) tableWS.setTitle("SCIPY Fit") - tableWS.addColumn(type='float', name="Spec Idx") + tableWS.addColumn(type="float", name="Spec Idx") for i in range(int(noOfMasses)): - tableWS.addColumn(type='float', name=f"Intensity {i}") - tableWS.addColumn(type='float', name=f"Width {i}") - tableWS.addColumn(type='float', name=f"Center {i}") - tableWS.addColumn(type='float', name="Norm Chi2") - tableWS.addColumn(type='float', name="No Iter") + tableWS.addColumn(type="float", name=f"Intensity {i}") + tableWS.addColumn(type="float", name=f"Width {i}") + tableWS.addColumn(type="float", name=f"Center {i}") + tableWS.addColumn(type="float", name="Norm Chi2") + tableWS.addColumn(type="float", name="No Iter") - for row in arrFitPars: # Pass array onto table ws + for row in arrFitPars: # Pass array onto table ws tableWS.addRow(row) return -def calculateNcpArr(ic, arrBestFitPars, resolutionPars, instrPars, kinematicArrays, ySpacesForEachMass): +def calculateNcpArr( + ic, arrBestFitPars, resolutionPars, instrPars, kinematicArrays, ySpacesForEachMass +): """Calculates the matrix of NCP from matrix of best fit parameters""" allNcpForEachMass = [] for i in range(len(arrBestFitPars)): - ncpForEachMass = calculateNcpRow( arrBestFitPars[i], ySpacesForEachMass[i], resolutionPars[i], instrPars[i], kinematicArrays[i], - ic - ) + ic, + ) allNcpForEachMass.append(ncpForEachMass) @@ -368,14 +454,18 @@ def calculateNcpArr(ic, arrBestFitPars, resolutionPars, instrPars, kinematicArra return allNcpForEachMass, allNcpTotal -def calculateNcpRow(initPars, ySpacesForEachMass, resolutionPars, instrPars, kinematicArrays, ic): +def calculateNcpRow( + initPars, ySpacesForEachMass, resolutionPars, instrPars, kinematicArrays, ic +): """input: all row shape - output: row shape with the ncpTotal for each mass""" + output: row shape with the ncpTotal for each mass""" - if np.all(initPars==0): + if np.all(initPars == 0): return np.zeros(ySpacesForEachMass.shape) - ncpForEachMass, ncpTotal = calculateNcpSpec(ic, initPars, ySpacesForEachMass, resolutionPars, instrPars, kinematicArrays) + ncpForEachMass, ncpTotal = calculateNcpSpec( + ic, initPars, ySpacesForEachMass, resolutionPars, instrPars, kinematicArrays + ) return ncpForEachMass @@ -386,19 +476,34 @@ def createNcpWorkspaces(ncpForEachMass, ncpTotal, ws, ic): ncpForEachMass = switchFirstTwoAxis(ncpForEachMass) # Use ws dataX to match with histogram data - dataX = ws.extractX()[:, :ncpTotal.shape[1]] # Make dataX match ncp shape automatically - assert ncpTotal.shape == dataX.shape, "DataX and DataY in ws need to be the same shape." - - ncpTotWS = createWS(dataX, ncpTotal, np.zeros(dataX.shape), ws.name()+"_TOF_Fitted_Profiles") + dataX = ws.extractX()[ + :, : ncpTotal.shape[1] + ] # Make dataX match ncp shape automatically + assert ( + ncpTotal.shape == dataX.shape + ), "DataX and DataY in ws need to be the same shape." + + ncpTotWS = createWS( + dataX, ncpTotal, np.zeros(dataX.shape), ws.name() + "_TOF_Fitted_Profiles" + ) MaskDetectors(Workspace=ncpTotWS, WorkspaceIndexList=ic.maskedDetectorIdx) - wsTotNCPSum = SumSpectra(InputWorkspace=ncpTotWS, OutputWorkspace=ncpTotWS.name()+"_Sum" ) + wsTotNCPSum = SumSpectra( + InputWorkspace=ncpTotWS, OutputWorkspace=ncpTotWS.name() + "_Sum" + ) # Individual ncp workspaces wsMNCPSum = [] for i, ncp_m in enumerate(ncpForEachMass): - ncpMWS = createWS(dataX, ncp_m, np.zeros(dataX.shape), ws.name()+"_TOF_Fitted_Profile_"+str(i)) + ncpMWS = createWS( + dataX, + ncp_m, + np.zeros(dataX.shape), + ws.name() + "_TOF_Fitted_Profile_" + str(i), + ) MaskDetectors(Workspace=ncpMWS, WorkspaceIndexList=ic.maskedDetectorIdx) - wsNCPSum = SumSpectra(InputWorkspace=ncpMWS, OutputWorkspace=ncpMWS.name()+"_Sum" ) + wsNCPSum = SumSpectra( + InputWorkspace=ncpMWS, OutputWorkspace=ncpMWS.name() + "_Sum" + ) wsMNCPSum.append(wsNCPSum) return wsTotNCPSum, wsMNCPSum @@ -410,19 +515,18 @@ def createWS(dataX, dataY, dataE, wsName): DataY=dataY.flatten(), DataE=dataE.flatten(), Nspec=len(dataY), - OutputWorkspace=wsName + OutputWorkspace=wsName, ) return ws def plotSumNCPFits(wsDataSum, wsTotNCPSum, wsMNCPSum, IC): - - if IC.runningSampleWS: # Skip saving figure if running bootstrap + if IC.runningSampleWS: # Skip saving figure if running bootstrap return lw = 2 - fig, ax = plt.subplots(subplot_kw={"projection":"mantid"}) + fig, ax = plt.subplots(subplot_kw={"projection": "mantid"}) ax.errorbar(wsDataSum, "k.", label="Spectra") ax.plot(wsTotNCPSum, "r-", label="Total NCP", linewidth=lw) @@ -434,7 +538,7 @@ def plotSumNCPFits(wsDataSum, wsTotNCPSum, wsMNCPSum, IC): ax.set_title("Sum of NCP fits") ax.legend() - fileName = wsDataSum.name()+"_NCP_Fits.pdf" + fileName = wsDataSum.name() + "_NCP_Fits.pdf" savePath = IC.figSavePath / fileName plt.savefig(savePath, bbox_inches="tight") plt.close(fig) @@ -451,30 +555,47 @@ def switchFirstTwoAxis(A): def extractMeans(wsName, IC): """Extract widths and intensities from tableWorkspace""" - fitParsTable = mtd[wsName+"_Best_Fit_NCP_Parameters"] + fitParsTable = mtd[wsName + "_Best_Fit_NCP_Parameters"] widths = np.zeros((IC.noOfMasses, fitParsTable.rowCount())) intensities = np.zeros(widths.shape) for i in range(IC.noOfMasses): widths[i] = fitParsTable.column(f"Width {i}") intensities[i] = fitParsTable.column(f"Intensity {i}") - meanWidths, stdWidths, meanIntensityRatios, stdIntensityRatios = calculateMeansAndStds(widths, intensities, IC) + ( + meanWidths, + stdWidths, + meanIntensityRatios, + stdIntensityRatios, + ) = calculateMeansAndStds(widths, intensities, IC) - assert len(widths) == IC.noOfMasses, "Widths and intensities must be in shape (noOfMasses, noOfSpec)" + assert ( + len(widths) == IC.noOfMasses + ), "Widths and intensities must be in shape (noOfMasses, noOfSpec)" return meanWidths, stdWidths, meanIntensityRatios, stdIntensityRatios -def createMeansAndStdTableWS(wsName, IC, meanWidths, stdWidths, meanIntensityRatios, stdIntensityRatios): - meansTableWS = CreateEmptyTableWorkspace(OutputWorkspace=wsName+"_Mean_Widths_And_Intensities") - meansTableWS.addColumn(type='float', name="Mass") - meansTableWS.addColumn(type='float', name="Mean Widths") - meansTableWS.addColumn(type='float', name="Std Widths") - meansTableWS.addColumn(type='float', name="Mean Intensities") - meansTableWS.addColumn(type='float', name="Std Intensities") +def createMeansAndStdTableWS( + wsName, IC, meanWidths, stdWidths, meanIntensityRatios, stdIntensityRatios +): + meansTableWS = CreateEmptyTableWorkspace( + OutputWorkspace=wsName + "_Mean_Widths_And_Intensities" + ) + meansTableWS.addColumn(type="float", name="Mass") + meansTableWS.addColumn(type="float", name="Mean Widths") + meansTableWS.addColumn(type="float", name="Std Widths") + meansTableWS.addColumn(type="float", name="Mean Intensities") + meansTableWS.addColumn(type="float", name="Std Intensities") print("\nCreated Table with means and std:") print("\nMass Mean \u00B1 Std Widths Mean \u00B1 Std Intensities\n") - for m, mw, stdw, mi, stdi in zip(IC.masses.astype(float), meanWidths, stdWidths, meanIntensityRatios, stdIntensityRatios): + for m, mw, stdw, mi, stdi in zip( + IC.masses.astype(float), + meanWidths, + stdWidths, + meanIntensityRatios, + stdIntensityRatios, + ): meansTableWS.addRow([m, mw, stdw, mi, stdi]) print(f"{m:5.2f} {mw:10.5f} \u00B1 {stdw:7.5f} {mi:10.5f} \u00B1 {stdi:7.5f}") print("\n") @@ -482,8 +603,9 @@ def createMeansAndStdTableWS(wsName, IC, meanWidths, stdWidths, meanIntensityRat def calculateMeansAndStds(widthsIn, intensitiesIn, IC): - - betterWidths, betterIntensities = filterWidthsAndIntensities(widthsIn, intensitiesIn, IC) + betterWidths, betterIntensities = filterWidthsAndIntensities( + widthsIn, intensitiesIn, IC + ) meanWidths = np.nanmean(betterWidths, axis=1) stdWidths = np.nanstd(betterWidths, axis=1) @@ -497,10 +619,12 @@ def calculateMeansAndStds(widthsIn, intensitiesIn, IC): def filterWidthsAndIntensities(widthsIn, intensitiesIn, IC): """Puts nans in places to be ignored""" - widths = widthsIn.copy() # Copy to avoid accidental changes in arrays + widths = widthsIn.copy() # Copy to avoid accidental changes in arrays intensities = intensitiesIn.copy() - zeroSpecs = np.all(widths==0, axis=0) # Catches all failed fits, not just masked spectra + zeroSpecs = np.all( + widths == 0, axis=0 + ) # Catches all failed fits, not just masked spectra widths[:, zeroSpecs] = np.nan intensities[:, zeroSpecs] = np.nan @@ -514,41 +638,58 @@ def filterWidthsAndIntensities(widthsIn, intensitiesIn, IC): betterWidths = np.where(filterMask, np.nan, widths) maskedIntensities = np.where(filterMask, np.nan, intensities) - betterIntensities = maskedIntensities / np.sum(maskedIntensities, axis=0) # Not nansum() + betterIntensities = maskedIntensities / np.sum( + maskedIntensities, axis=0 + ) # Not nansum() # When trying to estimate HToMassIdxRatio and normalization fails, skip normalization if np.all(np.isnan(betterIntensities)) & IC.runningPreliminary: - assert IC.noOfMSIterations == 0, "Calculation of mean intensities failed, cannot proceed with MS correction." \ - "Try to run again with noOfMSIterations=0." + assert IC.noOfMSIterations == 0, ( + "Calculation of mean intensities failed, cannot proceed with MS correction." + "Try to run again with noOfMSIterations=0." + ) betterIntensities = maskedIntensities else: pass - assert np.all(meanWidths!=np.nan), "At least one mean of widths is nan!" + assert np.all(meanWidths != np.nan), "At least one mean of widths is nan!" assert np.sum(filterMask) >= 1, "No widths survive filtering condition" - assert not(np.all(np.isnan(betterWidths))), "All filtered widths are nan" - assert not(np.all(np.isnan(betterIntensities))), "All filtered intensities are nan" - assert np.nanmax(betterWidths) != np.nanmin(betterWidths), f"All fitered widths have the same value: {np.nanmin(betterWidths)}" + assert not (np.all(np.isnan(betterWidths))), "All filtered widths are nan" + assert not (np.all(np.isnan(betterIntensities))), "All filtered intensities are nan" + assert np.nanmax(betterWidths) != np.nanmin( + betterWidths + ), f"All fitered widths have the same value: {np.nanmin(betterWidths)}" assert np.nanmax(betterIntensities) != np.nanmin( - betterIntensities), f"All fitered widths have the same value: {np.nanmin(betterIntensities)}" + betterIntensities + ), f"All fitered widths have the same value: {np.nanmin(betterIntensities)}" return betterWidths, betterIntensities -def fitNcpToSingleSpec(dataY, dataE, ySpacesForEachMass, resolutionPars, instrPars, kinematicArrays, ic): +def fitNcpToSingleSpec( + dataY, dataE, ySpacesForEachMass, resolutionPars, instrPars, kinematicArrays, ic +): """Fits the NCP and returns the best fit parameters for one spectrum""" - if np.all(dataY == 0) : - return np.zeros(len(ic.initPars)+3) + if np.all(dataY == 0): + return np.zeros(len(ic.initPars) + 3) result = optimize.minimize( errorFunction, ic.initPars, - args=(dataY, dataE, ySpacesForEachMass, resolutionPars, instrPars, kinematicArrays, ic), - method='SLSQP', - bounds = ic.bounds, - constraints=ic.constraints - ) + args=( + dataY, + dataE, + ySpacesForEachMass, + resolutionPars, + instrPars, + kinematicArrays, + ic, + ), + method="SLSQP", + bounds=ic.bounds, + constraints=ic.constraints, + ) fitPars = result["x"] @@ -557,48 +698,66 @@ def fitNcpToSingleSpec(dataY, dataE, ySpacesForEachMass, resolutionPars, instrPa return np.append(specFitPars, [result["fun"] / noDegreesOfFreedom, result["nit"]]) -def errorFunction(pars, dataY, dataE, ySpacesForEachMass, resolutionPars, instrPars, kinematicArrays, ic): +def errorFunction( + pars, + dataY, + dataE, + ySpacesForEachMass, + resolutionPars, + instrPars, + kinematicArrays, + ic, +): """Error function to be minimized, operates in TOF space""" - ncpForEachMass, ncpTotal = calculateNcpSpec(ic, pars, ySpacesForEachMass, resolutionPars, instrPars, kinematicArrays) + ncpForEachMass, ncpTotal = calculateNcpSpec( + ic, pars, ySpacesForEachMass, resolutionPars, instrPars, kinematicArrays + ) # Ignore any masked values from Jackknife or masked tof range - zerosMask = (dataY==0) + zerosMask = dataY == 0 ncpTotal = ncpTotal[~zerosMask] dataYf = dataY[~zerosMask] dataEf = dataE[~zerosMask] - if np.all(dataE==0): # When errors not present - return np.sum((ncpTotal - dataYf)**2) + if np.all(dataE == 0): # When errors not present + return np.sum((ncpTotal - dataYf) ** 2) - return np.sum((ncpTotal - dataYf)**2 / dataEf**2) + return np.sum((ncpTotal - dataYf) ** 2 / dataEf**2) -def calculateNcpSpec(ic, pars, ySpacesForEachMass, resolutionPars, instrPars, kinematicArrays): +def calculateNcpSpec( + ic, pars, ySpacesForEachMass, resolutionPars, instrPars, kinematicArrays +): """Creates a synthetic C(t) to be fitted to TOF values of a single spectrum, from J(y) and resolution functions - Shapes: datax (1, n), ySpacesForEachMass (4, n), res (4, 2), deltaQ (1, n), E0 (1,n), - where n is no of bins""" + Shapes: datax (1, n), ySpacesForEachMass (4, n), res (4, 2), deltaQ (1, n), E0 (1,n), + where n is no of bins""" masses, intensities, widths, centers = prepareArraysFromPars(ic, pars) v0, E0, deltaE, deltaQ = kinematicArrays gaussRes, lorzRes = caculateResolutionForEachMass( masses, ySpacesForEachMass, centers, resolutionPars, instrPars, kinematicArrays - ) + ) totalGaussWidth = np.sqrt(widths**2 + gaussRes**2) JOfY = pseudoVoigt(ySpacesForEachMass - centers, totalGaussWidth, lorzRes, ic) - FSE = - numericalThirdDerivative(ySpacesForEachMass, JOfY) * widths**4 / deltaQ * 0.72 + FSE = ( + -numericalThirdDerivative(ySpacesForEachMass, JOfY) + * widths**4 + / deltaQ + * 0.72 + ) - ncpForEachMass = intensities * (JOfY + FSE) * E0 * E0**(-0.92) * masses / deltaQ + ncpForEachMass = intensities * (JOfY + FSE) * E0 * E0 ** (-0.92) * masses / deltaQ ncpTotal = np.sum(ncpForEachMass, axis=0) return ncpForEachMass, ncpTotal def prepareArraysFromPars(ic, initPars): """Extracts the intensities, widths and centers from the fitting parameters - Reshapes all of the arrays to collumns, for the calculation of the ncp,""" + Reshapes all of the arrays to collumns, for the calculation of the ncp,""" masses = ic.masses[:, np.newaxis] intensities = initPars[::3].reshape(masses.shape) @@ -607,14 +766,22 @@ def prepareArraysFromPars(ic, initPars): return masses, intensities, widths, centers -def caculateResolutionForEachMass(masses, ySpacesForEachMass, centers, resolutionPars, instrPars, kinematicArrays): +def caculateResolutionForEachMass( + masses, ySpacesForEachMass, centers, resolutionPars, instrPars, kinematicArrays +): """Calculates the gaussian and lorentzian resolution output: two column vectors, each row corresponds to each mass""" - v0, E0, delta_E, delta_Q = kinematicsAtYCenters(ySpacesForEachMass, centers, kinematicArrays) + v0, E0, delta_E, delta_Q = kinematicsAtYCenters( + ySpacesForEachMass, centers, kinematicArrays + ) - gaussianResWidth = calcGaussianResolution(masses, v0, E0, delta_E, delta_Q, resolutionPars, instrPars) - lorentzianResWidth = calcLorentzianResolution(masses, v0, E0, delta_E, delta_Q, resolutionPars, instrPars) + gaussianResWidth = calcGaussianResolution( + masses, v0, E0, delta_E, delta_Q, resolutionPars, instrPars + ) + lorentzianResWidth = calcLorentzianResolution( + masses, v0, E0, delta_E, delta_Q, resolutionPars, instrPars + ) return gaussianResWidth, lorentzianResWidth @@ -643,40 +810,61 @@ def kinematicsAtYCenters(ySpacesForEachMass, centers, kinematicArrays): def calcGaussianResolution(masses, v0, E0, delta_E, delta_Q, resolutionPars, instrPars): # Currently the function that takes the most time in the fitting - assert masses.shape == (masses.size, 1), f"masses.shape: {masses.shape}. The shape of the masses array needs to be a collumn!" + assert masses.shape == ( + masses.size, + 1, + ), f"masses.shape: {masses.shape}. The shape of the masses array needs to be a collumn!" det, plick, angle, T0, L0, L1 = instrPars dE1, dTOF, dTheta, dL0, dL1, dE1_lorz = resolutionPars mN, Ef, en_to_vel, vf, hbar = loadConstants() - angle = angle * np.pi/180 + angle = angle * np.pi / 180 - dWdE1 = 1. + (E0 / Ef)**1.5 * (L1 / L0) - dWdTOF = 2. * E0 * v0 / L0 - dWdL1 = 2. * E0**1.5 / Ef**0.5 / L0 - dWdL0 = 2. * E0 / L0 + dWdE1 = 1.0 + (E0 / Ef) ** 1.5 * (L1 / L0) + dWdTOF = 2.0 * E0 * v0 / L0 + dWdL1 = 2.0 * E0**1.5 / Ef**0.5 / L0 + dWdL0 = 2.0 * E0 / L0 - dW2 = dWdE1**2*dE1**2 + dWdTOF**2*dTOF**2 + dWdL1**2*dL1**2 + dWdL0**2*dL0**2 + dW2 = ( + dWdE1**2 * dE1**2 + + dWdTOF**2 * dTOF**2 + + dWdL1**2 * dL1**2 + + dWdL0**2 * dL0**2 + ) # conversion from meV^2 to A^-2, dydW = (M/q)^2 - dW2 *= (masses / hbar**2 / delta_Q)**2 + dW2 *= (masses / hbar**2 / delta_Q) ** 2 - dQdE1 = 1. - (E0 / Ef)**1.5 * L1/L0 - np.cos(angle) * ((E0 / Ef)**0.5 - L1/L0 * E0/Ef) - dQdTOF = 2.*E0 * v0/L0 - dQdL1 = 2.*E0**1.5 / L0 / Ef**0.5 - dQdL0 = 2.*E0 / L0 - dQdTheta = 2. * np.sqrt(E0 * Ef) * np.sin(angle) - - dQ2 = dQdE1**2*dE1**2 + (dQdTOF**2*dTOF**2 + dQdL1**2*dL1**2 + dQdL0 - ** 2*dL0**2)*np.abs(Ef/E0*np.cos(angle)-1) + dQdTheta**2*dTheta**2 - dQ2 *= (mN / hbar**2 / delta_Q)**2 + dQdE1 = ( + 1.0 + - (E0 / Ef) ** 1.5 * L1 / L0 + - np.cos(angle) * ((E0 / Ef) ** 0.5 - L1 / L0 * E0 / Ef) + ) + dQdTOF = 2.0 * E0 * v0 / L0 + dQdL1 = 2.0 * E0**1.5 / L0 / Ef**0.5 + dQdL0 = 2.0 * E0 / L0 + dQdTheta = 2.0 * np.sqrt(E0 * Ef) * np.sin(angle) + + dQ2 = ( + dQdE1**2 * dE1**2 + + (dQdTOF**2 * dTOF**2 + dQdL1**2 * dL1**2 + dQdL0**2 * dL0**2) + * np.abs(Ef / E0 * np.cos(angle) - 1) + + dQdTheta**2 * dTheta**2 + ) + dQ2 *= (mN / hbar**2 / delta_Q) ** 2 # in A-1 #same as dy^2 = (dy/dw)^2*dw^2 + (dy/dq)^2*dq^2 gaussianResWidth = np.sqrt(dW2 + dQ2) return gaussianResWidth -def calcLorentzianResolution(masses, v0, E0, delta_E, delta_Q, resolutionPars, instrPars): - assert masses.shape == (masses.size, 1), "The shape of the masses array needs to be a collumn!" +def calcLorentzianResolution( + masses, v0, E0, delta_E, delta_Q, resolutionPars, instrPars +): + assert masses.shape == ( + masses.size, + 1, + ), "The shape of the masses array needs to be a collumn!" det, plick, angle, T0, L0, L1 = instrPars dE1, dTOF, dTheta, dL0, dL1, dE1_lorz = resolutionPars @@ -684,24 +872,27 @@ def calcLorentzianResolution(masses, v0, E0, delta_E, delta_Q, resolutionPars, i angle = angle * np.pi / 180 - dWdE1_lor = (1. + (E0/Ef)**1.5 * (L1/L0))**2 + dWdE1_lor = (1.0 + (E0 / Ef) ** 1.5 * (L1 / L0)) ** 2 # conversion from meV^2 to A^-2 - dWdE1_lor *= (masses / hbar**2 / delta_Q)**2 + dWdE1_lor *= (masses / hbar**2 / delta_Q) ** 2 - dQdE1_lor = (1. - (E0/Ef)**1.5 * L1/L0 - np.cos(angle) - * ((E0/Ef)**0.5 + L1/L0 * E0/Ef))**2 - dQdE1_lor *= (mN / hbar**2 / delta_Q)**2 + dQdE1_lor = ( + 1.0 + - (E0 / Ef) ** 1.5 * L1 / L0 + - np.cos(angle) * ((E0 / Ef) ** 0.5 + L1 / L0 * E0 / Ef) + ) ** 2 + dQdE1_lor *= (mN / hbar**2 / delta_Q) ** 2 - lorentzianResWidth = np.sqrt(dWdE1_lor + dQdE1_lor) * dE1_lorz # in A-1 + lorentzianResWidth = np.sqrt(dWdE1_lor + dQdE1_lor) * dE1_lorz # in A-1 return lorentzianResWidth def loadConstants(): """Output: the mass of the neutron, final energy of neutrons (selected by gold foil), factor to change energies into velocities, final velocity of neutron and hbar""" - mN=1.008 #a.m.u. - Ef=4906. # meV - en_to_vel = 4.3737 * 1.e-4 + mN = 1.008 # a.m.u. + Ef = 4906.0 # meV + en_to_vel = 4.3737 * 1.0e-4 vf = np.sqrt(Ef) * en_to_vel # m/us hbar = 2.0445 return mN, Ef, en_to_vel, vf, hbar @@ -709,36 +900,38 @@ def loadConstants(): def pseudoVoigt(x, sigma, gamma, IC): """Convolution between Gaussian with std sigma and Lorentzian with HWHM gamma""" - fg, fl = 2.*sigma*np.sqrt(2.*np.log(2.)), 2.*gamma - f = 0.5346 * fl + np.sqrt(0.2166*fl**2 + fg**2) - eta = 1.36603 * fl/f - 0.47719 * (fl/f)**2 + 0.11116 * (fl/f)**3 - sigma_v, gamma_v = f/(2.*np.sqrt(2.*np.log(2.))), f / 2. - pseudo_voigt = eta * lorentizian(x, gamma_v) + (1.-eta) * gaussian(x, sigma_v) - - norm = np.abs(np.trapz(pseudo_voigt, x, axis=1))[:, np.newaxis] if IC.normVoigt else 1 + fg, fl = 2.0 * sigma * np.sqrt(2.0 * np.log(2.0)), 2.0 * gamma + f = 0.5346 * fl + np.sqrt(0.2166 * fl**2 + fg**2) + eta = 1.36603 * fl / f - 0.47719 * (fl / f) ** 2 + 0.11116 * (fl / f) ** 3 + sigma_v, gamma_v = f / (2.0 * np.sqrt(2.0 * np.log(2.0))), f / 2.0 + pseudo_voigt = eta * lorentizian(x, gamma_v) + (1.0 - eta) * gaussian(x, sigma_v) + + norm = ( + np.abs(np.trapz(pseudo_voigt, x, axis=1))[:, np.newaxis] if IC.normVoigt else 1 + ) return pseudo_voigt / norm def gaussian(x, sigma): """Gaussian function centered at zero""" - gaussian = np.exp(-x**2/2/sigma**2) - gaussian /= np.sqrt(2.*np.pi)*sigma + gaussian = np.exp(-(x**2) / 2 / sigma**2) + gaussian /= np.sqrt(2.0 * np.pi) * sigma return gaussian def lorentizian(x, gamma): """Lorentzian centered at zero""" - lorentzian = gamma/np.pi / (x**2 + gamma**2) + lorentzian = gamma / np.pi / (x**2 + gamma**2) return lorentzian def numericalThirdDerivative(x, fun): - k6 = (- fun[:, 12:] + fun[:, :-12]) * 1 - k5 = (+ fun[:, 11:-1] - fun[:, 1:-11]) * 24 - k4 = (- fun[:, 10:-2] + fun[:, 2:-10]) * 192 - k3 = (+ fun[:, 9:-3] - fun[:, 3:-9]) * 488 - k2 = (+ fun[:, 8:-4] - fun[:, 4:-8]) * 387 - k1 = (- fun[:, 7:-5] + fun[:, 5:-7]) * 1584 + k6 = (-fun[:, 12:] + fun[:, :-12]) * 1 + k5 = (+fun[:, 11:-1] - fun[:, 1:-11]) * 24 + k4 = (-fun[:, 10:-2] + fun[:, 2:-10]) * 192 + k3 = (+fun[:, 9:-3] - fun[:, 3:-9]) * 488 + k2 = (+fun[:, 8:-4] - fun[:, 4:-8]) * 387 + k1 = (-fun[:, 7:-5] + fun[:, 5:-7]) * 1584 dev = k1 + k2 + k3 + k4 + k5 + k6 dev /= np.power(x[:, 7:-5] - x[:, 6:-6], 3) @@ -753,24 +946,38 @@ def numericalThirdDerivative(x, fun): def createWorkspacesForMSCorrection(ic, meanWidths, meanIntensityRatios, wsNCPM): """Creates _MulScattering and _TotScattering workspaces used for the MS correction""" - createSlabGeometry(ic, wsNCPM) # Sample properties for MS correction + createSlabGeometry(ic, wsNCPM) # Sample properties for MS correction - sampleProperties = calcMSCorrectionSampleProperties(ic, meanWidths, meanIntensityRatios) - print("\nThe sample properties for Multiple Scattering correction are:\n\n", - sampleProperties, "\n") + sampleProperties = calcMSCorrectionSampleProperties( + ic, meanWidths, meanIntensityRatios + ) + print( + "\nThe sample properties for Multiple Scattering correction are:\n\n", + sampleProperties, + "\n", + ) return createMulScatWorkspaces(ic, wsNCPM, sampleProperties) def createSlabGeometry(ic, wsNCPM): - half_height, half_width, half_thick = 0.5*ic.vertical_width, 0.5*ic.horizontal_width, 0.5*ic.thickness - xml_str = \ - " " \ - + " " % (half_width, -half_height, half_thick) \ - + " " % (half_width, half_height, half_thick) \ - + " " % (half_width, -half_height, -half_thick) \ - + " " % (-half_width, -half_height, half_thick) \ + half_height, half_width, half_thick = ( + 0.5 * ic.vertical_width, + 0.5 * ic.horizontal_width, + 0.5 * ic.thickness, + ) + xml_str = ( + ' ' + + ' ' + % (half_width, -half_height, half_thick) + + ' ' + % (half_width, half_height, half_thick) + + ' ' + % (half_width, -half_height, -half_thick) + + ' ' + % (-half_width, -half_height, half_thick) + "" + ) CreateSampleShape(wsNCPM, xml_str) @@ -779,8 +986,8 @@ def calcMSCorrectionSampleProperties(ic, meanWidths, meanIntensityRatios): masses = ic.masses.flatten() # If Backsscattering mode and H is present in the sample, add H to MS properties - if (ic.modeRunning == "BACKWARD"): - if (ic.HToMassIdxRatio is not None): # If H is present, ratio is a number + if ic.modeRunning == "BACKWARD": + if ic.HToMassIdxRatio is not None: # If H is present, ratio is a number masses = np.append(masses, 1.0079) meanWidths = np.append(meanWidths, 5.0) @@ -788,7 +995,7 @@ def calcMSCorrectionSampleProperties(ic, meanWidths, meanIntensityRatios): meanIntensityRatios = np.append(meanIntensityRatios, HIntensity) meanIntensityRatios /= np.sum(meanIntensityRatios) - MSProperties = np.zeros(3*len(masses)) + MSProperties = np.zeros(3 * len(masses)) MSProperties[::3] = masses MSProperties[1::3] = meanIntensityRatios MSProperties[2::3] = meanWidths @@ -807,8 +1014,11 @@ def createMulScatWorkspaces(ic, ws, sampleProperties): MS_amplitudes = sampleProperties[1::3] dens, trans = VesuvioThickness( - Masses=MS_masses, Amplitudes=MS_amplitudes, TransmissionGuess=ic.transmission_guess, Thickness=0.1 - ) + Masses=MS_masses, + Amplitudes=MS_amplitudes, + TransmissionGuess=ic.transmission_guess, + Thickness=0.1, + ) _TotScattering, _MulScattering = VesuvioCalculateMS( ws, @@ -817,22 +1027,30 @@ def createMulScatWorkspaces(ic, ws, sampleProperties): AtomicProperties=sampleProperties, BeamRadius=2.5, NumScatters=ic.multiple_scattering_order, - NumEventsPerRun=int(ic.number_of_events) - ) + NumEventsPerRun=int(ic.number_of_events), + ) data_normalisation = Integration(ws) simulation_normalisation = Integration("_TotScattering") for workspace in ("_MulScattering", "_TotScattering"): - Divide(LHSWorkspace=workspace, RHSWorkspace=simulation_normalisation, OutputWorkspace=workspace) - Multiply(LHSWorkspace=workspace, RHSWorkspace=data_normalisation, OutputWorkspace=workspace) - RenameWorkspace(InputWorkspace=workspace, OutputWorkspace=ws.name()+workspace) - SumSpectra(ws.name()+workspace, OutputWorkspace=ws.name()+workspace+"_Sum") - - DeleteWorkspaces( - [data_normalisation, simulation_normalisation, trans, dens] + Divide( + LHSWorkspace=workspace, + RHSWorkspace=simulation_normalisation, + OutputWorkspace=workspace, + ) + Multiply( + LHSWorkspace=workspace, + RHSWorkspace=data_normalisation, + OutputWorkspace=workspace, ) + RenameWorkspace(InputWorkspace=workspace, OutputWorkspace=ws.name() + workspace) + SumSpectra( + ws.name() + workspace, OutputWorkspace=ws.name() + workspace + "_Sum" + ) + + DeleteWorkspaces([data_normalisation, simulation_normalisation, trans, dens]) # The only remaining workspaces are the _MulScattering and _TotScattering - return mtd[ws.name()+"_MulScattering"] + return mtd[ws.name() + "_MulScattering"] def createWorkspacesForGammaCorrection(ic, meanWidths, meanIntensityRatios, wsNCPM): @@ -849,23 +1067,41 @@ def createWorkspacesForGammaCorrection(ic, meanWidths, meanIntensityRatios, wsNC ws = CloneWorkspace(InputWorkspace=inputWS, OutputWorkspace="tmpGC") for spec in range(ws.getNumberHistograms()): - background, corrected = VesuvioCalculateGammaBackground(InputWorkspace=inputWS, ComptonFunction=profiles, WorkspaceIndexList=spec) - ws.dataY(spec)[:], ws.dataE(spec)[:] = background.dataY(0)[:], background.dataE(0)[:] + background, corrected = VesuvioCalculateGammaBackground( + InputWorkspace=inputWS, ComptonFunction=profiles, WorkspaceIndexList=spec + ) + ws.dataY(spec)[:], ws.dataE(spec)[:] = ( + background.dataY(0)[:], + background.dataE(0)[:], + ) DeleteWorkspace(background) DeleteWorkspace(corrected) - RenameWorkspace(InputWorkspace="tmpGC", OutputWorkspace=inputWS+"_Gamma_Background") + RenameWorkspace( + InputWorkspace="tmpGC", OutputWorkspace=inputWS + "_Gamma_Background" + ) - Scale(InputWorkspace = inputWS+"_Gamma_Background", OutputWorkspace = inputWS+"_Gamma_Background", Factor=0.9, Operation="Multiply") - return mtd[inputWS+"_Gamma_Background"] + Scale( + InputWorkspace=inputWS + "_Gamma_Background", + OutputWorkspace=inputWS + "_Gamma_Background", + Factor=0.9, + Operation="Multiply", + ) + return mtd[inputWS + "_Gamma_Background"] def calcGammaCorrectionProfiles(masses, meanWidths, meanIntensityRatios): masses = masses.flatten() profiles = "" for mass, width, intensity in zip(masses, meanWidths, meanIntensityRatios): - profiles += "name=GaussianComptonProfile,Mass=" \ - + str(mass) + ",Width=" + str(width) \ - + ",Intensity=" + str(intensity) + ';' + profiles += ( + "name=GaussianComptonProfile,Mass=" + + str(mass) + + ",Width=" + + str(width) + + ",Intensity=" + + str(intensity) + + ";" + ) print("\n The sample properties for Gamma Correction are:\n", profiles) return profiles @@ -874,7 +1110,6 @@ class resultsObject: """Used to collect results from workspaces and store them in .npz files for testing.""" def __init__(self, ic): - allIterNcp = [] allFitWs = [] allTotNcp = [] @@ -883,21 +1118,21 @@ def __init__(self, ic): allMeanIntensities = [] allStdWidths = [] allStdIntensities = [] - j=0 + j = 0 while True: try: - wsIterName = ic.name+str(j) + wsIterName = ic.name + str(j) # Extract ws that were fitted ws = mtd[wsIterName] allFitWs.append(ws.extractY()) # Extract total ncp - totNcpWs = mtd[wsIterName+"_TOF_Fitted_Profiles"] + totNcpWs = mtd[wsIterName + "_TOF_Fitted_Profiles"] allTotNcp.append(totNcpWs.extractY()) # Extract best fit parameters - fitParTable = mtd[wsIterName+"_Best_Fit_NCP_Parameters"] + fitParTable = mtd[wsIterName + "_Best_Fit_NCP_Parameters"] bestFitPars = [] for key in fitParTable.keys(): bestFitPars.append(fitParTable.column(key)) @@ -906,9 +1141,11 @@ def __init__(self, ic): # Extract individual ncp allNCP = [] i = 0 - while True: # By default, looks for all ncp ws until it breaks + while True: # By default, looks for all ncp ws until it breaks try: - ncpWsToAppend = mtd[wsIterName+"_TOF_Fitted_Profile_"+str(i)] + ncpWsToAppend = mtd[ + wsIterName + "_TOF_Fitted_Profile_" + str(i) + ] allNCP.append(ncpWsToAppend.extractY()) i += 1 except KeyError: @@ -923,7 +1160,7 @@ def __init__(self, ic): allMeanIntensities.append(meansTable.column("Mean Intensities")) allStdIntensities.append(meansTable.column("Std Intensities")) - j+=1 + j += 1 except KeyError: break @@ -953,12 +1190,14 @@ def save(self): self.all_tot_ncp[:, self.maskedDetectorIdx, :] = np.nan savePath = self.resultsSavePath - np.savez(savePath, - all_fit_workspaces=self.all_fit_workspaces, - all_spec_best_par_chi_nit=self.all_spec_best_par_chi_nit, - all_mean_widths=self.all_mean_widths, - all_mean_intensities=self.all_mean_intensities, - all_std_widths=self.all_std_widths, - all_std_intensities=self.all_std_intensities, - all_tot_ncp=self.all_tot_ncp, - all_ncp_for_each_mass=self.all_ncp_for_each_mass) + np.savez( + savePath, + all_fit_workspaces=self.all_fit_workspaces, + all_spec_best_par_chi_nit=self.all_spec_best_par_chi_nit, + all_mean_widths=self.all_mean_widths, + all_mean_intensities=self.all_mean_intensities, + all_std_widths=self.all_std_widths, + all_std_intensities=self.all_std_intensities, + all_tot_ncp=self.all_tot_ncp, + all_ncp_for_each_mass=self.all_ncp_for_each_mass, + ) diff --git a/EVSVesuvio/vesuvio_analysis/bootstrap.py b/EVSVesuvio/vesuvio_analysis/bootstrap.py index 8677864f..f64465d1 100644 --- a/EVSVesuvio/vesuvio_analysis/bootstrap.py +++ b/EVSVesuvio/vesuvio_analysis/bootstrap.py @@ -1,6 +1,12 @@ from EVSVesuvio.vesuvio_analysis.fit_in_yspace import fitInYSpaceProcedure -from EVSVesuvio.vesuvio_analysis.procedures import runJointBackAndForwardProcedure, runIndependentIterativeProcedure -from EVSVesuvio.vesuvio_analysis.ICHelpers import buildFinalWSName, noOfHistsFromTOFBinning +from EVSVesuvio.vesuvio_analysis.procedures import ( + runJointBackAndForwardProcedure, + runIndependentIterativeProcedure, +) +from EVSVesuvio.vesuvio_analysis.ICHelpers import ( + buildFinalWSName, + noOfHistsFromTOFBinning, +) from mantid.api import AnalysisDataService, mtd from mantid.simpleapi import CloneWorkspace, SaveNexus, Load, SumSpectra from scipy import stats @@ -8,19 +14,21 @@ from pathlib import Path import time import matplotlib.pyplot as plt + plt.style.use("ggplot") currentPath = Path(__file__).parent.absolute() def runBootstrap(bckwdIC, fwdIC, bootIC, yFitIC): - checkValidInput(bootIC) - checkOutputDirExists(bckwdIC, fwdIC, bootIC) # Checks to see if those directories exits already + checkOutputDirExists( + bckwdIC, fwdIC, bootIC + ) # Checks to see if those directories exits already askUserConfirmation(bckwdIC, fwdIC, bootIC) AnalysisDataService.clear() - if bootIC.bootstrapType=="JACKKNIFE": + if bootIC.bootstrapType == "JACKKNIFE": return JackknifeProcedure(bckwdIC, fwdIC, bootIC, yFitIC) return bootstrapProcedure(bckwdIC, fwdIC, bootIC, yFitIC) @@ -28,8 +36,9 @@ def runBootstrap(bckwdIC, fwdIC, bootIC, yFitIC): def checkValidInput(bootIC): boot = bootIC.bootstrapType - assert (boot=="JACKKNIFE") | (boot=="BOOT_GAUSS_ERRS") | (boot=="BOOT_RESIDUALS"), \ - "bootstrapType not recognized. Options: 'JACKKNIFE', 'BOOT_GAUSS_ERRS', 'BOOT_RESIDUALS'" + assert ( + (boot == "JACKKNIFE") | (boot == "BOOT_GAUSS_ERRS") | (boot == "BOOT_RESIDUALS") + ), "bootstrapType not recognized. Options: 'JACKKNIFE', 'BOOT_GAUSS_ERRS', 'BOOT_RESIDUALS'" def checkOutputDirExists(bckwdIC, fwdIC, bootIC): @@ -37,19 +46,21 @@ def checkOutputDirExists(bckwdIC, fwdIC, bootIC): return proc = bootIC.procedure - if (proc=="BACKWARD") | (proc=="JOINT"): + if (proc == "BACKWARD") | (proc == "JOINT"): checkOutDirIC(bckwdIC) - if (proc=="FORWARD") | (proc=="JOINT"): + if (proc == "FORWARD") | (proc == "JOINT"): checkOutDirIC(fwdIC) return def checkOutDirIC(IC): if IC.bootSavePath.is_file() or IC.bootYFitSavePath.is_file(): - print(f"\nOutput data files were detected:" - f"\n{IC.bootSavePath.name}\n{IC.bootYFitSavePath.name}" - f"\nAborting Run of Bootstrap to prevent overwriting data." - f"\nTo avoid this issue you can change the number of samples to run.") + print( + f"\nOutput data files were detected:" + f"\n{IC.bootSavePath.name}\n{IC.bootYFitSavePath.name}" + f"\nAborting Run of Bootstrap to prevent overwriting data." + f"\nTo avoid this issue you can change the number of samples to run." + ) raise ValueError("Output data directories already exist. Aborted Bootstrap.") return @@ -58,25 +69,24 @@ def JackknifeProcedure(bckwdIC, fwdIC, bootIC, yFitIC): assert bootIC.procedure is not None proc = bootIC.procedure - if (proc=="FORWARD") | (proc=="BACKWARD"): + if (proc == "FORWARD") | (proc == "BACKWARD"): return bootstrapProcedure(bckwdIC, fwdIC, bootIC, yFitIC) - elif proc=="JOINT": # Do the Jackknife procedure separately - + elif proc == "JOINT": # Do the Jackknife procedure separately # Run original procedure to change fwdIC from running backward runOriginalBeforeBootstrap(bckwdIC, fwdIC, bootIC, yFitIC) - bootIC.procedure="BACKWARD" - bootIC.fitInYSpace="BACKWARD" + bootIC.procedure = "BACKWARD" + bootIC.fitInYSpace = "BACKWARD" bckwdJackRes = bootstrapProcedure(bckwdIC, fwdIC, bootIC, yFitIC) - bootIC.procedure="FORWARD" - bootIC.fitInYSpace="FORWARD" + bootIC.procedure = "FORWARD" + bootIC.fitInYSpace = "FORWARD" fwdJackRes = bootstrapProcedure(bckwdIC, fwdIC, bootIC, yFitIC) - return {**bckwdJackRes, **fwdJackRes} # For consistency + return {**bckwdJackRes, **fwdJackRes} # For consistency else: - raise ValueError ("Bootstrap procedure not recognized.") + raise ValueError("Bootstrap procedure not recognized.") def bootstrapProcedure(bckwdIC, fwdIC, bootIC, yFitIC): @@ -86,12 +96,16 @@ def bootstrapProcedure(bckwdIC, fwdIC, bootIC, yFitIC): Chooses fast or slow (correct) version of bootstrap depending on flag set in bootIC. Performs either independent or joint procedure depending of len(inputIC). """ - if bootIC.bootstrapType=="JACKKNIFE": - assert bootIC.procedure!='JOINT', "'JOINT' mode should not have reached Jackknife here." + if bootIC.bootstrapType == "JACKKNIFE": + assert ( + bootIC.procedure != "JOINT" + ), "'JOINT' mode should not have reached Jackknife here." AnalysisDataService.clear() - parentResults, parentWSnNCPs = runOriginalBeforeBootstrap(bckwdIC, fwdIC, bootIC, yFitIC) + parentResults, parentWSnNCPs = runOriginalBeforeBootstrap( + bckwdIC, fwdIC, bootIC, yFitIC + ) corrCoefs = autoCorrResiduals(parentWSnNCPs) nSamples = chooseNSamples(bootIC, parentWSnNCPs) @@ -105,20 +119,26 @@ def bootstrapProcedure(bckwdIC, fwdIC, bootIC, yFitIC): # Form each bootstrap workspace and run ncp fit with MS corrections for i in range(iStart, iEnd): AnalysisDataService.clear() - plt.close("all") # Not sure if previous step clears plt figures, so introduced this step to be safe + plt.close( + "all" + ) # Not sure if previous step clears plt figures, so introduced this step to be safe try: - sampleInputWS, parentWS = createSampleWS(parentWSNCPSavePaths, i, bootIC) # Creates ith sample + sampleInputWS, parentWS = createSampleWS( + parentWSNCPSavePaths, i, bootIC + ) # Creates ith sample except JackMaskCol: - continue # If Jackknife column already masked, skip to next column + continue # If Jackknife column already masked, skip to next column formSampleIC(bckwdIC, fwdIC, bootIC, sampleInputWS, parentWS) try: - iterResults = runMainProcedure(bckwdIC, fwdIC, bootIC, yFitIC) # Conversion to YSpace with masked column + iterResults = runMainProcedure( + bckwdIC, fwdIC, bootIC, yFitIC + ) # Conversion to YSpace with masked column except AssertionError: - continue # If the procedure fails, skip to next iteration + continue # If the procedure fails, skip to next iteration - storeBootIter(bootResults, i, iterResults) # Stores results for each iteration + storeBootIter(bootResults, i, iterResults) # Stores results for each iteration saveBootstrapResults(bootResults, bckwdIC, fwdIC) return bootResults @@ -126,24 +146,28 @@ def bootstrapProcedure(bckwdIC, fwdIC, bootIC, yFitIC): def askUserConfirmation(bckwdIC, fwdIC, bootIC): """Estimates running time for all samples and asks the user to confirm the run.""" - if not(bootIC.userConfirmation): # Skip user confirmation + if not (bootIC.userConfirmation): # Skip user confirmation return - tDict = storeRunnningTime(fwdIC, bckwdIC, bootIC) # Run times file path stores in bootIC + tDict = storeRunnningTime( + fwdIC, bckwdIC, bootIC + ) # Run times file path stores in bootIC proc = bootIC.procedure runTime = 0 - if (proc=="BACKWARD") | (proc=="JOINT"): + if (proc == "BACKWARD") | (proc == "JOINT"): runTime += calcRunTime(bckwdIC, tDict["tBackNoMS"], tDict["tBackPerMS"], bootIC) - if (proc=="FORWARD") | (proc=="JOINT"): + if (proc == "FORWARD") | (proc == "JOINT"): runTime += calcRunTime(fwdIC, tDict["tFowNoMS"], tDict["tFowPerMS"], bootIC) - userInput = input(f"\n\nEstimated time for Bootstrap procedure: {runTime/60:.1f} hours.\nProceed? (y/n): ") + userInput = input( + f"\n\nEstimated time for Bootstrap procedure: {runTime/60:.1f} hours.\nProceed? (y/n): " + ) if (userInput == "y") or (userInput == "Y"): return else: - raise KeyboardInterrupt ("Bootstrap procedure interrupted.") + raise KeyboardInterrupt("Bootstrap procedure interrupted.") def storeRunnningTime(fwdIC, bckwdIC, bootIC): @@ -151,18 +175,20 @@ def storeRunnningTime(fwdIC, bckwdIC, bootIC): savePath = bootIC.runTimesPath - if not(savePath.is_file()): + if not (savePath.is_file()): with open(savePath, "w") as txtFile: - txtFile.write("This file stores run times to estimate Bootstrap total run time.") + txtFile.write( + "This file stores run times to estimate Bootstrap total run time." + ) txtFile.write("\nTime in minutes.\n\n") resDict = {} with open(savePath, "r") as txtFile: for line in txtFile: - if line[0]=="{": # If line contains dictionary + if line[0] == "{": # If line contains dictionary resDict = eval(line) - if len(resDict)<4: + if len(resDict) < 4: # ans = input("Did not find necessary information to estimate runtime. Will run a short routine to store an # estimate. Please wait until this is finished. Press any key to continue.") resDict = buildRunTimes(fwdIC, bckwdIC) @@ -182,12 +208,12 @@ def buildRunTimes(fwdIC, bckwdIC): t0 = time.time() runIndependentIterativeProcedure(IC) t1 = time.time() - resDict["t"+mode+key] = (t1-t0) / 60 + resDict["t" + mode + key] = (t1 - t0) / 60 # Restore starting value IC.noOfMSIterations = oriMS # Correct times of only MS by subtacting time spend on fitting ncps - resDict["t"+mode+"PerMS"] -= 2 * resDict["t"+mode+"NoMS"] + resDict["t" + mode + "PerMS"] -= 2 * resDict["t" + mode + "NoMS"] return resDict @@ -196,20 +222,20 @@ def calcRunTime(IC, tNoMS, tPerMS, bootIC): if bootIC.skipMSIterations: timePerSample = tNoMS else: - timePerSample = tNoMS + (IC.noOfMSIterations) * (tNoMS+tPerMS) + timePerSample = tNoMS + (IC.noOfMSIterations) * (tNoMS + tPerMS) nSamples = bootIC.nSamples - if bootIC.bootstrapType=="JACKKNIFE": + if bootIC.bootstrapType == "JACKKNIFE": nSamples = 3 if bootIC.runningTest else noOfHistsFromTOFBinning(IC) - return nSamples * timePerSample + return nSamples * timePerSample def chooseLoopRange(bootIC, nSamples): iStart = 0 iEnd = nSamples - if bootIC.bootstrapType=="JACKKNIFE" and bootIC.runningTest: - iStart = int(nSamples/2) + if bootIC.bootstrapType == "JACKKNIFE" and bootIC.runningTest: + iStart = int(nSamples / 2) iEnd = iStart + 3 return iStart, iEnd @@ -230,19 +256,23 @@ def chooseNSamples(bootIC, parentWSnNCPs: dict): If Jackknife is running, no of samples is the number of bins in the workspace.""" nSamples = bootIC.nSamples - if bootIC.bootstrapType=="JACKKNIFE": - assert len(parentWSnNCPs) == 2, "Running Jackknife, supports only one IC at a time." - if bootIC.procedure=="FORWARD": + if bootIC.bootstrapType == "JACKKNIFE": + assert ( + len(parentWSnNCPs) == 2 + ), "Running Jackknife, supports only one IC at a time." + if bootIC.procedure == "FORWARD": key = "fwdNCP" - elif bootIC.procedure=="BACKWARD": + elif bootIC.procedure == "BACKWARD": key = "bckwdNCP" - nSamples = parentWSnNCPs[key].blocksize() # Number of cols from ncp workspace, accounts for missing last col or not + nSamples = parentWSnNCPs[ + key + ].blocksize() # Number of cols from ncp workspace, accounts for missing last col or not return nSamples def setICsToDefault(bckwdIC, fwdIC, yFitIC): - """Disables some features of yspace fit, makes sure the default """ + """Disables some features of yspace fit, makes sure the default""" # Disable Minos if yFitIC.runMinos: @@ -268,32 +298,38 @@ def runMainProcedure(bckwdIC, fwdIC, bootIC, yFitIC): resultsDict = {} - if (bootIC.procedure=="FORWARD") | (bootIC.procedure=="BACKWARD"): - - for mode, IC, key in zip(["FORWARD", "BACKWARD"], [fwdIC, bckwdIC], ["fwd", "bckwd"]): - - if bootIC.procedure==mode: - wsFinal, bckwdScatRes = runIndependentIterativeProcedure(IC, clearWS=False) - resultsDict[key+"Scat"] = bckwdScatRes - - if bootIC.bootstrapType=="JACKKNIFE": - yFitIC.maskTypeProcedure = 'NAN' # Enable NAN averaging in y-space fit + if (bootIC.procedure == "FORWARD") | (bootIC.procedure == "BACKWARD"): + for mode, IC, key in zip( + ["FORWARD", "BACKWARD"], [fwdIC, bckwdIC], ["fwd", "bckwd"] + ): + if bootIC.procedure == mode: + wsFinal, bckwdScatRes = runIndependentIterativeProcedure( + IC, clearWS=False + ) + resultsDict[key + "Scat"] = bckwdScatRes + + if bootIC.bootstrapType == "JACKKNIFE": + yFitIC.maskTypeProcedure = ( + "NAN" # Enable NAN averaging in y-space fit + ) bckwdYFitRes = fitInYSpaceProcedure(yFitIC, IC, wsFinal) - resultsDict[key+"YFit"] = bckwdYFitRes - - elif bootIC.procedure=="JOINT": + resultsDict[key + "YFit"] = bckwdYFitRes - ws, bckwdScatRes, fwdScatRes = runJointBackAndForwardProcedure(bckwdIC, fwdIC, clearWS=False) + elif bootIC.procedure == "JOINT": + ws, bckwdScatRes, fwdScatRes = runJointBackAndForwardProcedure( + bckwdIC, fwdIC, clearWS=False + ) resultsDict["bckwdScat"] = bckwdScatRes resultsDict["fwdScat"] = fwdScatRes - for mode, IC, key in zip(["FORWARD", "BACKWARD"], [fwdIC, bckwdIC], ["fwd", "bckwd"]): - - if (bootIC.fitInYSpace==mode) | (bootIC.fitInYSpace=="JOINT"): + for mode, IC, key in zip( + ["FORWARD", "BACKWARD"], [fwdIC, bckwdIC], ["fwd", "bckwd"] + ): + if (bootIC.fitInYSpace == mode) | (bootIC.fitInYSpace == "JOINT"): wsName = buildFinalWSName(IC.scriptName, mode, IC) fwdYFitRes = fitInYSpaceProcedure(yFitIC, IC, mtd[wsName]) - resultsDict[key+"YFit"] = fwdYFitRes + resultsDict[key + "YFit"] = fwdYFitRes else: raise ValueError("Bootstrap procedure not recognized.") @@ -307,17 +343,19 @@ def selectParentWorkspaces(bckwdIC, fwdIC, bootIC): """ parentWSnNCPsDict = {} - for mode, IC, key in zip(["FORWARD", "BACKWARD"], [fwdIC, bckwdIC], ["fwd", "bckwd"]): - - if (bootIC.procedure==mode) | (bootIC.procedure=="JOINT"): + for mode, IC, key in zip( + ["FORWARD", "BACKWARD"], [fwdIC, bckwdIC], ["fwd", "bckwd"] + ): + if (bootIC.procedure == mode) | (bootIC.procedure == "JOINT"): + wsIter = ( + str(IC.noOfMSIterations) if bootIC.skipMSIterations else "0" + ) # In case of skipping MS, select very last corrected ws - wsIter = str(IC.noOfMSIterations) if bootIC.skipMSIterations else "0" # In case of skipping MS, select very last corrected ws + parentWS = mtd[IC.name + wsIter] + parentNCP = mtd[parentWS.name() + "_TOF_Fitted_Profiles"] - parentWS = mtd[IC.name+wsIter] - parentNCP = mtd[parentWS.name()+"_TOF_Fitted_Profiles"] - - parentWSnNCPsDict[key+"WS"] = parentWS - parentWSnNCPsDict[key+"NCP"] = parentNCP + parentWSnNCPsDict[key + "WS"] = parentWS + parentWSnNCPsDict[key + "NCP"] = parentNCP return parentWSnNCPsDict @@ -328,23 +366,22 @@ def autoCorrResiduals(parentWSnNCP: dict): """ corrCoefs = {} for mode in ["bckwd", "fwd"]: - - try: # Look for workspaces in dictionary, skip if not present - parentWS = parentWSnNCP[mode+"WS"] - parentNCP = parentWSnNCP[mode+"NCP"] + try: # Look for workspaces in dictionary, skip if not present + parentWS = parentWSnNCP[mode + "WS"] + parentNCP = parentWSnNCP[mode + "NCP"] except KeyError: continue totNcp = parentNCP.extractY()[:, :] - dataY = parentWS.extractY()[:, :totNcp.shape[1]] # Missing last column or not + dataY = parentWS.extractY()[:, : totNcp.shape[1]] # Missing last column or not residuals = dataY - totNcp - lag = 1 # For lag-plot of self-correlation + lag = 1 # For lag-plot of self-correlation corr = np.zeros((len(residuals), 2)) for i, rowRes in enumerate(residuals): corr[i] = stats.pearsonr(rowRes[:-lag], rowRes[lag:]) - corrCoefs[mode+"Scat"] = corr + corrCoefs[mode + "Scat"] = corr return corrCoefs @@ -356,17 +393,19 @@ def initializeResults(parentResults: dict, nSamples, corrCoefs): bootResultObjs = {} for key in ["fwd", "bckwd"]: - - if key+"Scat" in parentResults: - bootResultObjs[key+"Scat"] = BootScattResults(parentResults[key+"Scat"], nSamples, corrCoefs[key+"Scat"]) - - if key+"YFit" in parentResults: - bootResultObjs[key+"YFit"] = BootYFitResults(parentResults[key+"YFit"], nSamples) + if key + "Scat" in parentResults: + bootResultObjs[key + "Scat"] = BootScattResults( + parentResults[key + "Scat"], nSamples, corrCoefs[key + "Scat"] + ) + + if key + "YFit" in parentResults: + bootResultObjs[key + "YFit"] = BootYFitResults( + parentResults[key + "YFit"], nSamples + ) return bootResultObjs class BootScattResults: - def __init__(self, parentResults, nSamples, corr): self.parentResult = parentResults.all_spec_best_par_chi_nit[-1] self.bootSamples = np.full((nSamples, *self.parentResult.shape), np.nan) @@ -376,16 +415,19 @@ def storeBootIterResults(self, j, bootResult): self.bootSamples[j] = bootResult.all_spec_best_par_chi_nit[-1] def saveResults(self, IC): - np.savez(IC.bootSavePath, boot_samples=self.bootSamples, - parent_result=self.parentResult, corr_residuals=self.corrResiduals) + np.savez( + IC.bootSavePath, + boot_samples=self.bootSamples, + parent_result=self.parentResult, + corr_residuals=self.corrResiduals, + ) def saveLog(self, IC): with open(IC.logFilePath, "a") as logFile: - logFile.write("\n"+IC.bootSavePathLog) + logFile.write("\n" + IC.bootSavePathLog) class BootYFitResults: - def __init__(self, parentResults, nSamples): self.parentPopt = parentResults.popt self.parentPerr = parentResults.perr @@ -395,12 +437,16 @@ def storeBootIterResults(self, j, bootResult): self.bootSamples[j] = bootResult.popt def saveResults(self, IC): - np.savez(IC.bootYFitSavePath, boot_samples=self.bootSamples, - parent_popt=self.parentPopt, parent_perr=self.parentPerr) + np.savez( + IC.bootYFitSavePath, + boot_samples=self.bootSamples, + parent_popt=self.parentPopt, + parent_perr=self.parentPerr, + ) def saveLog(self, IC): with open(IC.logFilePath, "a") as logFile: - logFile.write("\n"+IC.bootYFitSavePathLog) + logFile.write("\n" + IC.bootYFitSavePathLog) def storeBootIter(bootResultObjs: dict, j: int, bootIterResults: dict): @@ -412,16 +458,16 @@ def storeBootIter(bootResultObjs: dict, j: int, bootIterResults: dict): def saveBootstrapResults(bootResultObjs: dict, bckwdIC, fwdIC): for key, IC in zip(["bckwd", "fwd"], [bckwdIC, fwdIC]): for res in ["Scat", "YFit"]: - if key+res in bootResultObjs: - bootResultObjs[key+res].saveResults(IC) + if key + res in bootResultObjs: + bootResultObjs[key + res].saveResults(IC) return def saveBootstrapLogs(bootResultObjs: dict, bckwdIC, fwdIC): for key, IC in zip(["bckwd", "fwd"], [bckwdIC, fwdIC]): for res in ["Scat", "YFit"]: - if key+res in bootResultObjs: - bootResultObjs[key+res].saveLog(IC) + if key + res in bootResultObjs: + bootResultObjs[key + res].saveLog(IC) return @@ -451,17 +497,16 @@ def saveWorkspacesLocally(ws): def createSampleWS(parentWSNCPSavePaths: dict, j: int, bootIC): - boot = bootIC.bootstrapType - if boot=="JACKKNIFE": + if boot == "JACKKNIFE": return createJackknifeWS(parentWSNCPSavePaths, j) - elif boot=="BOOT_RESIDUALS": + elif boot == "BOOT_RESIDUALS": return createBootstrapWS(parentWSNCPSavePaths) - elif boot=="BOOT_GAUSS_ERRS": + elif boot == "BOOT_GAUSS_ERRS": return createBootstrapWS(parentWSNCPSavePaths, drawGauss=True) -def createBootstrapWS(parentWSNCPSavePaths:dict, drawGauss=False): +def createBootstrapWS(parentWSNCPSavePaths: dict, drawGauss=False): """ Creates bootstrap ws replica. Inputs: Experimental (parent) workspace and corresponding NCP total fit @@ -471,46 +516,54 @@ def createBootstrapWS(parentWSNCPSavePaths:dict, drawGauss=False): parentInputWS = {} for key in ["bckwd", "fwd"]: try: - parentWSPath = parentWSNCPSavePaths[key+"WS"] - totNcpWSPath = parentWSNCPSavePaths[key+"NCP"] + parentWSPath = parentWSNCPSavePaths[key + "WS"] + totNcpWSPath = parentWSNCPSavePaths[key + "NCP"] except KeyError: continue parentWS, totNcpWS = loadWorkspacesFromPath(parentWSPath, totNcpWSPath) totNcp = totNcpWS.extractY()[:, :] - dataY = parentWS.extractY()[:, :totNcp.shape[1]] # Missing last col or not - dataE = parentWS.extractE()[:, :totNcp.shape[1]] + dataY = parentWS.extractY()[:, : totNcp.shape[1]] # Missing last col or not + dataE = parentWS.extractE()[:, : totNcp.shape[1]] # Filter out masked columns - maskCols = np.all(dataY==0, axis=0) - dataY, totNcp, dataE = dataY[:, ~maskCols], totNcp[:, ~maskCols], dataE[:, ~maskCols] + maskCols = np.all(dataY == 0, axis=0) + dataY, totNcp, dataE = ( + dataY[:, ~maskCols], + totNcp[:, ~maskCols], + dataE[:, ~maskCols], + ) # Draw DataY from Gaussian distribution if drawGauss: bootDataY = np.random.normal(dataY, dataE) # Mean at dataY, width dataE - else: # Default, resample residuals + else: # Default, resample residuals residuals = dataY - totNcp bootRes = bootstrapResidualsSample(residuals) bootDataY = totNcp + bootRes # Add masked columns as in parent workspace fullBootDataY = np.zeros((len(bootDataY), len(maskCols))) - fullBootDataY[:, ~maskCols] = bootDataY # Set non-masked values + fullBootDataY[:, ~maskCols] = bootDataY # Set non-masked values # Pass dataY onto workspace - wsBoot = CloneWorkspace(parentWS, OutputWorkspace=parentWS.name()+"_Bootstrap") + wsBoot = CloneWorkspace( + parentWS, OutputWorkspace=parentWS.name() + "_Bootstrap" + ) for i, row in enumerate(fullBootDataY): - wsBoot.dataY(i)[:len(row)] = row # Last column will be ignored or not + wsBoot.dataY(i)[: len(row)] = row # Last column will be ignored or not if drawGauss: wsBoot.dataE(i)[:] = np.zeros(wsBoot.readE(i).size) - assert ~np.all(wsBoot.extractY() == parentWS.extractY()), "Bootstrap data not being correctly passed onto ws." + assert ~np.all( + wsBoot.extractY() == parentWS.extractY() + ), "Bootstrap data not being correctly passed onto ws." - bootInputWS[key+"WS"] = wsBoot - parentInputWS[key+"WS"] = parentWS - parentInputWS[key+"NCP"] = totNcpWS + bootInputWS[key + "WS"] = wsBoot + parentInputWS[key + "WS"] = parentWS + parentInputWS[key + "NCP"] = totNcpWS return bootInputWS, parentInputWS @@ -519,7 +572,7 @@ def bootstrapResidualsSample(residuals): bootRes = np.zeros(residuals.shape) for i, res in enumerate(residuals): - rowIdxs = np.random.randint(0, len(res), len(res)) # [low, high) + rowIdxs = np.random.randint(0, len(res), len(res)) # [low, high) bootRes[i] = res[rowIdxs] return bootRes @@ -534,11 +587,13 @@ def createJackknifeWS(parentWSNCPSavePaths: list, j: int): parentInputWS = {} # Jackknife does not have 'JOINT' option # Careful with this step if in future Jackknife allows for 'JOINT' internally - assert len(parentWSNCPSavePaths)==2, "Jackknife can only allow either forward or backward at a time." - for key in ["bckwd", "fwd"]: # Only one iteration is selected at a time + assert ( + len(parentWSNCPSavePaths) == 2 + ), "Jackknife can only allow either forward or backward at a time." + for key in ["bckwd", "fwd"]: # Only one iteration is selected at a time try: - parentWSPath = parentWSNCPSavePaths[key+"WS"] - totNcpWSPath = parentWSNCPSavePaths[key+"NCP"] + parentWSPath = parentWSNCPSavePaths[key + "WS"] + totNcpWSPath = parentWSNCPSavePaths[key + "NCP"] except KeyError: continue @@ -549,21 +604,25 @@ def createJackknifeWS(parentWSNCPSavePaths: list, j: int): jackDataY = dataY.copy() # Skip Jackknife procedure on columns that are already masked - if np.all(jackDataY[:, j]==0): + if np.all(jackDataY[:, j] == 0): raise JackMaskCol - jackDataY[:, j] = 0 # Masks j collumn with zeros + jackDataY[:, j] = 0 # Masks j collumn with zeros # DataE is not masked intentionally, to preserve errors that are used in the normalization of averaged NaN profile - wsJack = CloneWorkspace(parentWS, OutputWorkspace=parentWS.name()+"_Jackknife") + wsJack = CloneWorkspace( + parentWS, OutputWorkspace=parentWS.name() + "_Jackknife" + ) for i, yRow in enumerate(jackDataY): - wsJack.dataY(i)[:] = yRow # Last column will be ignored in ncp fit anyway + wsJack.dataY(i)[:] = yRow # Last column will be ignored in ncp fit anyway - assert np.all(wsJack.extractY() == jackDataY), "Bootstrap data not being correctly passed onto ws." + assert np.all( + wsJack.extractY() == jackDataY + ), "Bootstrap data not being correctly passed onto ws." - jackInputWS[key+"WS"] = wsJack - parentInputWS[key+"WS"] = parentWS - parentInputWS[key+"NCP"] = totNcpWS + jackInputWS[key + "WS"] = wsJack + parentInputWS[key + "WS"] = parentWS + parentInputWS[key + "NCP"] = totNcpWS return jackInputWS, parentInputWS @@ -572,6 +631,7 @@ class JackMaskCol(Exception): Custom exception used only to flag and skip a Jackknife iteration for a column that is already masked. """ + pass @@ -580,22 +640,23 @@ def loadWorkspacesFromPath(*savePaths): for path in savePaths: saveName = path.name.split(".")[0] ws = Load(str(path), OutputWorkspace=saveName) - SumSpectra(ws, OutputWorkspace=ws.name()+"_Sum") + SumSpectra(ws, OutputWorkspace=ws.name() + "_Sum") wsList.append(ws) return wsList -def formSampleIC(bckwdIC, fwdIC, bootIC, sampleInputWS:dict, parentWS:dict): +def formSampleIC(bckwdIC, fwdIC, bootIC, sampleInputWS: dict, parentWS: dict): """Adds atributes to initial conditions to start procedure with sample ws.""" - for mode, IC, key in zip(["FORWARD", "BACKWARD"], [fwdIC, bckwdIC], ["fwd", "bckwd"]): - - if (bootIC.procedure==mode) | (bootIC.procedure=="JOINT"): + for mode, IC, key in zip( + ["FORWARD", "BACKWARD"], [fwdIC, bckwdIC], ["fwd", "bckwd"] + ): + if (bootIC.procedure == mode) | (bootIC.procedure == "JOINT"): IC.runningSampleWS = True if bootIC.skipMSIterations: IC.noOfMSIterations = 0 - IC.sampleWS = sampleInputWS[key+"WS"] - IC.parentWS = parentWS[key+"WS"] + IC.sampleWS = sampleInputWS[key + "WS"] + IC.parentWS = parentWS[key + "WS"] diff --git a/EVSVesuvio/vesuvio_analysis/bootstrap_analysis.py b/EVSVesuvio/vesuvio_analysis/bootstrap_analysis.py index 531fe254..c39d07f3 100644 --- a/EVSVesuvio/vesuvio_analysis/bootstrap_analysis.py +++ b/EVSVesuvio/vesuvio_analysis/bootstrap_analysis.py @@ -1,38 +1,47 @@ from xml.dom import NotFoundErr -from EVSVesuvio.vesuvio_analysis.analysis_functions import calculateMeansAndStds, filterWidthsAndIntensities +from EVSVesuvio.vesuvio_analysis.analysis_functions import ( + calculateMeansAndStds, + filterWidthsAndIntensities, +) from EVSVesuvio.vesuvio_analysis.ICHelpers import setBootstrapDirs from EVSVesuvio.vesuvio_analysis.fit_in_yspace import selectModelAndPars import numpy as np -import matplotlib .pyplot as plt +import matplotlib.pyplot as plt from scipy import stats def runAnalysisOfStoredBootstrap(bckwdIC, fwdIC, yFitIC, bootIC, analysisIC, userCtr): - - if not(analysisIC.runAnalysis): + if not (analysisIC.runAnalysis): return - setBootstrapDirs(bckwdIC, fwdIC, bootIC, yFitIC) # Same function used to store data, to check below if dirs exist + setBootstrapDirs( + bckwdIC, fwdIC, bootIC, yFitIC + ) # Same function used to store data, to check below if dirs exist for IC in [bckwdIC, fwdIC]: - - if not(IC.bootSavePath.is_file()): + if not (IC.bootSavePath.is_file()): print("Bootstrap data files not found, unable to run analysis!") print(f"{IC.bootSavePath.name}") - continue # If main results are not present, assume ysapce results are also missing + continue # If main results are not present, assume ysapce results are also missing checkLogMatch(IC, isYFitFile=False) - bootParsRaw, parentParsRaw, nSamples, corrResiduals = readBootData(IC.bootSavePath) + bootParsRaw, parentParsRaw, nSamples, corrResiduals = readBootData( + IC.bootSavePath + ) checkResiduals(corrResiduals) - checkBootSamplesVSParent(bootParsRaw, parentParsRaw, IC) # Prints comparison + checkBootSamplesVSParent(bootParsRaw, parentParsRaw, IC) # Prints comparison - bootPars = bootParsRaw.copy() # By default do not filter means, copy to avoid accidental changes + bootPars = ( + bootParsRaw.copy() + ) # By default do not filter means, copy to avoid accidental changes if analysisIC.filterAvg: bootPars = filteredBootMeans(bootParsRaw.copy(), IC) try: print("\nCompare filtered parameters with parent:\n") - checkBootSamplesVSParent(bootPars, parentParsRaw, IC) # Prints comparison + checkBootSamplesVSParent( + bootPars, parentParsRaw, IC + ) # Prints comparison except AssertionError: print("\nUnable to calculate new means of filtered parameters.\n") @@ -40,17 +49,23 @@ def runAnalysisOfStoredBootstrap(bckwdIC, fwdIC, yFitIC, bootIC, analysisIC, use plotRawWidthsAndIntensities(analysisIC, IC, bootPars, parentParsRaw) # Calculate bootstrap histograms for mean widths and intensities - meanWidths, meanIntensities = calculateMeanWidthsIntensities(bootPars, IC, nSamples) + meanWidths, meanIntensities = calculateMeanWidthsIntensities( + bootPars, IC, nSamples + ) # If filer is on, check that it matches original procedure checkMeansProcedure(analysisIC, IC, meanWidths, meanIntensities, bootParsRaw) - plotMeanWidthsAndIntensities(analysisIC, IC, meanWidths, meanIntensities, parentParsRaw) + plotMeanWidthsAndIntensities( + analysisIC, IC, meanWidths, meanIntensities, parentParsRaw + ) plotMeansEvolution(analysisIC, meanWidths, meanIntensities) plot2DHistsWidthsAndIntensities(analysisIC, meanWidths, meanIntensities) - if not(IC.bootYFitSavePath.is_file()): - print("Bootstrap data file for y-space fit not found, unable to run analysis!") + if not (IC.bootYFitSavePath.is_file()): + print( + "Bootstrap data file for y-space fit not found, unable to run analysis!" + ) print(f"{IC.bootYFitSavePath.name}") continue @@ -74,8 +89,10 @@ def checkLogMatch(IC, isYFitFile): if line.split(" : ")[0] == currName: if line.strip("\n") == currentLog: return - raise NotFoundErr(currName+" found but corresponding log does not match.") - raise NotFoundErr(currName+" not found in logs file") + raise NotFoundErr( + currName + " found but corresponding log does not match." + ) + raise NotFoundErr(currName + " not found in logs file") def readBootData(dataPath): @@ -95,7 +112,9 @@ def readBootData(dataPath): print("\nCorrelation of coefficients not found!\n") failMask = np.all(np.isnan(bootParsRaw), axis=(1, 2)) - assert failMask.shape == (len(bootParsRaw),), f"Wrong shape of masking: {failMask.shape} != {bootParsRaw.shape} " + assert failMask.shape == ( + len(bootParsRaw), + ), f"Wrong shape of masking: {failMask.shape} != {bootParsRaw.shape} " if np.sum(failMask) > 0: print(f"\nNo of failed samples: {np.sum(failMask)}") print("\nUsing only good replicas ...\n") @@ -106,21 +125,25 @@ def readBootData(dataPath): print(f"Masked idxs with nans found: {maskedIdxs}") print(f"\nData files found:\n{dataPath.name}") print(f"\nNumber of samples in the file: {nSamples}") - assert ~np.all(bootParsRaw[-1] == parentParsRaw), "Error in Jackknife due to last column." + assert ~np.all( + bootParsRaw[-1] == parentParsRaw + ), "Error in Jackknife due to last column." return bootParsRaw, parentParsRaw, nSamples, corrResiduals def readYFitData(dataPath, yFitIC): """Resulting output has shape(no of pars, no of samples)""" - fitIdx = 0 # 0 to select Minuit, 1 to select LM + fitIdx = 0 # 0 to select Minuit, 1 to select LM bootYFitData = np.load(dataPath) try: bootYFitVals = bootYFitData["boot_samples"] except KeyError: - bootYFitVals = bootYFitData["boot_vals"] # To account for some previous samples - minFitVals = bootYFitVals[:, fitIdx, :-1].T # Select Minuit values and Discard last value chi2 + bootYFitVals = bootYFitData["boot_vals"] # To account for some previous samples + minFitVals = bootYFitVals[ + :, fitIdx, :-1 + ].T # Select Minuit values and Discard last value chi2 failMask = np.all(np.isnan(minFitVals), axis=0) if np.sum(failMask) > 0: @@ -142,7 +165,7 @@ def checkResiduals(corrRes): return corrCoef = corrRes[:, 0] - nCorrelated = np.sum(corrCoef>0.5) + nCorrelated = np.sum(corrCoef > 0.5) print(f"\nNumber of spectra with pearson r > 0.5: {nCorrelated}") return @@ -159,12 +182,16 @@ def checkBootSamplesVSParent(bestPars, parentPars, IC): meanBootWidths = np.mean(bootWidths, axis=0) meanBootIntensities = np.mean(bootIntensities, axis=0) - avgWidths, stdWidths, avgInt, stdInt = calculateMeansAndStds(meanBootWidths.T, meanBootIntensities.T, IC) + avgWidths, stdWidths, avgInt, stdInt = calculateMeansAndStds( + meanBootWidths.T, meanBootIntensities.T, IC + ) parentWidths = parentPars[:, 1::3] parentIntensities = parentPars[:, 0::3] - avgWidthsP, stdWidthsP, avgIntP, stdIntP = calculateMeansAndStds(parentWidths.T, parentIntensities.T, IC) + avgWidthsP, stdWidthsP, avgIntP, stdIntP = calculateMeansAndStds( + parentWidths.T, parentIntensities.T, IC + ) print("\nComparing Bootstrap means with parent means:\n") printResults(avgWidths, stdWidths, "Boot Widths") @@ -179,7 +206,9 @@ def printResults(arrM, arrE, mode): print(f"{mode} {i}: {m:>6.3f} \u00B1 {e:<6.3f}") -def filteredBootMeans(bestPars, IC): # Pass IC just to check flag for preliminary procedure +def filteredBootMeans( + bestPars, IC +): # Pass IC just to check flag for preliminary procedure """Use same filtering function used on original procedure""" # Extract Widths and Intensities from bootstrap samples @@ -188,7 +217,9 @@ def filteredBootMeans(bestPars, IC): # Pass IC just to check flag for prelimina # Perform the filter for i, (widths, intensities) in enumerate(zip(bootWidths, bootIntensities)): - filteredWidths, filteredIntensities = filterWidthsAndIntensities(widths.T, intensities.T, IC) + filteredWidths, filteredIntensities = filterWidthsAndIntensities( + widths.T, intensities.T, IC + ) bootWidths[i] = filteredWidths.T bootIntensities[i] = filteredIntensities.T @@ -206,7 +237,7 @@ def plotRawWidthsAndIntensities(analysisIC, IC, bootPars, parentPars): Plots histogram of means over spectra for each width or intensity. """ - if not(analysisIC.plotRawWidthsIntensities): + if not (analysisIC.plotRawWidthsIntensities): return parentWidths, parentIntensities = extractParentMeans(parentPars, IC) @@ -214,13 +245,16 @@ def plotRawWidthsAndIntensities(analysisIC, IC, bootPars, parentPars): fig, axs = plt.subplots(2, noOfMasses) - for axIdx, startIdx, kind, parentMeans in zip([0, 1], [1, 0], ["Width", "Intensity"], [parentWidths, parentIntensities]): - - for i, j in enumerate(range(startIdx, 3*noOfMasses, 3)): + for axIdx, startIdx, kind, parentMeans in zip( + [0, 1], [1, 0], ["Width", "Intensity"], [parentWidths, parentIntensities] + ): + for i, j in enumerate(range(startIdx, 3 * noOfMasses, 3)): axs[axIdx, i].set_title(f"{kind} {i}") idxSamples = selectRawSamplesPerIdx(bootPars, j) plotHists(axs[axIdx, i], idxSamples, disableCI=True, disableLeg=True) - axs[axIdx, i].axvline(parentMeans[i], 0.75, 0.97, color="b", ls="-", alpha=0.4) + axs[axIdx, i].axvline( + parentMeans[i], 0.75, 0.97, color="b", ls="-", alpha=0.4 + ) plt.show() return @@ -241,16 +275,15 @@ def calcMeansWithOriginalProc(bestPars, IC): bootWidths = bestPars[:, :, 1::3] bootIntensities = bestPars[:, :, 0::3] - bootMeanW = np.zeros((len(bootWidths[0,0,:]), len(bootWidths))) + bootMeanW = np.zeros((len(bootWidths[0, 0, :]), len(bootWidths))) # bootStdW = np.zeros(bootMeanW.shape) bootMeanI = np.zeros(bootMeanW.shape) # bootStdI = np.zeros(bootMeanW.shape) for j, (widths, intensities) in enumerate(zip(bootWidths, bootIntensities)): - meanW, stdW, meanI, stdI = calculateMeansAndStds(widths.T, intensities.T, IC) - bootMeanW[:, j] = meanW # Interested only in the means + bootMeanW[:, j] = meanW # Interested only in the means bootMeanI[:, j] = meanI return bootMeanW, bootMeanI @@ -264,12 +297,12 @@ def calculateMeanWidthsIntensities(bootPars, IC, nSamples): # Calculate bootstrap histograms for mean widths and intensities meanWidths = np.zeros((len(IC.masses), nSamples)) - for i, j in enumerate(range(1, 3*len(IC.masses), 3)): + for i, j in enumerate(range(1, 3 * len(IC.masses), 3)): idxSamples = selectRawSamplesPerIdx(bootPars, j) meanWidths[i, :] = np.nanmean(idxSamples, axis=0) meanIntensities = np.zeros((len(IC.masses), nSamples)) - for i, j in enumerate(range(0, 3*len(IC.masses), 3)): + for i, j in enumerate(range(0, 3 * len(IC.masses), 3)): idxSamples = selectRawSamplesPerIdx(bootPars, j) meanIntensities[i, :] = np.nanmean(idxSamples, axis=0) @@ -279,22 +312,26 @@ def calculateMeanWidthsIntensities(bootPars, IC, nSamples): def checkMeansProcedure(analysisIC, IC, meanWidths, meanIntensities, bootParsRaw): """Checks that filtering and averaging of Bootstrap samples follows the original procedure""" - if not(analysisIC.filterAvg): # When filtering not present, no comparison available + if not ( + analysisIC.filterAvg + ): # When filtering not present, no comparison available return - else: # Check that treatment of data matches original + else: # Check that treatment of data matches original meanWOri, meanIOri = calcMeansWithOriginalProc(bootParsRaw, IC) np.testing.assert_array_almost_equal(meanWOri, meanWidths) np.testing.assert_array_almost_equal(meanIOri, meanIntensities) return -def plotMeanWidthsAndIntensities(analysisIC, IC, meanWidths, meanIntensities, parentParsRaw): +def plotMeanWidthsAndIntensities( + analysisIC, IC, meanWidths, meanIntensities, parentParsRaw +): """ Most informative histograms, shows all mean widhts and intensities of Bootstrap samples """ - if not(analysisIC.plotMeanWidthsIntensities): + if not (analysisIC.plotMeanWidthsIntensities): return parentWidths, parentIntensities = extractParentMeans(parentParsRaw, IC) @@ -305,7 +342,9 @@ def plotMeanWidthsAndIntensities(analysisIC, IC, meanWidths, meanIntensities, pa axs[0].set_title("Histograms of mean Widths") axs[1].set_title("Histograms of mean Intensitiess") - for ax, means, parentMeans in zip(axs.flatten(), [meanWidths, meanIntensities], [parentWidths, parentIntensities]): + for ax, means, parentMeans in zip( + axs.flatten(), [meanWidths, meanIntensities], [parentWidths, parentIntensities] + ): plotHists(ax, means, disableAvg=True, disableCI=True) for pMean in parentMeans: ax.axvline(pMean, 0.75, 0.97, color="b", ls="-", alpha=0.4) @@ -317,7 +356,7 @@ def plotMeanWidthsAndIntensities(analysisIC, IC, meanWidths, meanIntensities, pa def plotMeansEvolution(IC, meanWidths, meanIntensities): """Shows how results of Bootstrap change depending on number of bootstrap samples""" - if not(IC.plotMeansEvolution): + if not (IC.plotMeansEvolution): return fig, axs = plt.subplots(2, 2) @@ -338,8 +377,7 @@ def plotMeansEvolution(IC, meanWidths, meanIntensities): def plotMeansEvolutionYFit(analysisIC, minuitFitVals): - - if not(analysisIC.plotMeansEvolution): + if not (analysisIC.plotMeansEvolution): return fig, ax = plt.subplots(2, 1) @@ -354,7 +392,6 @@ def plotMeansEvolutionYFit(analysisIC, minuitFitVals): def plotMeansOverNoSamples(ax, bootMeans, normFlag=False): - nSamples = len(bootMeans[0]) assert nSamples >= 10, "To plot evolution of means, need at least 10 samples!" noOfPoints = int(nSamples / 10) @@ -368,7 +405,7 @@ def plotMeansOverNoSamples(ax, bootMeans, normFlag=False): mean = np.mean(subSample, axis=1) - bounds = np.percentile(subSample, [16, 68+16], axis=1).T # 1 standard dev + bounds = np.percentile(subSample, [16, 68 + 16], axis=1).T # 1 standard dev assert bounds.shape == (len(subSample), 2), f"Wrong shape: {bounds.shape}" # errors = bounds - mean[:, np.newaxis] @@ -378,10 +415,18 @@ def plotMeansOverNoSamples(ax, bootMeans, normFlag=False): meansFinal = sampleMeans boundsFinal = sampleBounds ylabel = "Means and Errors Values" - if normFlag: # Rescale and normalize to last value, so all points are converging to zero + if ( + normFlag + ): # Rescale and normalize to last value, so all points are converging to zero lastValue = sampleMeans[:, -1][:, np.newaxis] - meansFinal = (sampleMeans - lastValue) / lastValue * 100 # Percent change to last value - boundsFinal = (sampleBounds - lastValue[:, np.newaxis, :]) / lastValue[:, np.newaxis, :] * 100 + meansFinal = ( + (sampleMeans - lastValue) / lastValue * 100 + ) # Percent change to last value + boundsFinal = ( + (sampleBounds - lastValue[:, np.newaxis, :]) + / lastValue[:, np.newaxis, :] + * 100 + ) ylabel = "Percent Change [%]" for i, (means, errors) in enumerate(zip(meansFinal, boundsFinal)): @@ -394,11 +439,12 @@ def plotMeansOverNoSamples(ax, bootMeans, normFlag=False): def plot2DHistsWidthsAndIntensities(IC, meanWidths, meanIntensities): - - if not(IC.plot2DHists): + if not (IC.plot2DHists): return - assert meanWidths.shape == meanIntensities.shape, "Widths and Intensities need to be the same shape." + assert ( + meanWidths.shape == meanIntensities.shape + ), "Widths and Intensities need to be the same shape." plot2DHists(meanWidths, "Widths") plot2DHists(meanIntensities, "Intensities") @@ -413,34 +459,35 @@ def plot2DHists(bootSamples, mode): for i in range(plotSize): for j in range(plotSize): - # axs[i, j].axis("off") - if j>i: + if j > i: axs[i, j].set_visible(False) - elif i==j: - if i>0: - orientation="horizontal" + elif i == j: + if i > 0: + orientation = "horizontal" else: - orientation="vertical" + orientation = "vertical" axs[i, j].hist(bootSamples[i], orientation=orientation) else: axs[i, j].hist2d(bootSamples[j], bootSamples[i]) - if j==0: + if j == 0: axs[i, j].set_ylabel(f"idx {i}") else: axs[i, j].get_yaxis().set_ticks([]) - if i==plotSize-1: + if i == plotSize - 1: axs[i, j].set_xlabel(f"idx{j}") else: axs[i, j].get_xaxis().set_ticks([]) - axs[i, j].set_title(f"r = {stats.pearsonr(bootSamples[i], bootSamples[j])[0]:.3f}") + axs[i, j].set_title( + f"r = {stats.pearsonr(bootSamples[i], bootSamples[j])[0]:.3f}" + ) fig.suptitle(f"1D and 2D Histograms of {mode}") @@ -448,10 +495,9 @@ def plot2DHists(bootSamples, mode): def printYFitParentPars(yFitIC, parentPopt, parentPerr): - model, defaultPars, sharedPars = selectModelAndPars(yFitIC.fitModel) sig = [p for p in defaultPars] - sig = ["y0"] + sig # Add intercept from outside function parameters + sig = ["y0"] + sig # Add intercept from outside function parameters print("\nParent parameters of y-sapce fit:\n") for p, m, e in zip(sig, parentPopt, parentPerr): @@ -461,16 +507,20 @@ def printYFitParentPars(yFitIC, parentPopt, parentPerr): def plotYFitHists(analysisIC, yFitIC, yFitHists): """Histogram for each parameter of model used for fit in y-space.""" - if not(analysisIC.plotYFitHists): + if not (analysisIC.plotYFitHists): return # Plot each parameter in an individual histogram - fig, axs = plt.subplots(2, int(np.ceil(len(yFitHists)/2)), figsize=(12, 7), tight_layout=True) + fig, axs = plt.subplots( + 2, int(np.ceil(len(yFitHists) / 2)), figsize=(12, 7), tight_layout=True + ) # To label each histogram, extract signature of function used for the fit model, defaultPars, sharedPars = selectModelAndPars(yFitIC.fitModel) - sig = [p for p in defaultPars] # Assumes defaultPars in same order as signature of the function - sig = ["y0"] + sig # Add intercept from outside function parameters + sig = [ + p for p in defaultPars + ] # Assumes defaultPars in same order as signature of the function + sig = ["y0"] + sig # Add intercept from outside function parameters for i, (ax, hist, par) in enumerate(zip(axs.flatten(), yFitHists, sig)): ax.set_title(f"Fit Parameter: {par}") @@ -478,7 +528,7 @@ def plotYFitHists(analysisIC, yFitIC, yFitHists): # Hide plots not in use: for ax in axs.flat: - if not ax.lines: # If empty list + if not ax.lines: # If empty list ax.set_visible(False) plt.show() @@ -486,8 +536,7 @@ def plotYFitHists(analysisIC, yFitIC, yFitHists): def plot2DHistsYFit(analysisIC, minuitFitVals): - - if not(analysisIC.plot2DHists): + if not (analysisIC.plot2DHists): return plot2DHists(minuitFitVals, "Y-space Fit parameters") @@ -499,28 +548,27 @@ def plotHists(ax, samples, disableCI=False, disableLeg=False, disableAvg=False): # ax.set_title(f"Histogram of {title}") for i, bootHist in enumerate(samples): - - if np.all(bootHist==0) or np.all(np.isnan(bootHist)): + if np.all(bootHist == 0) or np.all(np.isnan(bootHist)): continue mean = np.nanmean(bootHist) - bounds = np.percentile(bootHist, [16, 68+16]) # 1 std: 68%, 2 std: 95% + bounds = np.percentile(bootHist, [16, 68 + 16]) # 1 std: 68%, 2 std: 95% errors = bounds - mean leg = f"Row {i}: {mean:>6.3f} +{errors[1]:.3f} {errors[0]:.3f}" ax.hist(bootHist, histtype="step", label=leg, linewidth=1) ax.axvline(mean, 0.9, 0.97, color="k", ls="--", alpha=0.4) - if not(disableCI): + if not (disableCI): ax.axvspan(bounds[0], bounds[1], alpha=0.1, color="b") - if not(disableAvg): + if not (disableAvg): # Plot average over histograms avgHist = np.nanmean(samples, axis=0).flatten() ax.hist(avgHist, color="r", histtype="step", linewidth=2) ax.axvline(np.nanmean(avgHist), 0.75, 0.97, color="r", ls="--", linewidth=2) - if not(disableLeg): + if not (disableLeg): ax.legend(loc="upper center") @@ -542,7 +590,7 @@ def calcCorrWithScatAngle(samples, IC): firstSpec, lastSpec = IC.bootSavePath.name.split("_")[1].split("-") allSpec = ipMatrix[:, 0] - selectedSpec = (allSpec>=int(firstSpec)) & (allSpec>=int(lastSpec)) + selectedSpec = (allSpec >= int(firstSpec)) & (allSpec >= int(lastSpec)) thetas = ipMatrix[selectedSpec, 2] # Scattering angle on third column diff --git a/EVSVesuvio/vesuvio_analysis/deprecated_mantid_global_fit.py b/EVSVesuvio/vesuvio_analysis/deprecated_mantid_global_fit.py index 20135db3..b5d07428 100644 --- a/EVSVesuvio/vesuvio_analysis/deprecated_mantid_global_fit.py +++ b/EVSVesuvio/vesuvio_analysis/deprecated_mantid_global_fit.py @@ -18,25 +18,28 @@ def replaceNansWithZeros(ws): def artificialErrorsInUnphysicalBins(wsJoY): - wsGlobal = CloneWorkspace(InputWorkspace=wsJoY, OutputWorkspace=wsJoY.name()+'_Global') + wsGlobal = CloneWorkspace( + InputWorkspace=wsJoY, OutputWorkspace=wsJoY.name() + "_Global" + ) for j in range(wsGlobal.getNumberHistograms()): - wsGlobal.dataE(j)[wsGlobal.dataE(j)[:]==0] = 0.1 + wsGlobal.dataE(j)[wsGlobal.dataE(j)[:] == 0] = 0.1 - assert np.any(np.isnan(wsGlobal.extractE())) is False, "Nan present in input workspace need to be replaced by " \ - "zeros." + assert np.any(np.isnan(wsGlobal.extractE())) is False, ( + "Nan present in input workspace need to be replaced by " "zeros." + ) return wsGlobal def createOneOverQWs(wsQ): - wsInvQ = CloneWorkspace(InputWorkspace=wsQ, OutputWorkspace=wsQ.name()+"_Inverse") + wsInvQ = CloneWorkspace(InputWorkspace=wsQ, OutputWorkspace=wsQ.name() + "_Inverse") for j in range(wsInvQ.getNumberHistograms()): nonZeroFlag = wsInvQ.dataY(j)[:] != 0 wsInvQ.dataY(j)[nonZeroFlag] = 1 / wsInvQ.dataY(j)[nonZeroFlag] - ZeroIdxs = np.argwhere(wsInvQ.dataY(j)[:]==0) # Indxs of zero elements - if ZeroIdxs.size != 0: # When zeros are present - wsInvQ.dataY(j)[ZeroIdxs[0] - 1] = 0 # Put a zero before the first zero + ZeroIdxs = np.argwhere(wsInvQ.dataY(j)[:] == 0) # Indxs of zero elements + if ZeroIdxs.size != 0: # When zeros are present + wsInvQ.dataY(j)[ZeroIdxs[0] - 1] = 0 # Put a zero before the first zero return wsInvQ @@ -44,7 +47,7 @@ def createOneOverQWs(wsQ): def mantidGlobalFitProcedure(wsGlobal, wsQInv, wsRes, minimizer, yFitIC, wsFirstMass): """Original Implementation of Global Fit using Mantid""" - if yFitIC.fitModel=="SINGLE_GAUSSIAN": + if yFitIC.fitModel == "SINGLE_GAUSSIAN": convolution_template = """ (composite=Convolution,$domains=({0}); name=Resolution,Workspace={1},WorkspaceIndex={0}; @@ -61,7 +64,7 @@ def mantidGlobalFitProcedure(wsGlobal, wsQInv, wsRes, minimizer, yFitIC, wsFirst Sigma=6.0);ties=() ) )""" - elif yFitIC.fitModel=="GC_C4_C6": + elif yFitIC.fitModel == "GC_C4_C6": convolution_template = """ (composite=Convolution,$domains=({0}); name=Resolution,Workspace={1},WorkspaceIndex={0}; @@ -80,40 +83,49 @@ def mantidGlobalFitProcedure(wsGlobal, wsQInv, wsRes, minimizer, yFitIC, wsFirst ) )""" else: - raise ValueError(f"{yFitIC.fitModel} not supported for Global Mantid fit. Supported options: 'SINGLE_GAUSSIAN'," - f"'GC_C4_C6'") + raise ValueError( + f"{yFitIC.fitModel} not supported for Global Mantid fit. Supported options: 'SINGLE_GAUSSIAN'," + f"'GC_C4_C6'" + ) - print('\nGlobal fit in the West domain over 8 mixed banks\n') + print("\nGlobal fit in the West domain over 8 mixed banks\n") widths = [] for bank in range(8): - dets=[bank, bank+8, bank+16, bank+24] + dets = [bank, bank + 8, bank + 16, bank + 24] convolvedFunctionsList = [] ties = ["f0.f1.f1.f1.Sigma=f0.f1.f0.Sigma"] - datasets = {'InputWorkspace' : wsGlobal.name(), - 'WorkspaceIndex' : dets[0]} + datasets = {"InputWorkspace": wsGlobal.name(), "WorkspaceIndex": dets[0]} print("Detectors: ", dets) counter = 0 for i in dets: - print(f"Considering spectrum {wsGlobal.getSpectrumNumbers()[i]}") if wsGlobal.spectrumInfo().isMasked(i): print(f"Skipping masked spectrum {wsGlobal.getSpectrumNumbers()[i]}") continue - thisIterationFunction = convolution_template.format(counter, wsRes.name(), wsQInv.name()) + thisIterationFunction = convolution_template.format( + counter, wsRes.name(), wsQInv.name() + ) convolvedFunctionsList.append(thisIterationFunction) if counter > 0: - if yFitIC.fitModel == "SINGLE_GAUSSIAN": - ties.append('f{0}.f1.f0.Sigma= f{0}.f1.f1.f1.Sigma=f0.f1.f0.Sigma'.format(counter)) + ties.append( + "f{0}.f1.f0.Sigma= f{0}.f1.f1.f1.Sigma=f0.f1.f0.Sigma".format( + counter + ) + ) else: - ties.append('f{0}.f1.f0.Sigma= f{0}.f1.f1.f1.Sigma=f0.f1.f0.Sigma'.format(counter)) - ties.append('f{0}.f1.f0.c4=f0.f1.f0.c4'.format(counter)) - ties.append('f{0}.f1.f1.f1.c3=f0.f1.f1.f1.c3'.format(counter)) + ties.append( + "f{0}.f1.f0.Sigma= f{0}.f1.f1.f1.Sigma=f0.f1.f0.Sigma".format( + counter + ) + ) + ties.append("f{0}.f1.f0.c4=f0.f1.f0.c4".format(counter)) + ties.append("f{0}.f1.f1.f1.c3=f0.f1.f1.f1.c3".format(counter)) # Attach datasets datasets[f"InputWorkspace_{counter}"] = wsGlobal.name() @@ -121,18 +133,42 @@ def mantidGlobalFitProcedure(wsGlobal, wsQInv, wsRes, minimizer, yFitIC, wsFirst counter += 1 multifit_func = f"composite=MultiDomainFunction; {';'.join(convolvedFunctionsList)}; ties=({','.join(ties)})" - minimizer_string = f"{minimizer}, AbsError=0.00001, RealError=0.00001, MaxIterations=2000" + minimizer_string = ( + f"{minimizer}, AbsError=0.00001, RealError=0.00001, MaxIterations=2000" + ) # Unpack dictionary as arguments - Fit(multifit_func, Minimizer=minimizer_string, Output=wsFirstMass.name()+f'_Joy_Mixed_Banks_Bank_{str(bank)}_fit', **datasets) + Fit( + multifit_func, + Minimizer=minimizer_string, + Output=wsFirstMass.name() + f"_Joy_Mixed_Banks_Bank_{str(bank)}_fit", + **datasets, + ) # Select ws with fit results - ws=mtd[wsFirstMass.name()+f'_Joy_Mixed_Banks_Bank_{str(bank)}_fit_Parameters'] + ws = mtd[ + wsFirstMass.name() + f"_Joy_Mixed_Banks_Bank_{str(bank)}_fit_Parameters" + ] print(f"Bank: {str(bank)} -- sigma={ws.cell(2,1)} +/- {ws.cell(2,2)}") - widths.append(ws.cell(2,1)) - - DeleteWorkspace(wsFirstMass.name()+'_Joy_Mixed_Banks_Bank_'+str(bank)+'_fit_NormalisedCovarianceMatrix') - DeleteWorkspace(wsFirstMass.name()+'_Joy_Mixed_Banks_Bank_'+str(bank)+'_fit_Workspaces') - - print('\nAverage hydrogen standard deviation: ',np.mean(widths),' +/- ', np.std(widths)) + widths.append(ws.cell(2, 1)) + + DeleteWorkspace( + wsFirstMass.name() + + "_Joy_Mixed_Banks_Bank_" + + str(bank) + + "_fit_NormalisedCovarianceMatrix" + ) + DeleteWorkspace( + wsFirstMass.name() + + "_Joy_Mixed_Banks_Bank_" + + str(bank) + + "_fit_Workspaces" + ) + + print( + "\nAverage hydrogen standard deviation: ", + np.mean(widths), + " +/- ", + np.std(widths), + ) return widths diff --git a/EVSVesuvio/vesuvio_analysis/fit_in_yspace.py b/EVSVesuvio/vesuvio_analysis/fit_in_yspace.py index df12e3f1..c8ea26e1 100644 --- a/EVSVesuvio/vesuvio_analysis/fit_in_yspace.py +++ b/EVSVesuvio/vesuvio_analysis/fit_in_yspace.py @@ -2,7 +2,7 @@ import numpy as np from mantid.simpleapi import * from scipy import optimize -from scipy import signal +from scipy import signal from pathlib import Path from iminuit import Minuit, cost from iminuit.util import make_func_code, describe @@ -13,13 +13,14 @@ def fitInYSpaceProcedure(yFitIC, IC, wsTOF): - ncpForEachMass = extractNCPFromWorkspaces(wsTOF, IC) wsResSum, wsRes = calculateMantidResolutionFirstMass(IC, yFitIC, wsTOF) wsTOFMass0 = subtractAllMassesExceptFirst(IC, wsTOF, ncpForEachMass) - wsJoY, wsJoYAvg = ySpaceReduction(wsTOFMass0, IC.masses[0], yFitIC, ncpForEachMass[:, 0, :]) + wsJoY, wsJoYAvg = ySpaceReduction( + wsTOFMass0, IC.masses[0], yFitIC, ncpForEachMass[:, 0, :] + ) if yFitIC.symmetrisationFlag: wsJoYAvg = symmetrizeWs(wsJoYAvg) @@ -41,9 +42,13 @@ def fitInYSpaceProcedure(yFitIC, IC, wsTOF): def extractNCPFromWorkspaces(wsFinal, ic): """Extra function to extract ncps from loaded ws in mantid.""" - ncpForEachMass = mtd[wsFinal.name()+"_TOF_Fitted_Profile_0"].extractY()[np.newaxis, :, :] + ncpForEachMass = mtd[wsFinal.name() + "_TOF_Fitted_Profile_0"].extractY()[ + np.newaxis, :, : + ] for i in range(1, ic.noOfMasses): - ncpToAppend = mtd[wsFinal.name()+"_TOF_Fitted_Profile_" + str(i)].extractY()[np.newaxis, :, :] + ncpToAppend = mtd[wsFinal.name() + "_TOF_Fitted_Profile_" + str(i)].extractY()[ + np.newaxis, :, : + ] ncpForEachMass = np.append(ncpForEachMass, ncpToAppend, axis=0) # Ensure shape of ncp matches data @@ -51,7 +56,7 @@ def extractNCPFromWorkspaces(wsFinal, ic): assert shape[0] == ic.noOfMasses assert shape[1] == wsFinal.getNumberHistograms() # Final dimension can be missing last col or not - assert ((shape[2]==wsFinal.blocksize()) | (shape[2]==wsFinal.blocksize()-1)) + assert (shape[2] == wsFinal.blocksize()) | (shape[2] == wsFinal.blocksize() - 1) ncpForEachMass = switchFirstTwoAxis(ncpForEachMass) # Organizes ncp by spectra print("\nExtracted NCP profiles from workspaces.\n") @@ -61,18 +66,24 @@ def extractNCPFromWorkspaces(wsFinal, ic): def calculateMantidResolutionFirstMass(IC, yFitIC, ws): mass = IC.masses[0] - resName = ws.name()+"_Resolution" + resName = ws.name() + "_Resolution" for index in range(ws.getNumberHistograms()): - VesuvioResolution(Workspace=ws,WorkspaceIndex=index,Mass=mass,OutputWorkspaceYSpace="tmp") - Rebin(InputWorkspace="tmp", Params=yFitIC.rebinParametersForYSpaceFit, OutputWorkspace="tmp") - - if index == 0: # Ensures that workspace has desired units - RenameWorkspace("tmp", resName) + VesuvioResolution( + Workspace=ws, WorkspaceIndex=index, Mass=mass, OutputWorkspaceYSpace="tmp" + ) + Rebin( + InputWorkspace="tmp", + Params=yFitIC.rebinParametersForYSpaceFit, + OutputWorkspace="tmp", + ) + + if index == 0: # Ensures that workspace has desired units + RenameWorkspace("tmp", resName) else: AppendSpectra(resName, "tmp", OutputWorkspace=resName) MaskDetectors(resName, WorkspaceIndexList=IC.maskedDetectorIdx) - wsResSum = SumSpectra(InputWorkspace=resName, OutputWorkspace=resName+"_Sum") + wsResSum = SumSpectra(InputWorkspace=resName, OutputWorkspace=resName + "_Sum") normalise_workspace(wsResSum) DeleteWorkspace("tmp") @@ -95,16 +106,18 @@ def subtractAllMassesExceptFirst(ic, ws, ncpForEachMass): dataX, dataY, dataE = extractWS(ws) # Adjust for last column missing or not - dataY[:, :ncpTotalExceptFirst.shape[1]] -= ncpTotalExceptFirst + dataY[:, : ncpTotalExceptFirst.shape[1]] -= ncpTotalExceptFirst # Ignore any masked bins (columns) from initial ws - mask = np.all(ws.extractY()==0, axis=0) + mask = np.all(ws.extractY() == 0, axis=0) dataY[:, mask] = 0 - wsSubMass = CloneWorkspace(InputWorkspace=ws, OutputWorkspace=ws.name()+"_Mass0") + wsSubMass = CloneWorkspace(InputWorkspace=ws, OutputWorkspace=ws.name() + "_Mass0") passDataIntoWS(dataX, dataY, dataE, wsSubMass) MaskDetectors(Workspace=wsSubMass, WorkspaceIndexList=ic.maskedDetectorIdx) - SumSpectra(InputWorkspace=wsSubMass.name(), OutputWorkspace=wsSubMass.name()+"_Sum") + SumSpectra( + InputWorkspace=wsSubMass.name(), OutputWorkspace=wsSubMass.name() + "_Sum" + ) return wsSubMass @@ -120,13 +133,14 @@ def ySpaceReduction(wsTOF, mass0, yFitIC, ncp): rebinPars = yFitIC.rebinParametersForYSpaceFit - if np.any(np.all(wsTOF.extractY()==0, axis=0)): # Masked columns present - - if yFitIC.maskTypeProcedure=="NAN": + if np.any(np.all(wsTOF.extractY() == 0, axis=0)): # Masked columns present + if yFitIC.maskTypeProcedure == "NAN": # Build special workspace to store accumulated points wsJoY = convertToYSpace(wsTOF, mass0) xp = buildXRangeFromRebinPars(yFitIC) - wsJoYB = dataXBining(wsJoY, xp) # Unusual ws with several dataY points per each dataX point + wsJoYB = dataXBining( + wsJoY, xp + ) # Unusual ws with several dataY points per each dataX point # Need normalisation values from NCP masked workspace wsTOFNCP = replaceZerosWithNCP(wsTOF, ncp) @@ -134,18 +148,22 @@ def ySpaceReduction(wsTOF, mass0, yFitIC, ncp): wsJoYNCPN, wsJoYInt = rebinAndNorm(wsJoYNCP, rebinPars) # Normalize spectra of specieal workspace - wsJoYN = Divide(wsJoYB, wsJoYInt, OutputWorkspace=wsJoYB.name()+"_Normalised") + wsJoYN = Divide( + wsJoYB, wsJoYInt, OutputWorkspace=wsJoYB.name() + "_Normalised" + ) wsJoYAvg = weightedAvgXBins(wsJoYN, xp) return wsJoYN, wsJoYAvg - elif yFitIC.maskTypeProcedure=="NCP": + elif yFitIC.maskTypeProcedure == "NCP": wsTOF = replaceZerosWithNCP(wsTOF, ncp) else: - raise ValueError(""" + raise ValueError( + """ Masked TOF bins were found but no valid procedure in y-space fit was selected. Options: 'NAN', 'NCP' - """) + """ + ) wsJoY = convertToYSpace(wsTOF, mass0) wsJoYN, wsJoYI = rebinAndNorm(wsJoY, rebinPars) @@ -154,14 +172,19 @@ def ySpaceReduction(wsTOF, mass0, yFitIC, ncp): def convertToYSpace(wsTOF, mass0): - wsJoY = ConvertToYSpace(wsTOF, Mass=mass0, OutputWorkspace=wsTOF.name()+"_JoY") + wsJoY = ConvertToYSpace(wsTOF, Mass=mass0, OutputWorkspace=wsTOF.name() + "_JoY") return wsJoY def rebinAndNorm(wsJoY, rebinPars): - wsJoYR = Rebin(InputWorkspace=wsJoY, Params=rebinPars, FullBinsOnly=True, OutputWorkspace=wsJoY.name()+"_Rebinned") - wsJoYInt = Integration(wsJoYR, OutputWorkspace=wsJoYR.name()+"_Integrated") - wsJoYNorm = Divide(wsJoYR, wsJoYInt, OutputWorkspace=wsJoYR.name()+"_Normalised") + wsJoYR = Rebin( + InputWorkspace=wsJoY, + Params=rebinPars, + FullBinsOnly=True, + OutputWorkspace=wsJoY.name() + "_Rebinned", + ) + wsJoYInt = Integration(wsJoYR, OutputWorkspace=wsJoYR.name() + "_Integrated") + wsJoYNorm = Divide(wsJoYR, wsJoYInt, OutputWorkspace=wsJoYR.name() + "_Normalised") return wsJoYNorm, wsJoYInt @@ -170,20 +193,24 @@ def replaceZerosWithNCP(ws, ncp): Replaces columns of bins with zeros on dataY with NCP provided. """ dataX, dataY, dataE = extractWS(ws) - mask = np.all(dataY==0, axis=0) # Masked Cols + mask = np.all(dataY == 0, axis=0) # Masked Cols - dataY[:, mask] = ncp[:, mask[:ncp.shape[1]]] # mask of ncp adjusted for last col present or not + dataY[:, mask] = ncp[ + :, mask[: ncp.shape[1]] + ] # mask of ncp adjusted for last col present or not - wsMasked = CloneWorkspace(ws, OutputWorkspace=ws.name()+"_NCPMasked") + wsMasked = CloneWorkspace(ws, OutputWorkspace=ws.name() + "_NCPMasked") passDataIntoWS(dataX, dataY, dataE, wsMasked) - SumSpectra(wsMasked, OutputWorkspace=wsMasked.name()+"_Sum") + SumSpectra(wsMasked, OutputWorkspace=wsMasked.name() + "_Sum") return wsMasked def buildXRangeFromRebinPars(yFitIC): # Range used in case mask is set to NAN - first, step, last = [float(s) for s in yFitIC.rebinParametersForYSpaceFit.split(",")] - xp = np.arange(first, last, step) + step/2 # Correction to match Mantid range + first, step, last = [ + float(s) for s in yFitIC.rebinParametersForYSpaceFit.split(",") + ] + xp = np.arange(first, last, step) + step / 2 # Correction to match Mantid range return xp @@ -194,32 +221,35 @@ def dataXBining(ws, xp): Output ws has several dataY values per dataX point. """ - assert np.min(xp[:-1]-xp[1:]) == np.max(xp[:-1]-xp[1:]), "Bin widths need to be the same." - step = xp[1] - xp[0] # Calculate step from first two numbers + assert np.min(xp[:-1] - xp[1:]) == np.max( + xp[:-1] - xp[1:] + ), "Bin widths need to be the same." + step = xp[1] - xp[0] # Calculate step from first two numbers # Form bins with xp being the centers - bins = np.append(xp, [xp[-1]+step]) - step/2 + bins = np.append(xp, [xp[-1] + step]) - step / 2 dataX, dataY, dataE = extractWS(ws) # Loop below changes only the values of DataX for i, x in enumerate(dataX): - # Select only valid range xr - mask = (xnp.max(bins)) + mask = (x < np.min(bins)) | (x > np.max(bins)) xr = x[~mask] idxs = np.digitize(xr, bins) - newXR = np.array([xp[idx] for idx in idxs-1]) # Bin idx 1 refers to first bin ie idx 0 of centers + newXR = np.array( + [xp[idx] for idx in idxs - 1] + ) # Bin idx 1 refers to first bin ie idx 0 of centers # Pad invalid values with nans newX = x - newX[mask] = np.nan # Cannot use 0 as to not be confused with a dataX value + newX[mask] = np.nan # Cannot use 0 as to not be confused with a dataX value newX[~mask] = newXR - dataX[i] = newX # Update DataX + dataX[i] = newX # Update DataX # Mask DataE values in same places as DataY values - dataE[dataY==0] = 0 + dataE[dataY == 0] = 0 - wsXBins = CloneWorkspace(ws, OutputWorkspace=ws.name()+"_XBinned") + wsXBins = CloneWorkspace(ws, OutputWorkspace=ws.name() + "_XBinned") wsXBins = passDataIntoWS(dataX, dataY, dataE, wsXBins) return wsXBins @@ -230,7 +260,13 @@ def weightedAvgXBins(wsXBins, xp): meansY, meansE = weightedAvgXBinsArr(dataX, dataY, dataE, xp) - wsYSpaceAvg = CreateWorkspace(DataX=xp, DataY=meansY, DataE=meansE, NSpec=1, OutputWorkspace=wsXBins.name()+"_WeightedAvg") + wsYSpaceAvg = CreateWorkspace( + DataX=xp, + DataY=meansY, + DataE=meansE, + NSpec=1, + OutputWorkspace=wsXBins.name() + "_WeightedAvg", + ) return wsYSpaceAvg @@ -247,31 +283,31 @@ def weightedAvgXBinsArr(dataX, dataY, dataE, xp): for i in range(len(xp)): # Perform weighted average over all dataY and dataE values with the same xp[i] # Change shape to column to match weighted average function - pointMask = dataX==xp[i] + pointMask = dataX == xp[i] allY = dataY[pointMask][:, np.newaxis] allE = dataE[pointMask][:, np.newaxis] # If no points were found for a given abcissae - if (np.sum(pointMask)==0): + if np.sum(pointMask) == 0: mY, mE = 0, 0 # Mask with zeros # If one point was found, set to that point - elif (np.sum(pointMask)==1): + elif np.sum(pointMask) == 1: mY, mE = allY.flatten(), allE.flatten() # Weighted avg over all spectra and several points per spectra else: # Case of bootstrap replica with no errors - if np.all(dataE==0): + if np.all(dataE == 0): mY = avgArr(allY) mE = 0 # Default for most cases else: - mY, mE = weightedAvgArr(allY, allE) # Outputs masked values as zeros + mY, mE = weightedAvgArr(allY, allE) # Outputs masked values as zeros # DataY and DataE should never reach NaN, but safeguard in case they do - if (mE==np.nan) | (mY==np.nan): + if (mE == np.nan) | (mY == np.nan): mY, mE = 0, 0 meansY[i] = mY @@ -283,12 +319,18 @@ def weightedAvgXBinsArr(dataX, dataY, dataE, xp): def weightedAvgCols(wsYSpace): """Returns ws with weighted avg of columns of input ws""" dataX, dataY, dataE = extractWS(wsYSpace) - if np.all(dataE==0): # Bootstrap case where errors are not used + if np.all(dataE == 0): # Bootstrap case where errors are not used meanY = avgArr(dataY) meanE = np.zeros(meanY.shape) else: meanY, meanE = weightedAvgArr(dataY, dataE) - wsYSpaceAvg = CreateWorkspace(DataX=dataX[0, :], DataY=meanY, DataE=meanE, NSpec=1, OutputWorkspace=wsYSpace.name()+"_WeightedAvg") + wsYSpaceAvg = CreateWorkspace( + DataX=dataX[0, :], + DataY=meanY, + DataE=meanE, + NSpec=1, + OutputWorkspace=wsYSpace.name() + "_WeightedAvg", + ) return wsYSpaceAvg @@ -298,14 +340,16 @@ def avgArr(dataYO): Ignores any zero values as being masked. """ - assert len(dataYO)>1, "Averaging needs more than one element." + assert len(dataYO) > 1, "Averaging needs more than one element." dataY = dataYO.copy() - dataY[dataY==0] = np.nan + dataY[dataY == 0] = np.nan meanY = np.nanmean(dataY, axis=0) - meanY[meanY==np.nan] = 0 + meanY[meanY == np.nan] = 0 - assert np.all(np.all(dataYO==0, axis=0)==(meanY==0)), "Columns of zeros should give zero." + assert np.all( + np.all(dataYO == 0, axis=0) == (meanY == 0) + ), "Columns of zeros should give zero." return meanY @@ -316,31 +360,45 @@ def weightedAvgArr(dataYO, dataEO): """ # Run some tests - assert dataYO.shape==dataEO.shape, "Y and E arrays should have same shape for weighted average." - assert np.all((dataYO==0)==(dataEO==0)), f"Masked zeros should match in DataY and DataE: {np.argwhere((dataYO==0)!=(dataEO==0))}" - assert np.all(np.isnan(dataYO)==np.isnan(dataEO)), "Masked nans should match in DataY and DataE." - assert len(dataYO) > 1, "Weighted average needs more than one element to be performed." + assert ( + dataYO.shape == dataEO.shape + ), "Y and E arrays should have same shape for weighted average." + assert np.all( + (dataYO == 0) == (dataEO == 0) + ), f"Masked zeros should match in DataY and DataE: {np.argwhere((dataYO==0)!=(dataEO==0))}" + assert np.all( + np.isnan(dataYO) == np.isnan(dataEO) + ), "Masked nans should match in DataY and DataE." + assert ( + len(dataYO) > 1 + ), "Weighted average needs more than one element to be performed." dataY = dataYO.copy() # Copy arrays not to change original data dataE = dataEO.copy() # Ignore invalid data by changing zeros to nans # If data is already masked with nans, it remains unaltered - zerosMask = (dataY==0) + zerosMask = dataY == 0 dataY[zerosMask] = np.nan dataE[zerosMask] = np.nan - meanY = np.nansum(dataY/np.square(dataE), axis=0) / np.nansum(1/np.square(dataE), axis=0) - meanE = np.sqrt(1 / np.nansum(1/np.square(dataE), axis=0)) + meanY = np.nansum(dataY / np.square(dataE), axis=0) / np.nansum( + 1 / np.square(dataE), axis=0 + ) + meanE = np.sqrt(1 / np.nansum(1 / np.square(dataE), axis=0)) # Change invalid data back to original masking format with zeros - nanInfMask = (meanE==np.inf) | (meanE==np.nan) | (meanY==np.nan) + nanInfMask = (meanE == np.inf) | (meanE == np.nan) | (meanY == np.nan) meanY[nanInfMask] = 0 meanE[nanInfMask] = 0 # Test that columns of zeros are left unchanged - assert np.all((meanY==0)==(meanE==0)), "Weighted avg output should have masks in the same DataY and DataE." - assert np.all((np.all(dataYO==0, axis=0) | np.all(np.isnan(dataYO), axis=0)) == (meanY==0)), "Masked cols should be ignored." + assert np.all( + (meanY == 0) == (meanE == 0) + ), "Weighted avg output should have masks in the same DataY and DataE." + assert np.all( + (np.all(dataYO == 0, axis=0) | np.all(np.isnan(dataYO), axis=0)) == (meanY == 0) + ), "Masked cols should be ignored." return meanY, meanE @@ -348,7 +406,7 @@ def weightedAvgArr(dataYO, dataEO): def normalise_workspace(ws_name): """Updates workspace with the normalised version.""" tmp_norm = Integration(ws_name) - Divide(LHSWorkspace=ws_name,RHSWorkspace=tmp_norm,OutputWorkspace=ws_name) + Divide(LHSWorkspace=ws_name, RHSWorkspace=tmp_norm, OutputWorkspace=ws_name) DeleteWorkspace("tmp_norm") @@ -374,13 +432,13 @@ def symmetrizeWs(avgYSpace): dataX, dataY, dataE = extractWS(avgYSpace) - if np.all(dataE==0): + if np.all(dataE == 0): dataYS = symArr(dataY) dataES = np.zeros(dataYS.shape) else: dataYS, dataES = weightedSymArr(dataY, dataE) - wsSym = CloneWorkspace(avgYSpace, OutputWorkspace=avgYSpace.name()+"_Symmetrised") + wsSym = CloneWorkspace(avgYSpace, OutputWorkspace=avgYSpace.name() + "_Symmetrised") wsSym = passDataIntoWS(dataX, dataYS, dataES, wsSym) return wsSym @@ -394,15 +452,19 @@ def symArr(dataYO): assert len(dataYO.shape) == 2, "Symmetrization is written for 2D arrays." dataY = dataYO.copy() # Copy arrays not to risk changing original data - coMask = dataY==0 + coMask = dataY == 0 dataY[coMask] = np.nan yFlip = np.flip(dataY, axis=1) - dataYS = np.nanmean(np.stack((dataY, yFlip)), axis=0) # Normal avg between two numbers, cut-offs get ignored + dataYS = np.nanmean( + np.stack((dataY, yFlip)), axis=0 + ) # Normal avg between two numbers, cut-offs get ignored - dataYS[dataYS==np.nan] = 0 - np.testing.assert_array_equal(dataYS, np.flip(dataYS, axis=1)), f"Symmetrisation failed in {np.argwhere(dataYS!=np.flip(dataYS))}" + dataYS[dataYS == np.nan] = 0 + np.testing.assert_array_equal( + dataYS, np.flip(dataYS, axis=1) + ), f"Symmetrisation failed in {np.argwhere(dataYS!=np.flip(dataYS))}" np.testing.assert_allclose(dataYS[coMask], np.flip(dataYO, axis=1)[coMask]) return dataYS @@ -414,12 +476,14 @@ def weightedSymArr(dataYO, dataEO): the final value will be the valid point. """ assert len(dataYO.shape) == 2, "Symmetrization is written for 2D arrays." - assert np.all((dataYO==0)==(dataEO==0)), "Masked values should have zeros on both dataY and dataE." + assert np.all( + (dataYO == 0) == (dataEO == 0) + ), "Masked values should have zeros on both dataY and dataE." dataY = dataYO.copy() # Copy arrays not to risk changing original data dataE = dataEO.copy() - cutOffMask = dataY==0 + cutOffMask = dataY == 0 # Change values of yerr to leave cut-offs unchanged during symmetrisation dataE[cutOffMask] = np.full(np.sum(cutOffMask), np.inf) @@ -427,17 +491,23 @@ def weightedSymArr(dataYO, dataEO): eFlip = np.flip(dataE, axis=1) # Inverse variance weighting - dataYS = (dataY/dataE**2 + yFlip/eFlip**2) / (1/dataE**2 + 1/eFlip**2) - dataES = 1 / np.sqrt(1/dataE**2 + 1/eFlip**2) + dataYS = (dataY / dataE**2 + yFlip / eFlip**2) / ( + 1 / dataE**2 + 1 / eFlip**2 + ) + dataES = 1 / np.sqrt(1 / dataE**2 + 1 / eFlip**2) # Deal with effects from previously changing dataE=np.inf - nanInfMask = (dataES==np.inf) | (dataES==np.nan) | (dataYS==np.nan) + nanInfMask = (dataES == np.inf) | (dataES == np.nan) | (dataYS == np.nan) dataYS[nanInfMask] = 0 dataES[nanInfMask] = 0 # Test that arrays are symmetrised - np.testing.assert_array_equal(dataYS, np.flip(dataYS, axis=1)), f"Symmetrisation failed in {np.argwhere(dataYS!=np.flip(dataYS))}" - np.testing.assert_array_equal(dataES, np.flip(dataES, axis=1)), f"Symmetrisation failed in {np.argwhere(dataES!=np.flip(dataES))}" + np.testing.assert_array_equal( + dataYS, np.flip(dataYS, axis=1) + ), f"Symmetrisation failed in {np.argwhere(dataYS!=np.flip(dataYS))}" + np.testing.assert_array_equal( + dataES, np.flip(dataES, axis=1) + ), f"Symmetrisation failed in {np.argwhere(dataES!=np.flip(dataES))}" # Test that cut-offs were not included in the symmetrisation np.testing.assert_allclose(dataYS[cutOffMask], np.flip(dataYO, axis=1)[cutOffMask]) @@ -447,10 +517,9 @@ def weightedSymArr(dataYO, dataEO): def fitProfileMinuit(yFitIC, wsYSpaceSym, wsRes): - dataX, dataY, dataE = extractFirstSpectra(wsYSpaceSym) resX, resY, resE = extractFirstSpectra(wsRes) - assert np.all(dataX==resX), "Resolution should operate on the same range as DataX" + assert np.all(dataX == resX), "Resolution should operate on the same range as DataX" model, defaultPars, sharedPars = selectModelAndPars(yFitIC.fitModel) @@ -459,17 +528,17 @@ def fitProfileMinuit(yFitIC, wsYSpaceSym, wsRes): def convolvedModel(x, y0, *pars): return y0 + signal.convolve(model(x, *pars), resDense, mode="same") * xDelta - signature = describe(model)[:] # Build signature of convolved function - signature[1:1] = ["y0"] # Add intercept as first fitting parameter after range 'x' + signature = describe(model)[:] # Build signature of convolved function + signature[1:1] = ["y0"] # Add intercept as first fitting parameter after range 'x' convolvedModel.func_code = make_func_code(signature) - defaultPars["y0"] = 0 # Add initialization of parameter to dictionary + defaultPars["y0"] = 0 # Add initialization of parameter to dictionary # Fit only valid values, ignore cut-offs dataXNZ, dataYNZ, dataENZ = selectNonZeros(dataX, dataY, dataE) # Fit with Minuit - if np.all(dataE==0): # Choose fitting without weights + if np.all(dataE == 0): # Choose fitting without weights costFun = MyLeastSquares(dataXNZ, dataYNZ, convolvedModel) else: costFun = cost.LeastSquares(dataXNZ, dataYNZ, dataENZ, convolvedModel) @@ -477,19 +546,23 @@ def convolvedModel(x, y0, *pars): m = Minuit(costFun, **defaultPars) m.limits["A"] = (0, None) - if yFitIC.fitModel=="DOUBLE_WELL": + if yFitIC.fitModel == "DOUBLE_WELL": m.limits["d"] = (0, None) m.limits["R"] = (0, None) - if yFitIC.fitModel=="SINGLE_GAUSSIAN": + if yFitIC.fitModel == "SINGLE_GAUSSIAN": m.simplex() m.migrad() - def constrFunc()->None: # No constraint function for gaussian profile + def constrFunc() -> None: # No constraint function for gaussian profile return + else: - def constrFunc(*pars): # Constrain physical model before convolution - return model(dataXNZ, *pars[1:]) # First parameter is intercept, not part of model() + + def constrFunc(*pars): # Constrain physical model before convolution + return model( + dataXNZ, *pars[1:] + ) # First parameter is intercept, not part of model() m.simplex() m.scipy(constraints=optimize.NonlinearConstraint(constrFunc, 0, np.inf)) @@ -498,17 +571,21 @@ def constrFunc(*pars): # Constrain physical model before convolution m.hesse() # Weighted Chi2 - chi2 = m.fval / (len(dataXNZ)-m.nfit) + chi2 = m.fval / (len(dataXNZ) - m.nfit) # Best fit and confidence band # Calculated for the whole range of dataX, including where zero - dataYFit, dataYCov = jacobi.propagate(lambda pars: convolvedModel(dataX, *pars), m.values, m.covariance) + dataYFit, dataYCov = jacobi.propagate( + lambda pars: convolvedModel(dataX, *pars), m.values, m.covariance + ) dataYSigma = np.sqrt(np.diag(dataYCov)) - dataYSigma *= chi2 # Weight the confidence band + dataYSigma *= chi2 # Weight the confidence band Residuals = dataY - dataYFit # Create workspace to store best fit curve and errors on the fit - wsMinFit = createFitResultsWorkspace(wsYSpaceSym, dataX, dataY, dataE, dataYFit, dataYSigma, Residuals) + wsMinFit = createFitResultsWorkspace( + wsYSpaceSym, dataX, dataY, dataE, dataYFit, dataYSigma, Residuals + ) saveMinuitPlot(yFitIC, wsMinFit, m) # Calculate correlation matrix @@ -542,52 +619,117 @@ def selectModelAndPars(modelFlag): """ if modelFlag == "SINGLE_GAUSSIAN": + def model(x, A, x0, sigma): - return A / (2*np.pi)**0.5 / sigma * np.exp(-(x-x0)**2/2/sigma**2) + return ( + A + / (2 * np.pi) ** 0.5 + / sigma + * np.exp(-((x - x0) ** 2) / 2 / sigma**2) + ) - defaultPars = {"A":1, "x0":0, "sigma":5} - sharedPars = ["sigma"] # Used only in Global fit + defaultPars = {"A": 1, "x0": 0, "sigma": 5} + sharedPars = ["sigma"] # Used only in Global fit + + elif modelFlag == "GC_C4_C6": - elif (modelFlag=="GC_C4_C6"): def model(x, A, x0, sigma1, c4, c6): - return A * np.exp(-(x-x0)**2/2/sigma1**2) / (np.sqrt(2*np.pi*sigma1**2)) \ - *(1 + c4/32*(16*((x-x0)/np.sqrt(2)/sigma1)**4 - -48*((x-x0)/np.sqrt(2)/sigma1)**2+12) - +c6/384*(64*((x-x0)/np.sqrt(2)/sigma1)**6 - -480*((x-x0)/np.sqrt(2)/sigma1)**4 + 720*((x-x0)/np.sqrt(2)/sigma1)**2 - 120)) + return ( + A + * np.exp(-((x - x0) ** 2) / 2 / sigma1**2) + / (np.sqrt(2 * np.pi * sigma1**2)) + * ( + 1 + + c4 + / 32 + * ( + 16 * ((x - x0) / np.sqrt(2) / sigma1) ** 4 + - 48 * ((x - x0) / np.sqrt(2) / sigma1) ** 2 + + 12 + ) + + c6 + / 384 + * ( + 64 * ((x - x0) / np.sqrt(2) / sigma1) ** 6 + - 480 * ((x - x0) / np.sqrt(2) / sigma1) ** 4 + + 720 * ((x - x0) / np.sqrt(2) / sigma1) ** 2 + - 120 + ) + ) + ) - defaultPars = {"A":1, "x0":0, "sigma1":6, "c4":0, "c6":0} - sharedPars = ["sigma1", "c4", "c6"] # Used only in Global fit + 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=="GC_C4": def model(x, A, x0, sigma1, c4): - return A * np.exp(-(x-x0)**2/2/sigma1**2) / (np.sqrt(2*np.pi*sigma1**2)) \ - *(1 + c4/32*(16*((x-x0)/np.sqrt(2)/sigma1)**4 - -48*((x-x0)/np.sqrt(2)/sigma1)**2+12)) + return ( + A + * np.exp(-((x - x0) ** 2) / 2 / sigma1**2) + / (np.sqrt(2 * np.pi * sigma1**2)) + * ( + 1 + + c4 + / 32 + * ( + 16 * ((x - x0) / np.sqrt(2) / sigma1) ** 4 + - 48 * ((x - x0) / np.sqrt(2) / sigma1) ** 2 + + 12 + ) + ) + ) - defaultPars = {"A":1, "x0":0, "sigma1":6, "c4":0} - sharedPars = ["sigma1", "c4"] # Used only in Global fit + defaultPars = {"A": 1, "x0": 0, "sigma1": 6, "c4": 0} + sharedPars = ["sigma1", "c4"] # Used only in Global fit + + elif modelFlag == "GC_C6": - elif modelFlag=="GC_C6": def model(x, A, x0, sigma1, c6): - return A * np.exp(-(x-x0)**2/2/sigma1**2) / (np.sqrt(2*np.pi*sigma1**2)) \ - *(1 + +c6/384*(64*((x-x0)/np.sqrt(2)/sigma1)**6 - -480*((x-x0)/np.sqrt(2)/sigma1)**4 + 720*((x-x0)/np.sqrt(2)/sigma1)**2 - 120)) + return ( + A + * np.exp(-((x - x0) ** 2) / 2 / sigma1**2) + / (np.sqrt(2 * np.pi * sigma1**2)) + * ( + 1 + + +c6 + / 384 + * ( + 64 * ((x - x0) / np.sqrt(2) / sigma1) ** 6 + - 480 * ((x - x0) / np.sqrt(2) / sigma1) ** 4 + + 720 * ((x - x0) / np.sqrt(2) / sigma1) ** 2 + - 120 + ) + ) + ) - defaultPars = {"A":1, "x0":0, "sigma1":6, "c6":0} - sharedPars = ["sigma1", "c6"] # Used only in Global fit + defaultPars = {"A": 1, "x0": 0, "sigma1": 6, "c6": 0} + sharedPars = ["sigma1", "c6"] # Used only in Global fit + + elif modelFlag == "DOUBLE_WELL": - elif modelFlag=="DOUBLE_WELL": def model(x, A, d, R, sig1, sig2): # h = 2.04 - theta = np.linspace(0, np.pi, 300)[:, np.newaxis] # 300 points seem like a good estimate for ~10 examples + theta = np.linspace(0, np.pi, 300)[ + :, np.newaxis + ] # 300 points seem like a good estimate for ~10 examples y = x[np.newaxis, :] - sigTH = np.sqrt( sig1**2*np.cos(theta)**2 + sig2**2*np.sin(theta)**2 ) - alpha = 2*( d*sig2*sig1*np.sin(theta) / sigTH )**2 - beta = ( 2*sig1**2*d*np.cos(theta) / sigTH**2 ) * y - denom = 2.506628 * sigTH * (1 + R**2 + 2*R*np.exp(-2*d**2*sig1**2)) - jp = np.exp( -y**2/(2*sigTH**2)) * (1 + R**2 + 2*R*np.exp(-alpha)*np.cos(beta)) / denom + sigTH = np.sqrt( + sig1**2 * np.cos(theta) ** 2 + sig2**2 * np.sin(theta) ** 2 + ) + alpha = 2 * (d * sig2 * sig1 * np.sin(theta) / sigTH) ** 2 + beta = (2 * sig1**2 * d * np.cos(theta) / sigTH**2) * y + denom = ( + 2.506628 + * sigTH + * (1 + R**2 + 2 * R * np.exp(-2 * d**2 * sig1**2)) + ) + jp = ( + np.exp(-(y**2) / (2 * sigTH**2)) + * (1 + R**2 + 2 * R * np.exp(-alpha) * np.cos(beta)) + / denom + ) jp *= np.sin(theta) JBest = np.trapz(jp, x=theta, axis=0) @@ -595,18 +737,26 @@ def model(x, A, d, R, sig1, sig2): JBest *= A return JBest - defaultPars = {"A":1, "d":1, "R":1, "sig1":3, "sig2":5} # TODO: Starting parameters and bounds? - sharedPars = ["d", "R", "sig1", "sig2"] # Only varying parameter is amplitude A + defaultPars = { + "A": 1, + "d": 1, + "R": 1, + "sig1": 3, + "sig2": 5, + } # TODO: Starting parameters and bounds? + sharedPars = ["d", "R", "sig1", "sig2"] # Only varying parameter is amplitude A - elif modelFlag=="ANSIO_GAUSSIAN": + elif modelFlag == "ANSIO_GAUSSIAN": # Ansiotropic case def model(x, A, sig1, sig2): # h = 2.04 theta = np.linspace(0, np.pi, 300)[:, np.newaxis] y = x[np.newaxis, :] - sigTH = np.sqrt( sig1**2*np.cos(theta)**2 + sig2**2*np.sin(theta)**2 ) - jp = np.exp( -y**2/(2*sigTH**2)) / (2.506628*sigTH) + sigTH = np.sqrt( + sig1**2 * np.cos(theta) ** 2 + sig2**2 * np.sin(theta) ** 2 + ) + jp = np.exp(-(y**2) / (2 * sigTH**2)) / (2.506628 * sigTH) jp *= np.sin(theta) JBest = np.trapz(jp, x=theta, axis=0) @@ -614,28 +764,31 @@ def model(x, A, sig1, sig2): JBest *= A return JBest - defaultPars = {"A":1, "sig1":3, "sig2":5} + defaultPars = {"A": 1, "sig1": 3, "sig2": 5} sharedPars = ["sig1", "sig2"] - elif modelFlag=="Gaussian3D": - def model(x, A, sig_x, sig_y, sig_z): + elif modelFlag == "Gaussian3D": + def model(x, A, sig_x, sig_y, sig_z): y = x[:, np.newaxis, np.newaxis] - n_steps = 50 # Low number of integration steps because otherwise too slow + n_steps = 50 # Low number of integration steps because otherwise too slow theta = np.linspace(0, np.pi / 2, n_steps)[np.newaxis, :, np.newaxis] phi = np.linspace(0, np.pi / 2, n_steps)[np.newaxis, np.newaxis, :] + S2_inv = ( + np.sin(theta) ** 2 * np.cos(phi) ** 2 / sig_x**2 + + np.sin(theta) ** 2 * np.sin(phi) ** 2 / sig_y**2 + + np.cos(theta) ** 2 / sig_z**2 + ) - S2_inv = np.sin(theta)**2 * np.cos(phi)**2 / sig_x**2 \ - + np.sin(theta)**2 * np.sin(phi)**2 / sig_y**2 \ - + np.cos(theta)**2 / sig_z**2 - - J = np.sin(theta) / S2_inv * np.exp(- y**2 / 2 * S2_inv) + J = np.sin(theta) / S2_inv * np.exp(-(y**2) / 2 * S2_inv) - J = np.trapz(J, x=phi, axis=2)[:, :, np.newaxis] # Keep shape + J = np.trapz(J, x=phi, axis=2)[:, :, np.newaxis] # Keep shape J = np.trapz(J, x=theta, axis=1) - J *= A * 2 / np.pi * 1 / np.sqrt(2 * np.pi) * 1 / (sig_x * sig_y * sig_z) # Normalisation + J *= ( + A * 2 / np.pi * 1 / np.sqrt(2 * np.pi) * 1 / (sig_x * sig_y * sig_z) + ) # Normalisation J = J.squeeze() return J @@ -643,15 +796,22 @@ def model(x, A, sig_x, sig_y, sig_z): sharedPars = ["sig_x", "sig_y", "sig_z"] else: - raise ValueError("Fitting Model not recognized, available options: 'SINGLE_GAUSSIAN', 'GC_C4_C6', 'GC_C4'") + raise ValueError( + "Fitting Model not recognized, available options: 'SINGLE_GAUSSIAN', 'GC_C4_C6', 'GC_C4'" + ) print("\nShared Parameters: ", [key for key in sharedPars]) - print("\nUnshared Parameters: ", [key for key in defaultPars if key not in sharedPars]) + print( + "\nUnshared Parameters: ", [key for key in defaultPars if key not in sharedPars] + ) - assert all(isinstance(item, str) for item in sharedPars), "Parameters in list must be strings." - assert describe(model)[-len(sharedPars) - :]==sharedPars, "Function signature needs to have shared parameters at the end: " \ - "model(*unsharedPars, *sharedPars)" + assert all( + isinstance(item, str) for item in sharedPars + ), "Parameters in list must be strings." + assert describe(model)[-len(sharedPars) :] == sharedPars, ( + "Function signature needs to have shared parameters at the end: " + "model(*unsharedPars, *sharedPars)" + ) return model, defaultPars, sharedPars @@ -661,7 +821,7 @@ def selectNonZeros(dataX, dataY, dataE): Selects non zero points. Uses zeros in dataY becasue dataE can be all zeros in one of the bootstrap types. """ - zeroMask = dataY==0 + zeroMask = dataY == 0 dataXNZ = dataX[~zeroMask] dataYNZ = dataY[~zeroMask] @@ -675,7 +835,7 @@ class MyLeastSquares: This structure is required for high compatibility with Minuit. """ - errordef = Minuit.LEAST_SQUARES # for Minuit to compute errors correctly + errordef = Minuit.LEAST_SQUARES # for Minuit to compute errors correctly def __init__(self, x, y, model): self.model = model # model predicts y for given x @@ -692,14 +852,18 @@ def ndata(self): return len(self.x) -def createFitResultsWorkspace(wsYSpaceSym, dataX, dataY, dataE, dataYFit, dataYSigma, Residuals): +def createFitResultsWorkspace( + wsYSpaceSym, dataX, dataY, dataE, dataYFit, dataYSigma, Residuals +): """Creates workspace similar to the ones created by Mantid Fit.""" - wsMinFit = 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()+"_Fitted_Minuit") + wsMinFit = 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() + "_Fitted_Minuit", + ) return wsMinFit @@ -710,7 +874,7 @@ def saveMinuitPlot(yFitIC, wsMinuitFit, mObj): for p, v, e in zip(mObj.parameters, mObj.values, mObj.errors): leg += f"${p}={v:.2f} \\pm {e:.2f}$\n" - fig, ax = plt.subplots(subplot_kw={"projection":"mantid"}) + 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.set_xlabel("YSpace") @@ -718,7 +882,7 @@ def saveMinuitPlot(yFitIC, wsMinuitFit, mObj): ax.set_title("Minuit Fit") ax.legend() - fileName = wsMinuitFit.name()+".pdf" + fileName = wsMinuitFit.name() + ".pdf" savePath = yFitIC.figSavePath / fileName plt.savefig(savePath, bbox_inches="tight") plt.close(fig) @@ -726,11 +890,13 @@ def saveMinuitPlot(yFitIC, wsMinuitFit, mObj): def createCorrelationTableWorkspace(wsYSpaceSym, parameters, corrMatrix): - tableWS = CreateEmptyTableWorkspace(OutputWorkspace=wsYSpaceSym.name()+"_Fitted_Minuit_NormalizedCovarianceMatrix") + tableWS = CreateEmptyTableWorkspace( + OutputWorkspace=wsYSpaceSym.name() + "_Fitted_Minuit_NormalizedCovarianceMatrix" + ) tableWS.setTitle("Minuit Fit") - tableWS.addColumn(type='str',name="Name") + tableWS.addColumn(type="str", name="Name") for p in parameters: - tableWS.addColumn(type='float',name=p) + tableWS.addColumn(type="float", name=p) for p, arr in zip(parameters, corrMatrix): tableWS.addRow([p] + list(arr)) @@ -744,7 +910,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.runMinos): minosAutoErr = list(np.zeros((len(parameters), 2))) minosManErr = list(np.zeros((len(parameters), 2))) return parameters, values, errors, minosAutoErr, minosManErr @@ -755,7 +921,9 @@ def runMinos(mObj, yFitIC, constrFunc, wsName): bestFitVals[p] = v bestFitErrs[p] = e - if (yFitIC.fitModel=="SINGLE_GAUSSIAN"): # Case with no positivity constraint, can use automatic minos() + if ( + yFitIC.fitModel == "SINGLE_GAUSSIAN" + ): # Case with no positivity constraint, can use automatic minos() mObj.minos() me = mObj.merrors @@ -768,9 +936,11 @@ def runMinos(mObj, yFitIC, constrFunc, wsName): if yFitIC.showPlots: plotAutoMinos(mObj, wsName) - else: # Case with positivity constraint on function, use manual implementation + 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) + merrors, fig = runAndPlotManualMinos( + mObj, constrFunc, bestFitVals, bestFitErrs, yFitIC.showPlots + ) # Same as above, but the other way around minosManErr = [] @@ -779,10 +949,10 @@ def runMinos(mObj, yFitIC, constrFunc, wsName): minosAutoErr = list(np.zeros(np.array(minosManErr).shape)) if yFitIC.showPlots: - fig.canvas.manager.set_window_title(wsName+"_Manual_Implementation_MINOS") + fig.canvas.manager.set_window_title(wsName + "_Manual_Implementation_MINOS") fig.show() - return parameters, values, errors, minosAutoErr, minosManErr + return parameters, values, errors, minosAutoErr, minosManErr def runAndPlotManualMinos(minuitObj, constrFunc, bestFitVals, bestFitErrs, showPlots): @@ -796,33 +966,41 @@ def runAndPlotManualMinos(minuitObj, constrFunc, bestFitVals, bestFitErrs, showP # Set format of subplots height = 2 - width = int(np.ceil(len(minuitObj.parameters)/2)) + width = int(np.ceil(len(minuitObj.parameters) / 2)) figsize = (12, 7) # Output plot to Mantid - fig, axs = plt.subplots(height, width, tight_layout=True, figsize=figsize, subplot_kw={ - 'projection':'mantid'}) #subplot_kw={'projection':'mantid'} + fig, axs = plt.subplots( + height, + width, + tight_layout=True, + figsize=figsize, + subplot_kw={"projection": "mantid"}, + ) # subplot_kw={'projection':'mantid'} # fig.canvas.manager.set_window_title("Plot of Manual Implementation MINOS") merrors = {} for p, ax in zip(minuitObj.parameters, axs.flat): - lerr, uerr = runMinosForPar(minuitObj, constrFunc, p, 2, ax, bestFitVals, bestFitErrs, showPlots) + lerr, uerr = runMinosForPar( + minuitObj, constrFunc, p, 2, ax, bestFitVals, bestFitErrs, showPlots + ) 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 + if not ax.lines: # If empty list ax.set_visible(False) # 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() + fig.legend(handle, label, loc="lower right") + # fig.show() return merrors, fig -def runMinosForPar(minuitObj, constrFunc, var:str, bound:int, ax, bestFitVals, bestFitErrs, showPlots): - +def runMinosForPar( + minuitObj, constrFunc, var: str, bound: int, ax, bestFitVals, bestFitErrs, showPlots +): resetMinuit(minuitObj, bestFitVals, bestFitErrs) # Run Fitting procedures again to be on the safe side and reset to minimum minuitObj.scipy(constraints=optimize.NonlinearConstraint(constrFunc, 0, np.inf)) @@ -832,13 +1010,13 @@ def runMinosForPar(minuitObj, constrFunc, var:str, bound:int, ax, bestFitVals, b varVal = minuitObj.values[var] varErr = minuitObj.errors[var] # Store fval of best fit - fValsMin = minuitObj.fval # Used to calculate error bands at the end + fValsMin = minuitObj.fval # Used to calculate error bands at the end varSpace = buildVarRange(bound, varVal, varErr) # Split variable space into right and left side lhsVarSpace, rhsVarSpace = np.split(varSpace, 2) - lhsVarSpace = np.flip(lhsVarSpace) # Flip to start at minimum + lhsVarSpace = np.flip(lhsVarSpace) # Flip to start at minimum for minimizer in ("Scipy", "Migrad"): resetMinuit(minuitObj, bestFitVals, bestFitErrs) @@ -847,14 +1025,22 @@ def runMinosForPar(minuitObj, constrFunc, var:str, bound:int, ax, bestFitVals, b resetMinuit(minuitObj, bestFitVals, bestFitErrs) lhsMinos = runMinosOnRange(minuitObj, var, lhsVarSpace, minimizer, constrFunc) - wholeMinos = np.concatenate((np.flip(lhsMinos), rhsMinos), axis=None) # Flip left hand side again + wholeMinos = np.concatenate( + (np.flip(lhsMinos), rhsMinos), axis=None + ) # Flip left hand side again - if minimizer == "Scipy": # Calculate minos errors from constrained scipy - lerr, uerr = errsFromMinosCurve(varSpace, varVal, wholeMinos, fValsMin, dChi2=1) + if minimizer == "Scipy": # Calculate minos errors from constrained scipy + lerr, uerr = errsFromMinosCurve( + varSpace, varVal, wholeMinos, fValsMin, dChi2=1 + ) ax.plot(varSpace, wholeMinos, label="fVals Constr Scipy") - elif minimizer == "Migrad": # Plot migrad as well to see the difference between constrained and unconstrained - plotProfile(ax, var, varSpace, wholeMinos, lerr, uerr, fValsMin, varVal, varErr) + elif ( + minimizer == "Migrad" + ): # Plot migrad as well to see the difference between constrained and unconstrained + plotProfile( + ax, var, varSpace, wholeMinos, lerr, uerr, fValsMin, varVal, varErr + ) else: raise ValueError("Minimizer not recognized.") @@ -873,29 +1059,29 @@ def resetMinuit(minuitObj, bestFitVals, bestFitErrs): def buildVarRange(bound, varVal, varErr): """Range of points over which cost function is evaluated.""" # Create variable space more dense near the minima using a quadratic density - limit = (bound*varErr)**(1/2) # Square root is corrected below + limit = (bound * varErr) ** (1 / 2) # Square root is corrected below varSpace = np.linspace(-limit, limit, 30) varSpace = varSpace**2 * np.sign(varSpace) + varVal - assert len(varSpace)%2 == 0, "Number of points in Minos range needs to be even" + assert len(varSpace) % 2 == 0, "Number of points in Minos range needs to be even" return varSpace def runMinosOnRange(minuitObj, var, varRange, minimizer, constrFunc): - result = np.zeros(varRange.size) minuitObj.fixed[var] = True # Unconstrained fit over side range for i, value in enumerate(varRange): - - minuitObj.values[var] = value # Fix variable + minuitObj.values[var] = value # Fix variable if minimizer == "Migrad": - minuitObj.migrad() # Fit + minuitObj.migrad() # Fit elif minimizer == "Scipy": - minuitObj.scipy(constraints=optimize.NonlinearConstraint(constrFunc, 0, np.inf)) + minuitObj.scipy( + constraints=optimize.NonlinearConstraint(constrFunc, 0, np.inf) + ) - result[i] = minuitObj.fval # Store minimum + result[i] = minuitObj.fval # Store minimum minuitObj.fixed[var] = False return result @@ -908,12 +1094,14 @@ def errsFromMinosCurve(varSpace, varVal, fValsScipy, fValsMin, dChi2=1): # Calculate points of intersection with line delta fmin val = 1 idxErr = np.argwhere(np.diff(np.sign(fValsScipyDense - fValsMin - 1))) - if idxErr.size != 2: # Intersections not found, do not plot error range - lerr, uerr = 0., 0. + if idxErr.size != 2: # Intersections not found, do not plot error range + lerr, uerr = 0.0, 0.0 else: lerr, uerr = varSpaceDense[idxErr].flatten() - varVal - if lerr*uerr >= 0: # Case where we get either two positive or two negative errors, ill defined profile + if ( + lerr * uerr >= 0 + ): # Case where we get either two positive or two negative errors, ill defined profile lerr, uerr = 0, 0 return lerr, uerr @@ -922,11 +1110,17 @@ def errsFromMinosCurve(varSpace, varVal, fValsScipy, fValsMin, dChi2=1): def plotAutoMinos(minuitObj, wsName): # Set format of subplots height = 2 - width = int(np.ceil(len(minuitObj.parameters)/2)) + width = int(np.ceil(len(minuitObj.parameters) / 2)) figsize = (12, 7) # Output plot to Mantid - fig, axs = plt.subplots(height, width, tight_layout=True, figsize=figsize, subplot_kw={'projection':'mantid'}) - fig.canvas.manager.set_window_title(wsName+"_Plot_Automatic_MINOS") + fig, axs = plt.subplots( + height, + width, + tight_layout=True, + figsize=figsize, + subplot_kw={"projection": "mantid"}, + ) + fig.canvas.manager.set_window_title(wsName + "_Plot_Automatic_MINOS") for p, ax in zip(minuitObj.parameters, axs.flat): loc, fvals, status = minuitObj.mnprofile(p, bound=2) @@ -940,12 +1134,12 @@ def plotAutoMinos(minuitObj, wsName): # Hide plots not in use: for ax in axs.flat: - if not ax.lines: # If empty list + if not ax.lines: # If empty list ax.set_visible(False) # 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.legend(handle, label, loc="lower right") fig.show() @@ -956,31 +1150,45 @@ def plotProfile(ax, var, varSpace, fValsMigrad, lerr, uerr, fValsMin, varVal, va fValsMigrad : y axis """ - ax.set_title(var+f" = {varVal:.3f} {lerr:.3f} {uerr:+.3f}") + ax.set_title(var + f" = {varVal:.3f} {lerr:.3f} {uerr:+.3f}") ax.plot(varSpace, fValsMigrad, label="fVals Migrad") - ax.axvspan(lerr+varVal, uerr+varVal, alpha=0.2, color="red", label="Minos error") - ax.axvspan(varVal-varErr, varVal+varErr, alpha=0.2, color="green", label="Hessian Std error") + ax.axvspan( + lerr + varVal, uerr + varVal, alpha=0.2, color="red", label="Minos error" + ) + ax.axvspan( + varVal - varErr, + varVal + varErr, + alpha=0.2, + color="green", + label="Hessian Std error", + ) ax.axvline(varVal, 0.03, 0.97, color="k", ls="--") - ax.axhline(fValsMin+1, 0.03, 0.97, color="k") + ax.axhline(fValsMin + 1, 0.03, 0.97, color="k") ax.axhline(fValsMin, 0.03, 0.97, color="k") -def createFitParametersTableWorkspace(wsYSpaceSym, parameters, values, errors, minosAutoErr, minosManualErr, chi2): +def createFitParametersTableWorkspace( + wsYSpaceSym, parameters, values, errors, minosAutoErr, minosManualErr, chi2 +): # Create Parameters workspace - tableWS = CreateEmptyTableWorkspace(OutputWorkspace=wsYSpaceSym.name()+"_Fitted_Minuit_Parameters") + tableWS = CreateEmptyTableWorkspace( + OutputWorkspace=wsYSpaceSym.name() + "_Fitted_Minuit_Parameters" + ) tableWS.setTitle("Minuit Fit") - tableWS.addColumn(type='str', name="Name") - tableWS.addColumn(type='float', name="Value") - tableWS.addColumn(type='float', name="Error") - tableWS.addColumn(type='float', name="Auto Minos Error-") - tableWS.addColumn(type='float', name="Auto Minos Error+") - tableWS.addColumn(type='float', name="Manual Minos Error-") - tableWS.addColumn(type='float', name="Manual Minos Error+") - - for p, v, e, mae, mme in zip(parameters, values, errors, minosAutoErr, minosManualErr): + tableWS.addColumn(type="str", name="Name") + tableWS.addColumn(type="float", name="Value") + tableWS.addColumn(type="float", name="Error") + tableWS.addColumn(type="float", name="Auto Minos Error-") + tableWS.addColumn(type="float", name="Auto Minos Error+") + tableWS.addColumn(type="float", name="Manual Minos Error-") + tableWS.addColumn(type="float", name="Manual Minos Error+") + + for p, v, e, mae, mme in zip( + parameters, values, errors, minosAutoErr, minosManualErr + ): tableWS.addRow([p, v, e, mae[0], mae[1], mme[0], mme[1]]) tableWS.addRow(["Cost function", chi2, 0, 0, 0, 0, 0]) @@ -996,11 +1204,13 @@ def oddPointsRes(x, res): assert x.size == res.size, "x and res need to be the same size!" if res.size % 2 == 0: - dens = res.size+1 # If even change to odd + dens = res.size + 1 # If even change to odd else: - dens = res.size # If odd, keep being odd + dens = res.size # If odd, keep being odd - xDense = np.linspace(np.min(x), np.max(x), dens) # Make gridd with odd number of points - peak at center + xDense = np.linspace( + np.min(x), np.max(x), dens + ) # Make gridd with odd number of points - peak at center xDelta = xDense[1] - xDense[0] resDense = np.interp(xDense, x, res) @@ -1009,16 +1219,15 @@ 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": - function=f"""composite=Convolution,FixResolution=true,NumDeriv=true; + print("\nFitting on the sum of spectra in the West domain ...\n") + for minimizer in ["Levenberg-Marquardt", "Simplex"]: + if yFitIC.fitModel == "SINGLE_GAUSSIAN": + 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.fitModel == "GC_C4_C6": function = f""" composite=Convolution,FixResolution=true,NumDeriv=true; name=Resolution,Workspace={wsRes.name()},WorkspaceIndex=0,X=(),Y=(); @@ -1027,7 +1236,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=(020s}: "+" ".join([f"{elem:7.8s}" for elem in tableWS.column(key)])) + if key == "Name": + print( + f"{key:>20s}: " + + " ".join([f"{elem:7.8s}" for elem in tableWS.column(key)]) + ) else: - print(f"{key:>20s}: "+" ".join([f"{elem:7.4f}" for elem in tableWS.column(key)])) + print( + f"{key:>20s}: " + + " ".join([f"{elem:7.4f}" for elem in tableWS.column(key)]) + ) print("\n") class ResultsYFitObject: - def __init__(self, ic, yFitIC, wsFinalName, wsYSpaceAvgName): # Extract most relevant information from ws wsFinal = mtd[wsFinalName] @@ -1148,26 +1366,31 @@ def __init__(self, ic, yFitIC, wsFinalName, wsYSpaceAvgName): self.fitModel = yFitIC.fitModel def save(self): - np.savez(self.savePath, - YSpaceSymSumDataY=self.YSpaceSymSumDataY, - YSpaceSymSumDataE=self.YSpaceSymSumDataE, - resolution=self.resolution, - HdataY=self.HdataY, - finalRawDataY=self.finalRawDataY, - finalRawDataE=self.finalRawDataE, - popt=self.popt, - perr=self.perr) + np.savez( + self.savePath, + YSpaceSymSumDataY=self.YSpaceSymSumDataY, + YSpaceSymSumDataE=self.YSpaceSymSumDataE, + resolution=self.resolution, + HdataY=self.HdataY, + finalRawDataY=self.finalRawDataY, + finalRawDataE=self.finalRawDataE, + popt=self.popt, + perr=self.perr, + ) def runGlobalFit(wsYSpace, wsRes, IC, yFitIC): - print("\nRunning GLobal Fit ...\n") dataX, dataY, dataE, dataRes, instrPars = extractData(wsYSpace, wsRes, IC) - dataX, dataY, dataE, dataRes, instrPars = takeOutMaskedSpectra(dataX, dataY, dataE, dataRes, instrPars) + dataX, dataY, dataE, dataRes, instrPars = takeOutMaskedSpectra( + dataX, dataY, dataE, dataRes, instrPars + ) idxList = groupDetectors(instrPars, yFitIC) - dataX, dataY, dataE, dataRes = avgWeightDetGroups(dataX, dataY, dataE, dataRes, idxList, yFitIC) + dataX, dataY, dataE, dataRes = avgWeightDetGroups( + dataX, dataY, dataE, dataRes, idxList, yFitIC + ) if yFitIC.symmetrisationFlag: dataY, dataE = weightedSymArr(dataY, dataE) @@ -1178,10 +1401,11 @@ def runGlobalFit(wsYSpace, wsRes, IC, yFitIC): for i, (x, y, yerr, res) in enumerate(zip(dataX, dataY, dataE, dataRes)): totCost += calcCostFun(model, i, x, y, yerr, res, sharedPars) - defaultPars["y0"] = 0 # Introduce default parameter for convolved model + defaultPars["y0"] = 0 # Introduce default parameter for convolved model - assert len(describe(totCost)) == len(sharedPars) + len(dataY)*(len(defaultPars) - - len(sharedPars)), f"Wrong parameters for Global Fit:\n{describe(totCost)}" + assert len(describe(totCost)) == len(sharedPars) + len(dataY) * ( + len(defaultPars) - len(sharedPars) + ), f"Wrong parameters for Global Fit:\n{describe(totCost)}" # Minuit Fit with global cost function and local+global parameters initPars = minuitInitialParameters(defaultPars, sharedPars, len(dataY)) @@ -1189,22 +1413,22 @@ def runGlobalFit(wsYSpace, wsRes, IC, yFitIC): print("\nRunning Global Fit ...\n") m = Minuit(totCost, **initPars) - for i in range(len(dataY)): # Set limits for unshared parameters - m.limits["A"+str(i)] = (0, np.inf) + for i in range(len(dataY)): # Set limits for unshared parameters + m.limits["A" + str(i)] = (0, np.inf) - if yFitIC.fitModel=="DOUBLE_WELL": - m.limits["d"] = (0, np.inf) # Shared parameters + if yFitIC.fitModel == "DOUBLE_WELL": + m.limits["d"] = (0, np.inf) # Shared parameters m.limits["R"] = (0, np.inf) t0 = time.time() - if yFitIC.fitModel=="SINGLE_GAUSSIAN": + if yFitIC.fitModel == "SINGLE_GAUSSIAN": m.simplex() m.migrad() else: - totSig = describe(totCost) # This signature has 'x' already removed + totSig = describe(totCost) # This signature has 'x' already removed sharedIdxs = [totSig.index(shPar) for shPar in sharedPars] - nCostFunctions = len(totCost) # Number of individual cost functions + nCostFunctions = len(totCost) # Number of individual cost functions x = dataX[0] def constr(*pars): @@ -1215,14 +1439,19 @@ def constr(*pars): Builds array with all constraints from individual functions. """ - sharedPars = [pars[i] for i in sharedIdxs] # sigma1, c4, c6 in original GC + sharedPars = [pars[i] for i in sharedIdxs] # sigma1, c4, c6 in original GC unsharedPars = np.delete(pars, sharedIdxs, None) - unsharedParsSplit = np.split(unsharedPars, nCostFunctions) # Splits unshared parameters per individual cost fun + unsharedParsSplit = np.split( + unsharedPars, nCostFunctions + ) # Splits unshared parameters per individual cost fun joinedGC = np.zeros(nCostFunctions * x.size) for i, unshParsModel in enumerate( - unsharedParsSplit): # Attention to format of unshared and shared parameters when calling model - joinedGC[i*x.size : (i+1)*x.size] = model(x, *unshParsModel[1:], *sharedPars) # Intercept is first of unshared parameters + unsharedParsSplit + ): # Attention to format of unshared and shared parameters when calling model + joinedGC[i * x.size : (i + 1) * x.size] = model( + x, *unshParsModel[1:], *sharedPars + ) # Intercept is first of unshared parameters return joinedGC @@ -1235,7 +1464,9 @@ 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 + 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}") @@ -1247,7 +1478,9 @@ def constr(*pars): if yFitIC.showPlots: plotGlobalFit(dataX, dataY, dataE, m, totCost, wsYSpace.name()) - return np.array(m.values), np.array(m.errors) # Pass into array to store values in variable + return np.array(m.values), np.array( + m.errors + ) # Pass into array to store values in variable def extractData(ws, wsRes, ic): @@ -1256,7 +1489,9 @@ def extractData(ws, wsRes, ic): dataX = ws.extractX() dataRes = wsRes.extractY() instrPars = loadInstrParsFileIntoArray(ic) - assert len(instrPars) == len(dataY), "Load of IP file not working correctly, probable issue with indexing." + assert len(instrPars) == len( + dataY + ), "Load of IP file not working correctly, probable issue with indexing." return dataX, dataY, dataE, dataRes, instrPars @@ -1269,7 +1504,7 @@ def loadInstrParsFileIntoArray(ic): def takeOutMaskedSpectra(dataX, dataY, dataE, dataRes, instrPars): - zerosRowMask = np.all(dataY==0, axis=1) + zerosRowMask = np.all(dataY == 0, axis=1) dataY = dataY[~zerosRowMask] dataE = dataE[~zerosRowMask] dataX = dataX[~zerosRowMask] @@ -1277,6 +1512,7 @@ def takeOutMaskedSpectra(dataX, dataY, dataE, dataRes, instrPars): instrPars = instrPars[~zerosRowMask] return dataX, dataY, dataE, dataRes, instrPars + # ------- Groupings @@ -1298,22 +1534,24 @@ def groupDetectors(ipData, yFitIC): L1 /= np.sum(L1) theta /= np.sum(theta) - L1 *= 2 # Bigger weight to L1 + L1 *= 2 # Bigger weight to L1 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) - centers = points[startingIdxs, :] # Centers of cluster groups, NOT fitting parameter + startingIdxs = np.linspace(0, len(points) - 1, yFitIC.nGlobalFitGroups).astype(int) + centers = points[ + startingIdxs, : + ] # Centers of cluster groups, NOT fitting parameter - if False: # Set to True to investigate problems with groupings + if False: # Set to True to investigate problems with groupings plotDetsAndInitialCenters(L1, theta, centers) clusters = kMeansClustering(points, centers) idxList = formIdxList(clusters) if yFitIC.showPlots: - fig, ax = plt.subplots(tight_layout=True, subplot_kw={'projection':'mantid'}) + fig, ax = plt.subplots(tight_layout=True, subplot_kw={"projection": "mantid"}) fig.canvas.manager.set_window_title("Grouping of detectors") plotFinalGroups(ax, ipData, idxList) fig.show() @@ -1325,12 +1563,18 @@ def checkNGroupsValid(yFitIC, ipData): nSpectra = len(ipData) # Number of spectra in the workspace - if (yFitIC.nGlobalFitGroups=="ALL"): + if yFitIC.nGlobalFitGroups == "ALL": yFitIC.nGlobalFitGroups = nSpectra else: - assert isinstance(yFitIC.nGlobalFitGroups, int), "Number of global groups needs to be an integer." - assert yFitIC.nGlobalFitGroups<=nSpectra, "Number of global groups needs to be less or equal to the no of unmasked spectra." - assert yFitIC.nGlobalFitGroups>0, "NUmber of global groups needs to be bigger than zero" + assert isinstance( + yFitIC.nGlobalFitGroups, int + ), "Number of global groups needs to be an integer." + assert ( + yFitIC.nGlobalFitGroups <= nSpectra + ), "Number of global groups needs to be less or equal to the no of unmasked spectra." + assert ( + yFitIC.nGlobalFitGroups > 0 + ), "NUmber of global groups needs to be bigger than zero" return @@ -1343,15 +1587,21 @@ def kMeansClustering(points, centers): Detectors with the same i assigned belong to the same group. """ - prevCenters = centers # Starting centers - while True: - clusters = closestCenter(points, prevCenters) # Form groups by assigning points to their closest center - centers = calculateCenters(points, clusters) # Recalculate centers of new groups + prevCenters = centers # Starting centers + while True: + clusters = closestCenter( + points, prevCenters + ) # Form groups by assigning points to their closest center + centers = calculateCenters( + points, clusters + ) # Recalculate centers of new groups if np.all(centers == prevCenters): break - assert np.isfinite(centers).all(), f"Invalid centers found:\n{centers}\nTry a different number for the groupings." + assert np.isfinite( + centers + ).all(), f"Invalid centers found:\n{centers}\nTry a different number for the groupings." prevCenters = centers @@ -1367,15 +1617,13 @@ def closestCenter(points, centers): """ clusters = np.zeros(len(points)) - for p in range(len(points)): # Iterate over each point - - distMin = np.inf # To be replaced in first iteration + for p in range(len(points)): # Iterate over each point + distMin = np.inf # To be replaced in first iteration for i in range(len(centers)): # Assign closest center to point - dist = pairDistance(points[p], centers[i]) - if dist < distMin: # Store minimum found + if dist < distMin: # Store minimum found distMin = dist closeCenter = i @@ -1385,7 +1633,7 @@ def closestCenter(points, centers): def pairDistance(p1, p2): "Calculates the distance between two points." - return np.sqrt(np.sum(np.square(p1-p2))) + return np.sqrt(np.sum(np.square(p1 - p2))) def calculateCenters(points, clusters): @@ -1395,7 +1643,9 @@ def calculateCenters(points, clusters): centers = np.zeros((nGroups, 2)) for i in range(nGroups): - centers[i, :] = np.mean(points[clusters==i, :], axis=0) # If cluster i is not present, returns nan + centers[i, :] = np.mean( + points[clusters == i, :], axis=0 + ) # If cluster i is not present, returns nan return centers @@ -1404,7 +1654,7 @@ def formIdxList(clusters): idxList = [] for i in np.unique(clusters): - idxs = np.argwhere(clusters==i).flatten() + idxs = np.argwhere(clusters == i).flatten() idxList.append(list(idxs)) # Print groupings information @@ -1419,11 +1669,13 @@ def formIdxList(clusters): def plotDetsAndInitialCenters(L1, theta, centers): """Used in debugging.""" - fig, ax = plt.subplots(tight_layout=True, subplot_kw={'projection':'mantid'}) + fig, ax = plt.subplots(tight_layout=True, subplot_kw={"projection": "mantid"}) fig.canvas.manager.set_window_title("Starting centroids for groupings") ax.scatter(L1, theta, alpha=0.3, color="r", label="Detectors") ax.scatter(centers[:, 0], centers[:, 1], color="k", label="Starting centroids") - ax.axes.xaxis.set_ticks([]) # Numbers plotted do not correspond to real numbers, so hide them + ax.axes.xaxis.set_ticks( + [] + ) # Numbers plotted do not correspond to real numbers, so hide them ax.axes.yaxis.set_ticks([]) ax.set_xlabel("L1") ax.set_ylabel("Theta") @@ -1448,6 +1700,7 @@ def plotFinalGroups(ax, ipData, idxList): ax.legend() return + # --------- Weighted Avg of detectors @@ -1456,9 +1709,11 @@ def avgWeightDetGroups(dataX, dataY, dataE, dataRes, idxList, yFitIC): Performs weighted average on each detector group given by the index list. The imput arrays do not include masked spectra. """ - assert ~np.any(np.all(dataY==0, axis=1)), f"Input data should not include masked spectra at: {np.argwhere(np.all(dataY==0, axis=1))}" + assert ~np.any( + 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.maskTypeProcedure == "NAN": return avgGroupsWithBins(dataX, dataY, dataE, dataRes, idxList, yFitIC) # Use Default for unmasked or NCP masked @@ -1475,24 +1730,31 @@ def avgGroupsOverCols(dataX, dataY, dataE, dataRes, idxList): wDataX, wDataY, wDataE, wDataRes = initiateZeroArr((len(idxList), len(dataY[0]))) for i, idxs in enumerate(idxList): - groupX, groupY, groupE, groupRes = extractArrByIdx(dataX, dataY, dataE, dataRes, idxs) + groupX, groupY, groupE, groupRes = extractArrByIdx( + dataX, dataY, dataE, dataRes, idxs + ) assert len(groupY) > 0, "Group with zero detectors found, invalid." - if len(groupY) == 1: # Cannot use weight avg in single spec, wrong results + if len(groupY) == 1: # Cannot use weight avg in single spec, wrong results meanY, meanE = groupY, groupE meanRes = groupRes else: meanY, meanE = weightedAvgArr(groupY, groupE) - meanRes = np.nanmean(groupRes, axis=0) # Nans are not present but safeguard + meanRes = np.nanmean(groupRes, axis=0) # Nans are not present but safeguard - assert np.all(groupX[0] == np.mean(groupX, axis=0)), "X values should not change with groups" + assert np.all( + groupX[0] == np.mean(groupX, axis=0) + ), "X values should not change with groups" - for wsData, mean in zip([wDataX, wDataY, wDataE, wDataRes], [groupX[0], meanY, meanE, meanRes]): + for wsData, mean in zip( + [wDataX, wDataY, wDataE, wDataRes], [groupX[0], meanY, meanE, meanRes] + ): wsData[i] = mean - assert ~np.any(np.all(wDataY==0, axis=1) - ), f"Some avg weights in groups are not being performed:\n{np.argwhere(np.all(wDataY==0, axis=1))}" + assert ~np.any( + np.all(wDataY == 0, axis=1) + ), f"Some avg weights in groups are not being performed:\n{np.argwhere(np.all(wDataY==0, axis=1))}" return wDataX, wDataY, wDataE, wDataRes @@ -1509,13 +1771,17 @@ def avgGroupsWithBins(dataX, dataY, dataE, dataRes, idxList, yFitIC): wDataX, wDataY, wDataE, wDataRes = initiateZeroArr((len(idxList), len(meanX))) for i, idxs in enumerate(idxList): - groupX, groupY, groupE, groupRes = extractArrByIdx(dataX, dataY, dataE, dataRes, idxs) + groupX, groupY, groupE, groupRes = extractArrByIdx( + dataX, dataY, dataE, dataRes, idxs + ) meanY, meanE = weightedAvgXBinsArr(groupX, groupY, groupE, meanX) - meanRes = np.nanmean(groupRes, axis=0) # Nans are not present but safeguard + meanRes = np.nanmean(groupRes, axis=0) # Nans are not present but safeguard - for wsData, mean in zip([wDataX, wDataY, wDataE, wDataRes], [meanX, meanY, meanE, meanRes]): + for wsData, mean in zip( + [wDataX, wDataY, wDataE, wDataRes], [meanX, meanY, meanE, meanRes] + ): wsData[i] = mean return wDataX, wDataY, wDataE, wDataRes @@ -1526,7 +1792,7 @@ def initiateZeroArr(shape): wDataY = np.zeros(shape) wDataE = np.zeros(shape) wDataRes = np.zeros(shape) - return wDataX, wDataY, wDataE, wDataRes + return wDataX, wDataY, wDataE, wDataRes def extractArrByIdx(dataX, dataY, dataE, dataRes, idxs): @@ -1544,16 +1810,18 @@ def calcCostFun(model, i, x, y, yerr, res, sharedPars): def convolvedModel(xrange, y0, *pars): """Performs convolution first on high density grid and interpolates to desired x range""" - return y0 + signal.convolve(model(xrange, *pars), resDense, mode="same") * xDelta + return ( + y0 + signal.convolve(model(xrange, *pars), resDense, mode="same") * xDelta + ) signature = describe(model)[:] signature[1:1] = ["y0"] - costSig = [key if key in sharedPars else key+str(i) for key in signature] + costSig = [key if key in sharedPars else key + str(i) for key in signature] convolvedModel.func_code = make_func_code(costSig) # Select only valid data, i.e. when error is not 0 or nan or inf - nonZeros= (yerr!=0) & (yerr!=np.nan) & (yerr!=np.inf) & (y!=np.nan) + nonZeros = (yerr != 0) & (yerr != np.nan) & (yerr != np.inf) & (y != np.nan) xNZ = x[nonZeros] yNZ = y[nonZeros] yerrNZ = yerr[nonZeros] @@ -1573,12 +1841,11 @@ def minuitInitialParameters(defaultPars, sharedPars, nSpec): unsharedPars = [key for key in defaultPars if key not in sharedPars] for up in unsharedPars: for i in range(nSpec): - initPars[up+str(i)] = defaultPars[up] + initPars[up + str(i)] = defaultPars[up] return initPars def plotGlobalFit(dataX, dataY, dataE, mObj, totCost, wsName): - if len(dataY) > 10: print("\nToo many axes to show in figure, skipping the plot ...\n") return @@ -1586,12 +1853,12 @@ def plotGlobalFit(dataX, dataY, dataE, mObj, totCost, wsName): rows = 2 fig, axs = plt.subplots( rows, - int(np.ceil(len(dataY)/rows)), + int(np.ceil(len(dataY) / rows)), figsize=(15, 8), tight_layout=True, - subplot_kw={'projection':'mantid'} + subplot_kw={"projection": "mantid"}, ) - fig.canvas.manager.set_window_title(wsName+"_Plot_of_Global_Fit") + fig.canvas.manager.set_window_title(wsName + "_Plot_of_Global_Fit") # Data used in Global Fit for i, (x, y, yerr, ax) in enumerate(zip(dataX, dataY, dataE, axs.flat)): diff --git a/EVSVesuvio/vesuvio_analysis/procedures.py b/EVSVesuvio/vesuvio_analysis/procedures.py index 22c39d34..e50e412e 100644 --- a/EVSVesuvio/vesuvio_analysis/procedures.py +++ b/EVSVesuvio/vesuvio_analysis/procedures.py @@ -19,8 +19,12 @@ def runIndependentIterativeProcedure(IC, clearWS=True): def runJointBackAndForwardProcedure(bckwdIC, fwdIC, clearWS=True): - assert bckwdIC.modeRunning == "BACKWARD", "Missing backward IC, args usage: (bckwdIC, fwdIC)" - assert fwdIC.modeRunning == "FORWARD", "Missing forward IC, args usage: (bckwdIC, fwdIC)" + assert ( + bckwdIC.modeRunning == "BACKWARD" + ), "Missing backward IC, args usage: (bckwdIC, fwdIC)" + assert ( + fwdIC.modeRunning == "FORWARD" + ), "Missing forward IC, args usage: (bckwdIC, fwdIC)" # Clear worksapces before running one of the procedures below if clearWS: @@ -36,7 +40,9 @@ def runPreProcToEstHRatio(bckwdIC, fwdIC): Runs iterative procedure with alternating back and forward scattering. """ - assert bckwdIC.runningSampleWS is False, "Preliminary procedure not suitable for Bootstrap." + assert ( + bckwdIC.runningSampleWS is False + ), "Preliminary procedure not suitable for Bootstrap." fwdIC.runningPreliminary = True # Store original no of MS and set MS iterations to zero @@ -47,13 +53,12 @@ def runPreProcToEstHRatio(bckwdIC, fwdIC): nIter = askUserNoOfIterations() - HRatios = [] # List to store HRatios + HRatios = [] # List to store HRatios massIdxs = [] # Run preliminary forward with a good guess for the widths of non-H masses wsFinal, fwdScatResults = iterativeFitForDataReduction(fwdIC) - for i in range(int(nIter)): # Loop until convergence is achieved - - AnalysisDataService.clear() # Clears all Workspaces + for i in range(int(nIter)): # Loop until convergence is achieved + AnalysisDataService.clear() # Clears all Workspaces # Update H ratio massIdx, HRatio = calculateHToMassIdxRatio(fwdScatResults) @@ -67,7 +72,9 @@ def runPreProcToEstHRatio(bckwdIC, fwdIC): print(f"\nIdxs of masses for H ratio for each iteration: \n{massIdxs}") print(f"\nCorresponding H ratios: \n{HRatios}") - fwdIC.runningPreliminary = False # Change to default since end of preliminary procedure + fwdIC.runningPreliminary = ( + False # Change to default since end of preliminary procedure + ) # Set original number of MS iterations for IC, ori in zip([bckwdIC, fwdIC], oriMS): @@ -84,7 +91,9 @@ def runPreProcToEstHRatio(bckwdIC, fwdIC): def createTableWSHRatios(HRatios, massIdxs): - tableWS = CreateEmptyTableWorkspace(OutputWorkspace="H_Ratios_From_Preliminary_Procedure") + tableWS = CreateEmptyTableWorkspace( + OutputWorkspace="H_Ratios_From_Preliminary_Procedure" + ) tableWS.setTitle("H Ratios and Idxs at each iteration") tableWS.addColumn(type="int", name="iter") tableWS.addColumn(type="float", name="H Ratio") @@ -96,10 +105,14 @@ def createTableWSHRatios(HRatios, massIdxs): def askUserNoOfIterations(): print("\nH was detected but HToMassIdxRatio was not provided.") - print("\nSugested preliminary procedure:\n\nrun_forward\nfor n:\n estimate_HToMassIdxRatio\n run_backward\n" - " run_forward") - userInput = input("\n\nDo you wish to run preliminary procedure to estimate HToMassIdxRatio? (y/n)") - if not((userInput=="y") or (userInput=="Y")): + print( + "\nSugested preliminary procedure:\n\nrun_forward\nfor n:\n estimate_HToMassIdxRatio\n run_backward\n" + " run_forward" + ) + userInput = input( + "\n\nDo you wish to run preliminary procedure to estimate HToMassIdxRatio? (y/n)" + ) + if not ((userInput == "y") or (userInput == "Y")): raise KeyboardInterrupt("Preliminary procedure interrupted.") nIter = int(input("\nHow many iterations do you wish to run? n=")) @@ -116,8 +129,12 @@ def calculateHToMassIdxRatio(fwdScatResults): # To find idx of mass in backward scattering, take out first mass H fwdIntensitiesNoH = fwdMeanIntensityRatios[1:] - massIdx = np.argmax(fwdIntensitiesNoH) # Idex of forward inensities, which include H - assert fwdIntensitiesNoH[massIdx] != 0, "Cannot estimate H intensity since maximum peak from backscattering is zero." + massIdx = np.argmax( + fwdIntensitiesNoH + ) # Idex of forward inensities, which include H + assert ( + fwdIntensitiesNoH[massIdx] != 0 + ), "Cannot estimate H intensity since maximum peak from backscattering is zero." HRatio = fwdMeanIntensityRatios[0] / fwdIntensitiesNoH[massIdx] @@ -144,9 +161,9 @@ def setInitFwdParsFromBackResults(bckwdScatResults, bckwdIC, fwdIC): backMeanIntensityRatios = bckwdScatResults.all_mean_intensities[-1] if isHPresent(fwdIC.masses): - - assert len(backMeanWidths) == fwdIC.noOfMasses-1, "H Mass present, no of masses in frontneeds to be bigger" \ - "than back by 1." + assert len(backMeanWidths) == fwdIC.noOfMasses - 1, ( + "H Mass present, no of masses in frontneeds to be bigger" "than back by 1." + ) # Use H ratio to calculate intensity ratios HIntensity = bckwdIC.HToMassIdxRatio * backMeanIntensityRatios[bckwdIC.massIdx] @@ -160,34 +177,37 @@ def setInitFwdParsFromBackResults(bckwdScatResults, bckwdIC, fwdIC): # Set forward widths from backscattering fwdIC.initPars[4::3] = backMeanWidths # Fix all widths except for H, i.e. the first one - fwdIC.bounds[4::3] = backMeanWidths[:, np.newaxis] * np.ones((1,2)) - - else: # H mass not present anywhere + fwdIC.bounds[4::3] = backMeanWidths[:, np.newaxis] * np.ones((1, 2)) - assert len(backMeanWidths) == fwdIC.noOfMasses, "H Mass not present, no of masses needs to be the same for" \ - "front and back scattering." + else: # H mass not present anywhere + assert len(backMeanWidths) == fwdIC.noOfMasses, ( + "H Mass not present, no of masses needs to be the same for" + "front and back scattering." + ) # Set widths and intensity ratios fwdIC.initPars[1::3] = backMeanWidths fwdIC.initPars[0::3] = backMeanIntensityRatios - if len(backMeanWidths) > 1: # In the case of single mass, width is not fixed + if len(backMeanWidths) > 1: # In the case of single mass, width is not fixed # Fix all widhts except first - fwdIC.bounds[4::3] = backMeanWidths[1:][:, np.newaxis] * np.ones((1,2)) + fwdIC.bounds[4::3] = backMeanWidths[1:][:, np.newaxis] * np.ones((1, 2)) - print("\nChanged initial conditions of forward scattering according to mean widhts and intensity ratios from " - "backscattering.\n") + print( + "\nChanged initial conditions of forward scattering according to mean widhts and intensity ratios from " + "backscattering.\n" + ) return def isHPresent(masses) -> bool: + Hmask = np.abs(masses - 1) / 1 < 0.1 # H mass whithin 10% of 1 au - Hmask = np.abs(masses-1)/1 < 0.1 # H mass whithin 10% of 1 au - - if np.any(Hmask): # H present - + 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 ( + 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 diff --git a/EVSVesuvio/vesuvio_analysis/run_script.py b/EVSVesuvio/vesuvio_analysis/run_script.py index c6e51f7f..e62a8761 100644 --- a/EVSVesuvio/vesuvio_analysis/run_script.py +++ b/EVSVesuvio/vesuvio_analysis/run_script.py @@ -1,13 +1,32 @@ -from EVSVesuvio.vesuvio_analysis.ICHelpers import buildFinalWSName, completeICFromInputs, completeBootIC, completeYFitIC +from EVSVesuvio.vesuvio_analysis.ICHelpers import ( + buildFinalWSName, + completeICFromInputs, + completeBootIC, + completeYFitIC, +) from EVSVesuvio.vesuvio_analysis.bootstrap import runBootstrap from EVSVesuvio.vesuvio_analysis.fit_in_yspace import fitInYSpaceProcedure -from EVSVesuvio.vesuvio_analysis.procedures import runIndependentIterativeProcedure, runJointBackAndForwardProcedure, \ - runPreProcToEstHRatio, createTableWSHRatios, isHPresent +from EVSVesuvio.vesuvio_analysis.procedures import ( + runIndependentIterativeProcedure, + runJointBackAndForwardProcedure, + runPreProcToEstHRatio, + createTableWSHRatios, + isHPresent, +) from mantid.api import mtd -def runScript(userCtr, scriptName, wsBackIC, wsFrontIC, bckwdIC, fwdIC, yFitIC, bootIC, yes_to_all=False): - +def runScript( + userCtr, + scriptName, + wsBackIC, + wsFrontIC, + bckwdIC, + fwdIC, + yFitIC, + bootIC, + yes_to_all=False, +): # Set extra attributes from user attributes completeICFromInputs(fwdIC, scriptName, wsFrontIC) completeICFromInputs(bckwdIC, scriptName, wsBackIC) @@ -16,7 +35,9 @@ def runScript(userCtr, scriptName, wsBackIC, wsFrontIC, bckwdIC, fwdIC, yFitIC, checkInputs(userCtr) checkInputs(bootIC) - assert not(userCtr.runRoutine & bootIC.runBootstrap), "Main routine and bootstrap both set to run!" + assert not ( + userCtr.runRoutine & bootIC.runBootstrap + ), "Main routine and bootstrap both set to run!" def runProcedure(): proc = userCtr.procedure # Shorthad to make it easier to read @@ -25,18 +46,21 @@ def runProcedure(): return ranPreliminary = False - if (proc=="BACKWARD") | (proc=="JOINT"): + if (proc == "BACKWARD") | (proc == "JOINT"): if isHPresent(fwdIC.masses) & (bckwdIC.HToMassIdxRatio is None): - HRatios, massIdxs = runPreProcToEstHRatio(bckwdIC, fwdIC) # Sets H ratio to bckwdIC automatically + HRatios, massIdxs = runPreProcToEstHRatio( + bckwdIC, fwdIC + ) # Sets H ratio to bckwdIC automatically ranPreliminary = True - assert (isHPresent(fwdIC.masses) != (bckwdIC.HToMassIdxRatio is None) - ), "When H is not present, HToMassIdxRatio has to be set to None" + assert isHPresent(fwdIC.masses) != ( + bckwdIC.HToMassIdxRatio is None + ), "When H is not present, HToMassIdxRatio has to be set to None" - if (proc=="BACKWARD"): + if proc == "BACKWARD": res = runIndependentIterativeProcedure(bckwdIC) - if (proc=="FORWARD"): + if proc == "FORWARD": res = runIndependentIterativeProcedure(fwdIC) - if (proc=="JOINT"): + if proc == "JOINT": res = runJointBackAndForwardProcedure(bckwdIC, fwdIC) # If preliminary procedure ran, make TableWS with H ratios values @@ -48,39 +72,47 @@ def runProcedure(): wsNames = [] ICs = [] for mode, IC in zip(["BACKWARD", "FORWARD"], [bckwdIC, fwdIC]): - if (userCtr.fitInYSpace==mode) | (userCtr.fitInYSpace=="JOINT"): + if (userCtr.fitInYSpace == mode) | (userCtr.fitInYSpace == "JOINT"): wsNames.append(buildFinalWSName(scriptName, mode, IC)) ICs.append(IC) # If bootstrap is not None, run bootstrap procedure and finish if bootIC.runBootstrap: - assert (bootIC.procedure=="FORWARD") | (bootIC.procedure=="BACKWARD") | (bootIC.procedure=="JOINT"), "Invalid Bootstrap procedure." + assert ( + (bootIC.procedure == "FORWARD") + | (bootIC.procedure == "BACKWARD") + | (bootIC.procedure == "JOINT") + ), "Invalid Bootstrap procedure." return runBootstrap(bckwdIC, fwdIC, bootIC, yFitIC), None # Default workflow for procedure + fit in y space if userCtr.runRoutine: # Check if final ws are loaded: - wsInMtd = [ws in mtd for ws in wsNames] # Bool list - if (len(wsInMtd)>0) and all(wsInMtd): # When wsName is empty list, loop doesn't run + wsInMtd = [ws in mtd for ws in wsNames] # Bool list + if (len(wsInMtd) > 0) and all( + wsInMtd + ): # When wsName is empty list, loop doesn't run for wsName, IC in zip(wsNames, ICs): resYFit = fitInYSpaceProcedure(yFitIC, IC, mtd[wsName]) - return None, resYFit # To match return below. + return None, resYFit # To match return below. - checkUserClearWS(yes_to_all) # Check if user is OK with cleaning all workspaces + checkUserClearWS(yes_to_all) # Check if user is OK with cleaning all workspaces res = runProcedure() resYFit = None for wsName, IC in zip(wsNames, ICs): resYFit = fitInYSpaceProcedure(yFitIC, IC, mtd[wsName]) - return res, resYFit # Return results used only in tests + return res, resYFit # Return results used only in tests def checkUserClearWS(yes_to_all=False): """If any workspace is loaded, check if user is sure to start new procedure.""" if not yes_to_all and len(mtd) != 0: - userInput = input("This action will clean all current workspaces to start anew. Proceed? (y/n): ") + userInput = input( + "This action will clean all current workspaces to start anew. Proceed? (y/n): " + ) if (userInput == "y") | (userInput == "Y"): pass else: @@ -89,7 +121,6 @@ def checkUserClearWS(yes_to_all=False): def checkInputs(crtIC): - try: if ~crtIC.runRoutine: return @@ -98,7 +129,12 @@ def checkInputs(crtIC): 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 ( + (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 diff --git a/EVSVesuvio/vesuvio_analysis/tests/test_analysis_functions.py b/EVSVesuvio/vesuvio_analysis/tests/test_analysis_functions.py index ff40241c..b59bbdd8 100644 --- a/EVSVesuvio/vesuvio_analysis/tests/test_analysis_functions.py +++ b/EVSVesuvio/vesuvio_analysis/tests/test_analysis_functions.py @@ -21,5 +21,5 @@ def test_extract_ws_calls_extract_X_Y_and_E(self): self.ws.extractE.assert_called_once() -if __name__ == '__main__': +if __name__ == "__main__": unittest.main() diff --git a/config.flake8 b/config.flake8 index 9c61d7bb..87ff3f8a 100644 --- a/config.flake8 +++ b/config.flake8 @@ -4,9 +4,6 @@ exclude = .git, build, buildconfig, - EVSVesuvio/vesuvio_analysis/scribles, - EVSVesuvio/experiments/original_inputs, - EVSVesuvio/original_inputs.py, unpackaged max-complexity = 20 max-line-length = 140