| 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 | |
|---|
| 24 | import sys |
|---|
| 25 | import os |
|---|
| 26 | import platform |
|---|
| 27 | if '2' in platform.python_version_tuple()[0]: |
|---|
| 28 | import cPickle |
|---|
| 29 | else: |
|---|
| 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 | |
|---|
| 36 | import wx |
|---|
| 37 | import numpy as np |
|---|
| 38 | import matplotlib as mpl |
|---|
| 39 | try: |
|---|
| 40 | import OpenGL as ogl |
|---|
| 41 | except ImportError: |
|---|
| 42 | pass |
|---|
| 43 | import scipy as sp |
|---|
| 44 | |
|---|
| 45 | import GSASIIpath |
|---|
| 46 | GSASIIpath.SetVersionNumber("$Revision: 4364 $") |
|---|
| 47 | import GSASIIfiles as G2fil |
|---|
| 48 | import GSASIIplot as G2plt |
|---|
| 49 | import GSASIIctrlGUI as G2G |
|---|
| 50 | import GSASIIobj as G2obj |
|---|
| 51 | |
|---|
| 52 | __version__ = '0.0.1' |
|---|
| 53 | |
|---|
| 54 | # math to do F-test |
|---|
| 55 | def 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 | |
|---|
| 77 | def 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 | |
|---|
| 99 | def 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 | |
|---|
| 105 | def 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 | |
|---|
| 130 | class 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 | #====================================================================== |
|---|
| 744 | if __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() |
|---|
