Changeset 1893


Ignore:
Timestamp:
Jun 17, 2015 11:47:40 AM (7 years ago)
Author:
toby
Message:

fix background fitting by interpolation

Location:
trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIpwd.py

    r1889 r1893  
    13771377               
    13781378    def BackgroundPrint(Background,sigDict):
    1379         if Background[0][1]:
    1380             print 'Background coefficients for',Background[0][0],'function'
    1381             ptfmt = "%12.5f"
    1382             ptstr =  'values:'
    1383             sigstr = 'esds  :'
    1384             for i,back in enumerate(Background[0][3:]):
    1385                 ptstr += ptfmt % (back)
     1379        print 'Background coefficients for',Background[0][0],'function'
     1380        ptfmt = "%12.5f"
     1381        ptstr =  'value: '
     1382        sigstr = 'esd  : '
     1383        for i,back in enumerate(Background[0][3:]):
     1384            ptstr += ptfmt % (back)
     1385            if Background[0][1]:
    13861386                sigstr += ptfmt % (sigDict['Back;'+str(i)])
     1387            if len(ptstr) > 75:
     1388                print ptstr
     1389                if Background[0][1]: print sigstr
     1390                ptstr =  'value: '
     1391                sigstr = 'esd  : '
     1392        if len(ptstr) > 8:
    13871393            print ptstr
    1388             print sigstr
    1389         else:
    1390             print 'Background not refined'
     1394            if Background[0][1]: print sigstr
     1395
    13911396        if Background[1]['nDebye']:
    13921397            parms = ['DebyeA;','DebyeR;','DebyeU;']
     
    14021407                print line
    14031408        if Background[1]['nPeaks']:
    1404             parms = ['BkPkpos','BkPkint','BkPksig','BkPkgam']
    1405             print 'Peaks in background coefficients'
     1409            print 'Coefficients for Background Peaks'
    14061410            ptfmt = "%15.3f"
    1407             names =   'names :'
    1408             ptstr =  'values:'
    1409             sigstr = 'esds  :'
    1410             for item in sigDict:
    1411                 if 'BkPk' in item:
    1412                     names += '%15s'%(item)
    1413                     sigstr += ptfmt%(sigDict[item])
    1414                     parm,id = item.split(';')
    1415                     ip = parms.index(parm)
    1416                     ptstr += ptfmt%(Background[1]['peaksList'][int(id)][2*ip])
    1417             print names
    1418             print ptstr
    1419             print sigstr
     1411            for j,pl in enumerate(Background[1]['peaksList']):
     1412                names =  'peak %3d:'%(j+1)
     1413                ptstr =  'values  :'
     1414                sigstr = 'esds    :'
     1415                for i,lbl in enumerate(['BkPkpos','BkPkint','BkPksig','BkPkgam']):
     1416                    val = pl[2*i]
     1417                    prm = lbl+";"+str(j)
     1418                    names += '%15s'%(prm)
     1419                    ptstr += ptfmt%(val)
     1420                    if prm in sigDict:
     1421                        sigstr += ptfmt%(sigDict[prm])
     1422                    else:
     1423                        sigstr += " "*15
     1424                print names
     1425                print ptstr
     1426                print sigstr
    14201427                           
    14211428    def SetInstParms(Inst):
     
    15621569        return np.sqrt(weights)*getPeakProfileDerv(dataType,parmdict,xdata,varylist,bakType)
    15631570           
    1564     def errPeakProfile(values,xdata,ydata, weights,dataType,parmdict,varylist,bakType,dlg):       
     1571    def errPeakProfile(values,xdata,ydata,weights,dataType,parmdict,varylist,bakType,dlg):       
    15651572        parmdict.update(zip(varylist,values))
    15661573        M = np.sqrt(weights)*(getPeakProfile(dataType,parmdict,xdata,varylist,bakType)-ydata)
     
    16081615            try:
    16091616                result = so.leastsq(errPeakProfile,values,Dfun=devPeakProfile,full_output=True,ftol=Ftol,col_deriv=True,
    1610                     args=(x[xBeg:xFin],y[xBeg:xFin],w[xBeg:xFin],dataType,parmDict,varyList,bakType,dlg))
     1617                   args=(x[xBeg:xFin],y[xBeg:xFin],w[xBeg:xFin],dataType,parmDict,varyList,bakType,dlg))
    16111618                ncyc = int(result[2]['nfev']/2)
    16121619            finally:
  • trunk/GSASIIpwdGUI.py

    r1889 r1893  
    2828import random as ran
    2929import cPickle
     30import scipy.interpolate as si
    3031import GSASIIpath
    3132GSASIIpath.SetVersionNumber("$Revision$")
     
    814815        limits = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,PatternId, 'Limits'))[1]
    815816        inst,inst2 = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,PatternId, 'Instrument Parameters'))
    816        
    817         X = [x for x,y in background[1]['FixedPoints'] if limits[0] <= x <= limits[1]]
    818         Y = [y for x,y in background[1]['FixedPoints'] if limits[0] <= x <= limits[1]]
    819         W = [1]*len(X)
    820         Z = [0]*len(X)
     817        # sort the points for convenience and then separate them; extend the range if needed
     818        background[1]['FixedPoints'] = sorted(background[1]['FixedPoints'],key=lambda pair:pair[0])       
     819        X = [x for x,y in background[1]['FixedPoints']]
     820        Y = [y for x,y in background[1]['FixedPoints']]
     821        if X[0] > limits[0]:
     822            X = [limits[0]] + X
     823            Y = [Y[0]] + Y
     824        if X[-1] < limits[1]:
     825            X += [limits[1]]
     826            Y += [Y[-1]]
     827        # interpolate the fixed points onto the grid of data points within limits
     828        pwddata = G2frame.PatternTree.GetItemPyData(PatternId)[1]
     829        xBeg = np.searchsorted(pwddata[0],limits[0])
     830        xFin = np.searchsorted(pwddata[0],limits[1])
     831        xdata = pwddata[0][xBeg:xFin]
     832        ydata = si.interp1d(X,Y)(xdata)
     833        #GSASIIpath.IPyBreak()
     834        W = [1]*len(xdata)
     835        Z = [0]*len(xdata)
    821836
    822837        # load instrument and background params
     
    835850        try:
    836851            G2pwd.DoPeakFit('LSQ',[],background,limits,inst,inst2,
    837                             np.array((X,Y,W,Z,Z,Z)),bakVary,False,controls)
     852                            np.array((xdata,ydata,W,Z,Z,Z)),bakVary,False,controls)
    838853        finally:
    839854            wx.EndBusyCursor()
    840855        # compute the background values and plot them
    841         pwddata = G2frame.PatternTree.GetItemPyData(PatternId)[1]
    842         x = pwddata[0]
    843         xBeg = np.searchsorted(x,limits[0])
    844         xFin = np.searchsorted(x,limits[1])
    845856        parmDict = {}
    846857        bakType,bakDict,bakVary = G2pwd.SetBackgroundParms(background)
     
    850861        pwddata[5] *= 0
    851862        pwddata[4][xBeg:xFin] = G2pwd.getBackground(
    852             '',parmDict,bakType,dataType,x[xBeg:xFin])[0]
     863            '',parmDict,bakType,dataType,xdata)[0]
    853864        G2plt.PlotPatterns(G2frame,plotType='PWDR')
    854865        # show the updated background values
Note: See TracChangeset for help on using the changeset viewer.