Changeset 373


Ignore:
Timestamp:
Sep 16, 2011 3:51:06 PM (10 years ago)
Author:
vondreele
Message:

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

Location:
trunk
Files:
6 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASII.py

    r372 r373  
    547547            sub = self.PatternTree.AppendItem(parent=self.root,text='Controls')
    548548            self.PatternTree.SetItemPyData(sub,{})
     549            sub = self.PatternTree.AppendItem(parent=self.root,text='Covariance')
     550            self.PatternTree.SetItemPyData(sub,{})
     551           
    549552               
    550553    class CopyDialog(wx.Dialog):
     
    10491052                    self.dirname = dlg.GetDirectory()
    10501053                    G2IO.ProjFileOpen(self)
     1054                    #patch
     1055                    if not G2gd.GetPatternTreeItemId(self,self.root,'Covariance'):
     1056                        sub = self.PatternTree.AppendItem(parent=self.root,text='Covariance')
     1057                        self.PatternTree.SetItemPyData(sub,{})
     1058                    #end patch
    10511059                    self.PatternTree.SetItemText(self.root,'Loaded Data: '+self.GSASprojectfile)
    10521060                    self.PatternTree.Expand(self.root)
     
    14231431    class ViewParmDialog(wx.Dialog):
    14241432        def __init__(self,parent,title,parmDict):
    1425             wx.Dialog.__init__(self,parent,-1,title,size=(250,430),
     1433            wx.Dialog.__init__(self,parent,-1,title,size=(260,430),
    14261434                pos=wx.DefaultPosition,style=wx.DEFAULT_DIALOG_STYLE)
    1427             panel = wx.Panel(self,size=(250,430))
     1435            panel = wx.Panel(self,size=(260,430))
    14281436            parmNames = parmDict.keys()
    14291437            parmNames.sort()
     
    14321440                parmData = parmDict[name]
    14331441                try:
    1434                     parmText += ' %s \t%12.4g \n'%(name.ljust(20)+'\t'+parmData[1],parmData[0])
     1442                    parmText += ' %s \t%12.4g \n'%(name.ljust(19)+'\t'+parmData[1],parmData[0])
    14351443                except TypeError:
    14361444                    pass
    14371445            parmTable = wx.TextCtrl(panel,-1,parmText,
    1438                 style=wx.TE_MULTILINE|wx.TE_READONLY,size=(240,400))
     1446                style=wx.TE_MULTILINE|wx.TE_READONLY,size=(250,400))
    14391447            mainSizer = wx.BoxSizer(wx.VERTICAL)
    14401448            mainSizer.Add(parmTable)
  • trunk/GSASIIgrid.py

    r372 r373  
    712712            Id = GetPatternTreeItemId(self,self.root, 'Notebook')
    713713            if Id: self.PatternTree.SetItemPyData(Id,data)
     714        if self.dataFrame.GetLabel() == 'Covariance':
     715            data = [self.dataDisplay.GetValue()]
     716            self.dataDisplay.Clear()
     717            Id = GetPatternTreeItemId(self,self.root, 'Covariance')
     718            if Id: self.PatternTree.SetItemPyData(Id,data)
    714719        if 'Phase Data for' in self.dataFrame.GetLabel():
    715720            if self.dataDisplay:
     
    750755            self.Refine.Enable(True)
    751756            UpdateControls(self,data)
     757        elif self.PatternTree.GetItemText(item) == 'Covariance':
     758            data = self.PatternTree.GetItemPyData(item)
     759            G2plt.PlotCovariance(self)
    752760        elif 'IMG' in self.PatternTree.GetItemText(item):
    753761            self.Image = item
  • trunk/GSASIIplot.py

    r360 r373  
    13611361    Page.canvas.draw()
    13621362
     1363def PlotCovariance(self):
     1364    Data = self.PatternTree.GetItemPyData(
     1365        G2gd.GetPatternTreeItemId(self,self.root, 'Covariance'))
     1366    if not Data:
     1367        print 'No covariance matrix available'
     1368        return
     1369    varyList = Data['varyList']
     1370    Xmax = len(varyList)
     1371    covArray = Data['covariance']
     1372
     1373    def OnPlotKeyPress(event):
     1374        newPlot = False
     1375        if event.key == 's':
     1376            choice = [m for m in mpl.cm.datad.keys() if not m.endswith("_r")]
     1377            choice.sort()
     1378            dlg = wx.SingleChoiceDialog(self,'Select','Color scheme',choice)
     1379            if dlg.ShowModal() == wx.ID_OK:
     1380                sel = dlg.GetSelection()
     1381                self.ContourColor = choice[sel]
     1382            else:
     1383                self.ContourColor = 'Spectral'
     1384            dlg.Destroy()
     1385        PlotCovariance(self)
     1386
     1387    def OnMotion(event):
     1388        if event.xdata and event.ydata:                 #avoid out of frame errors
     1389            xpos = int(event.xdata+.5)
     1390            ypos = int(event.ydata+.5)
     1391            if -1 < xpos < len(varyList) and -1 < ypos < len(varyList):
     1392                self.G2plotNB.status.SetFields(['Key: s to change colors',
     1393                    '%s - %s: %5.3f'%(varyList[xpos].ljust(19),varyList[ypos].ljust(19),covArray[xpos][ypos])])
     1394    try:
     1395        plotNum = self.G2plotNB.plotList.index('Covariance')
     1396        Page = self.G2plotNB.nb.GetPage(plotNum)
     1397        Page.figure.clf()
     1398        Plot = Page.figure.gca()
     1399        if not Page.IsShown():
     1400            Page.Show()
     1401    except ValueError:
     1402        Plot = self.G2plotNB.addMpl('Covariance').gca()
     1403        plotNum = self.G2plotNB.plotList.index('Covariance')
     1404        Page = self.G2plotNB.nb.GetPage(plotNum)
     1405        Page.canvas.mpl_connect('motion_notify_event', OnMotion)
     1406        Page.canvas.mpl_connect('key_press_event', OnPlotKeyPress)
     1407
     1408    Page.SetFocus()
     1409    self.G2plotNB.status.SetFields(['',''])   
     1410    acolor = mpl.cm.get_cmap(self.ContourColor)
     1411    Img = Plot.imshow(covArray,aspect='equal',cmap=acolor,interpolation='nearest',origin='lower')
     1412    colorBar = Page.figure.colorbar(Img)
     1413    Plot.set_title('Variance-Covariance matrix from LS refinement')
     1414    Plot.set_xlabel('Variable number')
     1415    Plot.set_ylabel('Variable number')
     1416    Page.canvas.draw()
     1417
    13631418           
    13641419def PlotExposedImage(self,newPlot=False,event=None):
  • trunk/GSASIIpwd.py

    r372 r373  
    621621            elif bakType == 'cosine':
    622622                yb += parmDict[key]*npcosd(xdata*iBak)
    623     elif bakType in ['interpolate',]:
     623    elif bakType in ['lin interpolate','inv interpolate','log interpolate',]:
    624624        if nBak == 1:
    625625            yb = np.ones_like(xdata)*parmDict[pfx+'Back:0']
     
    630630            yb = parmDict[pfx+'Back:0']*T1+parmDict[pfx+'Back:1']*T2
    631631        else:
    632             bakPos = np.linspace(xdata[0],xdata[-1],nBak,True)
     632            if bakType == 'lin interpolate':
     633                bakPos = np.linspace(xdata[0],xdata[-1],nBak,True)
     634            elif bakType == 'inv interpolate':
     635                bakPos = 1./np.linspace(1./xdata[-1],1./xdata[0],nBak,True)
     636            elif bakType == 'log interpolate':
     637                bakPos = np.exp(np.linspace(np.log(xdata[0]),np.log(xdata[-1]),nBak,True))
     638            bakPos[0] = xdata[0]
     639            bakPos[-1] = xdata[-1]
    633640            bakVals = np.zeros(nBak)
    634641            for i in range(nBak):
     
    654661            elif bakType == 'cosine':
    655662                dydb[iBak] = npcosd(xdata*iBak)
    656     elif bakType in ['interpolate',]:
     663    elif bakType in ['lin interpolate','inv interpolate','log interpolate',]:
    657664        if nBak == 1:
    658665            dydb[0] = np.ones_like(xdata)
     
    663670            dydb = [T1,T2]
    664671        else:
    665             bakPos = np.linspace(xdata[0],xdata[-1],nBak,True)
     672            if bakType == 'lin interpolate':
     673                bakPos = np.linspace(xdata[0],xdata[-1],nBak,True)
     674            elif bakType == 'inv interpolate':
     675                bakPos = 1./np.linspace(1./xdata[-1],1./xdata[0],nBak,True)
     676            elif bakType == 'log interpolate':
     677                bakPos = np.exp(np.linspace(np.log(xdata[0]),np.log(xdata[-1]),nBak,True))
     678            bakPos[0] = xdata[0]
     679            bakPos[-1] = xdata[-1]
    666680            dx = bakPos[1]-bakPos[0]
    667681            for i,pos in enumerate(bakPos):
    668682                if i == 0:
    669                     dydb[0] = np.where(xdata<bakPos[1],(bakPos[1]-xdata)/dx,0.)
     683                    dydb[0] = np.where(xdata<bakPos[1],(bakPos[1]-xdata)/(bakPos[1]-bakPos[0]),0.)
    670684                elif i == len(bakPos)-1:
    671                     dydb[i] = np.where(xdata>bakPos[-2],(bakPos[-1]-xdata)/dx,0.)
     685                    dydb[i] = np.where(xdata>bakPos[-2],(bakPos[-1]-xdata)/(bakPos[-1]-bakPos[-2]),0.)
    672686                else:
    673687                    dydb[i] = np.where(xdata>bakPos[i],
    674                         np.where(xdata<bakPos[i+1],(bakPos[i+1]-xdata)/dx,0.),
    675                         np.where(xdata>bakPos[i-1],(xdata-bakPos[i-1])/dx,0.))
     688                        np.where(xdata<bakPos[i+1],(bakPos[i+1]-xdata)/(bakPos[i+1]-bakPos[i]),0.),
     689                        np.where(xdata>bakPos[i-1],(xdata-bakPos[i-1])/(bakPos[i]-bakPos[i-1]),0.))
    676690    return dydb
    677691
  • trunk/GSASIIpwdGUI.py

    r372 r373  
    286286    self.dataDisplay = wx.Panel(self.dataFrame)
    287287    BackId = G2gd.GetPatternTreeItemId(self,self.PatternId, 'Background')
    288     Choices = ['chebyschev','cosine','interpolate']
     288    Choices = ['chebyschev','cosine','lin interpolate','inv interpolate','log interpolate']
    289289    mainSizer = wx.BoxSizer(wx.VERTICAL)
    290290    topSizer = wx.BoxSizer(wx.HORIZONTAL)
     
    618618    elif data['Type'] == 'Bragg-Brentano':
    619619        parms += [['Shift',' Sample displacement(\xb5m): ','%.2f',],
    620             ['Transparency',' Sample transparency(1/\xb5eff,\xc5): ','%.4f'],]
     620            ['Transparency',' Sample transparency(1/\xb5eff,cm): ','%.4f'],]
    621621    parms.append(['Temperature',' Sample temperature(K): ','%.2f'])
    622622    parms.append(['Pressure',' Sample pressure(MPa): ','%.3f'])
  • trunk/GSASIIstruct.py

    r372 r373  
    162162    return Histograms,Phases
    163163   
    164 def SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases):
     164def SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases,CovData):
    165165    ''' Updates gpxfile from all histograms that are found in any phase
    166166    and any phase that used a histogram
     
    169169        Histograms = dictionary of histograms as {name:data,...}
    170170        Phases = dictionary of phases that use histograms
     171        CovData = [varyList,covariance matrix]
    171172    '''
    172173                       
     
    203204                    phaseName = data[iphase][0]
    204205                    data[iphase][1] = Phases[phaseName]
     206        elif datum[0] == 'Covariance':
     207            print data
     208            varyList,covMatrix = CovData
     209            data[0][1] = {'varyList':varyList,'covariance':covMatrix}
    205210        try:
    206211            histogram = Histograms[datum[0]]
     
    839844        print '\n Background function: ',Background[0],' Refine?',bool(Background[1])
    840845        line = ' Coefficients: '
    841         for back in Background[3:]:
     846        for i,back in enumerate(Background[3:]):
    842847            line += '%10.3f'%(back)
     848            if i and not i%10:
     849                line += '\n'+15*' '
    843850        print line
    844851       
     
    11941201        if 'Bragg' in calcControls[hfx+'instType']:
    11951202            pos -= const*(4.*parmDict[hfx+'Shift']*cosd(pos/2.0)+ \
    1196                 1.e-7*parmDict[hfx+'Transparency']*sind(pos))            #trans(=1/mueff) in Angstroms
     1203                parmDict[hfx+'Transparency']*sind(pos)*100.0)            #trans(=1/mueff) in cm
    11971204        else:               #Debye-Scherrer - simple but maybe not right
    11981205            pos -= const*(parmDict[hfx+'DisplaceX']*cosd(pos)+parmDict[hfx+'DisplaceY']*sind(pos))
     
    13481355        if 'Bragg' in calcControls[hfx+'instType']:
    13491356            dpdSh = -4.*const*cosd(pos/2.0)
    1350             dpdTr = -1.e-7*const*sind(pos)
     1357            dpdTr = -const*sind(pos)*100.0
    13511358            return dpdA,dpdw,dpdZ,dpdSh,dpdTr,0.,0.
    13521359        else:               #Debye-Scherrer - simple but maybe not right
     
    15101517                   
    15111518def Refine(GPXfile,dlg):
     1519    import cPickle
    15121520   
    15131521    def dervRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg):
     
    16281636        print 'Rwp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rwp,chisq,GOF)
    16291637        try:
    1630             sig = np.sqrt(np.diag(result[1])*GOF)
     1638            covMatrix = result[1]*GOF
     1639            sig = np.sqrt(np.diag(covMatrix))
     1640            xvar = np.outer(sig,np.ones_like(sig))
     1641            cov = np.divide(np.divide(covMatrix,xvar),xvar.T)
    16311642            if np.any(np.isnan(sig)):
    16321643                print '*** Least squares aborted - some invalid esds possible ***'
     
    16431654                    break
    16441655
     1656#    print 'dependentParmList: ',G2mv.dependentParmList
     1657#    print 'arrayList: ',G2mv.arrayList
     1658#    print 'invarrayList: ',G2mv.invarrayList
     1659#    print 'indParmList: ',G2mv.indParmList
     1660#    print 'fixedDict: ',G2mv.fixedDict
     1661    covFile = ospath.splitext(GPXfile)[0]+'.gpxcov'
     1662    file = open(covFile,'wb')
     1663    cPickle.dump(varyList,file,1)
     1664    cPickle.dump(result[0],file,1)
     1665    cPickle.dump(covMatrix,file,1)
     1666    file.close()
    16451667    sigDict = dict(zip(varyList,sig))
    16461668    SetPhaseData(parmDict,sigDict,Phases)
    16471669    SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms)
    16481670    SetHistogramData(parmDict,sigDict,Histograms)
    1649     SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases)
     1671    SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases,[varyList,cov])
    16501672#for testing purposes!!!
    1651     import cPickle
    1652     file = open('structTestdata.dat','wb')
    1653     cPickle.dump(parmDict,file,1)
    1654     cPickle.dump(varyList,file,1)
    1655     for histogram in Histograms:
    1656         if 'PWDR' in histogram[:4]:
    1657             Histogram = Histograms[histogram]
    1658     cPickle.dump(Histogram,file,1)
    1659     cPickle.dump(Phases,file,1)
    1660     cPickle.dump(calcControls,file,1)
    1661     cPickle.dump(pawleyLookup,file,1)
    1662     file.close()
     1673#    file = open('structTestdata.dat','wb')
     1674#    cPickle.dump(parmDict,file,1)
     1675#    cPickle.dump(varyList,file,1)
     1676#    for histogram in Histograms:
     1677#        if 'PWDR' in histogram[:4]:
     1678#            Histogram = Histograms[histogram]
     1679#    cPickle.dump(Histogram,file,1)
     1680#    cPickle.dump(Phases,file,1)
     1681#    cPickle.dump(calcControls,file,1)
     1682#    cPickle.dump(pawleyLookup,file,1)
     1683#    file.close()
    16631684
    16641685def main():
Note: See TracChangeset for help on using the changeset viewer.