From 341e02797b472fcf3e85694a401c868078159acf Mon Sep 17 00:00:00 2001 From: GuiMacielPereira Date: Tue, 10 Dec 2024 17:05:19 +0000 Subject: [PATCH] Rename variables for inputs of y space --- src/mvesuvio/analysis_fitting.py | 100 +++++++++--------- src/mvesuvio/config/analysis_inputs.py | 32 ++++-- tests/data/analysis/inputs/analysis_test.py | 18 ++-- .../data/analysis/inputs/yspace_gauss_test.py | 2 +- tests/data/analysis/inputs/yspace_gc_test.py | 4 +- 5 files changed, 84 insertions(+), 72 deletions(-) diff --git a/src/mvesuvio/analysis_fitting.py b/src/mvesuvio/analysis_fitting.py index 0b811d4..3002800 100644 --- a/src/mvesuvio/analysis_fitting.py +++ b/src/mvesuvio/analysis_fitting.py @@ -27,13 +27,13 @@ def fitInYSpaceProcedure(yFitIC, IC, wsTOF): wsTOFMass0 = subtract_profiles_except_lightest(IC, wsTOF) - if yFitIC.subtractFSE: + if yFitIC.subtract_calculated_fse_from_data: wsFSEMass0 = find_ws_name_fse_first_mass(IC) wsTOFMass0 = Minus(wsTOFMass0, wsFSEMass0, OutputWorkspace=wsTOFMass0.name() + "_fse") wsJoY, wsJoYAvg = ySpaceReduction(wsTOFMass0, IC.masses[0], yFitIC, IC) - if yFitIC.symmetrisationFlag: + if yFitIC.do_symmetrisation: wsJoYAvg = symmetrizeWs(wsJoYAvg) fitProfileMinuit(yFitIC, wsJoYAvg, wsResSum) @@ -44,7 +44,7 @@ def fitInYSpaceProcedure(yFitIC, IC, wsTOF): yfitResults = ResultsYFitObject(IC, yFitIC, wsTOF.name(), wsJoYAvg.name()) yfitResults.save() - if yFitIC.globalFit: + if yFitIC.do_global_fit: runGlobalFit(wsJoY, wsRes, IC, yFitIC) return yfitResults @@ -76,7 +76,7 @@ def calculateMantidResolutionFirstMass(IC, yFitIC, ws): ) Rebin( InputWorkspace="tmp", - Params=yFitIC.rebinParametersForYSpaceFit, + Params=yFitIC.range_for_rebinning_in_y_space, OutputWorkspace="tmp", ) @@ -129,10 +129,10 @@ def ySpaceReduction(wsTOF, mass0, yFitIC, ic): ws_name_lightest_mass = ic.name + '_' + str(ic.number_of_iterations_for_corrections) + '_' + str(min(ic.masses)) + '_ncp' ncp = mtd[ws_name_lightest_mass].extractY() - rebinPars = yFitIC.rebinParametersForYSpaceFit + rebinPars = yFitIC.range_for_rebinning_in_y_space if np.any(np.all(wsTOF.extractY() == 0, axis=0)): # Masked columns present - if yFitIC.maskTypeProcedure == "NAN": + if yFitIC.mask_zeros_with == "nan": # Build special workspace to store accumulated points wsJoY = convertToYSpace(wsTOF, mass0) xp = buildXRangeFromRebinPars(yFitIC) @@ -152,14 +152,14 @@ def ySpaceReduction(wsTOF, mass0, yFitIC, ic): wsJoYAvg = weightedAvgXBins(wsJoYN, xp) return wsJoYN, wsJoYAvg - elif yFitIC.maskTypeProcedure == "NCP": + elif yFitIC.mask_zeros_with == "ncp": wsTOF = replaceZerosWithNCP(wsTOF, ncp) else: raise ValueError( """ Masked TOF bins were found but no valid procedure in y-space fit was selected. - Options: 'NAN', 'NCP' + Options: 'nan', 'ncp' """ ) @@ -206,7 +206,7 @@ def replaceZerosWithNCP(ws, ncp): def buildXRangeFromRebinPars(yFitIC): # Range used in case mask is set to NAN first, step, last = [ - float(s) for s in yFitIC.rebinParametersForYSpaceFit.split(",") + float(s) for s in yFitIC.range_for_rebinning_in_y_space.split(",") ] xp = np.arange(first, last, step) + step / 2 # Correction to match Mantid range return xp @@ -519,7 +519,7 @@ def fitProfileMinuit(yFitIC, wsYSpaceSym, wsRes): resX, resY, resE = extractFirstSpectra(wsRes) assert np.all(dataX == resX), "Resolution should operate on the same range as DataX" - model, defaultPars, sharedPars = selectModelAndPars(yFitIC.fitModel) + model, defaultPars, sharedPars = selectModelAndPars(yFitIC.fitting_model) xDelta, resDense = oddPointsRes(resX, resY) @@ -544,11 +544,11 @@ def convolvedModel(x, y0, *pars): m = Minuit(costFun, **defaultPars) m.limits["A"] = (0, None) - if yFitIC.fitModel == "DOUBLE_WELL": + if yFitIC.fitting_model== "doublewell": m.limits["d"] = (0, None) m.limits["R"] = (0, None) - if yFitIC.fitModel == "SINGLE_GAUSSIAN": + if yFitIC.fitting_model== "gauss": m.simplex() m.migrad() @@ -616,7 +616,7 @@ def selectModelAndPars(modelFlag): The defaultPars should be in the same order as the signature of the function """ - if modelFlag == "SINGLE_GAUSSIAN": + if modelFlag == "gauss": def model(x, A, x0, sigma): return ( @@ -629,7 +629,7 @@ def model(x, A, x0, sigma): defaultPars = {"A": 1, "x0": 0, "sigma": 5} sharedPars = ["sigma"] # Used only in Global fit - elif modelFlag == "GC_C4_C6": + elif modelFlag == "gcc4c6": def model(x, A, x0, sigma1, c4, c6): return ( @@ -659,7 +659,7 @@ def model(x, A, x0, sigma1, c4, c6): defaultPars = {"A": 1, "x0": 0, "sigma1": 6, "c4": 0, "c6": 0} sharedPars = ["sigma1", "c4", "c6"] # Used only in Global fit - elif modelFlag == "GC_C4": + elif modelFlag == "gcc4": def model(x, A, x0, sigma1, c4): return ( @@ -681,7 +681,7 @@ def model(x, A, x0, sigma1, c4): defaultPars = {"A": 1, "x0": 0, "sigma1": 6, "c4": 0} sharedPars = ["sigma1", "c4"] # Used only in Global fit - elif modelFlag == "GC_C6": + elif modelFlag == "gcc6": def model(x, A, x0, sigma1, c6): return ( @@ -704,7 +704,7 @@ def model(x, A, x0, sigma1, c6): defaultPars = {"A": 1, "x0": 0, "sigma1": 6, "c6": 0} sharedPars = ["sigma1", "c6"] # Used only in Global fit - elif modelFlag == "DOUBLE_WELL": + elif modelFlag == "doublewell": def model(x, A, d, R, sig1, sig2): # h = 2.04 @@ -744,7 +744,7 @@ def model(x, A, d, R, sig1, sig2): } # TODO: Starting parameters and bounds? sharedPars = ["d", "R", "sig1", "sig2"] # Only varying parameter is amplitude A - elif modelFlag == "ANSIO_GAUSSIAN": + elif modelFlag == "ansiogauss": # Ansiotropic case def model(x, A, sig1, sig2): # h = 2.04 @@ -765,7 +765,7 @@ def model(x, A, sig1, sig2): defaultPars = {"A": 1, "sig1": 3, "sig2": 5} sharedPars = ["sig1", "sig2"] - elif modelFlag == "Gaussian3D": + elif modelFlag == "gauss3d": def model(x, A, sig_x, sig_y, sig_z): y = x[:, np.newaxis, np.newaxis] @@ -795,7 +795,7 @@ def model(x, A, sig_x, sig_y, sig_z): else: raise ValueError( - "Fitting Model not recognized, available options: 'SINGLE_GAUSSIAN', 'GC_C4_C6', 'GC_C4'" + "Fitting Model not recognized, available options: 'gauss', 'gcc4c6', 'gcc4', 'gcc6', 'ansiogauss' gauss3d'" ) print("\nShared Parameters: ", [key for key in sharedPars]) @@ -915,7 +915,7 @@ def runMinos(mObj, yFitIC, constrFunc, wsName): errors = list(mObj.errors) # If minos is set not to run, ouput columns with zeros on minos errors - if not (yFitIC.runMinos): + if not (yFitIC.run_minos): minosAutoErr = list(np.zeros((len(parameters), 2))) minosManErr = list(np.zeros((len(parameters), 2))) return parameters, values, errors, minosAutoErr, minosManErr @@ -927,7 +927,7 @@ def runMinos(mObj, yFitIC, constrFunc, wsName): bestFitErrs[p] = e if ( - yFitIC.fitModel == "SINGLE_GAUSSIAN" + yFitIC.fitting_model== "gauss" ): # Case with no positivity constraint, can use automatic minos() mObj.minos() me = mObj.merrors @@ -938,13 +938,13 @@ def runMinos(mObj, yFitIC, constrFunc, wsName): minosAutoErr.append([me[p].lower, me[p].upper]) minosManErr = list(np.zeros(np.array(minosAutoErr).shape)) - if yFitIC.showPlots: + if yFitIC.show_plots: plotAutoMinos(mObj, wsName) else: # Case with positivity constraint on function, use manual implementation # Changes values of minuit obj m, do not use m below this point merrors, fig = runAndPlotManualMinos( - mObj, constrFunc, bestFitVals, bestFitErrs, yFitIC.showPlots + mObj, constrFunc, bestFitVals, bestFitErrs, yFitIC.show_plots ) # Same as above, but the other way around @@ -953,7 +953,7 @@ def runMinos(mObj, yFitIC, constrFunc, wsName): minosManErr.append(merrors[p]) minosAutoErr = list(np.zeros(np.array(minosManErr).shape)) - if yFitIC.showPlots: + if yFitIC.show_plots: fig.canvas.manager.set_window_title(wsName + "_Manual_Implementation_MINOS") fig.show() @@ -990,7 +990,6 @@ def runAndPlotManualMinos(minuitObj, constrFunc, bestFitVals, bestFitErrs, showP ) merrors[p] = np.array([lerr, uerr]) - # if showPlots: # Hide plots not in use: for ax in axs.flat: if not ax.lines: # If empty list @@ -999,7 +998,6 @@ def runAndPlotManualMinos(minuitObj, constrFunc, bestFitVals, bestFitErrs, showP # ALl axes share same legend, so set figure legend to first axis handle, label = axs[0, 0].get_legend_handles_labels() fig.legend(handle, label, loc="lower right") - # fig.show() return merrors, fig @@ -1226,13 +1224,13 @@ def oddPointsRes(x, res): def fitProfileMantidFit(yFitIC, wsYSpaceSym, wsRes): print("\nFitting on the sum of spectra in the West domain ...\n") for minimizer in ["Levenberg-Marquardt", "Simplex"]: - if yFitIC.fitModel == "SINGLE_GAUSSIAN": + if yFitIC.fitting_model== "gauss": function = f"""composite=Convolution,FixResolution=true,NumDeriv=true; name=Resolution,Workspace={wsRes.name()},WorkspaceIndex=0; name=UserFunction,Formula=y0 + A*exp( -(x-x0)^2/2/sigma^2)/(2*3.1415*sigma^2)^0.5, y0=0,A=1,x0=0,sigma=5, ties=()""" - elif yFitIC.fitModel == "GC_C4_C6": + elif yFitIC.fitting_model== "gcc4c6": function = f""" composite=Convolution,FixResolution=true,NumDeriv=true; name=Resolution,Workspace={wsRes.name()},WorkspaceIndex=0,X=(),Y=(); @@ -1241,7 +1239,7 @@ def fitProfileMantidFit(yFitIC, wsYSpaceSym, wsRes): (64*((x-x0)/sqrt(2)/sigma1)^6 - 480*((x-x0)/sqrt(2)/sigma1)^4 + 720*((x-x0)/sqrt(2)/sigma1)^2 - 120)), y0=0, A=1,x0=0,sigma1=4.0,c4=0.0,c6=0.0,ties=(),constraints=(07s} = {v:>8.4f} \u00B1 {e:<8.4f}") print("\n") - if yFitIC.showPlots: + if yFitIC.show_plots: plotGlobalFit(dataX, dataY, dataE, m, totCost, wsYSpace.name()) return np.array(m.values), np.array( @@ -1533,7 +1533,7 @@ def groupDetectors(ipData, yFitIC): checkNGroupsValid(yFitIC, ipData) - print(f"\nNumber of gropus: {yFitIC.nGlobalFitGroups}") + print(f"\nNumber of gropus: {yFitIC.number_of_global_fit_groups}") L1 = ipData[:, -1].copy() theta = ipData[:, 2].copy() @@ -1547,7 +1547,7 @@ def groupDetectors(ipData, yFitIC): points = np.vstack((L1, theta)).T assert points.shape == (len(L1), 2), "Wrong shape." # Initial centers of groups - startingIdxs = np.linspace(0, len(points) - 1, yFitIC.nGlobalFitGroups).astype(int) + startingIdxs = np.linspace(0, len(points) - 1, yFitIC.number_of_global_fit_groups).astype(int) centers = points[ startingIdxs, : ] # Centers of cluster groups, NOT fitting parameter @@ -1558,7 +1558,7 @@ def groupDetectors(ipData, yFitIC): clusters = kMeansClustering(points, centers) idxList = formIdxList(clusters) - if yFitIC.showPlots: + if yFitIC.show_plots: fig, ax = plt.subplots(tight_layout=True, subplot_kw={"projection": "mantid"}) fig.canvas.manager.set_window_title("Grouping of detectors") plotFinalGroups(ax, ipData, idxList) @@ -1571,17 +1571,17 @@ def checkNGroupsValid(yFitIC, ipData): nSpectra = len(ipData) # Number of spectra in the workspace - if yFitIC.nGlobalFitGroups == "ALL": - yFitIC.nGlobalFitGroups = nSpectra + if yFitIC.number_of_global_fit_groups == "all": + yFitIC.number_of_global_fit_groups = nSpectra else: assert isinstance( - yFitIC.nGlobalFitGroups, int + yFitIC.number_of_global_fit_groups, int ), "Number of global groups needs to be an integer." assert ( - yFitIC.nGlobalFitGroups <= nSpectra + yFitIC.number_of_global_fit_groups <= nSpectra ), "Number of global groups needs to be less or equal to the no of unmasked spectra." assert ( - yFitIC.nGlobalFitGroups > 0 + yFitIC.number_of_global_fit_groups > 0 ), "NUmber of global groups needs to be bigger than zero" return @@ -1721,7 +1721,7 @@ def avgWeightDetGroups(dataX, dataY, dataE, dataRes, idxList, yFitIC): np.all(dataY == 0, axis=1) ), f"Input data should not include masked spectra at: {np.argwhere(np.all(dataY==0, axis=1))}" - if yFitIC.maskTypeProcedure == "NAN": + if yFitIC.mask_zeros_with == "nan": return avgGroupsWithBins(dataX, dataY, dataE, dataRes, idxList, yFitIC) # Use Default for unmasked or NCP masked diff --git a/src/mvesuvio/config/analysis_inputs.py b/src/mvesuvio/config/analysis_inputs.py index 06bacb1..0bc4039 100644 --- a/src/mvesuvio/config/analysis_inputs.py +++ b/src/mvesuvio/config/analysis_inputs.py @@ -11,7 +11,7 @@ class SampleParameters: @dataclass class BackwardAnalysisInputs(SampleParameters): - run_this_scattering_type = True + run_this_scattering_type = False fit_in_y_space = False runs = "43066-43076" @@ -83,15 +83,27 @@ class ForwardAnalysisInputs(SampleParameters): @dataclass class YSpaceFitInputs: - showPlots = True - symmetrisationFlag = False - subtractFSE = True - rebinParametersForYSpaceFit = "-25, 0.5, 25" # Needs to be symetric - fitModel = "SINGLE_GAUSSIAN" # Options: 'SINGLE_GAUSSIAN', 'GC_C4', 'GC_C6', 'GC_C4_C6', 'DOUBLE_WELL', 'ANSIO_GAUSSIAN', 'Gaussian3D' - runMinos = True - globalFit = True # Performs global fit with Minuit by default - nGlobalFitGroups = 4 # Number or string "ALL" - maskTypeProcedure = "NAN" # Options: 'NCP', 'NAN', None + show_plots = True + do_symmetrisation = False + subtract_calculated_fse_from_data = True + range_for_rebinning_in_y_space = "-25, 0.5, 25" # Needs to be symetric + # Fitting model options + # 'gauss': Single Gaussian + # 'gcc4': Gram-Charlier with C4 parameter + # 'gcc6': Gram-Charlier with C6 parameter + # 'doublewell': Double Well function + # 'gauss3D': 3-Dimensional Gaussian + fitting_model = "gauss" + run_minos = True + do_global_fit = True # Performs global fit with Minuit by default + # Number of groups of detectors to perform global (simultaneous) fit on + # Either an integer less than the number of detectors + # or option 'all', which does not form groups and fits all spectra simultaneously and individualy + number_of_global_fit_groups = 4 + # Type of masking + # 'nan': Zeros in workspace being fit are ignored + # 'ncp': Zeros in workspace being fit are replaced by the fitted neutron compton profile + mask_zeros_with = "nan" ######################## diff --git a/tests/data/analysis/inputs/analysis_test.py b/tests/data/analysis/inputs/analysis_test.py index 180be19..8cb9a91 100644 --- a/tests/data/analysis/inputs/analysis_test.py +++ b/tests/data/analysis/inputs/analysis_test.py @@ -83,15 +83,15 @@ class ForwardAnalysisInputs(SampleParameters): class YSpaceFitInputs: - showPlots = False - symmetrisationFlag = True - subtractFSE = False - rebinParametersForYSpaceFit = "-20, 0.5, 20" # Needs to be symetric - fitModel = "SINGLE_GAUSSIAN" - runMinos = True - globalFit = None - nGlobalFitGroups = 4 - maskTypeProcedure = None + show_plots = False + do_symmetrisation = True + subtract_calculated_fse_from_data = False + range_for_rebinning_in_y_space = "-20, 0.5, 20" # Needs to be symetric + fitting_model = "SINGLE_GAUSSIAN" + run_minos = True + do_global_fit = None + number_of_global_fit_groups = 4 + mask_zeros_with = None if (__name__ == "__main__") or (__name__ == "mantidqt.widgets.codeeditor.execution"): diff --git a/tests/data/analysis/inputs/yspace_gauss_test.py b/tests/data/analysis/inputs/yspace_gauss_test.py index 72cb1cc..2ec5c95 100644 --- a/tests/data/analysis/inputs/yspace_gauss_test.py +++ b/tests/data/analysis/inputs/yspace_gauss_test.py @@ -7,7 +7,7 @@ BackwardAnalysisInputs.fit_in_y_space = False ForwardAnalysisInputs.run_this_scattering_type = True BackwardAnalysisInputs.run_this_scattering_type = False -YSpaceFitInputs.fitModel = "SINGLE_GAUSSIAN" +YSpaceFitInputs.fitting_model = "gauss" if (__name__ == "__main__") or (__name__ == "mantidqt.widgets.codeeditor.execution"): import mvesuvio diff --git a/tests/data/analysis/inputs/yspace_gc_test.py b/tests/data/analysis/inputs/yspace_gc_test.py index d5fa66c..d2b2a7d 100644 --- a/tests/data/analysis/inputs/yspace_gc_test.py +++ b/tests/data/analysis/inputs/yspace_gc_test.py @@ -8,8 +8,8 @@ BackwardAnalysisInputs.fit_in_y_space = False ForwardAnalysisInputs.run_this_scattering_type = True BackwardAnalysisInputs.run_this_scattering_type = False -YSpaceFitInputs.fitModel = "GC_C4_C6" -YSpaceFitInputs.symmetrisationFlag = False +YSpaceFitInputs.fitting_model = "gcc4c6" +YSpaceFitInputs.do_symmetrisation = False if (__name__ == "__main__") or (__name__ == "mantidqt.widgets.codeeditor.execution"):