Changeset 4335

Timestamp:

Mar 1, 2020 5:26:57 PM (4 years ago)

Author:

toby

Message:

add autoInt routine, add f-test; updates for documentation

Location:

trunk

Files:

• G2AutoInt.py (added)
• G2compare.py (modified) (11 diffs)
• GSASIIautoInt.py (added)
• docs/source/GSASIIGUIr.rst (modified) (1 diff)
• docs/source/index.rst (modified) (2 diffs)

trunk/G2compare.py

@@ -10 +10 @@
 ########### SVN repository information ###################
 '''
+*G2compare: Tool for project comparison*
+---------------------------------------------
 '''
 
@@ -46 +48 @@
 
 __version__ = '0.0.1'
+
+# math to do F-test 
+def RC2Ftest(npts,RChiSq0,nvar0,RChiSq1,nvar1):
+    '''Compute the F-test probability that a model expanded with added 
+    parameters (relaxed model) is statistically more likely than the 
+    constrained (base) model
+    :param int npts: number of observed diffraction data points
+    :param float RChiSq0: Reduced Chi**2 for the base model
+    :param int nvar0: number of refined variables in the base model
+    :param float RChiSq0: Reduced Chi**2 for the relaxed model
+    :param int nvar1: number of refined variables in the relaxed model
+    '''
+    if nvar1 == nvar0:
+        raise Exception("parameter # agree, test is not valid")
+    elif nvar1 < nvar0:
+        print("Warning: RwFtest used with base/relaxed models reversed")
+        RChiSq0,nvar0,RChiSq1,nvar1 = RChiSq1,nvar1,RChiSq0,nvar0
+    ratio = RChiSq0 / RChiSq1
+    nu1 = float(nvar1 - nvar0)
+    nu2 = float(npts - nvar1)
+    F = (ratio - 1.) * nu2 / nu1
+    import scipy.stats
+    return scipy.stats.f.cdf(F,nu1,nu2)
+
+def RwFtest(npts,Rwp0,nvar0,Rwp1,nvar1):
+    '''Compute the F-test probability that a model expanded with added 
+    parameters (relaxed model) is statistically more likely than the 
+    constrained (base) model
+    :param int npts: number of observed diffraction data points
+    :param float Rwp0: Weighted profile R-factor or GOF for the base model
+    :param int nvar0: number of refined variables in the base model
+    :param float Rwp1: Weighted profile R-factor or GOF for the relaxed model
+    :param int nvar1: number of refined variables in the relaxed model
+    '''
+    if nvar1 == nvar0:
+        raise Exception("parameter # agree, test is not valid")
+    elif nvar1 < nvar0:
+        print("Warning: RwFtest used with base/relaxed models reversed")
+        Rwp0,nvar0,Rwp1,nvar1 = Rwp1,nvar1,Rwp0,nvar0
+    ratio = (Rwp0 / Rwp1)**2
+    nu1 = float(nvar1 - nvar0)
+    nu2 = float(npts - nvar1)
+    F = (ratio - 1.) * nu2 / nu1
+    import scipy.stats
+    return scipy.stats.f.cdf(F,nu1,nu2)
 
 def cPickleLoad(fp):
@@ -89 +136 @@
             style=wx.DEFAULT_FRAME_STYLE ^ wx.CLOSE_BOX)
         self.G2plotNB = G2plt.G2PlotNoteBook(self.plotFrame,G2frame=self)
-        self.plotFrame.Show()
+        #self.plotFrame.Show()
+        self.plotFrame.Show(False)  # until we have some graphics, hide the plot window
         # menus
         Frame = self.GetTopLevelParent() # same as self
-        Menu = wx.MenuBar()
+        self.MenuBar = wx.MenuBar()
         File = wx.Menu(title='')
-        Menu.Append(menu=File, title='&File')
+        self.MenuBar.Append(menu=File, title='&File')
         item = File.Append(wx.ID_ANY,'&Import project...\tCtrl+O','Open a GSAS-II project file (*.gpx)')            
         self.Bind(wx.EVT_MENU, self.onLoadGPX, id=item.GetId())
@@ -101 +149 @@
 
         self.Mode = wx.Menu(title='')
-        Menu.Append(menu=self.Mode, title='&Mode')
+        self.MenuBar.Append(menu=self.Mode, title='&Mode')
         self.wxID_Mode = {}
         for m in "Histogram","Phase","Project":
@@ -109 +157 @@
         item = self.Mode.Append(wx.ID_ANY,'Set histogram filter','Set a filter for histograms to display')
         self.Bind(wx.EVT_MENU, self.onHistFilter, id=item.GetId())
-        
-        Frame.SetMenuBar(Menu)
+        modeMenu = wx.Menu(title='')
+        self.MenuBar.Append(menu=modeMenu, title='TBD')
+        self.modeMenuPos = self.MenuBar.FindMenu('TBD')
+        #item = modeMenu.Append(wx.ID_ANY,'test') # menu can't be empty
+        
+        Frame.SetMenuBar(self.MenuBar)
         # status bar
         self.Status = self.CreateStatusBar()
@@ -156 +208 @@
 
         self.PWDRfilter = None
+        self.SetModeMenu()
         
     def SelectGPX(self):
@@ -192 +245 @@
         for fil,mode in self.fileList:
             self.loadFile(fil)
-        
+        self.SetModeMenu()
+            
+    def SetModeMenu(self):
+        '''Create the mode-specific menu and its contents
+        '''
+        modeMenu = wx.Menu(title='')
+        oldmenu = self.MenuBar.Replace(self.modeMenuPos,modeMenu,self.getMode())
+        wx.CallAfter(oldmenu.Destroy)
+        if self.getMode() == "Histogram":
+            pass
+        elif self.getMode() == "Phase":
+            pass
+        elif self.getMode() == "Project":
+            item = modeMenu.Append(wx.ID_ANY,'F-test')
+            self.Bind(wx.EVT_MENU, self.onProjFtest, id=item.GetId())
+
     def loadFile(self,fil):
         '''read or reread a file
@@ -377 +445 @@
                 if not data[0][0].startswith('Covariance'): continue
                 Covar = data[0]
-                GSASIIpath.IPyBreak_base()
+                #GSASIIpath.IPyBreak_base()
                 #if self.PWDRfilter is not None: # implement filter
                 #    if self.PWDRfilter not in data[0][0]: continue
@@ -424 +492 @@
     #     else:
     #         event.Skip(False)
-
+    
+    def onProjFtest(self,event):
+        '''Compare two projects (selected here if more than two are present)
+        using the statistical F-test (aka Hamilton R-factor test), see:
+
+            * Hamilton, R. W. (1965), Acta Crystallogr. 18, 502-510.
+            * Prince, E., Mathematical Techniques in Crystallography and Materials Science, Second ed. (Springer-Verlag, New York, 1994).
+        '''
+        items = []
+        item, cookie = self.GPXtree.GetFirstChild(self.root)
+        while item:
+            items.append(item)
+            item, cookie = self.GPXtree.GetNextChild(self.root, cookie)
+        if len(items) < 2:
+            G2G.G2MessageBox(self,'F-test requires two more more projects','Need more projects')
+            return
+        elif len(items) == 2:
+            s0 = items[0]
+            baseDict = self.GPXtree.GetItemPyData(s0)
+            s1 = items[1]            
+            relxDict = self.GPXtree.GetItemPyData(s1)
+            # sort out the dicts in order of number of refined variables
+            if len(baseDict['varyList']) > len(relxDict['varyList']):
+                s0,s1,baseDict,relxDict = s1,s0,relxDict,baseDict
+        else:
+            # need to make selection here
+            sel = []
+            for i in items:
+                sel.append(self.GPXtree.GetItemText(i))
+            dlg = G2G.G2SingleChoiceDialog(self,'Select constrained refinement',
+                                             'Choose refinement',sel)
+            if dlg.ShowModal() == wx.ID_OK:
+                s0 = dlg.GetSelection()
+                dlg.Destroy()
+            else:
+                dlg.Destroy()
+                return
+            inds = list(range(len(items)))
+            del sel[s0]
+            del inds[s0]
+            dlg = G2G.G2SingleChoiceDialog(self,'Select relaxed refinement',
+                                             'Choose refinement',sel)
+            if dlg.ShowModal() == wx.ID_OK:
+                s1 = dlg.GetSelection()
+                s1 = inds[s1]
+                dlg.Destroy()
+            else:
+                dlg.Destroy()
+                return
+            baseDict = self.GPXtree.GetItemPyData(items[s0])
+            relxDict = self.GPXtree.GetItemPyData(items[s1])
+            if len(baseDict['varyList']) > len(relxDict['varyList']):
+                G2G.G2MessageBox(self,
+                            'F-test warning: constrained refinement has more '+
+                            'variables ({}) than relaxed refinement ({}). Swapping'
+                            .format(len(baseDict['varyList']), len(relxDict['varyList'])),
+                            'Fits reversed?')
+                s0,s1,baseDict,relxDict = s1,s0,relxDict,baseDict
+                baseDict,relxDict = relxDict,baseDict
+        if len(baseDict['varyList']) == len(relxDict['varyList']):
+            G2G.G2MessageBox(self,'F-test requires differing numbers of variables','F-test not valid')
+            return
+        elif baseDict['Rvals']['Nobs'] != relxDict['Rvals']['Nobs']:
+            G2G.G2MessageBox(self,'F-test requires same number of observations in each refinement','F-test not valid')
+            return
+        prob = RwFtest(baseDict['Rvals']['Nobs'],
+                   baseDict['Rvals']['GOF'],len(baseDict['varyList']),
+                   relxDict['Rvals']['GOF'],len(relxDict['varyList']))
+        fmt = "{} model is \n*  {}\n*  {} variables and Reduced Chi**2 = {:.3f}"
+        msg = fmt.format('Constrained',self.GPXtree.GetItemText(s0)[:-11],
+                       len(baseDict['varyList']),
+                       baseDict['Rvals']['GOF']**2)
+        msg += '\n\n'
+        msg += fmt.format('Relaxed',self.GPXtree.GetItemText(s1)[:-11],
+                       len(relxDict['varyList']),
+                       relxDict['Rvals']['GOF']**2)
+        msg += '\n\nCumulative F-test probability {:.2f}%\n'.format(prob*100)
+        if prob > 0.95:
+            msg += "The relaxed model is statistically preferred to the constrained model."
+        elif prob > 0.75:
+            msg += "There is little ability to differentiate between the two models on a statistical basis."
+        else:
+            msg += "The constrained model is statistically preferred to the relaxed model."
+        G2G.G2MessageBox(self,msg,'F-test result')
+        #GSASIIpath.IPyBreak_base()
+    
 if __name__ == '__main__':
     #if sys.platform == "darwin": 
@@ -450 +603 @@
     Frame = main(application) # start the GUI
     argLoadlist = sys.argv[1:]
-    if len(argLoadlist) == 0:
-        argLoadlist = ['/Users/toby/Scratch/copy.gpx',
-                       '/Users/toby/Scratch/CW_YAG.gpx']
     for arg in argLoadlist:
         fil = os.path.splitext(arg)[0] + '.gpx'
@@ -460 +610 @@
         else:
             print('File {} not found. Skipping'.format(fil))
+
+    # debug code to select in initial mode
+    #self=Frame
+    #self.Mode.FindItemById(self.wxID_Mode['Project']).Check(True)
+    #self.onRefresh(None)
+    
     application.MainLoop()

trunk/docs/source/GSASIIGUIr.rst

@@ -1 @@
-*GSAS-II GUI Support Routines*
+*GSAS-II GUI Utility Modules*
 ==============================
 

trunk/docs/source/index.rst

@@ -8 +8 @@
 withing the GSAS-II framework, those planning to understand how
 GSAS-II works, or for people wishing to develop scripting applications
-using the API (:mod:`GSASIIscriptable`).
+using the API (:mod:`GSASIIscriptable`). Note that many data structures
+used in GSAS-II are defined in module :mod:`GSASIIobj`.
 
 For information on obtaining or learning to use GSAS-II, please see
@@ -38 +39 @@
   imports.rst 
   exports.rst
+  G2tools.rst