diff --git a/GSASII/GSASIIdataGUI.py b/GSASII/GSASIIdataGUI.py index e74aa6b0..d4e8a9ec 100644 --- a/GSASII/GSASIIdataGUI.py +++ b/GSASII/GSASIIdataGUI.py @@ -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) @@ -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. diff --git a/GSASII/GSASIIexprGUI.py b/GSASII/GSASIIexprGUI.py index afd7e3d2..3e221492 100644 --- a/GSASII/GSASIIexprGUI.py +++ b/GSASII/GSASIIexprGUI.py @@ -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] = [] @@ -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() diff --git a/GSASII/GSASIIobj.py b/GSASII/GSASIIobj.py index 29206665..07e50a9d 100644 --- a/GSASII/GSASIIobj.py +++ b/GSASII/GSASIIobj.py @@ -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 @@ -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