Skip to content
Snippets Groups Projects
Commit 3174077c authored by Federico Montesino Pouzols's avatar Federico Montesino Pouzols Committed by GitHub
Browse files

Merge pull request #17186 from mantidproject/17091_Remove_Quest

Indirect Bayes Quest - Remove Quest code
parents 6bd025cf 18570226
No related branches found
No related tags found
No related merge requests found
#pylint: disable=no-init
from mantid.api import PythonAlgorithm, AlgorithmFactory
from mantid.kernel import StringListValidator, StringMandatoryValidator
from mantid.simpleapi import *
from mantid import config, logger, mtd
import os
class Quest(PythonAlgorithm):
def category(self):
return "Workflow\\MIDAS"
def summary(self):
return "This is a variation of the stretched exponential option of Quasi."
def PyInit(self):
self.declareProperty(name='InputType',defaultValue='File',validator=StringListValidator(['File','Workspace']),
doc='Origin of data input - File (.nxs) or Workspace')
self.declareProperty(name='Instrument',defaultValue='iris',validator=StringListValidator(['irs','iris','osi','osiris']),
doc='Instrument')
self.declareProperty(name='Analyser',defaultValue='graphite002',validator=StringListValidator(['graphite002','graphite004']),
doc='Analyser & reflection')
self.declareProperty(name='SamNumber',defaultValue='',validator=StringMandatoryValidator(), doc='Sample run number')
self.declareProperty(name='ResInputType',defaultValue='File',validator=StringListValidator(['File','Workspace']),
doc='Origin of res input - File (_res.nxs) or Workspace')
self.declareProperty(name='ResNumber',defaultValue='',validator=StringMandatoryValidator(), doc='Resolution run number')
self.declareProperty(name='ResNormInputType',defaultValue='File',validator=StringListValidator(['File','Workspace']),
doc='Origin of ResNorm input - File (_red.nxs) or Workspace')
self.declareProperty(name='ResNormNumber',defaultValue='',validator=StringMandatoryValidator(), doc='ResNorm run number')
self.declareProperty(name='ElasticOption',defaultValue=True, doc='Include elastic peak in fit')
self.declareProperty(name='BackgroundOption',defaultValue='Sloping',validator=StringListValidator(['Sloping','Flat','Zero']),
doc='Form of background to fit')
self.declareProperty(name='EnergyMin', defaultValue=-0.5, doc='Minimum energy for fit. Default=-0.5')
self.declareProperty(name='EnergyMax', defaultValue=0.5, doc='Maximum energy for fit. Default=0.5')
self.declareProperty(name='SamBinning', defaultValue=1, doc='Binning value(integer) for sample. Default=1')
self.declareProperty(name='NumberSigma', defaultValue=50, doc='Number of sigma values. Default=50')
self.declareProperty(name='NumberBeta', defaultValue=30, doc='Number of beta values. Default=30')
self.declareProperty(name='Sequence',defaultValue=True, doc='Switch Sequence Off/On')
self.declareProperty(name='Plot',defaultValue='None',validator=StringListValidator(['None','Sigma','Beta','All']),
doc='Plot options')
self.declareProperty(name='Verbose',defaultValue=True, doc='Switch Verbose Off/On')
self.declareProperty(name='Save',defaultValue=False, doc='Switch Save result to nxs file Off/On')
#pylint: disable=too-many-locals,too-many-branches
def PyExec(self):
from IndirectImport import run_f2py_compatibility_test, is_supported_f2py_platform
if is_supported_f2py_platform():
import IndirectBayes as Main
run_f2py_compatibility_test()
self.log().information('Quest input')
inType = self.getPropertyValue('InputType')
prefix = self.getPropertyValue('Instrument')
ana = self.getPropertyValue('Analyser')
sam = self.getPropertyValue('SamNumber')
rinType = self.getPropertyValue('ResInputType')
res = self.getPropertyValue('ResNumber')
elastic = self.getProperty('ElasticOption').value
bgd = self.getPropertyValue('BackgroundOption')
emin = self.getPropertyValue('EnergyMin')
emax = self.getPropertyValue('EnergyMax')
nbin = self.getPropertyValue('SamBinning')
nbins = [nbin, 1]
nbet = self.getProperty('NumberBeta').value
nsig = self.getProperty('NumberSigma').value
nbs = [nbet, nsig]
sname = prefix+sam+'_'+ana + '_red'
rname = prefix+res+'_'+ana + '_res'
erange = [float(emin), float(emax)]
if elastic:
o_el = 1
else:
o_el = 0
fitOp = [o_el, bgd, 0, 0]
loopOp = self.getProperty('Sequence').value
verbOp = self.getProperty('Verbose').value
plotOp = self.getPropertyValue('Plot')
saveOp = self.getProperty('Save').value
workdir = config['defaultsave.directory']
if not os.path.isdir(workdir):
workdir = os.getcwd()
logger.information('Default Save directory is not set. Defaulting to current working Directory: ' + workdir)
if inType == 'File':
spath = os.path.join(workdir, sname+'.nxs') # path name for sample nxs file
LoadNexusProcessed(Filename=spath, OutputWorkspace=sname)
Smessage = 'Sample from File : '+spath
else:
Smessage = 'Sample from Workspace : '+sname
if rinType == 'File':
rpath = os.path.join(workdir, rname+'.nxs') # path name for res nxs file
LoadNexusProcessed(Filename=rpath, OutputWorkspace=rname)
Rmessage = 'Resolution from File : '+rpath
else:
Rmessage = 'Resolution from Workspace : '+rname
if verbOp:
logger.notice(Smessage)
logger.notice(Rmessage)
Main.QuestRun(sname,rname,nbs,erange,nbins,fitOp,loopOp,plotOp,saveOp)
AlgorithmFactory.subscribe(Quest) # Register algorithm with Mantid
......@@ -101,27 +101,30 @@ class QuestTest(stresstesting.MantidStressTest):
return not platform.system() == "Windows"
def runTest(self):
import IndirectBayes as Main
nbins = [1, 1]
nbs = [50, 30]
sname = 'irs26176_graphite002_red'
rname = 'irs26173_graphite002_res'
erange = [-0.5, 0.5]
fitOp = [True, 'Sloping', False, False] #elastic, background, width, resnorm
loopOp = False
plotOp = 'None'
saveOp = False
spath = sname+'.nxs' # path name for sample nxs file
LoadNexusProcessed(Filename=spath, OutputWorkspace=sname)
sample = LoadNexusProcessed(Filename=spath, OutputWorkspace=sname)
rpath = rname+'.nxs' # path name for res nxs file
LoadNexusProcessed(Filename=rpath, OutputWorkspace=rname)
Main.QuestRun(sname,rname,nbs,erange,nbins,fitOp,loopOp,plotOp,saveOp)
res = LoadNexusProcessed(Filename=rpath, OutputWorkspace=rname)
BayesStretch(SampleWorkspace=sample,
ResolutionWorkspace=res,
EMin=-0.5,
EMax=0.5,
SampleBins=1,
Elastic=True,
Background='Sloping',
NumberSigma=30,
NumberBeta=50,
Loop=False,
OutputWorkspaceFit='fit_group',
OutputWorkspaceContour='contour_group')
def validate(self):
self.tolerance = 1e-1
self.disableChecking.append('SpectraMap')
return 'irs26176_graphite002_Qst_Fit','ISISIndirectBayes_QuestTest.nxs'
return 'irs26176_graphite002_Stretch_Fit','ISISIndirectBayes_QuestTest.nxs'
def cleanup(self):
filenames = ['irs26176_graphite002_Qst.lpt','irs26176_graphite002_Qss.ql2',
......
.. algorithm::
.. summary::
.. alias::
.. properties::
Description
-----------
This is a variation of the stretched exponential option of
`Quasi <http://www.mantidproject.org/IndirectBayes:Quasi>`__. For each spectrum a fit is performed
for a grid of :math:`\beta` and :math:`\sigma` values. The distribution of goodness of fit values
is plotted.
This routine was originally part of the MODES package. Note that this algorithm
uses F2Py and is currently only supported on windows.
Usage
-----
**Example - a basic example using Quest to fit a reduced workspace.**
.. code-block:: python
sam = LoadNexusProcessed("irs26173_graphite002_red.nxs", OutputWorkspace='irs26173_graphite002_red')
res = LoadNexusProcessed("irs26173_graphite002_res.nxs", OutputWorkspace='irs26173_graphite002_res')
van = LoadNexusProcessed("irs26176_graphite002_red.nxs", OutputWorkspace='irs26176_graphite002_red')
ResNorm(VanNumber='26176', ResNumber='26173', InputType='File', ResInputType='File', Instrument='irs', Analyser='graphite002', Plot='None', Version=1)
Quest(SamNumber='26176', ResNumber='26173', ResNormNumber='26176', InputType='File', ResInputType='File', ResNormInputType='Workspace', Instrument='irs', Analyser='graphite002')
.. categories::
.. sourcelink::
......@@ -12,7 +12,6 @@ if is_supported_f2py_platform():
QLr = import_f2py("QLres")
QLd = import_f2py("QLdata")
Qse = import_f2py("QLse")
Que = import_f2py("Quest")
resnorm = import_f2py("ResNorm")
else:
unsupported_message()
......@@ -68,203 +67,6 @@ def GetXYE(inWS,n,array_len):
E=PadArray(Ein,array_len)
return N,X,Y,E
# Quest programs
def CheckBetSig(nbs):
Nsig = int(nbs[1])
if Nsig == 0:
raise ValueError('Number of sigma points is Zero')
if Nsig > 200:
raise ValueError('Max number of sigma points is 200')
Nbet = int(nbs[0])
if Nbet == 0:
raise ValueError('Number of beta points is Zero')
if Nbet > 200:
raise ValueError('Max number of beta points is 200')
return Nbet,Nsig
def QuestRun(samWS,resWS,nbs,erange,nbins,Fit,Loop,Plot,Save):
StartTime('Quest')
#expand fit options
elastic, background, width, res_norm = Fit
#convert true/false to 1/0 for fortran
o_el = 1 if elastic else 0
o_w1 = 1 if width else 0
o_res = 1 if res_norm else 0
#fortran code uses background choices defined using the following numbers
if background == 'Sloping':
o_bgd = 2
elif background == 'Flat':
o_bgd = 1
elif background == 'Zero':
o_bgd = 0
fitOp = [o_el, o_bgd, o_w1, o_res]
workdir = config['defaultsave.directory']
if not os.path.isdir(workdir):
workdir = os.getcwd()
logger.information('Default Save directory is not set. Defaulting to current working Directory: ' + workdir)
array_len = 4096 # length of array in Fortran
CheckXrange(erange,'Energy')
nbin,nrbin = nbins[0],nbins[1]
logger.information('Sample is ' + samWS)
logger.information('Resolution is ' + resWS)
CheckAnalysers(samWS,resWS)
nsam,ntc = CheckHistZero(samWS)
if Loop != True:
nsam = 1
efix = getEfixed(samWS)
theta,Q = GetThetaQ(samWS)
nres = CheckHistZero(resWS)[0]
if nres == 1:
prog = 'Qst' # res file
else:
raise ValueError('Stretched Exp ONLY works with RES file')
logger.information(' Number of spectra = '+str(nsam))
logger.information(' Erange : '+str(erange[0])+' to '+str(erange[1]))
fname = samWS[:-4] + '_'+ prog
wrks=os.path.join(workdir, samWS[:-4])
logger.information(' lptfile : ' + wrks +'_Qst.lpt')
lwrk=len(wrks)
wrks.ljust(140,' ')
wrkr=resWS
wrkr.ljust(140,' ')
Nbet,Nsig = nbs[0], nbs[1]
eBet0 = np.zeros(Nbet) # set errors to zero
eSig0 = np.zeros(Nsig) # set errors to zero
rscl = 1.0
Qaxis = ''
for m in range(0,nsam):
logger.information('Group ' +str(m)+ ' at angle '+ str(theta[m]))
nsp = m+1
nout,bnorm,Xdat,Xv,Yv,Ev = CalcErange(samWS,m,erange,nbin)
Ndat = nout[0]
Imin = nout[1]
Imax = nout[2]
Nb,Xb,Yb,_ = GetXYE(resWS,0,array_len)
numb = [nsam, nsp, ntc, Ndat, nbin, Imin, Imax, Nb, nrbin, Nbet, Nsig]
reals = [efix, theta[m], rscl, bnorm]
xsout,ysout,xbout,ybout,zpout=Que.quest(numb,Xv,Yv,Ev,reals,fitOp,\
Xdat,Xb,Yb,wrks,wrkr,lwrk)
dataXs = xsout[:Nsig] # reduce from fixed Fortran array
dataYs = ysout[:Nsig]
dataXb = xbout[:Nbet]
dataYb = ybout[:Nbet]
zpWS = fname + '_Zp' +str(m)
if m > 0:
Qaxis += ','
Qaxis += str(Q[m])
dataXz = []
dataYz = []
dataEz = []
for n in range(0,Nsig):
yfit_list = np.split(zpout[:Nsig*Nbet],Nsig)
dataYzp = yfit_list[n]
dataXz = np.append(dataXz,xbout[:Nbet])
dataYz = np.append(dataYz,dataYzp[:Nbet])
dataEz = np.append(dataEz,eBet0)
CreateWorkspace(OutputWorkspace=zpWS, DataX=dataXz, DataY=dataYz, DataE=dataEz,
Nspec=Nsig, UnitX='MomentumTransfer',
VerticalAxisUnit='MomentumTransfer', VerticalAxisValues=dataXs)
unitx = mtd[zpWS].getAxis(0).setUnit("Label")
unitx.setLabel('beta' , '')
unity = mtd[zpWS].getAxis(1).setUnit("Label")
unity.setLabel('sigma' , '')
if m == 0:
xSig = dataXs
ySig = dataYs
eSig = eSig0
xBet = dataXb
yBet = dataYb
eBet = eBet0
groupZ = zpWS
else:
xSig = np.append(xSig,dataXs)
ySig = np.append(ySig,dataYs)
eSig = np.append(eSig,eSig0)
xBet = np.append(xBet,dataXb)
yBet = np.append(yBet,dataYb)
eBet = np.append(eBet,eBet0)
groupZ = groupZ +','+ zpWS
#create workspaces for sigma and beta
CreateWorkspace(OutputWorkspace=fname+'_Sigma', DataX=xSig, DataY=ySig, DataE=eSig,\
Nspec=nsam, UnitX='', VerticalAxisUnit='MomentumTransfer', VerticalAxisValues=Qaxis)
unitx = mtd[fname+'_Sigma'].getAxis(0).setUnit("Label")
unitx.setLabel('sigma' , '')
CreateWorkspace(OutputWorkspace=fname+'_Beta', DataX=xBet, DataY=yBet, DataE=eBet,\
Nspec=nsam, UnitX='', VerticalAxisUnit='MomentumTransfer', VerticalAxisValues=Qaxis)
unitx = mtd[fname+'_Beta'].getAxis(0).setUnit("Label")
unitx.setLabel('beta' , '')
group = fname + '_Sigma,'+ fname + '_Beta'
fit_workspace = fname+'_Fit'
contour_workspace = fname+'_Contour'
GroupWorkspaces(InputWorkspaces=group,OutputWorkspace=fit_workspace)
GroupWorkspaces(InputWorkspaces=groupZ,OutputWorkspace=contour_workspace)
#add sample logs to the output workspaces
CopyLogs(InputWorkspace=samWS, OutputWorkspace=fit_workspace)
QuestAddSampleLogs(fit_workspace, resWS, background, elastic, erange, nbin, Nsig, Nbet)
CopyLogs(InputWorkspace=samWS, OutputWorkspace=contour_workspace)
QuestAddSampleLogs(contour_workspace, resWS, background, elastic, erange, nbin, Nsig, Nbet)
if Save:
fpath = os.path.join(workdir,fit_workspace+'.nxs')
SaveNexusProcessed(InputWorkspace=fit_workspace, Filename=fpath)
cpath = os.path.join(workdir,contour_workspace+'.nxs')
SaveNexusProcessed(InputWorkspace=contour_workspace, Filename=cpath)
logger.information('Output file for Fit : ' + fpath)
logger.information('Output file for Contours : ' + cpath)
if Plot != 'None' and Loop == True:
QuestPlot(fname,Plot)
EndTime('Quest')
def QuestAddSampleLogs(workspace, res_workspace, background, elastic_peak, e_range, sample_binning, sigma, beta):
energy_min, energy_max = e_range
AddSampleLog(Workspace=workspace, LogName="res_file",
LogType="String", LogText=res_workspace)
AddSampleLog(Workspace=workspace, LogName="background",
LogType="String", LogText=str(background))
AddSampleLog(Workspace=workspace, LogName="elastic_peak",
LogType="String", LogText=str(elastic_peak))
AddSampleLog(Workspace=workspace, LogName="energy_min",
LogType="Number", LogText=str(energy_min))
AddSampleLog(Workspace=workspace, LogName="energy_max",
LogType="Number", LogText=str(energy_max))
AddSampleLog(Workspace=workspace, LogName="sample_binning",
LogType="Number", LogText=str(sample_binning))
AddSampleLog(Workspace=workspace, LogName="sigma",
LogType="Number", LogText=str(sigma))
AddSampleLog(Workspace=workspace, LogName="beta",
LogType="Number", LogText=str(beta))
def QuestPlot(inputWS,Plot):
if Plot == 'Sigma' or Plot == 'All':
MTD_PLOT.importMatrixWorkspace(inputWS+'_Sigma').plotGraph2D()
if Plot == 'Beta' or Plot == 'All':
MTD_PLOT.importMatrixWorkspace(inputWS+'_Beta').plotGraph2D()
# ResNorm programs
def ResNormRun(vname,rname,erange,nbin,Plot='None',Save=False):
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment