source: trunk/G2compare.py @ 4345

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

Add more image settings to preferences; prevent save of config.py after "no save"

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 31.3 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-06 21:05:22 +0000 (Fri, 06 Mar 2020) $
6# $Author: toby $
7# $Revision: 4345 $
8# $URL: trunk/G2compare.py $
9# $Id: G2compare.py 4345 2020-03-06 21:05:22Z 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: 4345 $")
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            item = modeMenu.Append(wx.ID_ANY,'Prince test')
281            self.Bind(wx.EVT_MENU, self.onHistPrinceTest, id=item.GetId())
282        elif self.getMode() == "Phase":
283            pass
284        elif self.getMode() == "Project":
285            item = modeMenu.Append(wx.ID_ANY,'F-test')
286            self.Bind(wx.EVT_MENU, self.onProjFtest, id=item.GetId())
287
288    def loadFile(self,fil):
289        '''read or reread a file
290        '''
291        if self.getMode() == "Histogram":
292            self.LoadPwdr(fil)
293        elif self.getMode() == "Phase":
294            self.LoadPhase(fil)
295        elif self.getMode() == "Project":
296            self.LoadProject(fil)
297        else:
298            print("mode not implemented")
299            #raise Exception("mode not implemented")
300       
301    def onLoadGPX(self,event):
302        '''Initial load of GPX file in response to a menu command
303        '''
304        fil = self.SelectGPX()
305        if not fil: return
306        if not os.path.exists(fil): return
307        self.fileList.append([fil,'GPX'])
308        self.loadFile(fil)
309
310    def LoadPwdr(self,fil):
311        '''Load PWDR entries from a .GPX file to the tree.
312        see :func:`GSASIIIO.ProjFileOpen`
313        '''
314        G2frame = self
315        filep = open(fil,'rb')
316        shortname = os.path.splitext(os.path.split(fil)[1])[0]
317
318        wx.BeginBusyCursor()
319        histLoadList = []
320        try:
321            while True:
322                try:
323                    data = cPickleLoad(filep)
324                except EOFError:
325                    break
326                if not data[0][0].startswith('PWDR'): continue
327                if self.PWDRfilter is not None: # implement filter
328                    if self.PWDRfilter not in data[0][0]: continue
329                data[0][0] += ' ('
330                data[0][0] += shortname
331                data[0][0] += ')'
332                histLoadList.append(data)
333                       
334        except Exception as errmsg:
335            if GSASIIpath.GetConfigValue('debug'):
336                print('\nError reading GPX file:',errmsg)
337                import traceback
338                print (traceback.format_exc())
339            msg = wx.MessageDialog(G2frame,message="Error reading file "+
340                str(fil)+". This is not a current GSAS-II .gpx file",
341                caption="Load Error",style=wx.ICON_ERROR | wx.OK | wx.STAY_ON_TOP)
342            msg.ShowModal()
343        finally:
344            filep.close()
345            wx.EndBusyCursor()
346
347        datum = None
348        for i,data in enumerate(histLoadList):
349            datum = data[0]
350            datum[0] = G2obj.MakeUniqueLabel(datum[0],self.histList)
351            Id = G2frame.GPXtree.AppendItem(parent=G2frame.root,text=datum[0])
352            self.histList.append(datum[0])
353            # if 'ranId' not in datum[1][0]: # patch: add random Id if not present
354            #     datum[1][0]['ranId'] = ran.randint(0,sys.maxsize)
355            G2frame.GPXtree.SetItemPyData(Id,datum[1][:3])  #temp. trim off junk (patch?)
356            for datus in data[1:]:
357                sub = G2frame.GPXtree.AppendItem(Id,datus[0])
358    #patch
359                if datus[0] == 'Instrument Parameters' and len(datus[1]) == 1:
360                    if datum[0].startswith('PWDR'):
361                        datus[1] = [dict(zip(datus[1][3],zip(datus[1][0],datus[1][1],datus[1][2]))),{}]
362                    else:
363                        datus[1] = [dict(zip(datus[1][2],zip(datus[1][0],datus[1][1]))),{}]
364                    for item in datus[1][0]:               #zip makes tuples - now make lists!
365                        datus[1][0][item] = list(datus[1][0][item])
366    #end patch
367                G2frame.GPXtree.SetItemPyData(sub,datus[1])
368        if datum: # was anything loaded?
369            print('project load successful for {}'.format(datum[0]))
370    #        G2frame.Status.SetStatusText('Mouse RB drag/drop to reorder',0)
371    #    G2frame.SetTitleByGPX()
372        self.GPXtree.Expand(self.root)
373       
374    def onHistFilter(self,event):
375        'Load a filter string via a dialog in response to a menu event'
376        lbl = ''
377        if self.PWDRfilter is not None:
378            lbl = self.PWDRfilter
379        dlg = G2G.SingleStringDialog(self,'Set string',
380                                'Set a string that must be in histogram name',
381                                 lbl,size=(400,-1))
382        if dlg.Show():
383            if dlg.GetValue().strip() == '':
384                self.PWDRfilter = None
385            else:
386                self.PWDRfilter = dlg.GetValue()
387            dlg.Destroy()
388            self.onRefresh(event)
389        else:
390            dlg.Destroy()
391
392    def LoadPhase(self,fil):
393        '''Load Phase entries from a .GPX file to the tree.
394        see :func:`GSASIIIO.ProjFileOpen`
395        '''
396        G2frame = self
397        filep = open(fil,'rb')
398        shortname = os.path.splitext(os.path.split(fil)[1])[0]
399
400        wx.BeginBusyCursor()
401        Phases = None
402        try:
403            while True:
404                try:
405                    data = cPickleLoad(filep)
406                except EOFError:
407                    break
408                if not data[0][0].startswith('Phase'): continue
409                Phases = data
410                #if self.PWDRfilter is not None: # implement filter
411                #    if self.PWDRfilter not in data[0][0]: continue
412                data[0][0] += ' ('
413                if Phases:
414                    data[0][0] += shortname
415                    data[0][0] += ')'
416                else: 
417                    data[0][0] += shortname
418                    data[0][0] += 'has no phases)'
419                Phases = data
420                break
421               
422        except Exception as errmsg:
423            if GSASIIpath.GetConfigValue('debug'):
424                print('\nError reading GPX file:',errmsg)
425                import traceback
426                print (traceback.format_exc())
427            msg = wx.MessageDialog(G2frame,message="Error reading file "+
428                str(fil)+". This is not a current GSAS-II .gpx file",
429                caption="Load Error",style=wx.ICON_ERROR | wx.OK | wx.STAY_ON_TOP)
430            msg.ShowModal()
431        finally:
432            filep.close()
433            wx.EndBusyCursor()
434
435        datum = None
436        if Phases:
437            datum = data[0]
438            #datum[0] = G2obj.MakeUniqueLabel(datum[0],self.histList)
439            Id = G2frame.GPXtree.AppendItem(parent=G2frame.root,text=datum[0])
440            G2frame.GPXtree.SetItemPyData(Id,datum[1])
441            for datus in data[1:]:
442                #datus[0] += ' ('
443                #datus[0] += shortname
444                #datus[0] += ')'
445                sub = G2frame.GPXtree.AppendItem(Id,datus[0])
446                G2frame.GPXtree.SetItemPyData(sub,datus[1])
447        if datum: # was anything loaded?
448            self.GPXtree.Expand(Id)
449            print('project load successful for {}'.format(datum[0]))
450    #        G2frame.Status.SetStatusText('Mouse RB drag/drop to reorder',0)
451    #    G2frame.SetTitleByGPX()
452        self.GPXtree.Expand(self.root)
453
454    def LoadProject(self,fil):
455        '''Load the Covariance entry from a .GPX file to the tree.
456        see :func:`GSASIIIO.ProjFileOpen`
457        '''
458        G2frame = self
459        filep = open(fil,'rb')
460        shortname = os.path.splitext(os.path.split(fil)[1])[0]
461
462        wx.BeginBusyCursor()
463        Phases = None
464        try:
465            while True:
466                try:
467                    data = cPickleLoad(filep)
468                except EOFError:
469                    break
470                if not data[0][0].startswith('Covariance'): continue
471                Covar = data[0]
472                #GSASIIpath.IPyBreak_base()
473                #if self.PWDRfilter is not None: # implement filter
474                #    if self.PWDRfilter not in data[0][0]: continue
475                Covar[0] = shortname + ' Covariance'
476                Id = G2frame.GPXtree.AppendItem(parent=G2frame.root,text=Covar[0])
477                G2frame.GPXtree.SetItemPyData(Id,Covar[1])
478                break
479            else:
480                print("{} does not have refinement results".format(shortname))
481        except Exception as errmsg:
482            if GSASIIpath.GetConfigValue('debug'):
483                print('\nError reading GPX file:',errmsg)
484                import traceback
485                print (traceback.format_exc())
486            msg = wx.MessageDialog(G2frame,message="Error reading file "+
487                str(fil)+". This is not a current GSAS-II .gpx file",
488                caption="Load Error",style=wx.ICON_ERROR | wx.OK | wx.STAY_ON_TOP)
489            msg.ShowModal()
490        finally:
491            filep.close()
492            wx.EndBusyCursor()
493        self.GPXtree.Expand(self.root)
494
495    def OnDataTreeSelChanged(self,event):
496        def ClearData(self):
497            '''Initializes the contents of the dataWindow panel
498            '''
499            self.Unbind(wx.EVT_SIZE)
500            self.SetBackgroundColour(wx.Colour(240,240,240))
501            Sizer = self.GetSizer()
502            if Sizer:
503                try:
504                    Sizer.Clear(True)
505                except:
506                    pass
507        G2frame = self
508        item = event.GetItem()
509        #print('selected',item)
510        if self.getMode() == "Project":
511            ClearData(G2frame.dataWindow)
512            data = G2frame.GPXtree.GetItemPyData(item)
513            if data is None:
514                self.plotFrame.Show(False)
515                return
516            text = ''
517            if 'Rvals' in data:
518                Nvars = len(data['varyList'])
519                Rvals = data['Rvals']
520                text = ('Residuals after last refinement:\n'+
521                        '\twR = {:.3f}\n\tchi**2 = {:.1f}\n\tGOF = {:.2f}').format(
522                        Rvals['Rwp'],Rvals['chisq'],Rvals['GOF'])
523                text += '\n\tNobs = {}\n\tNvals = {}\n\tSVD zeros = {}'.format(
524                    Rvals['Nobs'],Nvars,Rvals.get('SVD0',0.))
525                text += '\n\tmax shift/esd = {:.3f}'.format(Rvals.get('Max shft/sig',0.0))
526                if 'lamMax' in Rvals:
527                    text += '\n\tlog10 MaxLambda = {:.1f}'.format(np.log10(Rvals['lamMax']))
528                G2frame.dataWindow.GetSizer().Add(
529                    wx.StaticText(G2frame.dataWindow,wx.ID_ANY,text)
530                )
531                self.plotFrame.Show(True)
532                G2plt.PlotCovariance(G2frame,data)
533        else:
534            self.plotFrame.Show(False)
535            ClearData(G2frame.dataWindow)
536       
537        #print(self.GPXtree._getTreeItemsList(item))
538        # pltNum = self.G2plotNB.nb.GetSelection()
539        # print(pltNum)
540        # if pltNum >= 0:                         #to avoid the startup with no plot!
541        #     self.G2plotNB.nb.GetPage(pltNum)
542        #     NewPlot = False
543        # else:
544        #     NewPlot = True
545        #if self.getMode() == "Histogram":
546        #self.PatternId = self.PickId  = item
547        #G2plt.PlotPatterns(self,plotType='PWDR',newPlot=NewPlot)
548           
549    # def OnGPXtreeItemExpanded(self,event):
550    #     item = event.GetItem()
551    #     print('expanded',item)
552    #     print(self.GPXtree._getTreeItemsList(item))
553    #     if item == self.root:
554    #         event.StopPropagation()
555    #     else:
556    #         event.Skip(False)
557   
558    def onProjFtest(self,event):
559        '''Compare two projects (selected here if more than two are present)
560        using the statistical F-test (aka Hamilton R-factor test), see:
561
562            * Hamilton, R. W. (1965), Acta Crystallogr. 18, 502-510.
563            * Prince, E., Mathematical Techniques in Crystallography and Materials Science, Second ed. (Springer-Verlag, New York, 1994).
564        '''
565        items = []
566        item, cookie = self.GPXtree.GetFirstChild(self.root)
567        while item:
568            items.append(item)
569            item, cookie = self.GPXtree.GetNextChild(self.root, cookie)
570        if len(items) < 2:
571            G2G.G2MessageBox(self,'F-test requires two projects','Need more projects')
572            return
573        elif len(items) == 2:
574            s0 = items[0]
575            baseDict = self.GPXtree.GetItemPyData(s0)
576            s1 = items[1]           
577            relxDict = self.GPXtree.GetItemPyData(s1)
578            # sort out the dicts in order of number of refined variables
579            if len(baseDict['varyList']) > len(relxDict['varyList']):
580                s0,s1,baseDict,relxDict = s1,s0,relxDict,baseDict
581        else:
582            # need to make selection here
583            sel = []
584            for i in items:
585                sel.append(self.GPXtree.GetItemText(i))
586            dlg = G2G.G2SingleChoiceDialog(self,'Select constrained refinement',
587                                             'Choose refinement',sel)
588            if dlg.ShowModal() == wx.ID_OK:
589                s0 = dlg.GetSelection()
590                dlg.Destroy()
591            else:
592                dlg.Destroy()
593                return
594            inds = list(range(len(items)))
595            del sel[s0]
596            del inds[s0]
597            dlg = G2G.G2SingleChoiceDialog(self,'Select relaxed refinement',
598                                             'Choose refinement',sel)
599            if dlg.ShowModal() == wx.ID_OK:
600                s1 = dlg.GetSelection()
601                s1 = inds[s1]
602                dlg.Destroy()
603            else:
604                dlg.Destroy()
605                return
606            baseDict = self.GPXtree.GetItemPyData(items[s0])
607            relxDict = self.GPXtree.GetItemPyData(items[s1])
608            if len(baseDict['varyList']) > len(relxDict['varyList']):
609                G2G.G2MessageBox(self,
610                            'F-test warning: constrained refinement has more '+
611                            'variables ({}) than relaxed refinement ({}). Swapping'
612                            .format(len(baseDict['varyList']), len(relxDict['varyList'])),
613                            'Fits reversed?')
614                s0,s1,baseDict,relxDict = s1,s0,relxDict,baseDict
615                baseDict,relxDict = relxDict,baseDict
616        if len(baseDict['varyList']) == len(relxDict['varyList']):
617            G2G.G2MessageBox(self,'F-test requires differing numbers of variables','F-test not valid')
618            return
619        elif baseDict['Rvals']['Nobs'] != relxDict['Rvals']['Nobs']:
620            G2G.G2MessageBox(self,'F-test requires same number of observations in each refinement','F-test not valid')
621            return
622        missingVars = []
623        for var in baseDict['varyList']: 
624            if var not in relxDict['varyList']: 
625                missingVars.append(var)
626        txt = ''
627        postmsg = ''
628        if missingVars:
629            txt = ('*** Possible invalid use of F-test: '+
630                'The F-test requires that the constrained model be a subset '+
631                'of the relaxed model. Review the parameters shown below to '+
632                'confirm missing var(s) have new names in the relaxed model. '+
633                '***\n\n')
634            postmsg = ('\n\nThese parameters are in the constrained '+
635              'fit and are not in the relaxed fit:\n*  ')
636            for i,var in enumerate(missingVars):
637                if i > 0: postmsg += ', '
638                postmsg += var
639            postmsg += ('\nThese parameters are in the relaxed fit and not'+
640                ' in the constrained fit:\n*  ')
641            i = 0
642            for var in relxDict['varyList']: 
643                if var not in baseDict['varyList']: 
644                    if i > 0: postmsg += ', '
645                    i += 1
646                    postmsg += var
647        #GSASIIpath.IPyBreak_base()
648        prob = RwFtest(baseDict['Rvals']['Nobs'],
649                   baseDict['Rvals']['GOF'],len(baseDict['varyList']),
650                   relxDict['Rvals']['GOF'],len(relxDict['varyList']))
651        fmt = "{} model is \n{}\n{} variables and Reduced Chi**2 = {:.3f}"
652        msg = txt
653        msg += fmt.format('Constrained',self.GPXtree.GetItemText(s0)[:-11],
654                       len(baseDict['varyList']),
655                       baseDict['Rvals']['GOF']**2)
656        msg += '\n\n'
657        msg += fmt.format('Relaxed',self.GPXtree.GetItemText(s1)[:-11],
658                       len(relxDict['varyList']),
659                       relxDict['Rvals']['GOF']**2)
660        msg += '\n\nCumulative F-test probability {:.2f}%\n'.format(prob*100)
661        if prob > 0.95:
662            msg += "The relaxed model is statistically preferred to the constrained model."
663        elif prob > 0.75:
664            msg += "There is little ability to differentiate between the two models on a statistical basis."
665        else:
666            msg += "The constrained model is statistically preferred to the relaxed model."
667        msg += postmsg
668        G2G.G2MessageBox(self,msg,'F-test result')
669       
670    def onHistPrinceTest(self,event):
671        '''Compare two histograms (selected here if more than two are present)
672        using the statistical test proposed by Ted Prince in
673        Acta Cryst. B35 1099-1100. (1982). Also see Int. Tables Vol. C
674        (1st Ed.) chapter 8.4, 618-621 (1995).
675        '''
676        items = []
677        item, cookie = self.GPXtree.GetFirstChild(self.root)
678        while item:
679            items.append(item)
680            item, cookie = self.GPXtree.GetNextChild(self.root, cookie)
681
682        if len(items) < 2:
683            G2G.G2MessageBox(self,'Prince test requires two histograms','Need more')
684            return
685        elif len(items) == 2:
686            s0,s1 = items
687        else:
688            # need to make selection here
689            sel = []
690            for i in items:
691                sel.append(self.GPXtree.GetItemText(i))
692            dlg = G2G.G2SingleChoiceDialog(self,'Select one refinement',
693                                             'Choose refinement',sel)
694            if dlg.ShowModal() == wx.ID_OK:
695                s0 = dlg.GetSelection()
696                dlg.Destroy()
697            else:
698                dlg.Destroy()
699                return
700            inds = list(range(len(items)))
701            del sel[s0]
702            del inds[s0]
703            dlg = G2G.G2SingleChoiceDialog(self,'Select comparison refinement',
704                                             'Choose refinement',sel)
705            if dlg.ShowModal() == wx.ID_OK:
706                s1 = dlg.GetSelection()
707                s1 = inds[s1]
708                dlg.Destroy()
709            else:
710                dlg.Destroy()
711                return
712        model0 = self.GPXtree.GetItemPyData(s0)
713        data0 = model0[1]
714        model1 = self.GPXtree.GetItemPyData(s1)
715        data1 = model1[1]
716        if len(data0[0]) != len(data1[0]):
717            G2G.G2MessageBox(self,'Unable to test: differing numbers of data points','Comparison not valid')
718            return
719        if max(abs((data0[0]-data1[0])/data0[0])) > 0.01:
720            G2G.G2MessageBox(self,'Unable to use test: "X" values differ','Comparison not valid')
721            return
722        X = data0[3] - data1[3]
723        #Z = np.sqrt(data0[3]) * (data0[1] - (data0[3] + data1[3])/2)
724        X = (data0[3] - data1[3]) / np.sqrt(data0[1])
725        Z = (data0[1] - (data0[3] + data1[3])/2) / np.sqrt(data0[1])
726        lam = np.sum(X*Z) / np.sum(X)
727        sig = np.sqrt(
728            (np.sum(Z*Z) - lam*lam*np.sum(X*X)) / 
729            ((len(data0[0]) - 1) * np.sum(X*X))
730            )
731           
732    0 the x-postions (two-theta in degrees),
733    1 the intensity values (Yobs),
734    2 the weights for each Yobs value
735    3 the computed intensity values (Ycalc)
736    4 the background values
737    5 Yobs-Ycalc
738       
739        GSASIIpath.IPyBreak_base()
740               
741#======================================================================   
742if __name__ == '__main__':
743    #if sys.platform == "darwin":
744    #    application = G2App(0) # create the GUI framework
745    #else:
746    application = wx.App(0) # create the GUI framework
747    try:
748        GSASIIpath.SetBinaryPath(True)
749    except:
750        print('Unable to run with current setup, do you want to update to the')
751        try:
752            if '2' in platform.python_version_tuple()[0]:           
753                ans = raw_input("latest GSAS-II version? Update ([Yes]/no): ")
754            else:
755                ans = input("latest GSAS-II version? Update ([Yes]/no): ")               
756        except:
757            ans = 'no'
758        if ans.strip().lower() == "no":
759            import sys
760            print('Exiting')
761            sys.exit()
762        print('Updating...')
763        GSASIIpath.svnUpdateProcess()
764    GSASIIpath.InvokeDebugOpts()
765    Frame = main(application) # start the GUI
766    argLoadlist = sys.argv[1:]
767    for arg in argLoadlist:
768        fil = os.path.splitext(arg)[0] + '.gpx'
769        if os.path.exists(fil):
770            Frame.fileList.append([fil,'GPX'])
771            Frame.loadFile(fil)
772        else:
773            print('File {} not found. Skipping'.format(fil))
774
775    # debug code to select in initial mode
776    #self=Frame
777    #self.Mode.FindItemById(self.wxID_Mode['Project']).Check(True)
778    #self.onRefresh(None)
779   
780    application.MainLoop()
Note: See TracBrowser for help on using the repository browser.