Changeset 373

Timestamp:

Sep 16, 2011 3:51:06 PM (12 years ago)

Author:

vondreele

Message:

new covariance item in tree
new covariance matrix plotting
save cov matrix in file

Location:

trunk

Files:

• GSASII.py (modified) (4 diffs)
• GSASIIgrid.py (modified) (2 diffs)
• GSASIIplot.py (modified) (1 diff)
• GSASIIpwd.py (modified) (4 diffs)
• GSASIIpwdGUI.py (modified) (2 diffs)
• GSASIIstruct.py (modified) (9 diffs)

trunk/GSASII.py

@@ -547 +547 @@
             sub = self.PatternTree.AppendItem(parent=self.root,text='Controls')
             self.PatternTree.SetItemPyData(sub,{})
+            sub = self.PatternTree.AppendItem(parent=self.root,text='Covariance')
+            self.PatternTree.SetItemPyData(sub,{})
+            
                 
     class CopyDialog(wx.Dialog):
@@ -1049 +1052 @@
                     self.dirname = dlg.GetDirectory()
                     G2IO.ProjFileOpen(self)
+                    #patch
+                    if not G2gd.GetPatternTreeItemId(self,self.root,'Covariance'):
+                        sub = self.PatternTree.AppendItem(parent=self.root,text='Covariance')
+                        self.PatternTree.SetItemPyData(sub,{})
+                    #end patch
                     self.PatternTree.SetItemText(self.root,'Loaded Data: '+self.GSASprojectfile)
                     self.PatternTree.Expand(self.root)
@@ -1423 +1431 @@
     class ViewParmDialog(wx.Dialog):
         def __init__(self,parent,title,parmDict):
-            wx.Dialog.__init__(self,parent,-1,title,size=(250,430),
+            wx.Dialog.__init__(self,parent,-1,title,size=(260,430),
                 pos=wx.DefaultPosition,style=wx.DEFAULT_DIALOG_STYLE)
-            panel = wx.Panel(self,size=(250,430))
+            panel = wx.Panel(self,size=(260,430))
             parmNames = parmDict.keys()
             parmNames.sort()
@@ -1432 +1440 @@
                 parmData = parmDict[name]
                 try:
-                    parmText += ' %s \t%12.4g \n'%(name.ljust(20)+'\t'+parmData[1],parmData[0])
+                    parmText += ' %s \t%12.4g \n'%(name.ljust(19)+'\t'+parmData[1],parmData[0])
                 except TypeError:
                     pass
             parmTable = wx.TextCtrl(panel,-1,parmText,
-                style=wx.TE_MULTILINE|wx.TE_READONLY,size=(240,400))
+                style=wx.TE_MULTILINE|wx.TE_READONLY,size=(250,400))
             mainSizer = wx.BoxSizer(wx.VERTICAL)
             mainSizer.Add(parmTable)

trunk/GSASIIgrid.py

@@ -712 +712 @@
             Id = GetPatternTreeItemId(self,self.root, 'Notebook')
             if Id: self.PatternTree.SetItemPyData(Id,data)
+        if self.dataFrame.GetLabel() == 'Covariance':
+            data = [self.dataDisplay.GetValue()]
+            self.dataDisplay.Clear() 
+            Id = GetPatternTreeItemId(self,self.root, 'Covariance')
+            if Id: self.PatternTree.SetItemPyData(Id,data)
         if 'Phase Data for' in self.dataFrame.GetLabel():
             if self.dataDisplay: 
@@ -750 +755 @@
             self.Refine.Enable(True)
             UpdateControls(self,data)
+        elif self.PatternTree.GetItemText(item) == 'Covariance':
+            data = self.PatternTree.GetItemPyData(item)
+            G2plt.PlotCovariance(self)
         elif 'IMG' in self.PatternTree.GetItemText(item):
             self.Image = item

trunk/GSASIIplot.py

@@ -1361 +1361 @@
     Page.canvas.draw()
 
+def PlotCovariance(self):
+    Data = self.PatternTree.GetItemPyData(
+        G2gd.GetPatternTreeItemId(self,self.root, 'Covariance'))
+    if not Data:
+        print 'No covariance matrix available'
+        return
+    varyList = Data['varyList']
+    Xmax = len(varyList)
+    covArray = Data['covariance']
+
+    def OnPlotKeyPress(event):
+        newPlot = False
+        if event.key == 's':
+            choice = [m for m in mpl.cm.datad.keys() if not m.endswith("_r")]
+            choice.sort()
+            dlg = wx.SingleChoiceDialog(self,'Select','Color scheme',choice)
+            if dlg.ShowModal() == wx.ID_OK:
+                sel = dlg.GetSelection()
+                self.ContourColor = choice[sel]
+            else:
+                self.ContourColor = 'Spectral'
+            dlg.Destroy()
+        PlotCovariance(self)
+
+    def OnMotion(event):
+        if event.xdata and event.ydata:                 #avoid out of frame errors
+            xpos = int(event.xdata+.5)
+            ypos = int(event.ydata+.5)
+            if -1 < xpos < len(varyList) and -1 < ypos < len(varyList):
+                self.G2plotNB.status.SetFields(['Key: s to change colors',
+                    '%s - %s: %5.3f'%(varyList[xpos].ljust(19),varyList[ypos].ljust(19),covArray[xpos][ypos])])
+    try:
+        plotNum = self.G2plotNB.plotList.index('Covariance')
+        Page = self.G2plotNB.nb.GetPage(plotNum)
+        Page.figure.clf()
+        Plot = Page.figure.gca()
+        if not Page.IsShown():
+            Page.Show()
+    except ValueError:
+        Plot = self.G2plotNB.addMpl('Covariance').gca()
+        plotNum = self.G2plotNB.plotList.index('Covariance')
+        Page = self.G2plotNB.nb.GetPage(plotNum)
+        Page.canvas.mpl_connect('motion_notify_event', OnMotion)
+        Page.canvas.mpl_connect('key_press_event', OnPlotKeyPress)
+
+    Page.SetFocus()
+    self.G2plotNB.status.SetFields(['',''])    
+    acolor = mpl.cm.get_cmap(self.ContourColor)
+    Img = Plot.imshow(covArray,aspect='equal',cmap=acolor,interpolation='nearest',origin='lower')
+    colorBar = Page.figure.colorbar(Img)
+    Plot.set_title('Variance-Covariance matrix from LS refinement')
+    Plot.set_xlabel('Variable number')
+    Plot.set_ylabel('Variable number')
+    Page.canvas.draw()
+
             
 def PlotExposedImage(self,newPlot=False,event=None):

trunk/GSASIIpwd.py

@@ -621 +621 @@
             elif bakType == 'cosine':
                 yb += parmDict[key]*npcosd(xdata*iBak)
-    elif bakType in ['interpolate',]:
+    elif bakType in ['lin interpolate','inv interpolate','log interpolate',]:
         if nBak == 1:
             yb = np.ones_like(xdata)*parmDict[pfx+'Back:0']
@@ -630 +630 @@
             yb = parmDict[pfx+'Back:0']*T1+parmDict[pfx+'Back:1']*T2
         else:
-            bakPos = np.linspace(xdata[0],xdata[-1],nBak,True)
+            if bakType == 'lin interpolate':
+                bakPos = np.linspace(xdata[0],xdata[-1],nBak,True)
+            elif bakType == 'inv interpolate':
+                bakPos = 1./np.linspace(1./xdata[-1],1./xdata[0],nBak,True)
+            elif bakType == 'log interpolate':
+                bakPos = np.exp(np.linspace(np.log(xdata[0]),np.log(xdata[-1]),nBak,True))
+            bakPos[0] = xdata[0]
+            bakPos[-1] = xdata[-1]
             bakVals = np.zeros(nBak)
             for i in range(nBak):
@@ -654 +661 @@
             elif bakType == 'cosine':
                 dydb[iBak] = npcosd(xdata*iBak)
-    elif bakType in ['interpolate',]:
+    elif bakType in ['lin interpolate','inv interpolate','log interpolate',]:
         if nBak == 1:
             dydb[0] = np.ones_like(xdata)
@@ -663 +670 @@
             dydb = [T1,T2]
         else:
-            bakPos = np.linspace(xdata[0],xdata[-1],nBak,True)
+            if bakType == 'lin interpolate':
+                bakPos = np.linspace(xdata[0],xdata[-1],nBak,True)
+            elif bakType == 'inv interpolate':
+                bakPos = 1./np.linspace(1./xdata[-1],1./xdata[0],nBak,True)
+            elif bakType == 'log interpolate':
+                bakPos = np.exp(np.linspace(np.log(xdata[0]),np.log(xdata[-1]),nBak,True))
+            bakPos[0] = xdata[0]
+            bakPos[-1] = xdata[-1]
             dx = bakPos[1]-bakPos[0]
             for i,pos in enumerate(bakPos):
                 if i == 0:
-                    dydb[0] = np.where(xdata<bakPos[1],(bakPos[1]-xdata)/dx,0.)
+                    dydb[0] = np.where(xdata<bakPos[1],(bakPos[1]-xdata)/(bakPos[1]-bakPos[0]),0.)
                 elif i == len(bakPos)-1:
-                    dydb[i] = np.where(xdata>bakPos[-2],(bakPos[-1]-xdata)/dx,0.)
+                    dydb[i] = np.where(xdata>bakPos[-2],(bakPos[-1]-xdata)/(bakPos[-1]-bakPos[-2]),0.)
                 else:
                     dydb[i] = np.where(xdata>bakPos[i],
-                        np.where(xdata<bakPos[i+1],(bakPos[i+1]-xdata)/dx,0.),
-                        np.where(xdata>bakPos[i-1],(xdata-bakPos[i-1])/dx,0.))
+                        np.where(xdata<bakPos[i+1],(bakPos[i+1]-xdata)/(bakPos[i+1]-bakPos[i]),0.),
+                        np.where(xdata>bakPos[i-1],(xdata-bakPos[i-1])/(bakPos[i]-bakPos[i-1]),0.))
     return dydb
 

trunk/GSASIIpwdGUI.py

@@ -286 +286 @@
     self.dataDisplay = wx.Panel(self.dataFrame)
     BackId = G2gd.GetPatternTreeItemId(self,self.PatternId, 'Background')
-    Choices = ['chebyschev','cosine','interpolate']
+    Choices = ['chebyschev','cosine','lin interpolate','inv interpolate','log interpolate']
     mainSizer = wx.BoxSizer(wx.VERTICAL)
     topSizer = wx.BoxSizer(wx.HORIZONTAL)
@@ -618 +618 @@
     elif data['Type'] == 'Bragg-Brentano':
         parms += [['Shift',' Sample displacement(\xb5m): ','%.2f',],
-            ['Transparency',' Sample transparency(1/\xb5eff,\xc5): ','%.4f'],]
+            ['Transparency',' Sample transparency(1/\xb5eff,cm): ','%.4f'],]
     parms.append(['Temperature',' Sample temperature(K): ','%.2f'])
     parms.append(['Pressure',' Sample pressure(MPa): ','%.3f'])

trunk/GSASIIstruct.py

@@ -162 +162 @@
     return Histograms,Phases
     
-def SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases):
+def SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases,CovData):
     ''' Updates gpxfile from all histograms that are found in any phase
     and any phase that used a histogram
@@ -169 +169 @@
         Histograms = dictionary of histograms as {name:data,...}
         Phases = dictionary of phases that use histograms
+        CovData = [varyList,covariance matrix]
     '''
                         
@@ -203 +204 @@
                     phaseName = data[iphase][0]
                     data[iphase][1] = Phases[phaseName]
+        elif datum[0] == 'Covariance':
+            print data
+            varyList,covMatrix = CovData
+            data[0][1] = {'varyList':varyList,'covariance':covMatrix}
         try:
             histogram = Histograms[datum[0]]
@@ -839 +844 @@
         print '\n Background function: ',Background[0],' Refine?',bool(Background[1])
         line = ' Coefficients: '
-        for back in Background[3:]:
+        for i,back in enumerate(Background[3:]):
             line += '%10.3f'%(back)
+            if i and not i%10:
+                line += '\n'+15*' '
         print line 
         
@@ -1194 +1201 @@
         if 'Bragg' in calcControls[hfx+'instType']:
             pos -= const*(4.*parmDict[hfx+'Shift']*cosd(pos/2.0)+ \
-                1.e-7*parmDict[hfx+'Transparency']*sind(pos))            #trans(=1/mueff) in Angstroms
+                parmDict[hfx+'Transparency']*sind(pos)*100.0)            #trans(=1/mueff) in cm
         else:               #Debye-Scherrer - simple but maybe not right
             pos -= const*(parmDict[hfx+'DisplaceX']*cosd(pos)+parmDict[hfx+'DisplaceY']*sind(pos))
@@ -1348 +1355 @@
         if 'Bragg' in calcControls[hfx+'instType']:
             dpdSh = -4.*const*cosd(pos/2.0)
-            dpdTr = -1.e-7*const*sind(pos)
+            dpdTr = -const*sind(pos)*100.0
             return dpdA,dpdw,dpdZ,dpdSh,dpdTr,0.,0.
         else:               #Debye-Scherrer - simple but maybe not right
@@ -1510 +1517 @@
                     
 def Refine(GPXfile,dlg):
+    import cPickle
     
     def dervRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg):
@@ -1628 +1636 @@
         print 'Rwp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rwp,chisq,GOF)
         try:
-            sig = np.sqrt(np.diag(result[1])*GOF)
+            covMatrix = result[1]*GOF
+            sig = np.sqrt(np.diag(covMatrix))
+            xvar = np.outer(sig,np.ones_like(sig))
+            cov = np.divide(np.divide(covMatrix,xvar),xvar.T)
             if np.any(np.isnan(sig)):
                 print '*** Least squares aborted - some invalid esds possible ***'
@@ -1643 +1654 @@
                     break
 
+#    print 'dependentParmList: ',G2mv.dependentParmList
+#    print 'arrayList: ',G2mv.arrayList
+#    print 'invarrayList: ',G2mv.invarrayList
+#    print 'indParmList: ',G2mv.indParmList
+#    print 'fixedDict: ',G2mv.fixedDict
+    covFile = ospath.splitext(GPXfile)[0]+'.gpxcov'
+    file = open(covFile,'wb')
+    cPickle.dump(varyList,file,1)
+    cPickle.dump(result[0],file,1)
+    cPickle.dump(covMatrix,file,1)
+    file.close()
     sigDict = dict(zip(varyList,sig))
     SetPhaseData(parmDict,sigDict,Phases)
     SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms)
     SetHistogramData(parmDict,sigDict,Histograms)
-    SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases)
+    SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases,[varyList,cov])
 #for testing purposes!!!
-    import cPickle
-    file = open('structTestdata.dat','wb')
-    cPickle.dump(parmDict,file,1)
-    cPickle.dump(varyList,file,1)
-    for histogram in Histograms:
-        if 'PWDR' in histogram[:4]:
-            Histogram = Histograms[histogram]
-    cPickle.dump(Histogram,file,1)
-    cPickle.dump(Phases,file,1)
-    cPickle.dump(calcControls,file,1)
-    cPickle.dump(pawleyLookup,file,1)
-    file.close()
+#    file = open('structTestdata.dat','wb')
+#    cPickle.dump(parmDict,file,1)
+#    cPickle.dump(varyList,file,1)
+#    for histogram in Histograms:
+#        if 'PWDR' in histogram[:4]:
+#            Histogram = Histograms[histogram]
+#    cPickle.dump(Histogram,file,1)
+#    cPickle.dump(Phases,file,1)
+#    cPickle.dump(calcControls,file,1)
+#    cPickle.dump(pawleyLookup,file,1)
+#    file.close()
 
 def main():