Changeset 1335 for trunk/GSASIIgrid.py


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

Redo pseudovar computation to handle constraint-generated vars

File:
1 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
Note: See TracChangeset for help on using the changeset viewer.