Changeset 296


Ignore:
Timestamp:
Jun 9, 2011 4:02:38 PM (12 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

Location:
trunk
Files:
7 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASIIIO.py

    r293 r296  
    804804                Id = item
    805805            item, cookie = self.PatternTree.GetNextChild(self.root, cookie)
    806         parms = ['PXC',data['wavelength'],0.0,0.99,1.0,-1.0,0.3,0.0,1.0,0.0,azm]    #set polarization for synchrotron radiation!
     806        parms = ['PXC',data['wavelength'],0.0,0.99,1.0,-0.10,0.4,0.30,1.0,0.0001,azm]    #set polarization for synchrotron radiation!
    807807        Y = self.Integrate[0][i]
    808808        W = 1./Y                    #probably not true
    … …  
    12321232   
    12331233    return Phase
     1234
     1235def ValEsd(value,esd=0,nTZ=False):                  #NOT complete - don't use
     1236    # returns value(esd) string; nTZ=True for no trailing zeros
     1237    # use esd < 0 for level of precision shown e.g. esd=-0.01 gives 2 places beyond decimal
     1238    #get the 2 significant digits in the esd
     1239    edig = lambda esd: int(round(10**(math.log10(esd) % 1+1)))
     1240    #get the number of digits to represent them
     1241    epl = lambda esd: 2+int(1.545-math.log10(10*edig(esd)))
     1242   
     1243    mdec = lambda esd: -int(math.log10(abs(esd)))
     1244    ndec = lambda esd: int(1.545-math.log10(abs(esd)))
     1245    if esd > 0:
     1246        fmt = '"%.'+str(ndec(esd))+'f(%d)"'
     1247        print fmt,ndec(esd),esd*10**(mdec(esd)+1)
     1248        return fmt%(value,int(esd*10**(mdec(esd)+1)))
     1249    elif esd < 0:
     1250         return str(round(value,mdec(esd)))
     1251    else:
     1252        text = "%F"%(value)
     1253        if nTZ:
     1254            return text.rstrip('0')
     1255        else:
     1256            return text
     1257
     1258def Fesd(value,esd=0,nTZ=False):
     1259#pythonized version of fortran routine in GSAS cifsubs directory - doesn't work correctly
     1260    nint = lambda x: int(round(x))
     1261    iExp = 0
     1262    if value == 0. and esd == 0.:
     1263        iDec = 1
     1264        iFld = 5
     1265    elif value == 0.:
     1266        iDec = max(0.,1.545-math.log10(abs(esd)))
     1267        iFld = 4+iDec
     1268    elif esd == 0.:
     1269        iDec = 5
     1270        iFld = max(1.,math.log10(abs(value)))+3+iDec
     1271    else:
     1272        iFld = math.log10(max(abs(esd),abs(value)))
     1273        if iFld < -4:
     1274            iExp = 1-iFld
     1275            iFld -= iExp
     1276        elif iFld > 8:
     1277            iExp = -iFld
     1278            iFld += iExp
     1279        if iExp:
     1280            value *= 10.**iExp
     1281            esd *= 10.**iExp
     1282        iDec = min(7,int(max(0.,1.545-math.log10(max(0.000001*abs(value),abs(esd))))))
     1283        iFld = max(1,iFld)+3+iDec
     1284    if esd <= 0.:
     1285        iSigw = 0
     1286    else:
     1287        iSig = nint(esd*(10.**iDec))
     1288        iSigw = 1
     1289        if iSig > 0:
     1290            iSigw = int(1.+math.log10(1.*iSig))
     1291    if iSigw > 2:
     1292        xmult = 10.**(iSigw-2)
     1293        value = xmult*nint(value/xmult)
     1294        iSig = xmult*nint(iSig/xmult)           
     1295        iSigw = int(1.+math.log10(1.*iSig))
     1296    if iSigw == 0:
     1297        fmt = '%.'+str(iDec)+'f'
     1298        string = fmt%(value)
     1299    elif iDec > 0:
     1300        fmt = '%.'+str(iDec)+'f(%d)'
     1301        string = fmt%(value,iSig)
     1302    else:
     1303        fmt = '%'+str(iFld)+'d(%d)'
     1304        string = fmt%(nint(value),iSig)
     1305    if iExp:
     1306        iDig = 1+math.log10(abs(1.*iExp))
     1307        if iExp > 0:
     1308            iDig += 1
     1309            string += str(-iExp)
     1310    return string
     1311   
     1312

  • trunk/GSASIIgrid.py

    r294 r296  
    4646] = [wx.NewId() for _init_coll_IndPeaks_Items in range(1)]
    4747
    48 [ wxID_UNDO,wxID_PEAKFIT,wxID_AUTOPEAKFIT,
    49 ] = [wx.NewId() for _init_coll_PEAK_Items in range(3)]
     48[ wxID_UNDO,wxID_LSQPEAKFIT,wxID_BFGSPEAKFIT,wxID_RESETSIGGAM,
     49] = [wx.NewId() for _init_coll_PEAK_Items in range(4)]
    5050
    5151[  wxID_INDEXPEAKS, wxID_REFINECELL, wxID_COPYCELL, wxID_MAKENEWPHASE,
    … …  
    187187        self.UnDo = parent.Append(help='Undo last least squares refinement',
    188188            id=wxID_UNDO, kind=wx.ITEM_NORMAL,text='UnDo')
    189         self.PeakFit = parent.Append(id=wxID_PEAKFIT, kind=wx.ITEM_NORMAL,text='PeakFit',
    190             help='Do single cycle of peak fitting least-squares refinement' )
    191         self.AutoPeakFit = parent.Append(id=wxID_AUTOPEAKFIT, kind=wx.ITEM_NORMAL,
    192             text='AutoPeakFit',help='Do peak fitting least-squares to convergence' )
     189        self.PeakFit = parent.Append(id=wxID_LSQPEAKFIT, kind=wx.ITEM_NORMAL,text='LSQ PeakFit',
     190            help='Peak fitting via least-squares' )
     191        self.PeakFit = parent.Append(id=wxID_BFGSPEAKFIT, kind=wx.ITEM_NORMAL,text='BFGS PeakFit',
     192            help='Peak fitting via BFGS algorithm' )
     193        self.ResetSigGam = parent.Append(id=wxID_RESETSIGGAM, kind=wx.ITEM_NORMAL,
     194            text='Reset sig and gam',help='Reset sigma and gamma to global fit' )
    193195           
    194196    def _init_coll_Index_Items(self,parent):
    … …  
    268270        self.UnDo.Enable(False)
    269271        self.PeakFit.Enable(False)
    270         self.AutoPeakFit.Enable(False)
    271272        self.IndexPeaks.Enable(False)
    272273        self.CopyCell.Enable(False)

  • trunk/GSASIIimgGUI.py

    r294 r296  
    550550    dataSizer.Add(wx.StaticText(parent=self.dataDisplay,label=' Tilt rotation'),0,
    551551        wx.ALIGN_CENTER_VERTICAL)
    552     rotSel = wx.TextCtrl(parent=self.dataDisplay,value=("%9.3f"%(data['rotation']-90.)),style=wx.TE_READONLY)
     552    rotSel = wx.TextCtrl(parent=self.dataDisplay,value=("%9.3f"%(data['rotation']-90.)),style=wx.TE_READONLY) #kluge to get rotation from vertical - see GSASIIimage
    553553    rotSel.SetBackgroundColour(VERY_LIGHT_GREY)
    554554    dataSizer.Add(rotSel,0,wx.ALIGN_CENTER_VERTICAL)

  • trunk/GSASIIlattice.py

    r264 r296  
    208208    return muT/Volume
    209209   
    210 def qRebin(x,y,nQbins):
    211 #Rebin data from x to qx
    212 #   input:
    213 #   X - array of q-values in increasing order
    214 #   y - array of intensities for the x array
    215 #   nQbins - number of q-bins desired
    216 #   returns:
    217 #   qx - array of q-values
    218 #   qy - array of rebinned intensities
    219     dQ = (x[-1]-x[0])/nQbins
    220     dQ2 = dQ/2.
    221     Qmin = x[0]
    222     yf = y[0]
    223     frac = 0.0
    224     jch = 0
    225     kch = 0
    226     qy = []
    227     for i in range(nQbins):
    228         Q = Qmin+dQ*i
    229         rom = Q
    230         rop = rom+dQ
    231         sdq = (1.-frac)*dQ
    232         qy.append(1-frac)*yf*dQ
    233        
    234 #      DELQ = (XT(NCHANS)-XT(1))/FLOAT(MCHANS)
    235 #      DELQ2 = DELQ/2.0
    236 #      QMIN = XT(1)
    237 #      JCH = 1
    238 #      KCH = 0
    239 #      YF = YT(1)
    240 #      FRAC = 0.0
    241 #      DQ = XT(2)-QMIN
    242 #      DO ICH = 1,MCHANS
    243 #        Q = QMIN+DELQ*FLOAT(ICH-1) !Front end interpolation
    244 #        ROM = Q
    245 #        ROP = ROM+DELQ
    246 #        SDQ = (1.0-FRAC)*DQ
    247 #        RY(ICH) = (1.0-FRAC)*YF*DQ !Sum intervening channels
    248 #        DO WHILE ( ROP.GT.XT(KCH+1) .AND. KCH.NE.NCHANS-1 )
    249 #          KCH = MIN(KCH+1,NCHANS-1)
    250 #        END DO
    251 #        DO LCH=JCH+1,KCH
    252 #          DQ = XT(LCH+1)-XT(LCH)
    253 #          SDQ = SDQ+DQ
    254 #          RY(ICH) = RY(ICH)+DQ*YT(LCH)
    255 #        END DO !Tail end interpolation
    256 #        XF = XT(MAX(1,KCH))
    257 #        YF = YT(MAX(1,KCH))
    258 #        DQ = XT(KCH+1)-XF
    259 #        DY = YT(KCH+1)-YF
    260 #        FRAC = (ROP-XF)/DQ
    261 #        SDQ = SDQ+DQ*FRAC
    262 #        RY(ICH) = RY(ICH)+YF*FRAC*DQ !Reset lower limit
    263 #        JCH = KCH
    264 #        RY(ICH) = RY(ICH)/SDQ
    265 #      END DO
    266 #      RETURN
    267 #      END
    268 
    269    
    270 
    271210#Permutations and Combinations
    272211# Four routines: combinations,uniqueCombinations, selections & permutations

  • trunk/GSASIIplot.py

    r289 r296  
    13201320                        Data['ring'].append([xpos,ypos])
    13211321                elif event.button == 3:
    1322                     print 'LB shift'
    13231322                    self.dataFrame.GetStatusBar().SetStatusText('Calibrating...')
    13241323                    if G2img.ImageCalibrate(self,Data):

  • 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()

  • trunk/GSASIIpwdGUI.py

    r289 r296  
    6262        file.close()
    6363        self.dataFrame.UnDo.Enable(True)
     64       
     65    def OnLSQPeakFit(event):
     66        OnPeakFit('LSQ')
     67       
     68    def OnBGFSPeakFit(event):
     69        OnPeakFit('BFGS')
    6470               
    65     def OnPeakFit(event):
     71    def OnPeakFit(FitPgm):
    6672        SaveState()
    67         print 'Peak Fitting - Do one cycle of peak fitting'
     73        print 'Peak Fitting:'
    6874        PatternId = self.PatternId
    6975        PickId = self.PickId
    … …  
    7682        inst = self.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(self,PatternId, 'Instrument Parameters'))
    7783        data = self.PatternTree.GetItemPyData(PatternId)[1]
    78         OK,smin,Rwp,runtime,GoOn = G2pwd.DoPeakFit(peaks,background,limits,inst,data)
     84        wx.BeginBusyCursor()
     85        try:
     86            G2pwd.DoPeakFit(FitPgm,peaks,background,limits,inst,data)
     87        finally:
     88            wx.EndBusyCursor()   
    7989        UpdatePeakGrid(self,peaks)
    8090        G2plt.PlotPatterns(self)
    81         if not OK:
    82             print 'Refinement failed'
    83             dlg = wx.MessageDialog(self, 'Do you want to reload now?', 'Refinement failed',  wx.YES_NO)
    84             try:
    85                 if dlg.ShowModal() == wx.ID_YES:
    86                     DoUnDo()
    87                     self.dataFrame.UnDo.Enable(False)
    88             finally:
    89                 dlg.Destroy()
    90         else:
    91             self.dataFrame.UnDo.Enable(True)
    92             print "%s%7.2f%s%12.6g" % ('Rwp = ',Rwp,'%, Smin = ',smin)
    93             print "%s%8.3f%s " % ('fitpeak time =',runtime,'s')
    94             print 'finished'
     91        print 'finished'
    9592        return
    9693       
    97     def OnAutoPeakFit(event):
    98         SaveState()
    99         print 'AutoPeak Fitting - run until minimized'
     94    def OnResetSigGam(event):
    10095        PatternId = self.PatternId
    10196        PickId = self.PickId
    10297        peaks = self.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(self,PatternId, 'Peak List'))
    10398        if not peaks:
    104             self.ErrorDialog('No peaks!','Nothing to fit!')
     99            self.ErrorDialog('No peaks!','Nothing to do!')
    105100            return
    106         background = self.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(self,PatternId, 'Background'))[0]
    107         limits = self.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(self,PatternId, 'Limits'))[1]
    108101        inst = self.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(self,PatternId, 'Instrument Parameters'))
    109         data = self.PatternTree.GetItemPyData(PatternId)[1]
    110         smin = 1.0e15
    111         GoOn = True
    112         while GoOn:
    113             osmin = smin
    114             OK,smin,Rwp,runtime,GoOn = G2pwd.DoPeakFit(peaks,background,limits,inst,data)
    115             UpdatePeakGrid(self,peaks)
    116             if not OK:
    117                 break
    118             G2plt.PlotPatterns(self)
    119             print "%s%7.2f%s%12.6g" % ('Rwp = ',Rwp,'%, Smin = ',smin)
    120             rat = (osmin-smin)/smin
    121             if rat < 1.0e-4: GoOn = False
    122         if not OK:
    123             print 'Refinement failed'
    124             dlg = wx.MessageDialog(self, 'Do you want to reload now?', 'Refinement failed',  wx.YES_NO)
    125             try:
    126                 if dlg.ShowModal() == wx.ID_YES:
    127                     DoUnDo()
    128                     self.dataFrame.UnDo.Enable(False)
    129             finally:
    130                 dlg.Destroy()
    131         else:
    132             self.dataFrame.UnDo.Enable(True)
    133             print "%s%8.3f%s " % ('fitpeak time =',runtime,'s per cycle')
    134             print 'finished'
    135         return       
    136 
     102        Inst = dict(zip(inst[3],inst[1]))
     103        print len(Inst['Type'])
     104        for peak in peaks:
     105            if Inst['Type'] in ['PXC','PNC']:
     106                peak[4] = Inst['U']*tand(peak[0]/2.0)**2+Inst['V']*tand(peak[0]/2.0)+Inst['W']
     107                peak[6] = Inst['X']/cosd(peak[0]/2.0)+Inst['Y']*tand(peak[0]/2.0)
     108        UpdatePeakGrid(self,peaks)
     109               
    137110    def RefreshPeakGrid(event):
    138111        r,c =  event.GetRow(),event.GetCol()
    … …  
    223196                if not len(self.PatternTree.GetItemPyData(self.PickId)):
    224197                    self.dataFrame.PeakFit.Enable(False)
    225                     self.dataFrame.AutoPeakFit.Enable(False)
    226198                       
    227199        elif colList:
    … …  
    249221        Status = self.dataFrame.CreateStatusBar()
    250222    self.Bind(wx.EVT_MENU, OnUnDo, id=G2gd.wxID_UNDO)
    251     self.Bind(wx.EVT_MENU, OnPeakFit, id=G2gd.wxID_PEAKFIT)
    252     self.Bind(wx.EVT_MENU, OnAutoPeakFit, id=G2gd.wxID_AUTOPEAKFIT)
     223    self.Bind(wx.EVT_MENU, OnLSQPeakFit, id=G2gd.wxID_LSQPEAKFIT)
     224    self.Bind(wx.EVT_MENU, OnBGFSPeakFit, id=G2gd.wxID_BFGSPEAKFIT)
     225    self.Bind(wx.EVT_MENU, OnResetSigGam, id=G2gd.wxID_RESETSIGGAM)
    253226    self.dataFrame.PeakFit.Enable(False)
    254     self.dataFrame.AutoPeakFit.Enable(False)
    255227    if data:
    256228        self.dataFrame.PeakFit.Enable(True)
    257         self.dataFrame.AutoPeakFit.Enable(True)
    258229    self.PickTable = []
    259230    rowLabels = []
Note: See TracChangeset for help on using the changeset viewer.