source: trunk/G2compare.py @ 4339

Last change on this file since 4339 was 4339, checked in by toby, 20 months ago

set svn flags; improve compare F-test & add info & covariance plot

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