Changeset 4335


Ignore:
Timestamp:
Mar 1, 2020 5:26:57 PM (20 months ago)
Author:
toby
Message:

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

Location:
trunk
Files:
2 added
3 edited

Legend:

Unmodified
Added
Removed
  • trunk/G2compare.py

    r4262 r4335  
    1010########### SVN repository information ###################
    1111'''
     12*G2compare: Tool for project comparison*
     13---------------------------------------------
    1214'''
    1315
     
    4648
    4749__version__ = '0.0.1'
     50
     51# math to do F-test
     52def RC2Ftest(npts,RChiSq0,nvar0,RChiSq1,nvar1):
     53    '''Compute the F-test probability that a model expanded with added
     54    parameters (relaxed model) is statistically more likely than the
     55    constrained (base) model
     56    :param int npts: number of observed diffraction data points
     57    :param float RChiSq0: Reduced Chi**2 for the base model
     58    :param int nvar0: number of refined variables in the base model
     59    :param float RChiSq0: Reduced Chi**2 for the relaxed model
     60    :param int nvar1: number of refined variables in the relaxed model
     61    '''
     62    if nvar1 == nvar0:
     63        raise Exception("parameter # agree, test is not valid")
     64    elif nvar1 < nvar0:
     65        print("Warning: RwFtest used with base/relaxed models reversed")
     66        RChiSq0,nvar0,RChiSq1,nvar1 = RChiSq1,nvar1,RChiSq0,nvar0
     67    ratio = RChiSq0 / RChiSq1
     68    nu1 = float(nvar1 - nvar0)
     69    nu2 = float(npts - nvar1)
     70    F = (ratio - 1.) * nu2 / nu1
     71    import scipy.stats
     72    return scipy.stats.f.cdf(F,nu1,nu2)
     73
     74def RwFtest(npts,Rwp0,nvar0,Rwp1,nvar1):
     75    '''Compute the F-test probability that a model expanded with added
     76    parameters (relaxed model) is statistically more likely than the
     77    constrained (base) model
     78    :param int npts: number of observed diffraction data points
     79    :param float Rwp0: Weighted profile R-factor or GOF for the base model
     80    :param int nvar0: number of refined variables in the base model
     81    :param float Rwp1: Weighted profile R-factor or GOF for the relaxed model
     82    :param int nvar1: number of refined variables in the relaxed model
     83    '''
     84    if nvar1 == nvar0:
     85        raise Exception("parameter # agree, test is not valid")
     86    elif nvar1 < nvar0:
     87        print("Warning: RwFtest used with base/relaxed models reversed")
     88        Rwp0,nvar0,Rwp1,nvar1 = Rwp1,nvar1,Rwp0,nvar0
     89    ratio = (Rwp0 / Rwp1)**2
     90    nu1 = float(nvar1 - nvar0)
     91    nu2 = float(npts - nvar1)
     92    F = (ratio - 1.) * nu2 / nu1
     93    import scipy.stats
     94    return scipy.stats.f.cdf(F,nu1,nu2)
    4895
    4996def cPickleLoad(fp):
     
    89136            style=wx.DEFAULT_FRAME_STYLE ^ wx.CLOSE_BOX)
    90137        self.G2plotNB = G2plt.G2PlotNoteBook(self.plotFrame,G2frame=self)
    91         self.plotFrame.Show()
     138        #self.plotFrame.Show()
     139        self.plotFrame.Show(False)  # until we have some graphics, hide the plot window
    92140        # menus
    93141        Frame = self.GetTopLevelParent() # same as self
    94         Menu = wx.MenuBar()
     142        self.MenuBar = wx.MenuBar()
    95143        File = wx.Menu(title='')
    96         Menu.Append(menu=File, title='&File')
     144        self.MenuBar.Append(menu=File, title='&File')
    97145        item = File.Append(wx.ID_ANY,'&Import project...\tCtrl+O','Open a GSAS-II project file (*.gpx)')           
    98146        self.Bind(wx.EVT_MENU, self.onLoadGPX, id=item.GetId())
     
    101149
    102150        self.Mode = wx.Menu(title='')
    103         Menu.Append(menu=self.Mode, title='&Mode')
     151        self.MenuBar.Append(menu=self.Mode, title='&Mode')
    104152        self.wxID_Mode = {}
    105153        for m in "Histogram","Phase","Project":
     
    109157        item = self.Mode.Append(wx.ID_ANY,'Set histogram filter','Set a filter for histograms to display')
    110158        self.Bind(wx.EVT_MENU, self.onHistFilter, id=item.GetId())
    111        
    112         Frame.SetMenuBar(Menu)
     159        modeMenu = wx.Menu(title='')
     160        self.MenuBar.Append(menu=modeMenu, title='TBD')
     161        self.modeMenuPos = self.MenuBar.FindMenu('TBD')
     162        #item = modeMenu.Append(wx.ID_ANY,'test') # menu can't be empty
     163       
     164        Frame.SetMenuBar(self.MenuBar)
    113165        # status bar
    114166        self.Status = self.CreateStatusBar()
     
    156208
    157209        self.PWDRfilter = None
     210        self.SetModeMenu()
    158211       
    159212    def SelectGPX(self):
     
    192245        for fil,mode in self.fileList:
    193246            self.loadFile(fil)
    194        
     247        self.SetModeMenu()
     248           
     249    def SetModeMenu(self):
     250        '''Create the mode-specific menu and its contents
     251        '''
     252        modeMenu = wx.Menu(title='')
     253        oldmenu = self.MenuBar.Replace(self.modeMenuPos,modeMenu,self.getMode())
     254        wx.CallAfter(oldmenu.Destroy)
     255        if self.getMode() == "Histogram":
     256            pass
     257        elif self.getMode() == "Phase":
     258            pass
     259        elif self.getMode() == "Project":
     260            item = modeMenu.Append(wx.ID_ANY,'F-test')
     261            self.Bind(wx.EVT_MENU, self.onProjFtest, id=item.GetId())
     262
    195263    def loadFile(self,fil):
    196264        '''read or reread a file
     
    377445                if not data[0][0].startswith('Covariance'): continue
    378446                Covar = data[0]
    379                 GSASIIpath.IPyBreak_base()
     447                #GSASIIpath.IPyBreak_base()
    380448                #if self.PWDRfilter is not None: # implement filter
    381449                #    if self.PWDRfilter not in data[0][0]: continue
     
    424492    #     else:
    425493    #         event.Skip(False)
    426 
     494   
     495    def onProjFtest(self,event):
     496        '''Compare two projects (selected here if more than two are present)
     497        using the statistical F-test (aka Hamilton R-factor test), see:
     498
     499            * Hamilton, R. W. (1965), Acta Crystallogr. 18, 502-510.
     500            * Prince, E., Mathematical Techniques in Crystallography and Materials Science, Second ed. (Springer-Verlag, New York, 1994).
     501        '''
     502        items = []
     503        item, cookie = self.GPXtree.GetFirstChild(self.root)
     504        while item:
     505            items.append(item)
     506            item, cookie = self.GPXtree.GetNextChild(self.root, cookie)
     507        if len(items) < 2:
     508            G2G.G2MessageBox(self,'F-test requires two more more projects','Need more projects')
     509            return
     510        elif len(items) == 2:
     511            s0 = items[0]
     512            baseDict = self.GPXtree.GetItemPyData(s0)
     513            s1 = items[1]           
     514            relxDict = self.GPXtree.GetItemPyData(s1)
     515            # sort out the dicts in order of number of refined variables
     516            if len(baseDict['varyList']) > len(relxDict['varyList']):
     517                s0,s1,baseDict,relxDict = s1,s0,relxDict,baseDict
     518        else:
     519            # need to make selection here
     520            sel = []
     521            for i in items:
     522                sel.append(self.GPXtree.GetItemText(i))
     523            dlg = G2G.G2SingleChoiceDialog(self,'Select constrained refinement',
     524                                             'Choose refinement',sel)
     525            if dlg.ShowModal() == wx.ID_OK:
     526                s0 = dlg.GetSelection()
     527                dlg.Destroy()
     528            else:
     529                dlg.Destroy()
     530                return
     531            inds = list(range(len(items)))
     532            del sel[s0]
     533            del inds[s0]
     534            dlg = G2G.G2SingleChoiceDialog(self,'Select relaxed refinement',
     535                                             'Choose refinement',sel)
     536            if dlg.ShowModal() == wx.ID_OK:
     537                s1 = dlg.GetSelection()
     538                s1 = inds[s1]
     539                dlg.Destroy()
     540            else:
     541                dlg.Destroy()
     542                return
     543            baseDict = self.GPXtree.GetItemPyData(items[s0])
     544            relxDict = self.GPXtree.GetItemPyData(items[s1])
     545            if len(baseDict['varyList']) > len(relxDict['varyList']):
     546                G2G.G2MessageBox(self,
     547                            'F-test warning: constrained refinement has more '+
     548                            'variables ({}) than relaxed refinement ({}). Swapping'
     549                            .format(len(baseDict['varyList']), len(relxDict['varyList'])),
     550                            'Fits reversed?')
     551                s0,s1,baseDict,relxDict = s1,s0,relxDict,baseDict
     552                baseDict,relxDict = relxDict,baseDict
     553        if len(baseDict['varyList']) == len(relxDict['varyList']):
     554            G2G.G2MessageBox(self,'F-test requires differing numbers of variables','F-test not valid')
     555            return
     556        elif baseDict['Rvals']['Nobs'] != relxDict['Rvals']['Nobs']:
     557            G2G.G2MessageBox(self,'F-test requires same number of observations in each refinement','F-test not valid')
     558            return
     559        prob = RwFtest(baseDict['Rvals']['Nobs'],
     560                   baseDict['Rvals']['GOF'],len(baseDict['varyList']),
     561                   relxDict['Rvals']['GOF'],len(relxDict['varyList']))
     562        fmt = "{} model is \n*  {}\n*  {} variables and Reduced Chi**2 = {:.3f}"
     563        msg = fmt.format('Constrained',self.GPXtree.GetItemText(s0)[:-11],
     564                       len(baseDict['varyList']),
     565                       baseDict['Rvals']['GOF']**2)
     566        msg += '\n\n'
     567        msg += fmt.format('Relaxed',self.GPXtree.GetItemText(s1)[:-11],
     568                       len(relxDict['varyList']),
     569                       relxDict['Rvals']['GOF']**2)
     570        msg += '\n\nCumulative F-test probability {:.2f}%\n'.format(prob*100)
     571        if prob > 0.95:
     572            msg += "The relaxed model is statistically preferred to the constrained model."
     573        elif prob > 0.75:
     574            msg += "There is little ability to differentiate between the two models on a statistical basis."
     575        else:
     576            msg += "The constrained model is statistically preferred to the relaxed model."
     577        G2G.G2MessageBox(self,msg,'F-test result')
     578        #GSASIIpath.IPyBreak_base()
     579   
    427580if __name__ == '__main__':
    428581    #if sys.platform == "darwin":
     
    450603    Frame = main(application) # start the GUI
    451604    argLoadlist = sys.argv[1:]
    452     if len(argLoadlist) == 0:
    453         argLoadlist = ['/Users/toby/Scratch/copy.gpx',
    454                        '/Users/toby/Scratch/CW_YAG.gpx']
    455605    for arg in argLoadlist:
    456606        fil = os.path.splitext(arg)[0] + '.gpx'
     
    460610        else:
    461611            print('File {} not found. Skipping'.format(fil))
     612
     613    # debug code to select in initial mode
     614    #self=Frame
     615    #self.Mode.FindItemById(self.wxID_Mode['Project']).Check(True)
     616    #self.onRefresh(None)
     617   
    462618    application.MainLoop()
  • trunk/docs/source/GSASIIGUIr.rst

    r3000 r4335  
    1 *GSAS-II GUI Support Routines*
     1*GSAS-II GUI Utility Modules*
    22==============================
    33
  • trunk/docs/source/index.rst

    r4249 r4335  
    88withing the GSAS-II framework, those planning to understand how
    99GSAS-II works, or for people wishing to develop scripting applications
    10 using the API (:mod:`GSASIIscriptable`).
     10using the API (:mod:`GSASIIscriptable`). Note that many data structures
     11used in GSAS-II are defined in module :mod:`GSASIIobj`.
    1112
    1213For information on obtaining or learning to use GSAS-II, please see
     
    3839  imports.rst
    3940  exports.rst
     41  G2tools.rst
    4042
Note: See TracChangeset for help on using the changeset viewer.