source: trunk/G2compare.py @ 4785

Last change on this file since 4785 was 4785, checked in by toby, 9 months ago

new table for text w/sorting

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 36.5 KB
Line 
1#!/usr/bin/env python
2# -*- coding: utf-8 -*-
3#GSAS-II Data/Model Comparison
4########### SVN repository information ###################
5# $Date: 2021-01-27 21:27:25 +0000 (Wed, 27 Jan 2021) $
6# $Author: toby $
7# $Revision: 4785 $
8# $URL: trunk/G2compare.py $
9# $Id: G2compare.py 4785 2021-01-27 21:27:25Z toby $
10########### SVN repository information ###################
11'''
12.. _G2compare:
13
14*G2compare: Tool for project comparison*
15---------------------------------------------
16'''
17
18#TODO:
19# Prince-test next
20# Make Phase unique? (need a phaselist)
21# more graphics
22# display more in datawindow
23
24import sys
25import os
26import platform
27import glob
28if '2' in platform.python_version_tuple()[0]:
29    import cPickle
30else:
31    try:
32        import _pickle as cPickle
33    except:
34        print('Warning: failed to import the optimized Py3 pickle (_pickle)')
35        import pickle as cPickle
36
37import wx
38import numpy as np
39import matplotlib as mpl
40try:
41    import OpenGL as ogl
42except ImportError:
43    pass
44import scipy as sp
45
46import GSASIIpath
47GSASIIpath.SetVersionNumber("$Revision: 4785 $")
48import GSASIIfiles as G2fil
49import GSASIIplot as G2plt
50import GSASIIdataGUI as G2gd
51import GSASIIctrlGUI as G2G
52import GSASIIobj as G2obj
53
54__version__ = '0.0.1'
55
56# math to do F-test
57def RC2Ftest(npts,RChiSq0,nvar0,RChiSq1,nvar1):
58    '''Compute the F-test probability that a model expanded with added
59    parameters (relaxed model) is statistically more likely than the
60    constrained (base) model
61    :param int npts: number of observed diffraction data points
62    :param float RChiSq0: Reduced Chi**2 for the base model
63    :param int nvar0: number of refined variables in the base model
64    :param float RChiSq0: Reduced Chi**2 for the relaxed model
65    :param int nvar1: number of refined variables in the relaxed model
66    '''
67    if nvar1 == nvar0:
68        raise Exception("parameter # agree, test is not valid")
69    elif nvar1 < nvar0:
70        print("Warning: RwFtest used with base/relaxed models reversed")
71        RChiSq0,nvar0,RChiSq1,nvar1 = RChiSq1,nvar1,RChiSq0,nvar0
72    ratio = RChiSq0 / RChiSq1
73    nu1 = float(nvar1 - nvar0)
74    nu2 = float(npts - nvar1)
75    F = (ratio - 1.) * nu2 / nu1
76    import scipy.stats
77    return scipy.stats.f.cdf(F,nu1,nu2)
78
79def RwFtest(npts,Rwp0,nvar0,Rwp1,nvar1):
80    '''Compute the F-test probability that a model expanded with added
81    parameters (relaxed model) is statistically more likely than the
82    constrained (base) model
83    :param int npts: number of observed diffraction data points
84    :param float Rwp0: Weighted profile R-factor or GOF for the base model
85    :param int nvar0: number of refined variables in the base model
86    :param float Rwp1: Weighted profile R-factor or GOF for the relaxed model
87    :param int nvar1: number of refined variables in the relaxed model
88    '''
89    if nvar1 == nvar0:
90        raise Exception("parameter # agree, test is not valid")
91    elif nvar1 < nvar0:
92        print("Warning: RwFtest used with base/relaxed models reversed")
93        Rwp0,nvar0,Rwp1,nvar1 = Rwp1,nvar1,Rwp0,nvar0
94    ratio = (Rwp0 / Rwp1)**2
95    nu1 = float(nvar1 - nvar0)
96    nu2 = float(npts - nvar1)
97    F = (ratio - 1.) * nu2 / nu1
98    import scipy.stats
99    return scipy.stats.f.cdf(F,nu1,nu2)
100
101def cPickleLoad(fp):
102    if '2' in platform.python_version_tuple()[0]:
103        return cPickle.load(fp)
104    else:
105       return cPickle.load(fp,encoding='latin-1')
106           
107def main(application):
108    '''Start up the GSAS-II GUI'''                       
109    knownVersions = ['3.6','3.7','3.8','3.9']
110    if platform.python_version()[:3] not in knownVersions: 
111        dlg = wx.MessageDialog(None, 
112                'GSAS-II requires Python 3.6+\n Yours is '+sys.version.split()[0],
113                'Python version error',  wx.OK)
114        try:
115            dlg.ShowModal()
116        finally:
117            dlg.Destroy()
118        sys.exit()
119           
120    application.main = MakeTopWindow(None)  # application.main is the main wx.Frame
121    application.SetTopWindow(application.main)
122    # save the current package versions
123    application.main.PackageVersions = G2fil.get_python_versions([wx, mpl, np, sp, ogl])
124    try:
125        application.SetAppDisplayName('GSAS-II Compare')
126    except:
127        pass
128    #application.GetTopWindow().SendSizeEvent()
129    application.GetTopWindow().Show(True)
130    return application.GetTopWindow()
131   
132class MakeTopWindow(wx.Frame):
133    '''Define the main frame and its associated menu items
134    '''
135    def __init__(self, parent):
136        size = wx.Size(700,450)
137        wx.Frame.__init__(self, name='dComp', parent=parent,
138            size=size,style=wx.DEFAULT_FRAME_STYLE, title='GSAS-II data/model comparison')
139        # misc vars
140        self.VcovColor = 'RdYlGn'       
141        # plot window
142        try:
143            size = GSASIIpath.GetConfigValue('Plot_Size')
144            if type(size) is tuple:
145                pass
146            elif type(size) is str:
147                size = eval(size)
148            else:
149                raise Exception
150        except:
151            size = wx.Size(700,600)               
152        self.plotFrame = wx.Frame(None,-1,'dComp Plots',size=size,
153            style=wx.DEFAULT_FRAME_STYLE ^ wx.CLOSE_BOX)
154        self.G2plotNB = G2plt.G2PlotNoteBook(self.plotFrame,G2frame=self)
155        #self.plotFrame.Show()
156        self.plotFrame.Show(False)  # until we have some graphics, hide the plot window
157        # menus
158        Frame = self.GetTopLevelParent() # same as self
159        self.MenuBar = wx.MenuBar()
160        File = wx.Menu(title='')
161        self.MenuBar.Append(menu=File, title='&File')
162        item = File.Append(wx.ID_ANY,'&Import single project...\tCtrl+O','Open a GSAS-II project file (*.gpx)')           
163        self.Bind(wx.EVT_MENU, self.onLoadGPX, id=item.GetId())
164        item = File.Append(wx.ID_ANY,'&Import multiple projects...\tCtrl+M','Open a GSAS-II project file (*.gpx)')           
165        self.Bind(wx.EVT_MENU, self.onLoadMultGPX, id=item.GetId())
166        item = File.Append(wx.ID_ANY,'&Wildcard import of projects...\tCtrl+W','Open a GSAS-II project file (*.gpx)')           
167        self.Bind(wx.EVT_MENU, self.onLoadWildGPX, id=item.GetId())
168#        item = File.Append(wx.ID_ANY,'&Import selected...','Open a GSAS-II project file (*.gpx)')
169#        self.Bind(wx.EVT_MENU, self.onLoadSel, id=item.GetId())
170
171        self.Mode = wx.Menu(title='')
172        self.MenuBar.Append(menu=self.Mode, title='&Mode')
173        self.wxID_Mode = {}
174        for i,m in enumerate(("Histogram","Phase","Project")):
175            self.wxID_Mode[m] = i+1
176            item = self.Mode.AppendRadioItem(i+1,m+'\tCtrl+{}'.format(i+1),
177                                                 'Display {}s'.format(m))
178            self.Bind(wx.EVT_MENU, self.onRefresh, id=item.GetId())
179        item = self.Mode.Append(wx.ID_ANY,'Set histogram filter','Set a filter for histograms to display')
180        self.Bind(wx.EVT_MENU, self.onHistFilter, id=item.GetId())
181        modeMenu = wx.Menu(title='')
182        self.MenuBar.Append(menu=modeMenu, title='TBD')
183        self.modeMenuPos = self.MenuBar.FindMenu('TBD')
184        #item = modeMenu.Append(wx.ID_ANY,'test') # menu can't be empty
185       
186        Frame.SetMenuBar(self.MenuBar)
187        # status bar
188        self.Status = self.CreateStatusBar()
189        self.Status.SetFieldsCount(2)
190        # split the frame and add the tree
191        self.mainPanel = wx.SplitterWindow(self, wx.ID_ANY, style=wx.SP_LIVE_UPDATE|wx.SP_3D)
192        self.mainPanel.SetMinimumPaneSize(100)
193        self.treePanel = wx.Panel(self.mainPanel, wx.ID_ANY,
194            style = wx.TAB_TRAVERSAL|wx.SUNKEN_BORDER)
195       
196#        self.dataWindow = G2DataWindow(self.mainPanel)
197        self.dataWindow = wx.ScrolledWindow(self.mainPanel)
198        dataSizer = wx.BoxSizer(wx.VERTICAL)
199        self.dataWindow.SetSizer(dataSizer)
200        self.mainPanel.SplitVertically(self.treePanel, self.dataWindow, 200)
201        self.Status.SetStatusWidths([200,-1])   # make these match?
202       
203        treeSizer = wx.BoxSizer(wx.VERTICAL)
204        self.treePanel.SetSizer(treeSizer)
205        self.GPXtree = G2G.G2TreeCtrl(id=wx.ID_ANY,
206            parent=self.treePanel, size=self.treePanel.GetClientSize(),style=wx.TR_DEFAULT_STYLE )
207        #TreeId = self.GPXtree.Id
208
209        treeSizer.Add(self.GPXtree,1,wx.EXPAND|wx.ALL,0)
210        #self.GPXtree.Bind(wx.EVT_TREE_SEL_CHANGED,self.OnDataTreeSelChanged)
211        self.GPXtree.Bind(wx.EVT_TREE_SEL_CHANGED,self.OnDataTreeSelChanged)
212        # self.GPXtree.Bind(wx.EVT_TREE_ITEM_RIGHT_CLICK,self.OnDataTreeSelChanged)
213        # self.GPXtree.Bind(wx.EVT_TREE_ITEM_COLLAPSED,
214        #     self.OnGPXtreeItemCollapsed, id=TreeId)
215        #self.GPXtree.Bind(wx.EVT_TREE_ITEM_EXPANDED,
216        #     self.OnGPXtreeItemExpanded, id=TreeId)
217        # self.GPXtree.Bind(wx.EVT_TREE_DELETE_ITEM,
218        #     self.OnGPXtreeItemDelete, id=TreeId)
219        # self.GPXtree.Bind(wx.EVT_TREE_KEY_DOWN,
220        #     self.OnGPXtreeKeyDown, id=TreeId)
221        # self.GPXtree.Bind(wx.EVT_TREE_BEGIN_RDRAG,
222        #     self.OnGPXtreeBeginRDrag, id=TreeId)       
223        # self.GPXtree.Bind(wx.EVT_TREE_END_DRAG,
224        #     self.OnGPXtreeEndDrag, id=TreeId)       
225        self.root = self.GPXtree.root       
226        self.Bind(wx.EVT_CLOSE, lambda event: sys.exit())
227
228        self.fileList = []  # list of files read for use in Reload
229        self.histList = []  # list of histograms loaded for unique naming
230        self.projList = []
231
232        self.PWDRfilter = None
233        for win,var in ((self.plotFrame,'Plot_Pos'),):
234        #for win,var in ((self,'Main_Pos'),(self.plotFrame,'Plot_Pos')):
235            try:
236                pos = GSASIIpath.GetConfigValue(var)
237                if type(pos) is str: pos = eval(pos)
238                win.SetPosition(pos)
239                if G2gd.GetDisplay(pos) is None: win.Center()
240            except:
241                if GSASIIpath.GetConfigValue(var):
242                    print('Value for config {} {} is invalid'.format(var,GSASIIpath.GetConfigValue(var)))
243                    win.Center()
244        self.SetModeMenu()
245       
246    def SelectGPX(self):
247        '''Select a .GPX file to be read
248        '''
249        dlg = wx.FileDialog(self, 'Choose GSAS-II project file', 
250                wildcard='GSAS-II project file (*.gpx)|*.gpx',style=wx.FD_OPEN)
251        try:
252            if dlg.ShowModal() != wx.ID_OK: return
253            fil = os.path.splitext(dlg.GetPath())[0]+'.gpx'
254        finally:
255            dlg.Destroy()
256        if os.path.exists(fil):
257            #self.fileList.append([fil,'GPX'])
258            return fil
259        else:
260            print('File {} not found, skipping'.format(fil))
261            return
262
263    def SelectMultGPX(self):
264        '''Select multiple .GPX files to be read
265        '''
266        dlg = wx.FileDialog(self, 'Choose GSAS-II project file', 
267                wildcard='GSAS-II project file (*.gpx)|*.gpx',
268                style=wx.FD_OPEN|wx.FD_MULTIPLE)
269        try:
270            if dlg.ShowModal() != wx.ID_OK: return
271            files = dlg.GetPaths()
272        finally:
273            dlg.Destroy()
274        fileList = []
275        for f in files:
276            fil = os.path.splitext(f)[0]+'.gpx'
277            if os.path.exists(fil):
278                if fil not in fileList:
279                    fileList.append(fil)
280            else:
281                print('File {} not found, skipping'.format(fil))
282        return fileList
283       
284    def getMode(self):
285        '''returns the display mode (one of "Histogram","Phase","Project").
286        Could return '?' in case of an error.
287        '''
288        for m in self.wxID_Mode:
289            if self.Mode.FindItemById(self.wxID_Mode[m]).IsChecked():
290                break
291        else:
292            m = '?'
293        return m
294   
295    def onRefresh(self,event):
296        '''reread all files, in response to a change in mode, etc.
297        '''
298        self.GPXtree.DeleteChildren(self.root)  # delete tree contents
299        self.histList = []  # clear list of loaded histograms
300        self.projList = []
301        for fil,mode in self.fileList:
302            self.loadFile(fil)
303        self.doneLoad()
304        self.SetModeMenu()
305           
306    def SetModeMenu(self):
307        '''Create the mode-specific menu and its contents
308        '''
309        modeMenu = wx.Menu(title='')
310        oldmenu = self.MenuBar.Replace(self.modeMenuPos,modeMenu,self.getMode())
311        wx.CallAfter(oldmenu.Destroy)
312        if self.getMode() == "Histogram":
313            item = modeMenu.Append(wx.ID_ANY,'Prince test')
314            self.Bind(wx.EVT_MENU, self.onHistPrinceTest, id=item.GetId())
315        elif self.getMode() == "Phase":
316            pass
317        elif self.getMode() == "Project":
318            item = modeMenu.Append(wx.ID_ANY,'F-test')
319            self.Bind(wx.EVT_MENU, self.onProjFtest, id=item.GetId())
320
321    def loadFile(self,fil):
322        '''read or reread a file
323        '''
324        if self.getMode() == "Histogram":
325            self.LoadPwdr(fil)
326        elif self.getMode() == "Phase":
327            self.LoadPhase(fil)
328        elif self.getMode() == "Project":
329            self.LoadProject(fil)
330        else:
331            print("mode not implemented")
332            #raise Exception("mode not implemented")
333
334    def doneLoad(self):
335        self.GPXtree.Expand(self.root)
336        if self.getMode() == "Project":
337            overId = self.GPXtree.InsertItem(pos=0,parent=self.root,text='Project Overview')
338            self.GPXtree.SelectItem(overId)
339
340    def onLoadGPX(self,event):
341        '''Initial load of GPX file in response to a menu command
342        '''
343        fil = self.SelectGPX()
344        if not fil: return
345        if not os.path.exists(fil): return
346        self.fileList.append([fil,'GPX'])
347        self.loadFile(fil)
348        self.doneLoad()
349
350    def onLoadMultGPX(self,event):
351        '''Initial load of multiple GPX files in response to a menu command
352        '''
353        for fil in self.SelectMultGPX():
354            if not os.path.exists(fil): continue
355            self.fileList.append([fil,'GPX'])
356            self.loadFile(fil)
357        self.doneLoad()
358
359    def onLoadWildGPX(self,event,wildcard=None):
360        '''Initial load of GPX file in response to a menu command
361        '''
362        home = os.path.abspath(os.getcwd())
363        hlp = '''Enter a wildcard version of a file name.
364The directory is assumed to be "{}" unless specified otherwise.
365Extensions will be set to .gpx and .bak files will be ignored unless
366the wildcard string includes "bak". For example, "*abc*" will match any
367.gpx file in that directory containing "abc". String "/tmp/A*" will match
368files in "/tmp" beginning with "A". Supplying two strings, "A*" and "B*bak*"
369will match files names beginning with "A" or "B", but ".bak" files will
370be included for the files beginning with "B" only.
371'''.format(home)
372        if wildcard is None:
373            dlg = G2G.MultiStringDialog(self, 'Enter wildcard file names', 
374                    ['wild-card 1'] , values=['*'], 
375                    lbl='Provide string(s) with "*" to find matching files',
376                    addRows=True, hlp=hlp)
377            res = dlg.Show()
378            wl = dlg.GetValues()
379            dlg.Destroy()
380            if not res: return
381        else:
382            wl = [wildcard]
383        for w in wl:
384            if not os.path.split(w)[0]:
385                w = os.path.join(home,w)
386            w = os.path.splitext(w)[0] + '.gpx'
387            for fil in glob.glob(w): 
388                if not os.path.exists(fil): continue
389                if '.bak' in fil and 'bak' not in w: continue
390                if fil in [i for i,j in self.fileList]: continue
391                self.fileList.append([fil,'GPX'])
392                self.loadFile(fil)
393        self.doneLoad()
394
395    def LoadPwdr(self,fil):
396        '''Load PWDR entries from a .GPX file to the tree.
397        see :func:`GSASIIIO.ProjFileOpen`
398        '''
399        G2frame = self
400        filep = open(fil,'rb')
401        shortname = os.path.splitext(os.path.split(fil)[1])[0]
402
403        wx.BeginBusyCursor()
404        histLoadList = []
405        try:
406            while True:
407                try:
408                    data = cPickleLoad(filep)
409                except EOFError:
410                    break
411                if not data[0][0].startswith('PWDR'): continue
412                if self.PWDRfilter is not None: # implement filter
413                    if self.PWDRfilter not in data[0][0]: continue
414                data[0][0] += ' ('
415                data[0][0] += shortname
416                data[0][0] += ')'
417                histLoadList.append(data)
418                       
419        except Exception as errmsg:
420            if GSASIIpath.GetConfigValue('debug'):
421                print('\nError reading GPX file:',errmsg)
422                import traceback
423                print (traceback.format_exc())
424            msg = wx.MessageDialog(G2frame,message="Error reading file "+
425                str(fil)+". This is not a current GSAS-II .gpx file",
426                caption="Load Error",style=wx.ICON_ERROR | wx.OK | wx.STAY_ON_TOP)
427            msg.ShowModal()
428        finally:
429            filep.close()
430            wx.EndBusyCursor()
431
432        datum = None
433        for i,data in enumerate(histLoadList):
434            datum = data[0]
435            datum[0] = G2obj.MakeUniqueLabel(datum[0],self.histList)
436            Id = G2frame.GPXtree.AppendItem(parent=G2frame.root,text=datum[0])
437            self.histList.append(datum[0])
438            # if 'ranId' not in datum[1][0]: # patch: add random Id if not present
439            #     datum[1][0]['ranId'] = ran.randint(0,sys.maxsize)
440            G2frame.GPXtree.SetItemPyData(Id,datum[1][:3])  #temp. trim off junk (patch?)
441            for datus in data[1:]:
442                sub = G2frame.GPXtree.AppendItem(Id,datus[0])
443    #patch
444                if datus[0] == 'Instrument Parameters' and len(datus[1]) == 1:
445                    if datum[0].startswith('PWDR'):
446                        datus[1] = [dict(zip(datus[1][3],zip(datus[1][0],datus[1][1],datus[1][2]))),{}]
447                    else:
448                        datus[1] = [dict(zip(datus[1][2],zip(datus[1][0],datus[1][1]))),{}]
449                    for item in datus[1][0]:               #zip makes tuples - now make lists!
450                        datus[1][0][item] = list(datus[1][0][item])
451    #end patch
452                G2frame.GPXtree.SetItemPyData(sub,datus[1])
453        if datum: # was anything loaded?
454            print('data load successful for {}'.format(datum[0]))
455    #        G2frame.Status.SetStatusText('Mouse RB drag/drop to reorder',0)
456    #    G2frame.SetTitleByGPX()
457        self.GPXtree.Expand(self.root)
458       
459    def onHistFilter(self,event):
460        'Load a filter string via a dialog in response to a menu event'
461        lbl = ''
462        if self.PWDRfilter is not None:
463            lbl = self.PWDRfilter
464        dlg = G2G.SingleStringDialog(self,'Set string',
465                                'Set a string that must be in histogram name',
466                                 lbl,size=(400,-1))
467        if dlg.Show():
468            if dlg.GetValue().strip() == '':
469                self.PWDRfilter = None
470            else:
471                self.PWDRfilter = dlg.GetValue()
472            dlg.Destroy()
473            self.onRefresh(event)
474        else:
475            dlg.Destroy()
476
477    def LoadPhase(self,fil):
478        '''Load Phase entries from a .GPX file to the tree.
479        see :func:`GSASIIIO.ProjFileOpen`
480        '''
481        G2frame = self
482        filep = open(fil,'rb')
483        shortname = os.path.splitext(os.path.split(fil)[1])[0]
484
485        wx.BeginBusyCursor()
486        Phases = None
487        try:
488            while True:
489                try:
490                    data = cPickleLoad(filep)
491                except EOFError:
492                    break
493                if not data[0][0].startswith('Phase'): continue
494                Phases = data
495                #if self.PWDRfilter is not None: # implement filter
496                #    if self.PWDRfilter not in data[0][0]: continue
497                data[0][0] += ' ('
498                if Phases:
499                    data[0][0] += shortname
500                    data[0][0] += ')'
501                else: 
502                    data[0][0] += shortname
503                    data[0][0] += 'has no phases)'
504                Phases = data
505                break
506               
507        except Exception as errmsg:
508            if GSASIIpath.GetConfigValue('debug'):
509                print('\nError reading GPX file:',errmsg)
510                import traceback
511                print (traceback.format_exc())
512            msg = wx.MessageDialog(G2frame,message="Error reading file "+
513                str(fil)+". This is not a current GSAS-II .gpx file",
514                caption="Load Error",style=wx.ICON_ERROR | wx.OK | wx.STAY_ON_TOP)
515            msg.ShowModal()
516        finally:
517            filep.close()
518            wx.EndBusyCursor()
519
520        datum = None
521        if Phases:
522            datum = data[0]
523            #datum[0] = G2obj.MakeUniqueLabel(datum[0],self.histList)
524            Id = G2frame.GPXtree.AppendItem(parent=G2frame.root,text=datum[0])
525            G2frame.GPXtree.SetItemPyData(Id,datum[1])
526            for datus in data[1:]:
527                #datus[0] += ' ('
528                #datus[0] += shortname
529                #datus[0] += ')'
530                sub = G2frame.GPXtree.AppendItem(Id,datus[0])
531                G2frame.GPXtree.SetItemPyData(sub,datus[1])
532        if datum: # was anything loaded?
533            self.GPXtree.Expand(Id)
534            print('Phase load successful for {}'.format(datum[0]))
535    #        G2frame.Status.SetStatusText('Mouse RB drag/drop to reorder',0)
536    #    G2frame.SetTitleByGPX()
537        self.GPXtree.Expand(self.root)
538
539    def LoadProject(self,fil):
540        '''Load the Covariance entry from a .GPX file to the tree.
541        see :func:`GSASIIIO.ProjFileOpen`
542        '''
543        import datetime
544        G2frame = self
545        filep = open(fil,'rb')
546        saved = datetime.datetime.fromtimestamp(os.path.getmtime(fil)).strftime("%Y-%h-%d %H:%M")
547        shortname = os.path.splitext(os.path.split(fil)[1])[0]
548        projInfo = [shortname,saved]
549        wx.BeginBusyCursor()
550        #Phases = None
551        #G2frame.GPXtree.SetItemPyData(Id,Covar[1])
552        try:
553            while True:
554                try:
555                    data = cPickleLoad(filep)
556                except EOFError:
557                    break
558                #if data[0][0] == 'Controls' and 'LastSavedUsing' in data[0][1]:
559                if not data[0][0].startswith('Covariance'): continue
560                Covar = data[0]
561                f = '{:d}'
562                if 'varyList' in data[0][1]:
563                    projInfo += [f.format(len(data[0][1]['varyList']))]
564                else:
565                    projInfo += ['?']
566                for v in 'Nobs','GOF':
567                    if 'Rvals' in data[0][1] and v in data[0][1]['Rvals']:
568                        projInfo += [f.format(data[0][1]['Rvals'][v])]
569                    else:
570                        projInfo += ['?']
571                    f = '{:6.2f}'
572                #GSASIIpath.IPyBreak_base()
573                #if self.PWDRfilter is not None: # implement filter
574                #    if self.PWDRfilter not in data[0][0]: continue
575                Covar[0] = shortname # + ' Covariance'
576                Id = G2frame.GPXtree.AppendItem(parent=G2frame.root,text=Covar[0])
577                G2frame.GPXtree.SetItemPyData(Id,Covar[1])
578                self.projList.append(projInfo)
579                break
580            else:
581                print("{} does not have refinement results".format(shortname))
582        except Exception as errmsg:
583            if GSASIIpath.GetConfigValue('debug'):
584                print('\nError reading GPX file:',errmsg)
585                import traceback
586                print (traceback.format_exc())
587            msg = wx.MessageDialog(G2frame,message="Error reading file "+
588                str(fil)+". This is not a current GSAS-II .gpx file",
589                caption="Load Error",style=wx.ICON_ERROR | wx.OK | wx.STAY_ON_TOP)
590            msg.ShowModal()
591        finally:
592            filep.close()
593            wx.EndBusyCursor()
594
595    def OnDataTreeSelChanged(self,event):
596        def ClearData(self):
597            '''Initializes the contents of the dataWindow panel
598            '''
599            self.Unbind(wx.EVT_SIZE)
600            self.SetBackgroundColour(wx.Colour(240,240,240))
601            Sizer = self.GetSizer()
602            if Sizer:
603                try:
604                    Sizer.Clear(True)
605                except:
606                    pass
607        G2frame = self
608        item = event.GetItem()
609        #print('selected',item)
610        lbl = G2frame.GPXtree.GetItemText(item)
611        if self.getMode() == "Project" and lbl == 'Project Overview':
612            ClearData(G2frame.dataWindow)
613            #import imp
614            #imp.reload(G2G)
615            pnl = G2G.SortableLstCtrl(G2frame.dataWindow)
616            h = ["File", "last saved", "vars", "Nobs", "GOF"]
617            j = [ 0,       0,           1,      1,      1]
618            pnl.PopulateHeader(h,j)
619            for i,line in enumerate(self.projList):
620                pnl.PopulateLine(i,line)
621            G2frame.dataWindow.GetSizer().Add(pnl,1,wx.EXPAND)
622            pnl.SetColWidth(0,maxwidth=170)
623            for i in range(1,len(h)):
624                pnl.SetColWidth(i,minwidth=50)
625            G2frame.dataWindow.SendSizeEvent()
626        elif self.getMode() == "Project":
627            ClearData(G2frame.dataWindow)
628            data = G2frame.GPXtree.GetItemPyData(item)
629            if data is None:
630                self.plotFrame.Show(False)
631                return
632            text = ''
633            if 'Rvals' in data:
634                Nvars = len(data['varyList'])
635                Rvals = data['Rvals']
636                text = ('Residuals after last refinement:\n'+
637                        '\twR = {:.3f}\n\tchi**2 = {:.1f}\n\tGOF = {:.2f}').format(
638                        Rvals['Rwp'],Rvals['chisq'],Rvals['GOF'])
639                text += '\n\tNobs = {}\n\tNvals = {}\n\tSVD zeros = {}'.format(
640                    Rvals['Nobs'],Nvars,Rvals.get('SVD0',0.))
641                text += '\n\tmax shift/esd = {:.3f}'.format(Rvals.get('Max shft/sig',0.0))
642                if 'lamMax' in Rvals:
643                    text += '\n\tlog10 MaxLambda = {:.1f}'.format(np.log10(Rvals['lamMax']))
644                text += '\n\tReduced χ**2 = {:.2f}'.format(Rvals['GOF']**2)
645                G2frame.dataWindow.GetSizer().Add(
646                    wx.StaticText(G2frame.dataWindow,wx.ID_ANY,text)
647                )
648                self.plotFrame.Show(True)
649                G2plt.PlotCovariance(G2frame,data)
650        else:
651            self.plotFrame.Show(False)
652            ClearData(G2frame.dataWindow)
653       
654        #print(self.GPXtree._getTreeItemsList(item))
655        # pltNum = self.G2plotNB.nb.GetSelection()
656        # print(pltNum)
657        # if pltNum >= 0:                         #to avoid the startup with no plot!
658        #     self.G2plotNB.nb.GetPage(pltNum)
659        #     NewPlot = False
660        # else:
661        #     NewPlot = True
662        #if self.getMode() == "Histogram":
663        #self.PatternId = self.PickId  = item
664        #G2plt.PlotPatterns(self,plotType='PWDR',newPlot=NewPlot)
665           
666    # def OnGPXtreeItemExpanded(self,event):
667    #     item = event.GetItem()
668    #     print('expanded',item)
669    #     print(self.GPXtree._getTreeItemsList(item))
670    #     if item == self.root:
671    #         event.StopPropagation()
672    #     else:
673    #         event.Skip(False)
674   
675    def onProjFtest(self,event):
676        '''Compare two projects (selected here if more than two are present)
677        using the statistical F-test (aka Hamilton R-factor test), see:
678
679            * Hamilton, R. W. (1965), Acta Crystallogr. 18, 502-510.
680            * Prince, E., Mathematical Techniques in Crystallography and Materials Science, Second ed. (Springer-Verlag, New York, 1994).
681        '''
682        items = []
683        item, cookie = self.GPXtree.GetFirstChild(self.root)
684        while item:
685            items.append(item)
686            item, cookie = self.GPXtree.GetNextChild(self.root, cookie)
687        if len(items) < 2:
688            G2G.G2MessageBox(self,'F-test requires two projects','Need more projects')
689            return
690        elif len(items) == 2:
691            s0 = items[0]
692            baseDict = self.GPXtree.GetItemPyData(s0)
693            s1 = items[1]           
694            relxDict = self.GPXtree.GetItemPyData(s1)
695            # sort out the dicts in order of number of refined variables
696            if len(baseDict['varyList']) > len(relxDict['varyList']):
697                s0,s1,baseDict,relxDict = s1,s0,relxDict,baseDict
698        else:
699            # need to make selection here
700            sel = []
701            for i in items:
702                sel.append(self.GPXtree.GetItemText(i))
703            dlg = G2G.G2SingleChoiceDialog(self,'Select constrained refinement',
704                                             'Choose refinement',sel)
705            if dlg.ShowModal() == wx.ID_OK:
706                s0 = dlg.GetSelection()
707                dlg.Destroy()
708            else:
709                dlg.Destroy()
710                return
711            inds = list(range(len(items)))
712            del sel[s0]
713            del inds[s0]
714            dlg = G2G.G2SingleChoiceDialog(self,'Select relaxed refinement',
715                                             'Choose refinement',sel)
716            if dlg.ShowModal() == wx.ID_OK:
717                s1 = dlg.GetSelection()
718                s1 = inds[s1]
719                dlg.Destroy()
720            else:
721                dlg.Destroy()
722                return
723            baseDict = self.GPXtree.GetItemPyData(items[s0])
724            relxDict = self.GPXtree.GetItemPyData(items[s1])
725            if len(baseDict['varyList']) > len(relxDict['varyList']):
726                G2G.G2MessageBox(self,
727                            'F-test warning: constrained refinement has more '+
728                            'variables ({}) than relaxed refinement ({}). Swapping'
729                            .format(len(baseDict['varyList']), len(relxDict['varyList'])),
730                            'Fits reversed?')
731                s0,s1,baseDict,relxDict = s1,s0,relxDict,baseDict
732                baseDict,relxDict = relxDict,baseDict
733        if len(baseDict['varyList']) == len(relxDict['varyList']):
734            G2G.G2MessageBox(self,'F-test requires differing numbers of variables','F-test not valid')
735            return
736        elif baseDict['Rvals']['Nobs'] != relxDict['Rvals']['Nobs']:
737            G2G.G2MessageBox(self,'F-test requires same number of observations in each refinement','F-test not valid')
738            return
739        missingVars = []
740        for var in baseDict['varyList']: 
741            if var not in relxDict['varyList']: 
742                missingVars.append(var)
743        txt = ''
744        postmsg = ''
745        if missingVars:
746            txt = ('*** Possible invalid use of F-test: '+
747                'The F-test requires that the constrained model be a subset '+
748                'of the relaxed model. Review the parameters shown below to '+
749                'confirm missing var(s) have new names in the relaxed model. '+
750                '***\n\n')
751            postmsg = ('\n\nThese parameters are in the constrained '+
752              'fit and are not in the relaxed fit:\n*  ')
753            for i,var in enumerate(missingVars):
754                if i > 0: postmsg += ', '
755                postmsg += var
756            postmsg += ('\nThese parameters are in the relaxed fit and not'+
757                ' in the constrained fit:\n*  ')
758            i = 0
759            for var in relxDict['varyList']: 
760                if var not in baseDict['varyList']: 
761                    if i > 0: postmsg += ', '
762                    i += 1
763                    postmsg += var
764        #GSASIIpath.IPyBreak_base()
765        prob = RwFtest(baseDict['Rvals']['Nobs'],
766                   baseDict['Rvals']['GOF'],len(baseDict['varyList']),
767                   relxDict['Rvals']['GOF'],len(relxDict['varyList']))
768        fmt = "{} model is \n{}\n{} variables and Reduced Chi**2 = {:.3f}"
769        msg = txt
770        msg += fmt.format('Constrained',self.GPXtree.GetItemText(s0)[:-11],
771                       len(baseDict['varyList']),
772                       baseDict['Rvals']['GOF']**2)
773        msg += '\n\n'
774        msg += fmt.format('Relaxed',self.GPXtree.GetItemText(s1)[:-11],
775                       len(relxDict['varyList']),
776                       relxDict['Rvals']['GOF']**2)
777        msg += '\n\nCumulative F-test probability {:.2f}%\n'.format(prob*100)
778        if prob > 0.95:
779            msg += "The relaxed model is statistically preferred to the constrained model."
780        elif prob > 0.75:
781            msg += "There is little ability to differentiate between the two models on a statistical basis."
782        else:
783            msg += "The constrained model is statistically preferred to the relaxed model."
784        msg += postmsg
785        G2G.G2MessageBox(self,msg,'F-test result')
786       
787    def onHistPrinceTest(self,event):
788        '''Compare two histograms (selected here if more than two are present)
789        using the statistical test proposed by Ted Prince in
790        Acta Cryst. B35 1099-1100. (1982). Also see Int. Tables Vol. C
791        (1st Ed.) chapter 8.4, 618-621 (1995).
792        '''
793        items = []
794        item, cookie = self.GPXtree.GetFirstChild(self.root)
795        while item:
796            items.append(item)
797            item, cookie = self.GPXtree.GetNextChild(self.root, cookie)
798
799        if len(items) < 2:
800            G2G.G2MessageBox(self,'Prince test requires two histograms','Need more')
801            return
802        elif len(items) == 2:
803            s0,s1 = items
804        else:
805            # need to make selection here
806            sel = []
807            for i in items:
808                sel.append(self.GPXtree.GetItemText(i))
809            dlg = G2G.G2SingleChoiceDialog(self,'Select one refinement',
810                                             'Choose refinement',sel)
811            if dlg.ShowModal() == wx.ID_OK:
812                s0 = dlg.GetSelection()
813                dlg.Destroy()
814            else:
815                dlg.Destroy()
816                return
817            inds = list(range(len(items)))
818            del sel[s0]
819            del inds[s0]
820            dlg = G2G.G2SingleChoiceDialog(self,'Select comparison refinement',
821                                             'Choose refinement',sel)
822            if dlg.ShowModal() == wx.ID_OK:
823                s1 = dlg.GetSelection()
824                s1 = inds[s1]
825                dlg.Destroy()
826            else:
827                dlg.Destroy()
828                return
829        model0 = self.GPXtree.GetItemPyData(s0)
830        data0 = model0[1]
831        model1 = self.GPXtree.GetItemPyData(s1)
832        data1 = model1[1]
833        if len(data0[0]) != len(data1[0]):
834            G2G.G2MessageBox(self,'Unable to test: differing numbers of data points','Comparison not valid')
835            return
836        if max(abs((data0[0]-data1[0])/data0[0])) > 0.01:
837            G2G.G2MessageBox(self,'Unable to use test: "X" values differ','Comparison not valid')
838            return
839        # X = data0[3] - data1[3]
840        # #Z = np.sqrt(data0[3]) * (data0[1] - (data0[3] + data1[3])/2)
841        # X = (data0[3] - data1[3]) / np.sqrt(data0[1])
842        # Z = (data0[1] - (data0[3] + data1[3])/2) / np.sqrt(data0[1])
843        # lam = np.sum(X*Z) / np.sum(X)
844        # sig = np.sqrt(
845        #     (np.sum(Z*Z) - lam*lam*np.sum(X*X)) /
846        #     ((len(data0[0]) - 1) * np.sum(X*X))
847        #     )
848           
849#    0 the x-postions (two-theta in degrees),
850#    1 the intensity values (Yobs),
851#    2 the weights for each Yobs value
852#    3 the computed intensity values (Ycalc)
853#    4 the background values
854#    5 Yobs-Ycalc
855       
856        GSASIIpath.IPyBreak_base()
857               
858#======================================================================   
859if __name__ == '__main__':
860    #if sys.platform == "darwin":
861    #    application = G2App(0) # create the GUI framework
862    #else:
863    application = wx.App(0) # create the GUI framework
864    try:
865        GSASIIpath.SetBinaryPath(True)
866    except:
867        print('Unable to run with current setup, do you want to update to the')
868        try:
869#            if '2' in platform.python_version_tuple()[0]:           
870#                ans = raw_input("latest GSAS-II version? Update ([Yes]/no): ")
871#            else:
872                ans = input("latest GSAS-II version? Update ([Yes]/no): ")               
873        except:
874            ans = 'no'
875        if ans.strip().lower() == "no":
876            print('Exiting')
877            sys.exit()
878        print('Updating...')
879        GSASIIpath.svnUpdateProcess()
880    GSASIIpath.InvokeDebugOpts()
881    Frame = main(application) # start the GUI
882    argLoadlist = sys.argv[1:]
883    for arg in argLoadlist:
884        fil = os.path.splitext(arg)[0] + '.gpx'
885        if os.path.exists(fil):
886            Frame.fileList.append([fil,'GPX'])
887            Frame.loadFile(fil)
888        else:
889            print('File {} not found. Skipping'.format(fil))
890    Frame.doneLoad()
891    # debug code to select Project in initial mode
892    #Frame.onLoadWildGPX(None,wildcard='*')
893    #Frame.Mode.FindItemById(Frame.wxID_Mode['Project']).Check(True)
894    #Frame.onRefresh(None)
895   
896    application.MainLoop()
Note: See TracBrowser for help on using the repository browser.