source: trunk/G2compare.py @ 4364

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

minor docs changes; remove remaining .xor. reference

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