Changeset 353 for trunk/GSASIIpwd.py


Ignore:
Timestamp:
Aug 24, 2011 4:07:29 PM (11 years ago)
Author:
vondreele
Message:

modify psvfcj.for to use sh/l only
add derivative routine to pypowder.for
GSASIIgrid.py - new LS controls for derivative type
GSASIIlattice.py - work on docs as per P Jemian
GSASIIpwd.py - major mods to use FCJpsvoight fortran code & derivatives
GSASIIpwdGUI.py - mod to use LS controls

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASIIpwd.py

    r350 r353  
    492492fcjde = fcjde_gen(name='fcjde',shapes='t,s,dx')
    493493               
    494 def getBackground(pfx,parmDict,bakType,xdata):
    495     yb = np.zeros_like(xdata)
    496     if bakType == 'chebyschev':
    497         iBak = 0
    498         while True:
    499             key = pfx+'Back:'+str(iBak)
    500             try:
    501                 yb += parmDict[key]*(xdata-xdata[0])**iBak
    502                 iBak += 1
    503             except KeyError:
    504                 break
    505     return yb
    506 
    507494def getWidths(pos,sig,gam,shl):
    508495    widths = [np.sqrt(sig)/100.,gam/200.]
    … …  
    535522    return intens*Df(xdata)*DX/dx
    536523
     524def getBackground(pfx,parmDict,bakType,xdata):
     525    yb = np.zeros_like(xdata)
     526    if bakType == 'chebyschev':
     527        iBak = 0
     528        while True:
     529            key = pfx+'Back:'+str(iBak)
     530            try:
     531                yb += parmDict[key]*(xdata-xdata[0])**iBak
     532                iBak += 1
     533            except KeyError:
     534                break
     535    return yb
     536   
     537def getBackgroundDerv(pfx,parmDict,bakType,xdata):
     538    dydb = []
     539    if bakType == 'chebyschev':
     540        iBak = 0
     541        while True:
     542            if pfx+'Back:'+str(iBak) in parmDict:
     543                dydb.append((xdata-xdata[0])**iBak)
     544                iBak += 1
     545            else:
     546                break
     547    return dydb
     548
    537549#use old fortran routine
    538 def getFCJVoigt3(pos,intens,sig,gam,shl,xdata):
     550def getFCJVoigt3(pos,sig,gam,shl,xdata):
    539551   
    540552    Df = pyd.pypsvfcj(len(xdata),xdata-pos,pos,sig,gam,shl)
    541553    Df /= np.sum(Df)
    542     return intens*Df
     554    return Df
     555
     556def getdFCJVoigt3(pos,sig,gam,shl,xdata):
     557   
     558    Df,dFdp,dFds,dFdg,dFdsh = pyd.pydpsvfcj(len(xdata),xdata-pos,pos,sig,gam,shl) #might have to make these numpy arrays?
     559    sumDf = np.sum(Df)
     560    return Df/sumDf,dFdp,dFds,dFdg,dFdsh
     561   
    543562
    544563def getPeakProfile(parmDict,xdata,varyList,bakType):
    … …  
    552571    X = parmDict['X']
    553572    Y = parmDict['Y']
    554     shl = max(parmDict['SH/L'],0.0005)
     573    shl = max(parmDict['SH/L'],0.002)
    555574    Ka2 = False
    556575    if 'Lam1' in parmDict.keys():
    … …  
    589608            elif not iBeg-iFin:     #peak above high limit
    590609                return yb+yc
    591             yc[iBeg:iFin] += getFCJVoigt3(pos,intens,sig,gam,shl,xdata[iBeg:iFin])
     610            yc[iBeg:iFin] += intens*getFCJVoigt3(pos,sig,gam,shl,xdata[iBeg:iFin])
    592611            if Ka2:
    593612                pos2 = pos+lamRatio*tand(pos/2.0)       # + 360/pi * Dlam/lam * tan(th)
    … …  
    596615                iFin = min(lenX,iFin+kdelt)
    597616                if iBeg-iFin:
    598                     yc[iBeg:iFin] += getFCJVoigt3(pos2,intens*kRatio,sig,gam,shl,xdata[iBeg:iFin])
     617                    yc[iBeg:iFin] += intens*kRatio*getFCJVoigt3(pos2,sig,gam,shl,xdata[iBeg:iFin])
    599618            iPeak += 1
    600619        except KeyError:        #no more peaks to process
    601620            return yb+yc
    602 
     621           
     622def getPeakProfileDerv(parmDict,xdata,varyList,bakType):
     623# needs to return np.array([dMdx1,dMdx2,...]) in same order as varylist = backVary,insVary,peakVary order
     624    dMdv = np.zeros(shape=(len(varyList),len(xdata)))
     625    if 'Back:0' in varyList:            #background derivs are in front if present
     626        dMdb = getBackgroundDerv('',parmDict,bakType,xdata)
     627        dMdv[0:len(dMdb)] = dMdb
     628       
     629    dx = xdata[1]-xdata[0]
     630    U = parmDict['U']
     631    V = parmDict['V']
     632    W = parmDict['W']
     633    X = parmDict['X']
     634    Y = parmDict['Y']
     635    shl = max(parmDict['SH/L'],0.002)
     636    Ka2 = False
     637    if 'Lam1' in parmDict.keys():
     638        Ka2 = True
     639        lamRatio = 360*(parmDict['Lam2']-parmDict['Lam1'])/(np.pi*parmDict['Lam1'])
     640        kRatio = parmDict['I(L2)/I(L1)']
     641    iPeak = 0
     642    while True:
     643        try:
     644            pos = parmDict['pos'+str(iPeak)]
     645            intens = parmDict['int'+str(iPeak)]
     646            sigName = 'sig'+str(iPeak)
     647            tanth = tand(pos/2.0)
     648            costh = cosd(pos/2.0)
     649            if sigName in varyList:
     650                sig = parmDict[sigName]
     651            else:
     652                sig = U*tanth**2+V*tanth+W
     653                dsdU = tanth**2
     654                dsdV = tanth
     655                dsdW = 1.0
     656            sig = max(sig,0.001)          #avoid neg sigma
     657            gamName = 'gam'+str(iPeak)
     658            if gamName in varyList:
     659                gam = parmDict[gamName]
     660            else:
     661                gam = X/costh+Y*tanth
     662                dgdX = 1.0/costh
     663                dgdY = tanth
     664            gam = max(gam,0.001)             #avoid neg gamma
     665            Wd,fmin,fmax = getWidths(pos,sig,gam,shl)
     666            iBeg = np.searchsorted(xdata,pos-fmin)
     667            lenX = len(xdata)
     668            if not iBeg:
     669                iFin = np.searchsorted(xdata,pos+fmin)
     670            elif iBeg == lenX:
     671                iFin = iBeg
     672            else:
     673                iFin = min(lenX,iBeg+int((fmin+fmax)/dx))
     674            if not iBeg+iFin:       #peak below low limit
     675                iPeak += 1
     676                continue
     677            elif not iBeg-iFin:     #peak above high limit
     678                break
     679            dMdpk = np.zeros(shape=(6,len(xdata)))
     680            dMdipk = getdFCJVoigt3(pos,sig,gam,shl,xdata[iBeg:iFin])
     681            for i in range(5):
     682                dMdpk[i][iBeg:iFin] = 100.*dx*intens*dMdipk[i]
     683            dMdpk[0][iBeg:iFin] = dMdipk[0]
     684            dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4]}
     685            if Ka2:
     686                pos2 = pos+lamRatio*tand(pos/2.0)       # + 360/pi * Dlam/lam * tan(th)
     687                kdelt = int((pos2-pos)/dx)               
     688                iBeg = min(lenX,iBeg+kdelt)
     689                iFin = min(lenX,iFin+kdelt)
     690                if iBeg-iFin:
     691                    dMdipk = getdFCJVoigt3(pos2,sig,gam,shl,xdata[iBeg:iFin])
     692                    for i in range(5):
     693                        dMdpk[i][iBeg:iFin] += 100.*dx*intens*kRatio*dMdipk[i]
     694                    dMdpk[0][iBeg:iFin] += kRatio*dMdipk[0]
     695                    dMdpk[5][iBeg:iFin] += dMdipk[0]
     696                    dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':dMdpk[5]*intens}
     697            for parmName in ['pos','int','sig','gam']:
     698                try:
     699                    idx = varyList.index(parmName+str(iPeak))
     700                    dMdv[idx] = dervDict[parmName]
     701                except ValueError:
     702                    pass
     703            if 'U' in varyList:
     704                dMdv[varyList.index('U')] += dsdU*dervDict['sig']
     705            if 'V' in varyList:
     706                dMdv[varyList.index('V')] += dsdV*dervDict['sig']
     707            if 'W' in varyList:
     708                dMdv[varyList.index('W')] += dsdW*dervDict['sig']
     709            if 'X' in varyList:
     710                dMdv[varyList.index('X')] += dgdX*dervDict['gam']
     711            if 'Y' in varyList:
     712                dMdv[varyList.index('Y')] += dgdY*dervDict['gam']
     713            if 'SH/L' in varyList:
     714                dMdv[varyList.index('SH/L')] += dervDict['shl']         #problem here
     715            if 'I(L2)/I(L1)' in varyList:
     716                dMdv[varyList.index('I(L2)/I(L1)')] += dervDict['L1/L2']
     717            iPeak += 1
     718        except KeyError:        #no more peaks to process
     719            break
     720    return dMdv
     721       
    603722def Dict2Values(parmdict, varylist):
    604723    '''Use before call to leastsq to setup list of values for the parameters
    … …  
    611730    parmdict.update(zip(varylist,values))
    612731   
    613 def DoPeakFit(FitPgm,Peaks,Background,Limits,Inst,data,oneCycle=False):
     732def DoPeakFit(FitPgm,Peaks,Background,Limits,Inst,data,oneCycle=False,controls=None):
    614733   
    615734    def SetBackgroundParms(Background):
    … …  
    651770        insVary = []
    652771        for i,flag in enumerate(insFlags):
    653             if flag:
     772            if flag and insNames[i] in ['U','V','W','X','Y','SH/L','I(L2)/I(L1)']:
    654773                insVary.append(insNames[i])
    655774        instDict = dict(zip(insNames,insVals))
    656775        instDict['X'] = max(instDict['X'],0.01)
    657776        instDict['Y'] = max(instDict['Y'],0.01)
    658         instDict['SH/L'] = max(instDict['SH/L'],0.0001)
     777        instDict['SH/L'] = max(instDict['SH/L'],0.002)
    659778        return dataType,instDict,insVary
    660779       
    … …  
    743862                    ptstr += 12*' '
    744863            print '%s'%(('Peak'+str(i+1)).center(8)),ptstr
     864               
     865    def devPeakProfile(values, xdata, ydata, weights, parmdict, varylist,bakType,dlg):
     866        parmdict.update(zip(varylist,values))
     867        return np.sqrt(weights)*getPeakProfileDerv(parmdict,xdata,varylist,bakType)
    745868           
    746869    def errPeakProfile(values, xdata, ydata, weights, parmdict, varylist,bakType,dlg):       
    747870        parmdict.update(zip(varylist,values))
    748         M = np.sqrt(weights)*(ydata-getPeakProfile(parmdict,xdata,varylist,bakType))
     871        M = np.sqrt(weights)*(getPeakProfile(parmdict,xdata,varylist,bakType)-ydata)
    749872        Rwp = min(100.,np.sqrt(np.sum(M**2)/np.sum(weights*ydata**2))*100.)
    750873        if dlg:
    … …  
    754877        return M
    755878       
    756     Ftol = 0.0001
     879    if controls:
     880        Ftol = controls['min dM/M']
     881        derivType = controls['deriv type']
     882    else:
     883        Ftol = 0.0001
     884        derivType = 'analytic'
    757885    if oneCycle:
    758886        Ftol = 1.0
    … …  
    778906            dlg.SetPosition(wx.Point(screenSize[2]-Size[0]-305,screenSize[1]+5))
    779907            try:
    780                 result = so.leastsq(errPeakProfile,values,full_output=True,epsfcn=1.e-8,ftol=Ftol,
    781                     args=(x[xBeg:xFin],y[xBeg:xFin],w[xBeg:xFin],parmDict,varyList,bakType,dlg))
     908                if derivType == 'analytic':
     909                    result = so.leastsq(errPeakProfile,values,Dfun=devPeakProfile,full_output=True,ftol=Ftol,col_deriv=True,
     910                        args=(x[xBeg:xFin],y[xBeg:xFin],w[xBeg:xFin],parmDict,varyList,bakType,dlg))
     911                    ncyc = int(result[2]['nfev']/2)
     912                else:
     913                    result = so.leastsq(errPeakProfile,values,full_output=True,ftol=Ftol,epsfcn=1.e-8,
     914                        args=(x[xBeg:xFin],y[xBeg:xFin],w[xBeg:xFin],parmDict,varyList,bakType,dlg))
     915                    ncyc = int(result[2]['nfev']/len(varyList))
    782916            finally:
    783917                dlg.Destroy()
    784918            runtime = time.time()-begin   
    785919            chisq = np.sum(result[2]['fvec']**2)
    786             ncyc = int(result[2]['nfev']/len(varyList))
    787920            Values2Dict(parmDict, varyList, result[0])
    788921            Rwp = np.sqrt(chisq/np.sum(w[xBeg:xFin]*y[xBeg:xFin]**2))*100.      #to %
    … …  
    9291062        'Lam1':1.540500,'Lam2':1.544300,'I(L2)/I(L1)':0.5,
    9301063        }
     1064    global parmDict2
     1065    parmDict2 = {
     1066        'pos0':5.7,'int0':10.0,'sig0':0.5,'gam0':1.0,
     1067        'U':1.,'V':-1.,'W':0.1,'X':0.0,'Y':2.,'SH/L':0.004,
     1068        'Back0':5.,'Back1':-0.02,'Back2':.004,
     1069        'Lam1':1.540500,'Lam2':1.544300,'I(L2)/I(L1)':0.5,
     1070        }
    9311071    global varyList
    9321072    varyList = []
    … …  
    9491089    print '100+6*Ka1-2 peaks=1200 peaks',time.time()-time0
    9501090   
     1091def test2(name,delt):
     1092    if NeedTestData: TestData()
     1093    varyList = [name,]
     1094    xdata = np.linspace(5.6,5.8,800)
     1095    hplot = plotter.add('derivatives test for '+name).gca()
     1096    hplot.plot(xdata,getPeakProfileDerv(parmDict2,xdata,varyList,bakType)[0])
     1097    y0 = getPeakProfile(parmDict2,xdata,varyList,bakType)
     1098    parmDict2[name] += delt
     1099    y1 = getPeakProfile(parmDict2,xdata,varyList,bakType)
     1100    hplot.plot(xdata,(y1-y0)/delt,'r+')
     1101   
    9511102   
    9521103
    … …  
    9551106    global plotter
    9561107    plotter = plot.PlotNotebook()
    957     test0()
     1108#    test0()
     1109    test2('I(L2)/I(L1)',.0001)
    9581110    print "OK"
    9591111    plotter.StartEventLoop()
Note: See TracChangeset for help on using the changeset viewer.