Skip to content

Commit

Permalink
Split yspace class into back and front
Browse files Browse the repository at this point in the history
This allows backward and forward fitting to have their own fitting parameters.
Updated rest of files to reflect this change.
  • Loading branch information
GuiMacielPereira committed Jan 15, 2025
1 parent ade7d8f commit 35a22de
Show file tree
Hide file tree
Showing 3 changed files with 60 additions and 38 deletions.
63 changes: 32 additions & 31 deletions src/mvesuvio/analysis_fitting.py
Original file line number Diff line number Diff line change
Expand Up @@ -15,39 +15,39 @@
repoPath = Path(__file__).absolute().parent # Path to the repository


def fitInYSpaceProcedure(yFitIC, IC, wsTOF):
def fitInYSpaceProcedure(IC, wsTOF):

try:
wsTOF = mtd[wsTOF]
except KeyError:
print(f"Workspace to fit {wsTOF} not found.")
return

wsResSum, wsRes = calculateMantidResolutionFirstMass(IC, yFitIC, wsTOF)
wsResSum, wsRes = calculateMantidResolutionFirstMass(IC, wsTOF)

wsTOFMass0 = subtract_profiles_except_lightest(IC, wsTOF)

if yFitIC.subtract_calculated_fse_from_data:
if IC.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)
wsJoY, wsJoYAvg = ySpaceReduction(wsTOFMass0, IC)

if yFitIC.do_symmetrisation:
if IC.do_symmetrisation:
wsJoYAvg = symmetrizeWs(wsJoYAvg)

fitProfileMinuit(yFitIC, wsJoYAvg, wsResSum)
fitProfileMantidFit(yFitIC, wsJoYAvg, wsResSum)
fitProfileMinuit(IC, wsJoYAvg, wsResSum)
fitProfileMantidFit(IC, wsJoYAvg, wsResSum)

printYSpaceFitResults(wsJoYAvg.name())

yfitResults = ResultsYFitObject(IC, yFitIC, wsTOF.name(), wsJoYAvg.name())
yfitResults = ResultsYFitObject(IC, wsTOF.name(), wsJoYAvg.name())
yfitResults.save()

if yFitIC.do_global_fit:
runGlobalFit(wsJoY, wsRes, IC, yFitIC)
if IC.do_global_fit:
runGlobalFit(wsJoY, wsRes, IC)

save_workspaces(yFitIC)
save_workspaces(IC)
return yfitResults


Expand All @@ -67,7 +67,7 @@ def find_ws_name_fse_first_mass(ic):
return ws_names_fse[np.argmin(ws_masses)]


def calculateMantidResolutionFirstMass(IC, yFitIC, ws):
def calculateMantidResolutionFirstMass(IC, ws):
mass = IC.masses[0]

resName = ws.name() + "_Resolution"
Expand All @@ -77,7 +77,7 @@ def calculateMantidResolutionFirstMass(IC, yFitIC, ws):
)
Rebin(
InputWorkspace="tmp",
Params=yFitIC.range_for_rebinning_in_y_space,
Params=IC.range_for_rebinning_in_y_space,
OutputWorkspace="tmp",
)

Expand Down Expand Up @@ -124,19 +124,20 @@ def switchFirstTwoAxis(A):
return np.stack(np.split(A, len(A), axis=0), axis=2)[0]


def ySpaceReduction(wsTOF, mass0, yFitIC, ic):
def ySpaceReduction(wsTOF, ic):
"""Seperate procedures depending on masking specified."""

mass0 = ic.masses[0]
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.range_for_rebinning_in_y_space
rebinPars = ic.range_for_rebinning_in_y_space

if np.any(np.all(wsTOF.extractY() == 0, axis=0)): # Masked columns present
if yFitIC.mask_zeros_with == "nan":
if ic.mask_zeros_with == "nan":
# Build special workspace to store accumulated points
wsJoY = convertToYSpace(wsTOF, mass0)
xp = buildXRangeFromRebinPars(yFitIC)
xp = buildXRangeFromRebinPars(ic)
wsJoYB = dataXBining(
wsJoY, xp
) # Unusual ws with several dataY points per each dataX point
Expand Down Expand Up @@ -1319,7 +1320,7 @@ def printYSpaceFitResults(wsJoYName):


class ResultsYFitObject:
def __init__(self, ic, yFitIC, wsFinalName, wsYSpaceAvgName):
def __init__(self, ic, wsFinalName, wsYSpaceAvgName):
# Extract most relevant information from ws
wsFinal = mtd[wsFinalName]
wsResSum = mtd[wsFinalName + "_Resolution_Sum"]
Expand All @@ -1339,19 +1340,19 @@ def __init__(self, ic, yFitIC, wsFinalName, wsYSpaceAvgName):
poptList = []
perrList = []
try:
wsFitMinuit = mtd[wsJoYAvg.name() + "_minuit_" + yFitIC.fitting_model + "_Parameters"]
wsFitMinuit = mtd[wsJoYAvg.name() + "_minuit_" + ic.fitting_model + "_Parameters"]
poptList.append(wsFitMinuit.column("Value"))
perrList.append(wsFitMinuit.column("Error"))
except:
pass
try:
wsFitLM = mtd[wsJoYAvg.name() + "_lm_" + yFitIC.fitting_model + "_Parameters"]
wsFitLM = mtd[wsJoYAvg.name() + "_lm_" + ic.fitting_model + "_Parameters"]
poptList.append(wsFitLM.column("Value"))
perrList.append(wsFitLM.column("Error"))
except:
pass
try:
wsFitSimplex = mtd[wsJoYAvg.name() + "_simplex_" + yFitIC.fitting_model + "_Parameters"]
wsFitSimplex = mtd[wsJoYAvg.name() + "_simplex_" + ic.fitting_model + "_Parameters"]
poptList.append(wsFitSimplex.column("Value"))
perrList.append(wsFitSimplex.column("Error"))
except:
Expand All @@ -1371,7 +1372,7 @@ def __init__(self, ic, yFitIC, wsFinalName, wsYSpaceAvgName):
self.perr = perr

self.savePath = ic.ySpaceFitSavePath
self.fitting_model= yFitIC.fitting_model
self.fitting_model= ic.fitting_model

def save(self):
np.savez(
Expand All @@ -1387,23 +1388,23 @@ def save(self):
)


def runGlobalFit(wsYSpace, wsRes, IC, yFitIC):
def runGlobalFit(wsYSpace, wsRes, IC):
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
)

idxList = groupDetectors(instrPars, yFitIC)
idxList = groupDetectors(instrPars, IC)
dataX, dataY, dataE, dataRes = avgWeightDetGroups(
dataX, dataY, dataE, dataRes, idxList, yFitIC
dataX, dataY, dataE, dataRes, idxList, IC
)

if yFitIC.do_symmetrisation:
if IC.do_symmetrisation:
dataY, dataE = weightedSymArr(dataY, dataE)

model, defaultPars, sharedPars = selectModelAndPars(yFitIC.fitting_model)
model, defaultPars, sharedPars = selectModelAndPars(IC.fitting_model)

totCost = 0
for i, (x, y, yerr, res) in enumerate(zip(dataX, dataY, dataE, dataRes)):
Expand All @@ -1424,12 +1425,12 @@ def runGlobalFit(wsYSpace, wsRes, IC, yFitIC):
for i in range(len(dataY)): # Set limits for unshared parameters
m.limits["A" + str(i)] = (0, np.inf)

if yFitIC.fitting_model== "doublewell":
if IC.fitting_model== "doublewell":
m.limits["d"] = (0, np.inf) # Shared parameters
m.limits["R"] = (0, np.inf)

t0 = time.time()
if yFitIC.fitting_model== "gauss":
if IC.fitting_model== "gauss":
m.simplex()
m.migrad()

Expand Down Expand Up @@ -1477,8 +1478,8 @@ def constr(*pars):

create_table_for_global_fit_parameters(wsYSpace.name(), m, chi2)

if yFitIC.show_plots:
plotGlobalFit(dataX, dataY, dataE, m, totCost, wsYSpace.name(), yFitIC)
if IC.show_plots:
plotGlobalFit(dataX, dataY, dataE, m, totCost, wsYSpace.name(), IC)

# Pass into array to store values in variable
return np.array(m.values), np.array(m.errors)
Expand Down
29 changes: 25 additions & 4 deletions src/mvesuvio/config/analysis_inputs.py
Original file line number Diff line number Diff line change
Expand Up @@ -11,8 +11,8 @@ class SampleParameters:

@dataclass
class BackwardAnalysisInputs(SampleParameters):
run_this_scattering_type = False
fit_in_y_space = False
run_this_scattering_type = True
fit_in_y_space = True

runs = "43066-43076"
empty_runs = "41876-41923"
Expand Down Expand Up @@ -45,6 +45,29 @@ class BackwardAnalysisInputs(SampleParameters):
do_gamma_correction = False


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"


@dataclass
class ForwardAnalysisInputs(SampleParameters):
run_this_scattering_type = True
Expand Down Expand Up @@ -81,8 +104,6 @@ class ForwardAnalysisInputs(SampleParameters):
do_gamma_correction = True


@dataclass
class YSpaceFitInputs:
show_plots = True
do_symmetrisation = False
subtract_calculated_fse_from_data = True
Expand Down
6 changes: 3 additions & 3 deletions src/mvesuvio/main/run_routine.py
Original file line number Diff line number Diff line change
Expand Up @@ -32,7 +32,6 @@ def setup(self):

self.bckwd_ai = ai.BackwardAnalysisInputs
self.fwd_ai = ai.ForwardAnalysisInputs
self.yFitIC = ai.YSpaceFitInputs

# Names of workspaces to check if they exist to skip analysis
self.ws_to_fit_y_space = []
Expand Down Expand Up @@ -63,7 +62,8 @@ def setup(self):
# TODO: sort out yfit inputs
figSavePath = self.experiment_path / "figures"
figSavePath.mkdir(exist_ok=True)
self.yFitIC.figSavePath = figSavePath
self.fwd_ai.figSavePath = figSavePath
self.bckwd_ai.figSavePath = figSavePath


def import_from_inputs(self):
Expand Down Expand Up @@ -95,7 +95,7 @@ def run(self):

def runAnalysisFitting(self):
for wsName, i_cls in zip(self.ws_to_fit_y_space, self.classes_to_fit_y_space):
self.fitting_result = fitInYSpaceProcedure(self.yFitIC, i_cls, wsName)
self.fitting_result = fitInYSpaceProcedure(i_cls, wsName)
return


Expand Down

0 comments on commit 35a22de

Please sign in to comment.