Changeset 1893

Timestamp:

Jun 17, 2015 11:47:40 AM (7 years ago)

Author:

toby

Message:

fix background fitting by interpolation

Location:

trunk

Files:

• GSASIIpwd.py (modified) (4 diffs)
• GSASIIpwdGUI.py (modified) (4 diffs)

trunk/GSASIIpwd.py

@@ -1377 +1377 @@
                 
     def BackgroundPrint(Background,sigDict):
-        if Background[0][1]:
-            print 'Background coefficients for',Background[0][0],'function'
-            ptfmt = "%12.5f"
-            ptstr =  'values:'
-            sigstr = 'esds  :'
-            for i,back in enumerate(Background[0][3:]):
-                ptstr += ptfmt % (back)
+        print 'Background coefficients for',Background[0][0],'function'
+        ptfmt = "%12.5f"
+        ptstr =  'value: '
+        sigstr = 'esd  : '
+        for i,back in enumerate(Background[0][3:]):
+            ptstr += ptfmt % (back)
+            if Background[0][1]:
                 sigstr += ptfmt % (sigDict['Back;'+str(i)])
+            if len(ptstr) > 75:
+                print ptstr
+                if Background[0][1]: print sigstr
+                ptstr =  'value: '
+                sigstr = 'esd  : '
+        if len(ptstr) > 8:
             print ptstr
-            print sigstr
-        else:
-            print 'Background not refined'
+            if Background[0][1]: print sigstr
+
         if Background[1]['nDebye']:
             parms = ['DebyeA;','DebyeR;','DebyeU;']
@@ -1402 +1407 @@
                 print line
         if Background[1]['nPeaks']:
-            parms = ['BkPkpos','BkPkint','BkPksig','BkPkgam']
-            print 'Peaks in background coefficients'
+            print 'Coefficients for Background Peaks'
             ptfmt = "%15.3f"
-            names =   'names :'
-            ptstr =  'values:'
-            sigstr = 'esds  :'
-            for item in sigDict:
-                if 'BkPk' in item:
-                    names += '%15s'%(item)
-                    sigstr += ptfmt%(sigDict[item])
-                    parm,id = item.split(';')
-                    ip = parms.index(parm)
-                    ptstr += ptfmt%(Background[1]['peaksList'][int(id)][2*ip])
-            print names
-            print ptstr
-            print sigstr
+            for j,pl in enumerate(Background[1]['peaksList']):
+                names =  'peak %3d:'%(j+1)
+                ptstr =  'values  :'
+                sigstr = 'esds    :'
+                for i,lbl in enumerate(['BkPkpos','BkPkint','BkPksig','BkPkgam']):
+                    val = pl[2*i]
+                    prm = lbl+";"+str(j)
+                    names += '%15s'%(prm)
+                    ptstr += ptfmt%(val)
+                    if prm in sigDict:
+                        sigstr += ptfmt%(sigDict[prm])
+                    else:
+                        sigstr += " "*15
+                print names
+                print ptstr
+                print sigstr
                             
     def SetInstParms(Inst):
@@ -1562 +1569 @@
         return np.sqrt(weights)*getPeakProfileDerv(dataType,parmdict,xdata,varylist,bakType)
             
-    def errPeakProfile(values,xdata,ydata, weights,dataType,parmdict,varylist,bakType,dlg):        
+    def errPeakProfile(values,xdata,ydata,weights,dataType,parmdict,varylist,bakType,dlg):        
         parmdict.update(zip(varylist,values))
         M = np.sqrt(weights)*(getPeakProfile(dataType,parmdict,xdata,varylist,bakType)-ydata)
@@ -1608 +1615 @@
             try:
                 result = so.leastsq(errPeakProfile,values,Dfun=devPeakProfile,full_output=True,ftol=Ftol,col_deriv=True,
-                    args=(x[xBeg:xFin],y[xBeg:xFin],w[xBeg:xFin],dataType,parmDict,varyList,bakType,dlg))
+                   args=(x[xBeg:xFin],y[xBeg:xFin],w[xBeg:xFin],dataType,parmDict,varyList,bakType,dlg))
                 ncyc = int(result[2]['nfev']/2)
             finally:

trunk/GSASIIpwdGUI.py

@@ -28 +28 @@
 import random as ran
 import cPickle
+import scipy.interpolate as si
 import GSASIIpath
 GSASIIpath.SetVersionNumber("$Revision$")
@@ -814 +815 @@
         limits = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,PatternId, 'Limits'))[1]
         inst,inst2 = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,PatternId, 'Instrument Parameters'))
-        
-        X = [x for x,y in background[1]['FixedPoints'] if limits[0] <= x <= limits[1]]
-        Y = [y for x,y in background[1]['FixedPoints'] if limits[0] <= x <= limits[1]]
-        W = [1]*len(X)
-        Z = [0]*len(X)
+        # sort the points for convenience and then separate them; extend the range if needed
+        background[1]['FixedPoints'] = sorted(background[1]['FixedPoints'],key=lambda pair:pair[0])        
+        X = [x for x,y in background[1]['FixedPoints']]
+        Y = [y for x,y in background[1]['FixedPoints']]
+        if X[0] > limits[0]:
+            X = [limits[0]] + X
+            Y = [Y[0]] + Y
+        if X[-1] < limits[1]:
+            X += [limits[1]]
+            Y += [Y[-1]]
+        # interpolate the fixed points onto the grid of data points within limits
+        pwddata = G2frame.PatternTree.GetItemPyData(PatternId)[1]
+        xBeg = np.searchsorted(pwddata[0],limits[0])
+        xFin = np.searchsorted(pwddata[0],limits[1])
+        xdata = pwddata[0][xBeg:xFin]
+        ydata = si.interp1d(X,Y)(xdata)
+        #GSASIIpath.IPyBreak()
+        W = [1]*len(xdata)
+        Z = [0]*len(xdata)
 
         # load instrument and background params
@@ -835 +850 @@
         try:
             G2pwd.DoPeakFit('LSQ',[],background,limits,inst,inst2,
-                            np.array((X,Y,W,Z,Z,Z)),bakVary,False,controls)
+                            np.array((xdata,ydata,W,Z,Z,Z)),bakVary,False,controls)
         finally:
             wx.EndBusyCursor()
         # compute the background values and plot them
-        pwddata = G2frame.PatternTree.GetItemPyData(PatternId)[1]
-        x = pwddata[0]
-        xBeg = np.searchsorted(x,limits[0])
-        xFin = np.searchsorted(x,limits[1])
         parmDict = {}
         bakType,bakDict,bakVary = G2pwd.SetBackgroundParms(background)
@@ -850 +861 @@
         pwddata[5] *= 0
         pwddata[4][xBeg:xFin] = G2pwd.getBackground(
-            '',parmDict,bakType,dataType,x[xBeg:xFin])[0]
+            '',parmDict,bakType,dataType,xdata)[0]
         G2plt.PlotPatterns(G2frame,plotType='PWDR')
         # show the updated background values