Changeset 1335


Ignore:
Timestamp:
May 7, 2014 2:38:30 PM (8 years ago)
Author:
toby
Message:

Redo pseudovar computation to handle constraint-generated vars

Location:
trunk
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIgrid.py

    r1323 r1335  
    39193919        return parmDict,modVaryList,data[name]['sig']
    39203920
    3921     def EvalPSvarDeriv(calcobj,parmDict,var,step):
     3921    def UpdateParmDict(parmDict):
     3922        '''generate the atom positions and the direct & reciprocal cell values,
     3923        because they might be needed to evaluate the pseudovar
     3924        '''
     3925        Ddict = dict(zip(['D11','D22','D33','D12','D13','D23'],
     3926                         ['A'+str(i) for i in range(6)])
     3927                     )
     3928        delList = []
     3929        phaselist = []
     3930        for item in parmDict:
     3931            if ':' not in item: continue
     3932            key = item.split(':')
     3933            if len(key) < 3: continue
     3934            # remove the dA[xyz] terms, they would only bring confusion
     3935            if key[2].startswith('dA'):
     3936                delList.append(item)
     3937            # compute and update the corrected reciprocal cell terms using the Dij values
     3938            elif key[2] in Ddict:
     3939                if key[0] not in phaselist: phaselist.append(key[0])
     3940                akey = key[0]+'::'+Ddict[key[2]]
     3941                parmDict[akey] += parmDict[item]
     3942                delList.append(item)
     3943        for item in delList:
     3944            del parmDict[item]               
     3945        for i in phaselist:
     3946            pId = int(i)
     3947            # apply cell symmetry
     3948            A,zeros = G2stIO.cellFill(str(pId)+'::',SGdata[pId],parmDict,zeroDict[pId])
     3949            # convert to direct cell & add the unique terms to the dictionary
     3950            for i,val in enumerate(G2lat.A2cell(A)):
     3951                if i in uniqCellIndx[pId]:
     3952                    lbl = str(pId)+'::'+cellUlbl[i]
     3953                    parmDict[lbl] = val
     3954            lbl = str(pId)+'::'+'vol'
     3955            parmDict[lbl] = G2lat.calc_V(A)
     3956        return parmDict
     3957
     3958    def EvalPSvarDeriv(calcobj,parmDict,var,ESD):
    39223959        '''Evaluate an expression derivative with respect to a
    3923         GSAS-II variable name.
     3960        GSAS-II variable name.
     3961
     3962        Note this likely could be faster if the loop over calcobjs were done
     3963        inside after the Dict was created.
    39243964        '''
    3925         varList = [var]
    3926         origVals = [parmDict[var]]
    3927         minusVals = [parmDict[var]-step]
    3928         plusVals = [parmDict[var]+step]
    3929         if var in Rcelldict:
    3930             pId = var.split(':')[0] # get phase number
    3931             pfx = pId + '::'
    3932             pId = int(pId)
    3933             Rcell = Rcelldict.copy()
    3934             Rcell.update({v:parmDict[v] for v in Rcelldict if v in parmDict})
    3935             Rcell[var] = parmDict[var] - step
    3936             Aminus,zeros = G2stIO.cellFill(pfx,SGdata[pId],Rcell,zeroDict[pId])
    3937             Rcell[var] =  parmDict[var] + step
    3938             Aplus,zeros = G2stIO.cellFill(pfx,SGdata[pId],Rcell,zeroDict[pId])
    3939             for i,(mval,pval) in enumerate(zip(G2lat.A2cell(Aminus),G2lat.A2cell(Aplus))):
    3940                 if i in uniqCellIndx[pId]:
    3941                     lbl = pfx+cellUlbl[i]
    3942                     varList.append(lbl)
    3943                     origVals.append(parmDict[lbl])
    3944                     minusVals.append(mval)
    3945                     plusVals.append(pval)
    3946             lbl = pfx+'vol'
    3947             varList.append(lbl)
    3948             origVals.append(parmDict[lbl])
    3949             minusVals.append(G2lat.calc_V(Aminus))
    3950             plusVals.append(G2lat.calc_V(Aplus))
    3951 
    3952         calcobj.UpdateVars(varList,plusVals)
    3953         pVal = calcobj.EvalExpression()
    3954         calcobj.UpdateVars(varList,minusVals)
    3955         mVal = calcobj.EvalExpression()
    3956         calcobj.UpdateVars(varList,origVals)
    3957         return (pVal - mVal) / (2.*step)
    3958 
     3965        step = ESD/10
     3966        Ddict = dict(zip(['D11','D22','D33','D12','D13','D23'],
     3967                         ['A'+str(i) for i in range(6)])
     3968                     )
     3969        results = []
     3970        phaselist = []
     3971        for incr in step,-step:
     3972            VparmDict = parmDict.copy()           
     3973            # as saved, the parmDict has updated 'A[xyz]' values, but 'dA[xyz]'
     3974            # values are not zeroed: fix that!
     3975            VparmDict.update({item:0.0 for item in parmDict if 'dA' in item})
     3976            VparmDict[var] += incr
     3977            G2mv.Dict2Map(VparmDict,[]) # apply constraints
     3978            # generate the atom positions and the direct & reciprocal cell values now, because they might
     3979            # needed to evaluate the pseudovar
     3980            for item in VparmDict:
     3981                if ':' not in item: continue
     3982                key = item.split(':')
     3983                if len(key) < 3: continue
     3984                # apply any new shifts to atom positions
     3985                if key[2].startswith('dA'):
     3986                    VparmDict[''.join(item.split('d'))] += VparmDict[item]
     3987                    VparmDict[item] = 0.0
     3988                # compute and update the corrected reciprocal cell terms using the Dij values
     3989                if key[2] in Ddict:
     3990                    if key[0] not in phaselist: phaselist.append(key[0])
     3991                    akey = key[0]+'::'+Ddict[key[2]]
     3992                    VparmDict[akey] += VparmDict[item]
     3993            for i in phaselist:
     3994                pId = int(i)
     3995                # apply cell symmetry
     3996                A,zeros = G2stIO.cellFill(str(pId)+'::',SGdata[pId],VparmDict,zeroDict[pId])
     3997                # convert to direct cell & add the unique terms to the dictionary
     3998                for i,val in enumerate(G2lat.A2cell(A)):
     3999                    if i in uniqCellIndx[pId]:
     4000                        lbl = str(pId)+'::'+cellUlbl[i]
     4001                        VparmDict[lbl] = val
     4002                lbl = str(pId)+'::'+'vol'
     4003                VparmDict[lbl] = G2lat.calc_V(A)
     4004            # dict should be fully updated, use it & calculate
     4005            calcobj.SetupCalc(VparmDict)
     4006            results.append(calcobj.EvalExpression())
     4007        return (results[0] - results[1]) / (2.*step)
     4008       
    39594009    def EnableParFitEqMenus():
    39604010        'Enables or disables the Parametric Fit menu items that require existing defs'
     
    43004350        esdList = []
    43014351        for seqnum,name in enumerate(histNames):
    4302             parmDict,modVaryL,ESDL = CreatePSvarDict(seqnum,name)
    4303             calcobj.SetupCalc(parmDict)
     4352            sigs = data[name]['sig']
     4353            G2mv.InitVars()
     4354            parmDict = data[name].get('parmDict')
     4355            if parmDict:
     4356                constraintInfo = data[name].get('constraintInfo')
     4357                groups,parmlist,constrDict,fixedList,ihst = constraintInfo
     4358                varyList = data[name]['varyList']
     4359                parmDict = data[name]['parmDict']
     4360                G2mv.GenerateConstraints(groups,parmlist,varyList,constrDict,fixedList,parmDict,SeqHist=ihst)
     4361                derivs = np.array(
     4362                    [EvalPSvarDeriv(calcobj,parmDict.copy(),var,ESD)
     4363                     for var,ESD in zip(varyList,sigs)]
     4364                    )
     4365                esdList.append(np.sqrt(
     4366                    np.inner(derivs,np.inner(data[name]['covMatrix'],derivs.T))
     4367                    ))
     4368                PSvarDict = parmDict.copy()
     4369                UpdateParmDict(PSvarDict)
     4370                #calcobj.SetupCalc(PSvarDict)
     4371                calcobj.UpdateDict(PSvarDict)
     4372            else:
     4373                PSvarDict,unused,unused = CreatePSvarDict(0,histNames[0])
     4374                calcobj.SetupCalc(PSvarDict)
    43044375            valList.append(calcobj.EvalExpression())
    4305             derivs = np.array(
    4306                 [EvalPSvarDeriv(calcobj,parmDict,var,ESD/10.) for var,ESD in zip(modVaryL,ESDL)]
    4307                 )
    4308             esdList.append(np.sqrt(
    4309                 np.inner(derivs,np.inner(data[name]['covMatrix'],derivs.T))
    4310                 ))
     4376        if not esdList:
     4377            esdList = None
    43114378        colList += [valList]
    43124379        colSigs += [esdList]
     
    43164383
    43174384    # Make dict needed for creating & editing pseudovars (PSvarDict).
    4318     PSvarDict,unused,unused = CreatePSvarDict(0,histNames[0])
     4385    name = histNames[0]
     4386    parmDict = data[name].get('parmDict')
     4387    if parmDict:
     4388        PSvarDict = parmDict.copy()
     4389        UpdateParmDict(PSvarDict)
     4390    else:
     4391        print 'Sequential refinement needs to be rerun to obtain ESDs for PseudoVariables'
     4392        PSvarDict,unused,unused = CreatePSvarDict(0,histNames[0])
    43194393    # Also dicts of dependent (depVarDict) & independent vars (indepVarDict)
    43204394    # for Parametric fitting from the data table
  • trunk/GSASIIobj.py

    r1332 r1335  
    18481848        not evaluation of expression.)
    18491849        '''
    1850        
     1850        # Patch: for old-style expressions with a (now removed step size)
     1851        for v in self.eObj.assgnVars:
     1852            if not isinstance(self.eObj.assgnVars[v], basestring):
     1853                self.eObj.assgnVars[v] = self.eObj.assgnVars[v][0]
     1854
    18511855    def SetupCalc(self,parmDict):
    18521856        '''Do all preparations to use the expression for computation.
     
    19151919        for var,val in zip(varList,valList):
    19161920            self.exprDict[self.lblLookup.get(var,'undefined: '+var)] = val
     1921
     1922    def UpdateDict(self,parmDict):
     1923        '''Update the dict for the expression with values in a dict
     1924        :param list parmDict: a dict of values some of which may be in use here
     1925        '''
     1926        for var in parmDict:
     1927            if var in self.lblLookup:
     1928                self.exprDict[self.lblLookup[var]] = parmDict[var]
    19171929           
    19181930    def EvalExpression(self):
  • trunk/GSASIIstrMain.py

    r1299 r1335  
    293293            groups,parmlist = G2mv.GroupConstraints(constrDict)
    294294            G2mv.GenerateConstraints(groups,parmlist,varyList,constrDict,fixedList,parmDict,SeqHist=ihst)
     295            constraintInfo = (groups,parmlist,constrDict,fixedList,ihst)
    295296        except:
    296297            print ' *** ERROR - your constraints are internally inconsistent ***'
     
    366367        newCellDict = copy.deepcopy(G2stMth.GetNewCellParms(parmDict,varyList))
    367368        newAtomDict = copy.deepcopy(G2stMth.ApplyXYZshifts(parmDict,varyList))
    368         covData = {'variables':result[0],'varyList':varyList,'sig':sig,'Rvals':Rvals,
    369                    'varyListStart':varyListStart,
    370                    'covMatrix':covMatrix,'title':histogram,'newAtomDict':newAtomDict,
    371                    'newCellDict':newCellDict,'depParmDict':depParmDict}
     369        histRefData = {
     370            'variables':result[0],'varyList':varyList,'sig':sig,'Rvals':Rvals,
     371            'varyListStart':varyListStart,
     372            'covMatrix':covMatrix,'title':histogram,'newAtomDict':newAtomDict,
     373            'newCellDict':newCellDict,'depParmDict':depParmDict,
     374            'constraintInfo':constraintInfo,
     375            'parmDict':parmDict}
     376        SeqResult[histogram] = histRefData
    372377        G2stMth.ApplyRBModels(parmDict,Phases,rigidbodyDict,True)
    373378#        G2stIO.SetRigidBodyModels(parmDict,sigDict,rigidbodyDict,printFile)
    374379        G2stIO.SetHistogramPhaseData(parmDict,sigDict,Phases,Histo,ifPrint,printFile)
    375380        G2stIO.SetHistogramData(parmDict,sigDict,Histo,ifPrint,printFile)
    376         SeqResult[histogram] = covData
    377         G2stIO.SetUsedHistogramsAndPhases(GPXfile,Histo,Phases,rigidbodyDict,covData,makeBack)
     381        G2stIO.SetUsedHistogramsAndPhases(GPXfile,Histo,Phases,rigidbodyDict,histRefData,makeBack)
    378382        makeBack = False
    379383        NewparmDict = {}
  • trunk/GSASIIstrMath.py

    r1322 r1335  
    882882   
    883883    :param dict parmDict: parameter dictionary
    884     :param list varyList: list of variables
     884    :param list varyList: list of variables (not used!)
    885885    :returns: newAtomDict - dictionary of new atomic coordinate names & values; key is parameter shift name
    886886
Note: See TracChangeset for help on using the changeset viewer.