Changeset 296 for trunk/GSASIIpwd.py


Ignore:
Timestamp:
Jun 9, 2011 4:02:38 PM (11 years ago)
Author:
vondreele
Message:

major change to peak fitting in GSASIIpwd.py & GSASIIpwdGUI.py
"Auto peakfit" removed from GSASIIgrid.py
ValSig? & Fesd moved from GSASIIimage.py to GSASIIIO.py
Change listing output for image calibration - now shows esds for parameters
GSASIIlattice.py - rebin removed - not to be used ever

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIpwd.py

    r292 r296  
     1#/usr/bin/env python
     2# -*- coding: utf-8 -*-
    13#GSASII powder calculation module
    24########### SVN repository information ###################
     
    1113import wx
    1214import time
     15
    1316import numpy as np
    1417import scipy as sp
    1518import numpy.linalg as nl
     19from numpy.fft import ifft, fft, fftshift
    1620import scipy.interpolate as si
    1721import scipy.stats as st
     22import scipy.optimize as so
     23
    1824import GSASIIpath
    19 import pypowder as pyp              #assumes path has been amended to include correctr bin directory
     25#import pypowder as pyp
    2026import GSASIIplot as G2plt
    2127import GSASIIlattice as G2lat
    2228import GSASIIElem as G2elem
    2329import GSASIIgrid as G2gd
     30import GSASIIIO as G2IO
    2431
    2532# trig functions in degrees
     
    4350npT2q = lambda tth,wave: 2.0*np.pi*npT2stl(tth,wave)
    4451   
    45 np.seterr(divide='ignore')      #this is for the FCJ functions
    46 
    47 #Peak shape definitions
    48 # Finger-Cox_Jephcoat D(2phi,2th) function for S/L = H/L
    49 
    50 class fcjde_gen(st.rv_continuous):
    51    
    52     """
    53     Finger-Cox-Jephcoat D(2phi,2th) function for S/L = H/L
    54     Ref: J. Appl. Cryst. (1994) 27, 892-900.
    55     Parameters
    56     -----------------------------------------
    57     x: array like 2-theta steps (-1 to 1)
    58     t: 2-theta position of peak
    59     s: sum(S/L,H/L); S: sample height, H: detector opening,
    60         L: sample to detector opening distance
    61     dx: 2-theta step size in deg
    62     Result for fcj.pdf
    63     -----------------------------------------
    64     if t < 90:
    65         X = dx*x+t
    66     else:
    67         X = 180.-dx*x-t
    68     if X < t & s = S/L+H/L:
    69         fcj.pdf = [1/sqrt({cos(x)**2/cos(t)**2}-1) - 1/s]/|cos(X)|   
    70     if X >= t:
    71         fcj.pdf = 0   
    72     """
    73     def _pdf(self,x,t,s,dx):
    74         T = np.where(t<=90.,dx*x+t,180.-dx*x-t)
    75         ax = npcosd(T)**2
    76         bx = npcosd(t)**2
    77         bx = np.where(ax>bx,bx,ax)
    78         fx = np.where(ax>bx,(np.sqrt(bx/(ax-bx))-1./s)/np.sqrt(ax),0.0)
    79         fx = np.where(fx > 0,fx,0.0)
    80         return fx
    81 #    def _cdf(self, x):
    82 #    def _ppf(self, q):
    83 #    def _sf(self, x):
    84 #    def _isf(self, q):
    85 #    def _stats(self):
    86 #    def _entropy(self):
    87        
    88 fcjde = fcjde_gen(name='fcjde')
    89                
    90        
    91 # Finger-Cox_Jephcoat D(2phi,2th) function for S/L != H/L
    92 
    93 class fcjd_gen(st.rv_continuous):
    94     """
    95     Finger-Cox-Jephcoat D(2phi,2th) function for S/L != H/L
    96     Ref: J. Appl. Cryst. (1994) 27, 892-900.
    97     Parameters
    98     -----------------------------------------
    99     x: array like 2-theta step numbers (-1 to 1)
    100     t: 2-theta position of peak
    101     h: detector opening height/sample to detector opening distance
    102     s: sample height/sample to detector opening distance
    103     dx: 2-theta step size in deg
    104     Result for fcj2.pdf
    105     -----------------------------------------
    106     if t < 90:
    107         X = dx*x+t
    108     else:
    109         X = 180.-dx*x-t
    110     infl = acos(cos(t)*sqrt((h-s)**2+1))
    111     if X < infl:   
    112         fcj.pdf = [1/sqrt({cos(X)**2/cos(t)**2}-1) - 1/shl]/cos(X)   
    113     for X < t & s = S/L+H/L   
    114     fcj.pdf(x,t,s) = 0   
    115     for 2phi >= 2th
    116     """
    117     def _pdf(self,x,t,h,s,dx):
    118         T = np.where(t<=90.,dx*x+t,180.-dx*x-t)
    119         a = np.abs(npcosd(t))*(np.sqrt((h-s)**2+1.))
    120         infl = np.where((a <= 1.),np.abs(npacosd(a)),T)
    121         ax = npcosd(T)**2
    122         bx = npcosd(t)**2
    123         bx = np.where(ax>bx,bx,ax)
    124         H = np.where(ax>bx,np.sqrt((ax/bx)-1.),0.0)
    125         W1 = h+s-H
    126         W2 = np.where ((h > s),2.*s,2.*h)
    127         fx = 2.*h*np.sqrt((ax/bx)-1.)*np.sqrt(ax)
    128         fx = np.where(fx>0.0,1./fx,0.0)
    129         fx = np.where((T < infl),fx*W1,fx*W2)
    130         return np.where((fx > 0.),fx,0.0)
    131 #    def _cdf(self, x):
    132 #    def _ppf(self, q):
    133 #    def _sf(self, x):
    134 #    def _isf(self, q):
    135 #    def _stats(self):
    136 #    def _entropy(self):
    137        
    138 fcjd = fcjd_gen(name='fcjd')
    139                
    140 # Finger-Cox_Jephcoat D(2phi,2th) function for S/L != H/L using sum & difference
    141 
    142 class fcjdsd_gen(st.rv_continuous):
    143     """
    144     Finger-Cox-Jephcoat D(2phi,2th) function for S/L != H/L using sum & difference
    145     Ref: J. Appl. Cryst. (1994) 27, 892-900.
    146     Parameters
    147     -----------------------------------------
    148     x: array like 2-theta positions
    149     t: 2-theta position of peak
    150     s: sum(S/L,H/L); S: sample height, H: detector opening,
    151         L: sample to detector opening distance
    152     d: difference(S/L,H/L)
    153     dx: 2-theta step size in deg
    154     Result for fcj2.pdf
    155     -----------------------------------------
    156     if t < 90:
    157         X = dx*x+t
    158     else:
    159         X = 180.-dx*x-t
    160    
    161     fcj.pdf(x,t,s,d) = [1/sqrt({cos(X)**2/cos(t)**2}-1) - 1/shl]/cos(2phi)
    162    
    163     for 2phi < 2tth & shl = S/L+H/L
    164    
    165     fcj.pdf(x,tth,shl) = 0
    166    
    167     for 2phi >= 2th
    168     """
    169     def _argcheck(self,t,s,d,dx):
    170         return (t > 0)&(s > 0)&(abs(d) < s)
    171     def _pdf(self,x,t,s,d,dx):
    172         T = np.where(t<=90.,dx*x+t,180.-dx*x-t)
    173         a = np.abs(npcosd(t))*np.sqrt(d**2+1.)
    174         infl = np.where((a < 1.),np.abs(npacosd(a)),T)
    175         ax = npcosd(T)**2
    176         bx = npcosd(t)**2
    177         bx = np.where(ax>bx,bx,ax)
    178         H = np.where(ax>bx,np.sqrt((ax/bx)-1.),0.0)
    179         W1 = s-H
    180         W2 = np.where ((d > 0),s-d,s+d)
    181         fx = np.where(ax>bx,1./((s+d)*np.sqrt((ax/bx)-1.)*np.sqrt(ax)),0.0)
    182         fx = np.where((T < infl),fx*W1,fx*W2)
    183         return np.where((fx > 0.),fx,0.0)
    184 #    def _cdf(self, x):
    185 #    def _ppf(self, q):
    186 #    def _sf(self, x):
    187 #    def _isf(self, q):
    188 #    def _stats(self):
    189 #    def _entropy(self):
    190        
    191 fcjdsd = fcjdsd_gen(name='fcjdsd')
    192                
     52       
    19353def factorize(num):
    19454    ''' Provide prime number factors for integer num
     
    23393    return plist
    23494
     95def ValuesOut(parmdict, varylist):
     96    '''Use before call to leastsq to setup list of values for the parameters
     97    in parmdict, as selected by key in varylist'''
     98    return [parmdict[key] for key in varylist]
     99   
     100def ValuesIn(parmdict, varylist, values):
     101    ''' Use after call to leastsq to update the parameter dictionary with
     102    values corresponding to keys in varylist'''
     103    parmdict.update(zip(varylist,values))
     104   
    235105def Transmission(Geometry,Abs,Diam):
    236106#Calculate sample transmission
     
    547417        Amu += el['FormulaNo']*el['mu']
    548418       
    549 
    550 def ValEsd(value,esd=0,nTZ=False):                  #NOT complete - don't use
    551     # returns value(esd) string; nTZ=True for no trailing zeros
    552     # use esd < 0 for level of precision shown e.g. esd=-0.01 gives 2 places beyond decimal
    553     #get the 2 significant digits in the esd
    554     edig = lambda esd: int(round(10**(math.log10(esd) % 1+1)))
    555     #get the number of digits to represent them
    556     epl = lambda esd: 2+int(1.545-math.log10(10*edig(esd)))
    557    
    558     mdec = lambda esd: -int(math.log10(abs(esd)))
    559     ndec = lambda esd: int(1.545-math.log10(abs(esd)))
    560     if esd > 0:
    561         fmt = '"%.'+str(ndec(esd))+'f(%d)"'
    562         print fmt,ndec(esd),esd*10**(mdec(esd)+1)
    563         return fmt%(value,int(esd*10**(mdec(esd)+1)))
    564     elif esd < 0:
    565          return str(round(value,mdec(esd)))
     419#GSASII peak fitting routine: Finger, Cox & Jephcoat model       
     420
     421np.seterr(divide='ignore')
     422
     423# Voigt function = convolution (norm,cauchy)
     424
     425class voigt_gen(st.rv_continuous):
     426    '''
     427    Voigt function
     428    Parameters
     429    -----------------------------------------
     430    '''
     431    def _pdf(self,x,s,g):
     432        z = np.empty([2,len(x)])
     433        z[0] = st.norm.pdf(x,scale=s)
     434        z[1] = st.cauchy.pdf(x,scale=g)
     435        Z = fft(z)
     436        return fftshift(ifft(Z.prod(axis=0))).real
     437#    def _cdf(self, x):
     438#    def _ppf(self, q):
     439#    def _sf(self, x):
     440#    def _isf(self, q):
     441#    def _stats(self):
     442#    def _entropy(self):
     443
     444voigt = voigt_gen(name='voigt',shapes='ts,g')
     445       
     446# Finger-Cox_Jephcoat D(2phi,2th) function for S/L = H/L
     447
     448class fcjde_gen(st.rv_continuous):
     449    """
     450    Finger-Cox-Jephcoat D(2phi,2th) function for S/L = H/L
     451    Ref: J. Appl. Cryst. (1994) 27, 892-900.
     452    Parameters
     453    -----------------------------------------
     454    x: array -1 to 1
     455    t: 2-theta position of peak
     456    s: sum(S/L,H/L); S: sample height, H: detector opening,
     457        L: sample to detector opening distance
     458    dx: 2-theta step size in deg
     459    Result for fcj.pdf
     460    -----------------------------------------
     461    if x < t & s = S/L+H/L:
     462        fcj.pdf = [1/sqrt({cos(x)**2/cos(t)**2}-1) - 1/s]/cos(x)   
     463    if x >= t:
     464        fcj.pdf = 0   
     465    """
     466    def _pdf(self,x,t,s,dx):
     467#        T = np.where(t<=90.,dx*x+t,180.-dx*x-t)
     468        T = dx*x+t
     469        ax = npcosd(T)**2
     470        bx = npcosd(t)**2
     471        bx = np.where(ax>bx,bx,ax)
     472        fx = np.where(ax>bx,(np.sqrt(bx/(ax-bx))-1./s)/np.sqrt(ax),0.0)
     473        fx = np.where(fx > 0.,fx,0.0)
     474        return fx
     475#    def _cdf(self, x):
     476#    def _ppf(self, q):
     477#    def _sf(self, x):
     478#    def _isf(self, q):
     479#    def _stats(self):
     480#    def _entropy(self):
     481       
     482fcjde = fcjde_gen(name='fcjde',shapes='t,s,dx')
     483               
     484def getBackground(parmDict,bakType,xdata):
     485    if bakType == 'chebyschev':
     486        yb = np.zeros_like(xdata)
     487        iBak = 0
     488        while True:
     489            key = 'Back'+str(iBak)
     490            try:
     491                yb += parmDict[key]*(xdata-xdata[0])**iBak
     492                iBak += 1
     493            except KeyError:
     494                return yb
     495               
     496def getPeakProfile(parmDict,xdata,varyList,bakType):
     497   
     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)
     531    yc = np.zeros_like(yb)
     532    U = parmDict['U']
     533    V = parmDict['V']
     534    W = parmDict['W']
     535    X = parmDict['X']
     536    Y = parmDict['Y']
     537    shl = max(parmDict['SH/L'],0.001)
     538    Ka2 = False
     539    if 'Lam1' in parmDict.keys():
     540        Ka2 = True
     541        lamRatio = parmDict['Lam2']/parmDict['Lam1']
     542        kRatio = parmDict['I(L2)/I(L1)']
     543    iPeak = 0
     544    while True:
     545        try:
     546            pos = parmDict['pos'+str(iPeak)]
     547            intens = parmDict['int'+str(iPeak)]
     548            sigName = 'sig'+str(iPeak)
     549            if sigName in varyList:
     550                sig = parmDict[sigName]
     551            else:
     552                sig = U*tand(pos/2.0)**2+V*tand(pos/2.0)+W
     553            sig = max(sig,0.001)          #avoid neg sigma
     554            gamName = 'gam'+str(iPeak)
     555            if gamName in varyList:
     556                gam = parmDict[gamName]
     557            else:
     558                gam = X/cosd(pos/2.0)+Y*tand(pos/2.0)
     559            gam = max(gam,0.01)             #avoid neg gamma
     560            Wd,fmin,fmax = getWidths(pos,sig,gam,shl)
     561            if pos > 90:
     562                fmin,fmax = [fmax,fmin]         
     563            iBeg = np.searchsorted(xdata,pos-fmin)
     564            iFin = np.searchsorted(xdata,pos+fmax)
     565            if not iBeg+iFin:       #peak below low limit
     566                iPeak += 1
     567                continue
     568            elif not iBeg-iFin:     #peak above high limit
     569                return yb+yc
     570            yc[iBeg:iFin] += getFCJVoigt(pos,intens,sig,gam,shl,xdata[iBeg:iFin])
     571            if Ka2:
     572                pos2 = 2.0*asind(lamRatio*sind(pos/2.0))
     573                Wd,fmin,fmax = getWidths(pos2,sig,gam,shl)
     574                if pos > 90:
     575                    fmin,fmax = [fmax,fmin]         
     576                iBeg = np.searchsorted(xdata,pos2-fmin)
     577                iFin = np.searchsorted(xdata,pos2+fmax)
     578                yc[iBeg:iFin] += getFCJVoigt(pos2,intens*kRatio,sig,gam,shl,xdata[iBeg:iFin])
     579            iPeak += 1
     580        except KeyError:        #no more peaks to process
     581            return yb+yc
     582       
     583def getWidths(pos,sig,gam,shl):
     584    if shl:
     585        widths = [np.sqrt(sig)/100.,gam/200.,.1]
    566586    else:
    567         text = "%F"%(value)
    568         if nTZ:
    569             return text.rstrip('0')
    570         else:
    571             return text
    572 
    573        
    574 #GSASII peak fitting routine: Thompson, Cox & Hastings; Finger, Cox & Jephcoat model       
    575 
    576 def DoPeakFit(peaks,background,limits,inst,data):
    577    
    578     def backgroundPrint(background,sigback):
    579         if background[1]:
    580             print 'Background coefficients for',background[0],'function'
     587        widths = [np.sqrt(sig)/100.,gam/200.]
     588    fwhm = 2.355*widths[0]+2.*widths[1]
     589    fmin = 10.*(fwhm+shl*abs(npcosd(pos)))
     590    fmax = 15.0*fwhm
     591    return widths,fmin,fmax
     592               
     593def DoPeakFit(FitPgm,Peaks,Background,Limits,Inst,data):
     594   
     595    def SetBackgroundParms(Background):
     596        bakType,bakFlag = Background[:2]
     597        backVals = Background[3:]
     598        backNames = ['Back'+str(i) for i in range(len(backVals))]
     599        if bakFlag: #returns backNames as varyList = backNames
     600            return bakType,dict(zip(backNames,backVals)),backNames
     601        else:       #no background varied; varyList = []
     602            return bakType,dict(zip(backNames,backVals)),[]
     603       
     604    def GetBackgroundParms(parmList,Background):
     605        iBak = 0
     606        while True:
     607            try:
     608                bakName = 'Back'+str(iBak)
     609                Background[iBak+3] = parmList[bakName]
     610                iBak += 1
     611            except KeyError:
     612                break
     613               
     614    def BackgroundPrint(Background,sigDict):
     615        if Background[1]:
     616            print 'Background coefficients for',Background[0],'function'
    581617            ptfmt = "%12.5f"
    582618            ptstr =  'values:'
    583619            sigstr = 'esds  :'
    584             for i,back in enumerate(background[3:]):
     620            for i,back in enumerate(Background[3:]):
    585621                ptstr += ptfmt % (back)
    586                 sigstr += ptfmt % (sigback[i+3])
     622                sigstr += ptfmt % (sigDict['Back'+str(i)])
    587623            print ptstr
    588624            print sigstr
    589625        else:
    590626            print 'Background not refined'
    591    
    592     def instPrint(instVal,siginst,insLabels):
     627           
     628    def SetInstParms(Inst):
     629        insVals,insFlags,insNames = Inst[1:4]
     630        dataType = insVals[0]
     631        insVary = []
     632        for i,flag in enumerate(insFlags):
     633            if flag:
     634                insVary.append(insNames[i])
     635        instDict = dict(zip(insNames,insVals))
     636        instDict['X'] = max(instDict['X'],0.1)
     637        instDict['Y'] = max(instDict['Y'],0.1)
     638        instDict['SH/L'] = max(instDict['SH/L'],0.001)
     639        return dataType,instDict,insVary
     640       
     641    def GetInstParms(parmDict,Inst,varyList,Peaks):
     642        instNames = Inst[3]
     643        for i,name in enumerate(instNames):
     644            Inst[1][i] = parmDict[name]
     645        iPeak = 0
     646        while True:
     647            try:
     648                sigName = 'sig'+str(iPeak)
     649                pos = parmDict['pos'+str(iPeak)]
     650                if sigName not in varyList:
     651                    parmDict[sigName] = parmDict['U']*tand(pos/2.0)**2+parmDict['V']*tand(pos/2.0)+parmDict['W']
     652                gamName = 'gam'+str(iPeak)
     653                if gamName not in varyList:
     654                    parmDict[gamName] = parmDict['X']/cosd(pos/2.0)+parmDict['Y']*tand(pos/2.0)
     655                iPeak += 1
     656            except KeyError:
     657                break
     658       
     659    def InstPrint(Inst,sigDict):
    593660        print 'Instrument Parameters:'
    594661        ptfmt = "%12.6f"
     
    596663        ptstr =  'values:'
    597664        sigstr = 'esds  :'
    598         for i,value in enumerate(instVal):
    599             ptlbls += "%s" % (insLabels[i].center(12))
    600             ptstr += ptfmt % (value)
    601             if siginst[i]:
    602                 sigstr += ptfmt % (siginst[i])
     665        instNames = Inst[3][1:]
     666        for i,name in enumerate(instNames):
     667            ptlbls += "%s" % (name.center(12))
     668            ptstr += ptfmt % (Inst[1][i+1])
     669            if name in sigDict:
     670                sigstr += ptfmt % (sigDict[name])
    603671            else:
    604672                sigstr += 12*' '
     
    606674        print ptstr
    607675        print sigstr
    608    
    609     def peaksPrint(peaks,sigpeaks):
    610         print 'Peak list='
    611 
    612     begin = time.time()
    613     Np = len(peaks[0])
    614     DataType = inst[1][0]
    615     instVal = inst[1][1:]
    616     Insref = inst[2][1:]
    617     insLabels = inst[3][1:]
    618     Ka2 = False
    619     Ioff = 3
    620     if len(instVal) == 12:
    621         lamratio = instVal[1]/instVal[0]
    622         Ka2 = True
    623         Ioff = 5
    624     insref = Insref[len(Insref)-7:-1]               #just U,V,W,X,Y,SH/L
    625     for peak in peaks:
    626         dip = []
    627         dip.append(tand(peak[0]/2.0)**2)
    628         dip.append(tand(peak[0]/2.0))
    629         dip.append(1.0/cosd(peak[0]/2.0))
    630         dip.append(tand(peak[0]/2.0))
    631         peak.append(dip)
    632     B = background[2]
    633     bcof = background[3:3+B]
    634     Bv = 0
    635     if background[1]:
    636         Bv = B
     676
     677    def setPeaksParms(Peaks):
     678        peakNames = []
     679        peakVary = []
     680        peakVals = []
     681        names = ['pos','int','sig','gam']
     682        for i,peak in enumerate(Peaks):
     683            for j in range(4):
     684                peakVals.append(peak[2*j])
     685                parName = names[j]+str(i)
     686                peakNames.append(parName)
     687                if peak[2*j+1]:
     688                    peakVary.append(parName)
     689        return dict(zip(peakNames,peakVals)),peakVary
     690               
     691    def GetPeaksParms(parmDict,Peaks,varyList):
     692        names = ['pos','int','sig','gam']
     693        for i,peak in enumerate(Peaks):
     694            for j in range(4):
     695                pos = parmDict['pos'+str(i)]
     696                parName = names[j]+str(i)
     697                if parName in varyList:
     698                    peak[2*j] = parmDict[parName]
     699                elif 'sig' in parName:
     700                    peak[2*j] = parmDict['U']*tand(pos/2.0)**2+parmDict['V']*tand(pos/2.0)+parmDict['W']
     701                elif 'gam' in parName:
     702                    peak[2*j] = parmDict['X']/cosd(pos/2.0)+parmDict['Y']*tand(pos/2.0)
     703                       
     704    def PeaksPrint(parmDict,sigDict,varyList):
     705        print 'Peak coefficients:'
     706        names = ['pos','int','sig','gam']
     707        head = 15*' '
     708        for name in names:
     709            head += name.center(12)+'esd'.center(12)
     710        print head
     711        ptfmt = {'pos':"%12.5f",'int':"%12.1f",'sig':"%12.3f",'gam':"%12.3f"}
     712        for i,peak in enumerate(Peaks):
     713            ptstr =  ':'
     714            for j in range(4):
     715                name = names[j]
     716                parName = name+str(i)
     717                ptstr += ptfmt[name] % (parmDict[parName])
     718                if parName in varyList:
     719#                    ptstr += G2IO.ValEsd(parmDict[parName],sigDict[parName])
     720                    ptstr += ptfmt[name] % (sigDict[parName])
     721                else:
     722#                    ptstr += G2IO.ValEsd(parmDict[parName],0.0)
     723                    ptstr += 12*' '
     724            print '%s'%(('Peak'+str(i+1)).center(8)),ptstr
     725           
     726    def errPeakProfile(values, xdata, ydata, weights, parmdict, varylist,bakType):       
     727        parmdict.update(zip(varylist,values))
     728        return np.sqrt(weights)*(ydata-getPeakProfile(parmdict,xdata,varylist,bakType))
     729       
    637730    x,y,w,yc,yb,yd = data               #these are numpy arrays!
    638     V = []
    639     A = []
    640     swobs = 0.0
    641     smin = 0.0
    642     first = True
    643     GoOn = True
    644     Go = True
    645     dlg = wx.ProgressDialog("Elapsed time","Fitting peaks to pattern",len(x), \
    646         style = wx.PD_ELAPSED_TIME|wx.PD_AUTO_HIDE|wx.PD_REMAINING_TIME|wx.PD_CAN_ABORT)
    647     screenSize = wx.DisplaySize()
    648     Size = dlg.GetSize()
    649     dlg.SetPosition(wx.Point(screenSize[0]-Size[0]-300,0))
    650     try:
    651         i = 0
    652         for xi in x :
    653             Go = dlg.Update(i)[0]
    654             if GoOn:
    655                 GoOn = Go
    656             if limits[0] <= xi <= limits[1]:
    657                 yb[i] = 0.0
    658                 dp = []
    659                 for j in range(B):
    660                     t = (xi-limits[0])**j
    661                     yb[i] += t*bcof[j]
    662                     if background[1]:
    663                         dp.append(t)
    664                 yc[i] = yb[i]
    665                 Iv = 0
    666                 for j in range(6):
    667                     if insref[j]:
    668                         dp.append(0.0)
    669                         Iv += 1
    670                 for peak in peaks:
    671                     dip = peak[-1]
    672                     f = pyp.pypsvfcj(peak[2],xi-peak[0],peak[0],peak[4],peak[6],instVal[-2],0.0)
    673                     yc[i] += f[0]*peak[2]
    674                     if f[0] > 0.0:
    675                         j = 0
    676                         if insref[0]:              #U
    677                             dp[Bv+j] += f[3]*dip[0]
    678                             j += 1
    679                         if insref[1]:              #V
    680                             dp[Bv+j] += f[3]*dip[1]
    681                             j += 1
    682                         if insref[2]:              #W
    683                             dp[Bv+j] += f[3]
    684                             j += 1
    685                         if insref[3]:              #X
    686                             dp[Bv+j] += f[4]*dip[2]
    687                             j += 1
    688                         if insref[4]:              #Y
    689                             dp[Bv+j] += f[4]*dip[3]
    690                             j += 1
    691                         if insref[5]:              #SH/L
    692                             dp[Bv+j] += f[5]
    693                     if Ka2:
    694                        pos2 = 2.0*asind(lamratio*sind(peak[0]/2.0))
    695                        f2 = pyp.pypsvfcj(peak[2],xi-pos2,peak[0],peak[4],peak[6],instVal[-2],0.0)
    696                        yc[i] += f2[0]*peak[2]*instVal[3]
    697                        if f[0] > 0.0:
    698                            j = 0
    699                            if insref[0]:              #U
    700                                dp[Bv+j] += f2[3]*dip[0]*instVal[3]
    701                                j += 1
    702                            if insref[1]:              #V
    703                                dp[Bv+j] += f2[3]*dip[1]*instVal[3]
    704                                j += 1
    705                            if insref[2]:              #W
    706                                dp[Bv+j] += f2[3]*instVal[3]
    707                                j += 1
    708                            if insref[3]:              #X
    709                                dp[Bv+j] += f2[4]*dip[2]*instVal[3]
    710                                j += 1
    711                            if insref[4]:              #Y
    712                                dp[Bv+j] += f2[4]*dip[3]*instVal[3]
    713                                j += 1
    714                            if insref[5]:              #SH/L
    715                                dp[Bv+j] += f2[5]*instVal[3]                       
    716                     for j in range(0,Np,2):
    717                         if peak[j+1]: dp.append(f[j/2+1])
    718                 yd[i] = y[i]-yc[i]
    719                 swobs += w[i]*y[i]**2
    720                 t2 = w[i]*yd[i]
    721                 smin += t2*yd[i]
    722                 if first:
    723                     first = False
    724                     M = len(dp)
    725                     A = np.zeros(shape=(M,M))
    726                     V = np.zeros(shape=(M))
    727                 A,V = pyp.buildmv(t2,w[i],M,dp,A,V)
    728             i += 1
    729     finally:
    730         dlg.Destroy()
    731     Rwp = smin/swobs
    732     Rwp = math.sqrt(Rwp)*100.0
    733     norm = np.diag(A)
    734     for i,elm in enumerate(norm):
    735         if elm <= 0.0:
    736             print norm
    737             return False,0,0,0,False
    738     for i in xrange(len(V)):
    739         norm[i] = 1.0/math.sqrt(norm[i])
    740         V[i] *= norm[i]
    741         a = A[i]
    742         for j in xrange(len(V)):
    743             a[j] *= norm[i]
    744     A = np.transpose(A)
    745     for i in xrange(len(V)):
    746         a = A[i]
    747         for j in xrange(len(V)):
    748             a[j] *= norm[i]
    749     b = nl.solve(A,V)
    750     A = nl.inv(A)
    751     sig = np.diag(A)
    752     for i in xrange(len(V)):
    753         b[i] *= norm[i]
    754         sig[i] *= norm[i]*norm[i]
    755         sig[i] = math.sqrt(abs(sig[i]))
    756     sigback = [0,0,0]
    757     for j in range(Bv):
    758         background[j+3] += b[j]
    759         sigback.append(sig[j])
    760     backgroundPrint(background,sigback)
    761     k = 0
    762     delt = []
    763     if Ka2:
    764         siginst = [0,0,0,0,0]
    765     else:
    766         siginst = [0,0,0]
    767     for j in range(6):
    768         if insref[j]:
    769             instVal[j+Ioff] += b[Bv+k]
    770             siginst.append(sig[Bv+k])
    771             delt.append(b[Bv+k])
    772             k += 1
    773         else:
    774             delt.append(0.0)
    775             siginst.append(0.0)
    776     delt.append(0.0)                    #dummies for azm
    777     siginst.append(0.0)
    778     instPrint(instVal,siginst,insLabels)
    779     inst[1] = [DataType,]
    780     for val in instVal:
    781         inst[1].append(val)
    782     B = Bv+Iv
    783     for peak in peaks:
    784         del peak[-1]                        # remove dip from end
    785         delsig = delt[0]*tand(peak[0]/2.0)**2+delt[1]*tand(peak[0]/2.0)+delt[2]
    786         delgam = delt[3]/cosd(peak[0]/2.0)+delt[4]*tand(peak[0]/2.0)
    787         for j in range(0,len(peak[:-1]),2):
    788             if peak[j+1]:
    789                 peak[j] += b[B]
    790                 B += 1
    791         peak[4] += delsig
    792         if peak[4] < 0.0:
    793             print 'ERROR - negative sigma'
    794             return False,0,0,0,False           
    795         peak[6] += delgam
    796         if peak[6] < 0.0:
    797             print 'ERROR - negative gamma'
    798             return False,0,0,0,False
    799     runtime = time.time()-begin   
    800     data = [x,y,w,yc,yb,yd]
    801     return True,smin,Rwp,runtime,GoOn
    802 
     731    xBeg = np.searchsorted(x,Limits[0])
     732    xFin = np.searchsorted(x,Limits[1])
     733    bakType,bakDict,bakVary = SetBackgroundParms(Background)
     734    dataType,insDict,insVary = SetInstParms(Inst)
     735    peakDict,peakVary = setPeaksParms(Peaks)
     736    parmDict = {}
     737    parmDict.update(bakDict)
     738    parmDict.update(insDict)
     739    parmDict.update(peakDict)
     740    varyList = bakVary+insVary+peakVary
     741    while True:
     742        begin = time.time()
     743        values =  np.array(ValuesOut(parmDict, varyList))
     744        if FitPgm == 'LSQ':
     745            result = so.leastsq(errPeakProfile,values,
     746                args=(x[xBeg:xFin],y[xBeg:xFin],w[xBeg:xFin],parmDict,varyList,bakType),full_output=True)
     747            runtime = time.time()-begin   
     748            print 'Number of function calls:',result[2]['nfev'],' Number of observations: ',xFin-xBeg,' Number of parameters: ',len(varyList)
     749            print "%s%8.3f%s " % ('fitpeak time =',runtime,'s')
     750            ValuesIn(parmDict, varyList, result[0])
     751            chisq = np.sum(errPeakProfile(result[0],x[xBeg:xFin],y[xBeg:xFin],w[xBeg:xFin],parmDict,varyList,bakType)**2)
     752            Rwp = np.sqrt(chisq/np.sum(w[xBeg:xFin]*y[xBeg:xFin]**2))*100.      #to %
     753            GOF = chisq/(xFin-xBeg-len(varyList))
     754            print "%s%7.2f%s%12.6g%s%6.2f" % ('Rwp = ',Rwp,'%, chi**2 = ',chisq,' reduced chi**2 = ',GOF)
     755            try:
     756                sig = np.sqrt(np.diag(result[1])*chisq)
     757                break                   #refinement succeeded - finish up!
     758            except ValueError:          #result[1] is None on singular matrix
     759                print 'Refinement failed - singular matrix'
     760                Ipvt = result[2]['ipvt']
     761                for i,ipvt in enumerate(Ipvt):
     762                    if not np.sum(result[2]['fjac'],axis=1)[i]:
     763                        print 'Removing parameter: ',varyList[ipvt-1]
     764                        del(varyList[ipvt-1])
     765                        break
     766        elif FitPgm == 'BFGS':
     767            print 'Other program here'
     768            return
     769       
     770    sigDict = dict(zip(varyList,sig))
     771    yb[xBeg:xFin] = getBackground(parmDict,bakType,x[xBeg:xFin])
     772    yc[xBeg:xFin] = getPeakProfile(parmDict,x[xBeg:xFin],varyList,bakType)
     773    yd[xBeg:xFin] = y[xBeg:xFin]-yc[xBeg:xFin]
     774    GetBackgroundParms(parmDict,Background)
     775    BackgroundPrint(Background,sigDict)
     776    GetInstParms(parmDict,Inst,varyList,Peaks)
     777    InstPrint(Inst,sigDict)
     778    GetPeaksParms(parmDict,Peaks,varyList)   
     779    PeaksPrint(parmDict,sigDict,varyList)
     780   
    803781def CalcPDF(data,inst,xydata):
    804782    auxPlot = []
     
    881859    return auxPlot
    882860       
     861#testing data
     862import plot
     863NeedTestData = True
     864def TestData():
     865#    global NeedTestData
     866    NeedTestData = False
     867    global bakType
     868    bakType = 'chebyschev'
     869    global xdata
     870    xdata = np.linspace(4.0,40.0,36000)
     871    global parmDict0
     872    parmDict0 = {
     873        'pos0':5.6964,'int0':8835.8,'sig0':1.0,'gam0':1.0,
     874        'pos1':11.4074,'int1':3922.3,'sig1':1.0,'gam1':1.0,
     875        'pos2':20.6426,'int2':1573.7,'sig2':1.0,'gam2':1.0,
     876        'pos3':26.9568,'int3':925.1,'sig3':1.0,'gam3':1.0,
     877        'U':1.163,'V':-0.605,'W':0.093,'X':0.0,'Y':2.183,'SH/L':0.002,
     878        'Back0':5.384,'Back1':-0.015,'Back2':.004,
     879        }
     880    global parmDict1
     881    parmDict1 = {
     882        'pos0':13.4924,'int0':48697.6,'sig0':1.0,'gam0':1.0,
     883        'pos1':23.4360,'int1':43685.5,'sig1':1.0,'gam1':1.0,
     884        'pos2':27.1152,'int2':123712.6,'sig2':1.0,'gam2':1.0,
     885        'pos3':33.7196,'int3':65349.4,'sig3':1.0,'gam3':1.0,
     886        'pos4':36.1119,'int4':115829.8,'sig4':1.0,'gam4':1.0,
     887        'pos5':39.0122,'int5':6916.9,'sig5':1.0,'gam5':1.0,
     888        'U':22.75,'V':-17.596,'W':10.594,'X':1.577,'Y':5.778,'SH/L':0.002,
     889        'Back0':36.897,'Back1':-0.508,'Back2':.006,
     890        'Lam1':1.540500,'Lam2':1.544300,'I(L2)/I(L1)':0.5,
     891        }
     892    global varyList
     893    varyList = []
     894
     895def test0():
     896    if NeedTestData: TestData()
     897    msg = 'test '
     898    gplot = plotter.add('FCJ-Voigt, 11BM').gca()
     899    gplot.plot(xdata,getBackground(parmDict0,bakType,xdata))   
     900    gplot.plot(xdata,getPeakProfile(parmDict0,xdata,varyList,bakType))
     901    fplot = plotter.add('FCJ-Voigt, Ka1+2').gca()
     902    fplot.plot(xdata,getBackground(parmDict1,bakType,xdata))   
     903    fplot.plot(xdata,getPeakProfile(parmDict1,xdata,varyList,bakType))
     904    time0 = time.time()
     905    for i in range(100):
     906        y = getPeakProfile(parmDict1,xdata,varyList,bakType)
     907    print '100+6*Ka1-2 peaks=3600 peaks',time.time()-time0
     908   
     909
     910if __name__ == '__main__':
     911    global plotter
     912    plotter = plot.PlotNotebook()
     913    test0()
     914    print "OK"
     915    plotter.StartEventLoop()
Note: See TracChangeset for help on using the changeset viewer.