Changeset 342 for trunk/GSASIIpwd.py


Ignore:
Timestamp:
Aug 5, 2011 2:35:43 PM (12 years ago)
Author:
vondreele
Message:

major modifications to get 1st "working" version of Refine
Add GSASIImapvars.py
fix cell refinement in indexing window
single cycle for peak refinement

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASIIpwd.py

    r335 r342  
    9292    return plist
    9393
    94 def ValuesOut(parmdict, varylist):
    95     '''Use before call to leastsq to setup list of values for the parameters
    96     in parmdict, as selected by key in varylist'''
    97     return [parmdict[key] for key in varylist]
    98    
    99 def ValuesIn(parmdict, varylist, values):
    100     ''' Use after call to leastsq to update the parameter dictionary with
    101     values corresponding to keys in varylist'''
    102     parmdict.update(zip(varylist,values))
    103    
    10494def Transmission(Geometry,Abs,Diam):
    10595#Calculate sample transmission
    … …  
    482472fcjde = fcjde_gen(name='fcjde',shapes='t,s,dx')
    483473               
    484 def getBackground(parmDict,bakType,xdata):
     474def getBackground(pfx,parmDict,bakType,xdata):
     475    yb = np.zeros_like(xdata)
    485476    if bakType == 'chebyschev':
    486         yb = np.zeros_like(xdata)
    487477        iBak = 0
    488478        while True:
    489             key = 'Back'+str(iBak)
     479            key = pfx+'Back:'+str(iBak)
    490480            try:
    491481                yb += parmDict[key]*(xdata-xdata[0])**iBak
    492482                iBak += 1
    493483            except KeyError:
    494                 return yb
    495                
     484                break
     485    return yb
     486
     487def calcPeakFFT(x,fxns,widths,pos,args):
     488    dx = x[1]-x[0]
     489    z = np.empty([len(fxns),len(x)])
     490    for i,(fxn,width,arg) in enumerate(zip(fxns,widths,args)):
     491        z[i] = fxn.pdf(x,loc=pos,scale=width,*arg)*dx
     492    Z = fft(z)
     493    if len(fxns)%2:
     494        return ifft(Z.prod(axis=0)).real
     495    else:
     496        return fftshift(ifft(Z.prod(axis=0))).real
     497                   
     498def getFCJVoigt(pos,intens,sig,gam,shl,xdata):
     499   
     500    DX = xdata[1]-xdata[0]
     501    widths,fmin,fmax = getWidths(pos,sig,gam,shl)
     502    x = np.linspace(pos-fmin,pos+fmin,1024)
     503    if pos > 90:
     504        fmin,fmax = [fmax,fmin]         
     505    dx = x[1]-x[0]
     506    arg = [pos,shl/10.,dx,]
     507    Df = fcjde.pdf(x,*arg,loc=pos,scale=1.0)
     508    if len(np.nonzero(Df)[0])>5:
     509        fxns = [st.norm,st.cauchy,fcjde]
     510        args = [[],[],arg]
     511    else:
     512        fxns = [st.norm,st.cauchy]
     513        args = [[],[]]
     514    Df = calcPeakFFT(x,fxns,widths,pos,args)
     515    Df /= np.sum(Df)
     516    Df = si.interp1d(x,Df,bounds_error=False,fill_value=0.0)
     517    return intens*Df(xdata)*DX/dx       #*10 to get close to old fxn???
     518                                   
    496519def getPeakProfile(parmDict,xdata,varyList,bakType):
    497520   
    498     def calcPeakFFT(x,fxns,widths,pos,args):
    499         dx = x[1]-x[0]
    500         z = np.empty([len(fxns),len(x)])
    501         for i,(fxn,width,arg) in enumerate(zip(fxns,widths,args)):
    502             z[i] = fxn.pdf(x,loc=pos,scale=width,*arg)*dx
    503         Z = fft(z)
    504         if len(fxns)%2:
    505             return ifft(Z.prod(axis=0)).real
    506         else:
    507             return fftshift(ifft(Z.prod(axis=0))).real
    508                        
    509     def getFCJVoigt(pos,intens,sig,gam,shl,xdata):
    510        
    511         DX = xdata[1]-xdata[0]
    512         widths,fmin,fmax = getWidths(pos,sig,gam,shl)
    513         x = np.linspace(pos-fmin,pos+fmin,1024)
    514         if pos > 90:
    515             fmin,fmax = [fmax,fmin]         
    516         dx = x[1]-x[0]
    517         arg = [pos,shl/10.,dx,]
    518         Df = fcjde.pdf(x,*arg,loc=pos,scale=1.0)
    519         if len(np.nonzero(Df)[0])>5:
    520             fxns = [st.norm,st.cauchy,fcjde]
    521             args = [[],[],arg]
    522         else:
    523             fxns = [st.norm,st.cauchy]
    524             args = [[],[]]
    525         Df = calcPeakFFT(x,fxns,widths,pos,args)
    526         Df /= np.sum(Df)
    527         Df = si.interp1d(x,Df,bounds_error=False,fill_value=0.0)
    528         return intens*Df(xdata)*DX/dx       #*10 to get close to old fxn???
    529                    
    530     yb = getBackground(parmDict,bakType,xdata)
     521    yb = getBackground('',parmDict,bakType,xdata)
    531522    yc = np.zeros_like(yb)
    532523    U = parmDict['U']
    … …  
    591582    return widths,fmin,fmax
    592583               
    593 def DoPeakFit(FitPgm,Peaks,Background,Limits,Inst,data):
     584def Dict2Values(parmdict, varylist):
     585    '''Use before call to leastsq to setup list of values for the parameters
     586    in parmdict, as selected by key in varylist'''
     587    return [parmdict[key] for key in varylist]
     588   
     589def Values2Dict(parmdict, varylist, values):
     590    ''' Use after call to leastsq to update the parameter dictionary with
     591    values corresponding to keys in varylist'''
     592    parmdict.update(zip(varylist,values))
     593   
     594def DoPeakFit(FitPgm,Peaks,Background,Limits,Inst,data,oneCycle=False):
    594595   
    595596    def SetBackgroundParms(Background):
    596597        bakType,bakFlag = Background[:2]
    597598        backVals = Background[3:]
    598         backNames = ['Back'+str(i) for i in range(len(backVals))]
     599        backNames = ['Back:'+str(i) for i in range(len(backVals))]
    599600        if bakFlag: #returns backNames as varyList = backNames
    600601            return bakType,dict(zip(backNames,backVals)),backNames
    … …  
    606607        while True:
    607608            try:
    608                 bakName = 'Back'+str(iBak)
     609                bakName = 'Back:'+str(iBak)
    609610                Background[iBak+3] = parmList[bakName]
    610611                iBak += 1
    … …  
    620621            for i,back in enumerate(Background[3:]):
    621622                ptstr += ptfmt % (back)
    622                 sigstr += ptfmt % (sigDict['Back'+str(i)])
     623                sigstr += ptfmt % (sigDict['Back:'+str(i)])
    623624            print ptstr
    624625            print sigstr
    … …  
    675676        print sigstr
    676677
    677     def setPeaksParms(Peaks):
     678    def SetPeaksParms(Peaks):
    678679        peakNames = []
    679680        peakVary = []
    … …  
    734735        return M
    735736       
     737    Ftol = 0.01
     738    if oneCycle:
     739        Ftol = 1.0
    736740    x,y,w,yc,yb,yd = data               #these are numpy arrays!
    737741    xBeg = np.searchsorted(x,Limits[0])
    … …  
    739743    bakType,bakDict,bakVary = SetBackgroundParms(Background)
    740744    dataType,insDict,insVary = SetInstParms(Inst)
    741     peakDict,peakVary = setPeaksParms(Peaks)
     745    peakDict,peakVary = SetPeaksParms(Peaks)
    742746    parmDict = {}
    743747    parmDict.update(bakDict)
    … …  
    747751    while True:
    748752        begin = time.time()
    749         values =  np.array(ValuesOut(parmDict, varyList))
     753        values =  np.array(Dict2Values(parmDict, varyList))
    750754        if FitPgm == 'LSQ':
    751755            dlg = wx.ProgressDialog('Residual','Peak fit Rwp = ',101.0,
    … …  
    755759            dlg.SetPosition(wx.Point(screenSize[2]-Size[0]-305,screenSize[1]+5))
    756760            try:
    757                 result = so.leastsq(errPeakProfile,values,
    758                     args=(x[xBeg:xFin],y[xBeg:xFin],w[xBeg:xFin],parmDict,varyList,bakType,dlg),full_output=True)
     761                result = so.leastsq(errPeakProfile,values,full_output=True,ftol=Ftol,
     762                    args=(x[xBeg:xFin],y[xBeg:xFin],w[xBeg:xFin],parmDict,varyList,bakType,dlg))
    759763            finally:
    760764                dlg.Destroy()
    … …  
    762766            chisq = np.sum(result[2]['fvec']**2)
    763767            ncyc = int(result[2]['nfev']/len(varyList))
    764             ValuesIn(parmDict, varyList, result[0])
     768            Values2Dict(parmDict, varyList, result[0])
    765769            Rwp = np.sqrt(chisq/np.sum(w[xBeg:xFin]*y[xBeg:xFin]**2))*100.      #to %
    766770            GOF = chisq/(xFin-xBeg-len(varyList))
    … …  
    786790       
    787791    sigDict = dict(zip(varyList,sig))
    788     yb[xBeg:xFin] = getBackground(parmDict,bakType,x[xBeg:xFin])
     792    yb[xBeg:xFin] = getBackground('',parmDict,bakType,x[xBeg:xFin])
    789793    yc[xBeg:xFin] = getPeakProfile(parmDict,x[xBeg:xFin],varyList,bakType)
    790794    yd[xBeg:xFin] = y[xBeg:xFin]-yc[xBeg:xFin]
    … …  
    913917    msg = 'test '
    914918    gplot = plotter.add('FCJ-Voigt, 11BM').gca()
    915     gplot.plot(xdata,getBackground(parmDict0,bakType,xdata))   
     919    gplot.plot(xdata,getBackground('',parmDict0,bakType,xdata))   
    916920    gplot.plot(xdata,getPeakProfile(parmDict0,xdata,varyList,bakType))
    917921    fplot = plotter.add('FCJ-Voigt, Ka1+2').gca()
    918     fplot.plot(xdata,getBackground(parmDict1,bakType,xdata))   
     922    fplot.plot(xdata,getBackground('',parmDict1,bakType,xdata))   
    919923    fplot.plot(xdata,getPeakProfile(parmDict1,xdata,varyList,bakType))
     924   
     925def test1():
     926    if NeedTestData: TestData()
    920927    time0 = time.time()
    921928    for i in range(100):
    922929        y = getPeakProfile(parmDict1,xdata,varyList,bakType)
    923     print '100+6*Ka1-2 peaks=3600 peaks',time.time()-time0
     930    print '100+6*Ka1-2 peaks=1200 peaks',time.time()-time0
     931   
    924932   
    925933
Note: See TracChangeset for help on using the changeset viewer.