Skip to content

Commit

Permalink
Uncertainty for Arbitrary Expressions (AdvancedPhotonSource#64)
Browse files Browse the repository at this point in the history 
* start work on entering an expression to compute uncertainties from GSAS-II parameters

* finish computation of user-supplied expression w/s.u.; minor bug fix in expression editor
  • Loading branch information
@briantoby
briantoby authored Sep 6, 2024
1 parent 4514980 commit ad46b21
Show file tree
Hide file tree
Showing 3 changed files with 136 additions and 13 deletions.
131 changes: 131 additions & 0 deletions GSASII/GSASIIdataGUI.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -802,6 802,8 @@ def _Add_CalculateMenuItems(self,parent):
self.Bind(wx.EVT_MENU, self.OnRefinePartials, id=item.GetId())
item = parent.Append(wx.ID_ANY,'&Parameter Impact\tCTRL I','Perform a derivative calculation')
self.Bind(wx.EVT_MENU, self.OnDerivCalc, id=item.GetId())
item = parent.Append(wx.ID_ANY,'Evaluate expression and s.u.','Perform uncertainty analysis on an expression of GSAS-II parameters')
self.Bind(wx.EVT_MENU, self.OnExpressionCalc, id=item.GetId())

item = parent.Append(wx.ID_ANY,'Save partials as csv','Save the computed partials as a csv file')
self.Refine.append(item)
Expand Down Expand Up @@ -5360,6 5362,135 @@ def OnDerivCalc(self,event):
G2G.G2ScrolledGrid(self,'Parameter Impact Results','Impact Results',tbl,colLbls,colTypes,
maxSize=(700,400),comment=' Cite: B.H. Toby, IUCrJ, to be published')

def OnExpressionCalc(self,event):
'''Compute an arbitrary expression (supplied by user) as well as the
(statistical) standard uncertainty on that expression.
Uses the :class:`GSASIIexprGUI.ExpressionDialog` to obtain
an expression which is evaluated using the
:class:`GSASIIobj.ExpressionObj` capability. Then the derivative of
the expression is computed numerically for every parameter in the
covariance matrix. Finally the derivative list is used to find
the s.u. on the expression using Ted Prince's method.
'''
import GSASIIexprGUI as G2exG
import GSASIIstrMath as G2stMth
def extendChanges():
'''Propagate changes due to constraint and rigid bodies
from varied parameters to dependent parameters
'''
# apply constraints
G2mv.Dict2Map(prms)
# apply rigid body constraints
G2stMth.ApplyRBModels(prms,Phases,rigidbodyDict)
# apply shifts to atoms
for dk in prms:
if not '::dA' in dk: continue
if prms[dk] == 0: continue
k = dk.replace('::dA','::A')
prms[k] = prms[dk]
prms[dk] = 0
# TODO: debug: find what has changed
#print(var,[k for k in prms if prms[k] != parmValDict[k]])
#print(var,[prms[k] for k in prms if prms[k] != parmValDict[k]])
# end debug
def striphist(var,insChar=''):
'strip a histogram number from a var name'
sv = var.split(':')
if len(sv) <= 1: return var
if sv[1]:
sv[1] = insChar
return ':'.join(sv)


name = histNames[0]
data = seqDict

parmDict,varyList = self.MakeLSParmDict()
Histograms,Phases = self.GetUsedHistogramsAndPhasesfromTree()
if not len(Phases) or not len(Histograms):
print('Expression computation not possible without definition of phases and histograms')
return
Constraints = self.GPXtree.GetItemPyData(GetGPXtreeItemId(self,self.root,'Constraints'))
errmsg,warnmsg = G2cnstG.CheckConstraints(self,Phases,Histograms,
Constraints,[],copy.copy(varyList))
item = GetGPXtreeItemId(self,self.root,'Covariance')
covData = self.GPXtree.GetItemPyData(item)
covData['covMatrix'] = covData.get('covMatrix',[])
parmValDict = {i:parmDict[i][0] for i in parmDict} # dict w/parm values only
G2mv.Map2Dict(parmValDict,G2mv.saveVaryList)
rigidbodyDict = self.GPXtree.GetItemPyData(GetGPXtreeItemId(self,self.root,'Rigid bodies'))

if self.testSeqRefineMode():
seqDict = self.GPXtree.GetItemPyData(GetGPXtreeItemId(self,self.root,'Sequential results'))
histNames = [h for h in self.testSeqRefineMode() if h in seqDict]
if len(histNames) == 0:
print('no histograms')
return
# wildcard the histogram values
for k in sorted(parmValDict,reverse=True):
ks = striphist(k,'*')
if ks != k:
parmValDict[ks] = parmValDict[k]
del parmValDict[k]

dlg = G2exG.ExpressionDialog(self,parmValDict,
header="Evaluate an expression of GSAS-II parameters",
VarLabel = "Expression",
fit=False,wildCard=self.testSeqRefineMode())
exprobj = dlg.Show(True)
if not exprobj: return

if not self.testSeqRefineMode():
histNames = [0]
for h in histNames:
prfx = ''
if self.testSeqRefineMode():
parmValDict = seqDict[h]['parmDict']
covMatrix = seqDict[h]['covMatrix']
CvaryList = seqDict[h]['varyList']
Csig = seqDict[h]['sig']
prfx = h ': '
else:
covMatrix = covData['covMatrix']
CvaryList = covData['varyList']
Csig = covData['sig']
# evaluate the expression
try:
calcobj = G2obj.ExpressionCalcObj(exprobj)
calcobj.SetupCalc(parmValDict)
value = calcobj.EvalExpression()
except:
continue
# evaluate the derivatives of the expression with respect to
# the varied parameters
if len(covMatrix) > 0:
derivs = []
for var,sig in zip(CvaryList,Csig):
prms = copy.copy(parmValDict)
if sig == 0:
derivs.append(0.0)
continue
prms[var] = sig
extendChanges()
calcobj.SetupCalc(prms)
valP = calcobj.EvalExpression()
prms[var] -= 2*sig
extendChanges()
calcobj.SetupCalc(prms)
valM = calcobj.EvalExpression()
derivs.append((valP-valM)/(2*sig))
Avec = np.array(derivs)
sig = np.sqrt(np.inner(Avec.T,np.inner(covMatrix,Avec)))
else:
sig = -1
if sig > 0:
print(f'{prfx}{exprobj.expression!r} value = {G2mth.ValEsd(value,sig)}')
elif sig == 0:
print(f'{prfx}{exprobj.expression!r} value = {G2mth.ValEsd(value,-abs(value)/10000)} (0 uncertainty)')
else:
print(f'{prfx}{exprobj.expression!r} value = {G2mth.ValEsd(value,-abs(value)/10000)} (no s.u.: no covariance available')

def OnRefine(self,event):
'''Perform a single refinement or a sequential refinement (depending on controls setting)
Called from the Calculate/Refine menu.
Expand Down
6 changes: 3 additions & 3 deletions GSASII/GSASIIexprGUI.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -38,7 38,7 @@ def IndexParmDict(parmDict,wildcard):
2 is histogram/phase parms, 3 is histogram parms and 4 are global parameters
'''
prex = re.compile('[0-9] ::.*')
hrex = re.compile(':[0-9] :.*')
hrex = re.compile(':[0-9\*] :.*')
parmLists = {}
for i in (1,2,3,4):
parmLists[i] = []
Expand Down Expand Up @@ -674,9 674,9 @@ def Repaint(self,exprObj):
self.depVarDict.get(self.dependentVar,'?')
)
if self.wildCard:
msg = " with first defined values"
msg = " with first defined value(s)"
else:
msg = " using parameter values"
msg = " using parameter value(s)"
self.setEvalResult(f"Expression evaluates to: {s}{depVal}{msg}")
self.OKbtn.Enable()
if self.ExtraBtn: self.ExtraBtn.Enable()
Expand Down
12 changes: 2 additions & 10 deletions GSASII/GSASIIobj.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1913,12 1913,8 @@ def __init__(self,exprObj):
'''Standard error evaluation where supplied by the evaluator
'''
# Patch: for old-style expressions with a (now removed step size)
if '2' in platform.python_version_tuple()[0]:
basestr = basestring
else:
basestr = str
for v in self.eObj.assgnVars:
if not isinstance(self.eObj.assgnVars[v], basestr):
if not isinstance(self.eObj.assgnVars[v], str):
self.eObj.assgnVars[v] = self.eObj.assgnVars[v][0]
self.parmDict = {}
'''A copy of the parameter dictionary, for distance and angle computation
Expand All @@ -1936,13 1932,9 @@ def SetupCalc(self,parmDict):

# look at first value in parmDict to determine its type
parmsInList = True
if '2' in platform.python_version_tuple()[0]:
basestr = basestring
else:
basestr = str
for key in parmDict:
val = parmDict[key]
if isinstance(val, basestr):
if isinstance(val, str):
parmsInList = False
break
try: # check if values are in lists
Expand Down

0 comments on commit ad46b21

Please sign in to comment.