1 | #!/usr/bin/env python |
---|
2 | # -*- coding: utf-8 -*- |
---|
3 | #GSAS-II Data/Model Comparison |
---|
4 | ########### SVN repository information ################### |
---|
5 | # $Date: $ |
---|
6 | # $Author: toby $ |
---|
7 | # $Revision: $ |
---|
8 | # $URL: $ |
---|
9 | # $Id: $ |
---|
10 | ########### SVN repository information ################### |
---|
11 | ''' |
---|
12 | *G2compare: Tool for project comparison* |
---|
13 | --------------------------------------------- |
---|
14 | ''' |
---|
15 | |
---|
16 | #TODO: in OnDataTreeSelChanged want to plot patterns |
---|
17 | # Make PWDR unique (use histlist) |
---|
18 | # add graphics |
---|
19 | # implement project |
---|
20 | |
---|
21 | import sys |
---|
22 | import os |
---|
23 | import platform |
---|
24 | if '2' in platform.python_version_tuple()[0]: |
---|
25 | import cPickle |
---|
26 | else: |
---|
27 | try: |
---|
28 | import _pickle as cPickle |
---|
29 | except: |
---|
30 | print('Warning: failed to import the optimized Py3 pickle (_pickle)') |
---|
31 | import pickle as cPickle |
---|
32 | |
---|
33 | import wx |
---|
34 | import numpy as np |
---|
35 | import matplotlib as mpl |
---|
36 | try: |
---|
37 | import OpenGL as ogl |
---|
38 | except ImportError: |
---|
39 | pass |
---|
40 | import scipy as sp |
---|
41 | |
---|
42 | import GSASIIpath |
---|
43 | GSASIIpath.SetVersionNumber("$Revision: 4154 $") |
---|
44 | import GSASIIfiles as G2fil |
---|
45 | import GSASIIplot as G2plt |
---|
46 | import GSASIIctrlGUI as G2G |
---|
47 | import GSASIIobj as G2obj |
---|
48 | |
---|
49 | __version__ = '0.0.1' |
---|
50 | |
---|
51 | # math to do F-test |
---|
52 | def 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 | |
---|
74 | def 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) |
---|
95 | |
---|
96 | def cPickleLoad(fp): |
---|
97 | if '2' in platform.python_version_tuple()[0]: |
---|
98 | return cPickle.load(fp) |
---|
99 | else: |
---|
100 | return cPickle.load(fp,encoding='latin-1') |
---|
101 | |
---|
102 | def main(application): |
---|
103 | '''Start up the GSAS-II GUI''' |
---|
104 | knownVersions = ['2.7','3.6','3.7','3.8'] |
---|
105 | if platform.python_version()[:3] not in knownVersions: |
---|
106 | dlg = wx.MessageDialog(None, |
---|
107 | 'GSAS-II requires Python 2.7.x or 3.6+\n Yours is '+sys.version.split()[0], |
---|
108 | 'Python version error', wx.OK) |
---|
109 | try: |
---|
110 | dlg.ShowModal() |
---|
111 | finally: |
---|
112 | dlg.Destroy() |
---|
113 | sys.exit() |
---|
114 | |
---|
115 | application.main = MakeTopWindow(None) # application.main is the main wx.Frame |
---|
116 | application.SetTopWindow(application.main) |
---|
117 | # save the current package versions |
---|
118 | application.main.PackageVersions = G2fil.get_python_versions([wx, mpl, np, sp, ogl]) |
---|
119 | try: |
---|
120 | application.SetAppDisplayName('GSAS-II Compare') |
---|
121 | except: |
---|
122 | pass |
---|
123 | #application.GetTopWindow().SendSizeEvent() |
---|
124 | application.GetTopWindow().Show(True) |
---|
125 | return application.GetTopWindow() |
---|
126 | |
---|
127 | class MakeTopWindow(wx.Frame): |
---|
128 | '''Define the main frame and its associated menu items |
---|
129 | ''' |
---|
130 | def __init__(self, parent): |
---|
131 | size = wx.Size(700,450) |
---|
132 | wx.Frame.__init__(self, name='dComp', parent=parent, |
---|
133 | size=size,style=wx.DEFAULT_FRAME_STYLE, title='GSAS-II data/model comparison') |
---|
134 | # plot window |
---|
135 | self.plotFrame = wx.Frame(None,-1,'dComp Plots',size=size, |
---|
136 | style=wx.DEFAULT_FRAME_STYLE ^ wx.CLOSE_BOX) |
---|
137 | self.G2plotNB = G2plt.G2PlotNoteBook(self.plotFrame,G2frame=self) |
---|
138 | #self.plotFrame.Show() |
---|
139 | self.plotFrame.Show(False) # until we have some graphics, hide the plot window |
---|
140 | # menus |
---|
141 | Frame = self.GetTopLevelParent() # same as self |
---|
142 | self.MenuBar = wx.MenuBar() |
---|
143 | File = wx.Menu(title='') |
---|
144 | self.MenuBar.Append(menu=File, title='&File') |
---|
145 | item = File.Append(wx.ID_ANY,'&Import project...\tCtrl+O','Open a GSAS-II project file (*.gpx)') |
---|
146 | self.Bind(wx.EVT_MENU, self.onLoadGPX, id=item.GetId()) |
---|
147 | # item = File.Append(wx.ID_ANY,'&Import selected...','Open a GSAS-II project file (*.gpx)') |
---|
148 | # self.Bind(wx.EVT_MENU, self.onLoadSel, id=item.GetId()) |
---|
149 | |
---|
150 | self.Mode = wx.Menu(title='') |
---|
151 | self.MenuBar.Append(menu=self.Mode, title='&Mode') |
---|
152 | self.wxID_Mode = {} |
---|
153 | for m in "Histogram","Phase","Project": |
---|
154 | i = self.wxID_Mode[m] = wx.NewId() |
---|
155 | item = self.Mode.AppendRadioItem(i,m,'Display {}s'.format(m)) |
---|
156 | self.Bind(wx.EVT_MENU, self.onRefresh, id=item.GetId()) |
---|
157 | item = self.Mode.Append(wx.ID_ANY,'Set histogram filter','Set a filter for histograms to display') |
---|
158 | self.Bind(wx.EVT_MENU, self.onHistFilter, id=item.GetId()) |
---|
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) |
---|
165 | # status bar |
---|
166 | self.Status = self.CreateStatusBar() |
---|
167 | self.Status.SetFieldsCount(2) |
---|
168 | # split the frame and add the tree |
---|
169 | self.mainPanel = wx.SplitterWindow(self, wx.ID_ANY, style=wx.SP_LIVE_UPDATE|wx.SP_3D) |
---|
170 | self.mainPanel.SetMinimumPaneSize(100) |
---|
171 | self.treePanel = wx.Panel(self.mainPanel, wx.ID_ANY, |
---|
172 | style = wx.TAB_TRAVERSAL|wx.SUNKEN_BORDER) |
---|
173 | |
---|
174 | # self.dataWindow = G2DataWindow(self.mainPanel) |
---|
175 | self.dataWindow = wx.ScrolledWindow(self.mainPanel) |
---|
176 | dataSizer = wx.BoxSizer(wx.VERTICAL) |
---|
177 | self.dataWindow.SetSizer(dataSizer) |
---|
178 | self.mainPanel.SplitVertically(self.treePanel, self.dataWindow, 200) |
---|
179 | self.Status.SetStatusWidths([200,-1]) # make these match? |
---|
180 | |
---|
181 | treeSizer = wx.BoxSizer(wx.VERTICAL) |
---|
182 | self.treePanel.SetSizer(treeSizer) |
---|
183 | self.GPXtree = G2G.G2TreeCtrl(id=wx.ID_ANY, |
---|
184 | parent=self.treePanel, size=self.treePanel.GetClientSize(),style=wx.TR_DEFAULT_STYLE ) |
---|
185 | TreeId = self.GPXtree.Id |
---|
186 | |
---|
187 | treeSizer.Add(self.GPXtree,1,wx.EXPAND|wx.ALL,0) |
---|
188 | #self.GPXtree.Bind(wx.EVT_TREE_SEL_CHANGED,self.OnDataTreeSelChanged) |
---|
189 | self.GPXtree.Bind(wx.EVT_TREE_SEL_CHANGED,self.OnDataTreeSelChanged) |
---|
190 | # self.GPXtree.Bind(wx.EVT_TREE_ITEM_RIGHT_CLICK,self.OnDataTreeSelChanged) |
---|
191 | # self.GPXtree.Bind(wx.EVT_TREE_ITEM_COLLAPSED, |
---|
192 | # self.OnGPXtreeItemCollapsed, id=TreeId) |
---|
193 | #self.GPXtree.Bind(wx.EVT_TREE_ITEM_EXPANDED, |
---|
194 | # self.OnGPXtreeItemExpanded, id=TreeId) |
---|
195 | # self.GPXtree.Bind(wx.EVT_TREE_DELETE_ITEM, |
---|
196 | # self.OnGPXtreeItemDelete, id=TreeId) |
---|
197 | # self.GPXtree.Bind(wx.EVT_TREE_KEY_DOWN, |
---|
198 | # self.OnGPXtreeKeyDown, id=TreeId) |
---|
199 | # self.GPXtree.Bind(wx.EVT_TREE_BEGIN_RDRAG, |
---|
200 | # self.OnGPXtreeBeginRDrag, id=TreeId) |
---|
201 | # self.GPXtree.Bind(wx.EVT_TREE_END_DRAG, |
---|
202 | # self.OnGPXtreeEndDrag, id=TreeId) |
---|
203 | self.root = self.GPXtree.root |
---|
204 | self.Bind(wx.EVT_CLOSE, lambda event: sys.exit()) |
---|
205 | |
---|
206 | self.fileList = [] # list of files read for use in Reload |
---|
207 | self.histList = [] # list of histograms loaded for unique naming |
---|
208 | |
---|
209 | self.PWDRfilter = None |
---|
210 | self.SetModeMenu() |
---|
211 | |
---|
212 | def SelectGPX(self): |
---|
213 | '''Select a .GPX file to be read |
---|
214 | ''' |
---|
215 | dlg = wx.FileDialog(self, 'Choose GSAS-II project file', |
---|
216 | wildcard='GSAS-II project file (*.gpx)|*.gpx',style=wx.FD_OPEN) |
---|
217 | try: |
---|
218 | if dlg.ShowModal() != wx.ID_OK: return |
---|
219 | fil = os.path.splitext(dlg.GetPath())[0]+'.gpx' |
---|
220 | finally: |
---|
221 | dlg.Destroy() |
---|
222 | if os.path.exists(fil): |
---|
223 | self.fileList.append([fil,'GPX']) |
---|
224 | return fil |
---|
225 | else: |
---|
226 | print('File {} not found, skipping'.format(fil)) |
---|
227 | return |
---|
228 | |
---|
229 | def getMode(self): |
---|
230 | '''returns the display mode (one of "Histogram","Phase","Project"). |
---|
231 | Could return '?' in case of an error. |
---|
232 | ''' |
---|
233 | for m in self.wxID_Mode: |
---|
234 | if self.Mode.FindItemById(self.wxID_Mode[m]).IsChecked(): |
---|
235 | break |
---|
236 | else: |
---|
237 | m = '?' |
---|
238 | return m |
---|
239 | |
---|
240 | def onRefresh(self,event): |
---|
241 | '''reread all files, in response to a change in mode, etc. |
---|
242 | ''' |
---|
243 | self.GPXtree.DeleteChildren(self.root) # delete tree contents |
---|
244 | self.histList = [] # clear list of loaded histograms |
---|
245 | for fil,mode in self.fileList: |
---|
246 | self.loadFile(fil) |
---|
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 | |
---|
263 | def loadFile(self,fil): |
---|
264 | '''read or reread a file |
---|
265 | ''' |
---|
266 | if self.getMode() == "Histogram": |
---|
267 | self.LoadPwdr(fil) |
---|
268 | elif self.getMode() == "Phase": |
---|
269 | self.LoadPhase(fil) |
---|
270 | elif self.getMode() == "Project": |
---|
271 | self.LoadProject(fil) |
---|
272 | else: |
---|
273 | print("mode not implemented") |
---|
274 | #raise Exception("mode not implemented") |
---|
275 | |
---|
276 | def onLoadGPX(self,event): |
---|
277 | '''Initial load of GPX file in response to a menu command |
---|
278 | ''' |
---|
279 | fil = self.SelectGPX() |
---|
280 | if not fil: return |
---|
281 | if not os.path.exists(fil): return |
---|
282 | self.fileList.append([fil,'GPX']) |
---|
283 | self.loadFile(fil) |
---|
284 | |
---|
285 | def LoadPwdr(self,fil): |
---|
286 | '''Load PWDR entries from a .GPX file to the tree. |
---|
287 | see :func:`GSASIIIO.ProjFileOpen` |
---|
288 | ''' |
---|
289 | G2frame = self |
---|
290 | filep = open(fil,'rb') |
---|
291 | shortname = os.path.splitext(os.path.split(fil)[1])[0] |
---|
292 | |
---|
293 | wx.BeginBusyCursor() |
---|
294 | histLoadList = [] |
---|
295 | try: |
---|
296 | while True: |
---|
297 | try: |
---|
298 | data = cPickleLoad(filep) |
---|
299 | except EOFError: |
---|
300 | break |
---|
301 | if not data[0][0].startswith('PWDR'): continue |
---|
302 | if self.PWDRfilter is not None: # implement filter |
---|
303 | if self.PWDRfilter not in data[0][0]: continue |
---|
304 | data[0][0] += ' (' |
---|
305 | data[0][0] += shortname |
---|
306 | data[0][0] += ')' |
---|
307 | histLoadList.append(data) |
---|
308 | |
---|
309 | except Exception as errmsg: |
---|
310 | if GSASIIpath.GetConfigValue('debug'): |
---|
311 | print('\nError reading GPX file:',errmsg) |
---|
312 | import traceback |
---|
313 | print (traceback.format_exc()) |
---|
314 | msg = wx.MessageDialog(G2frame,message="Error reading file "+ |
---|
315 | str(fil)+". This is not a current GSAS-II .gpx file", |
---|
316 | caption="Load Error",style=wx.ICON_ERROR | wx.OK | wx.STAY_ON_TOP) |
---|
317 | msg.ShowModal() |
---|
318 | finally: |
---|
319 | filep.close() |
---|
320 | wx.EndBusyCursor() |
---|
321 | |
---|
322 | datum = None |
---|
323 | for i,data in enumerate(histLoadList): |
---|
324 | datum = data[0] |
---|
325 | datum[0] = G2obj.MakeUniqueLabel(datum[0],self.histList) |
---|
326 | Id = G2frame.GPXtree.AppendItem(parent=G2frame.root,text=datum[0]) |
---|
327 | self.histList.append(datum[0]) |
---|
328 | # if 'ranId' not in datum[1][0]: # patch: add random Id if not present |
---|
329 | # datum[1][0]['ranId'] = ran.randint(0,sys.maxsize) |
---|
330 | G2frame.GPXtree.SetItemPyData(Id,datum[1][:3]) #temp. trim off junk (patch?) |
---|
331 | for datus in data[1:]: |
---|
332 | sub = G2frame.GPXtree.AppendItem(Id,datus[0]) |
---|
333 | #patch |
---|
334 | if datus[0] == 'Instrument Parameters' and len(datus[1]) == 1: |
---|
335 | if datum[0].startswith('PWDR'): |
---|
336 | datus[1] = [dict(zip(datus[1][3],zip(datus[1][0],datus[1][1],datus[1][2]))),{}] |
---|
337 | else: |
---|
338 | datus[1] = [dict(zip(datus[1][2],zip(datus[1][0],datus[1][1]))),{}] |
---|
339 | for item in datus[1][0]: #zip makes tuples - now make lists! |
---|
340 | datus[1][0][item] = list(datus[1][0][item]) |
---|
341 | #end patch |
---|
342 | G2frame.GPXtree.SetItemPyData(sub,datus[1]) |
---|
343 | if datum: # was anything loaded? |
---|
344 | print('project load successful for {}'.format(datum[0])) |
---|
345 | # G2frame.Status.SetStatusText('Mouse RB drag/drop to reorder',0) |
---|
346 | # G2frame.SetTitleByGPX() |
---|
347 | self.GPXtree.Expand(self.root) |
---|
348 | |
---|
349 | def onHistFilter(self,event): |
---|
350 | 'Load a filter string via a dialog in response to a menu event' |
---|
351 | lbl = '' |
---|
352 | if self.PWDRfilter is not None: |
---|
353 | lbl = self.PWDRfilter |
---|
354 | dlg = G2G.SingleStringDialog(self,'Set string', |
---|
355 | 'Set a string that must be in histogram name', |
---|
356 | lbl,size=(400,-1)) |
---|
357 | if dlg.Show(): |
---|
358 | if dlg.GetValue().strip() == '': |
---|
359 | self.PWDRfilter = None |
---|
360 | else: |
---|
361 | self.PWDRfilter = dlg.GetValue() |
---|
362 | dlg.Destroy() |
---|
363 | self.onRefresh(event) |
---|
364 | else: |
---|
365 | dlg.Destroy() |
---|
366 | |
---|
367 | def LoadPhase(self,fil): |
---|
368 | '''Load Phase entries from a .GPX file to the tree. |
---|
369 | see :func:`GSASIIIO.ProjFileOpen` |
---|
370 | ''' |
---|
371 | G2frame = self |
---|
372 | filep = open(fil,'rb') |
---|
373 | shortname = os.path.splitext(os.path.split(fil)[1])[0] |
---|
374 | |
---|
375 | wx.BeginBusyCursor() |
---|
376 | Phases = None |
---|
377 | try: |
---|
378 | while True: |
---|
379 | try: |
---|
380 | data = cPickleLoad(filep) |
---|
381 | except EOFError: |
---|
382 | break |
---|
383 | if not data[0][0].startswith('Phase'): continue |
---|
384 | Phases = data |
---|
385 | #if self.PWDRfilter is not None: # implement filter |
---|
386 | # if self.PWDRfilter not in data[0][0]: continue |
---|
387 | data[0][0] += ' (' |
---|
388 | if Phases: |
---|
389 | data[0][0] += shortname |
---|
390 | data[0][0] += ')' |
---|
391 | else: |
---|
392 | data[0][0] += shortname |
---|
393 | data[0][0] += 'has no phases)' |
---|
394 | Phases = data |
---|
395 | break |
---|
396 | |
---|
397 | except Exception as errmsg: |
---|
398 | if GSASIIpath.GetConfigValue('debug'): |
---|
399 | print('\nError reading GPX file:',errmsg) |
---|
400 | import traceback |
---|
401 | print (traceback.format_exc()) |
---|
402 | msg = wx.MessageDialog(G2frame,message="Error reading file "+ |
---|
403 | str(fil)+". This is not a current GSAS-II .gpx file", |
---|
404 | caption="Load Error",style=wx.ICON_ERROR | wx.OK | wx.STAY_ON_TOP) |
---|
405 | msg.ShowModal() |
---|
406 | finally: |
---|
407 | filep.close() |
---|
408 | wx.EndBusyCursor() |
---|
409 | |
---|
410 | datum = None |
---|
411 | if Phases: |
---|
412 | datum = data[0] |
---|
413 | #datum[0] = G2obj.MakeUniqueLabel(datum[0],self.histList) |
---|
414 | Id = G2frame.GPXtree.AppendItem(parent=G2frame.root,text=datum[0]) |
---|
415 | G2frame.GPXtree.SetItemPyData(Id,datum[1]) |
---|
416 | for datus in data[1:]: |
---|
417 | #datus[0] += ' (' |
---|
418 | #datus[0] += shortname |
---|
419 | #datus[0] += ')' |
---|
420 | sub = G2frame.GPXtree.AppendItem(Id,datus[0]) |
---|
421 | G2frame.GPXtree.SetItemPyData(sub,datus[1]) |
---|
422 | if datum: # was anything loaded? |
---|
423 | self.GPXtree.Expand(Id) |
---|
424 | print('project load successful for {}'.format(datum[0])) |
---|
425 | # G2frame.Status.SetStatusText('Mouse RB drag/drop to reorder',0) |
---|
426 | # G2frame.SetTitleByGPX() |
---|
427 | self.GPXtree.Expand(self.root) |
---|
428 | |
---|
429 | def LoadProject(self,fil): |
---|
430 | '''Load the Covariance entry from a .GPX file to the tree. |
---|
431 | see :func:`GSASIIIO.ProjFileOpen` |
---|
432 | ''' |
---|
433 | G2frame = self |
---|
434 | filep = open(fil,'rb') |
---|
435 | shortname = os.path.splitext(os.path.split(fil)[1])[0] |
---|
436 | |
---|
437 | wx.BeginBusyCursor() |
---|
438 | Phases = None |
---|
439 | try: |
---|
440 | while True: |
---|
441 | try: |
---|
442 | data = cPickleLoad(filep) |
---|
443 | except EOFError: |
---|
444 | break |
---|
445 | if not data[0][0].startswith('Covariance'): continue |
---|
446 | Covar = data[0] |
---|
447 | #GSASIIpath.IPyBreak_base() |
---|
448 | #if self.PWDRfilter is not None: # implement filter |
---|
449 | # if self.PWDRfilter not in data[0][0]: continue |
---|
450 | Covar[0] = shortname + ' Covariance' |
---|
451 | Id = G2frame.GPXtree.AppendItem(parent=G2frame.root,text=Covar[0]) |
---|
452 | G2frame.GPXtree.SetItemPyData(Id,Covar[1]) |
---|
453 | break |
---|
454 | else: |
---|
455 | print("{} does not have refinement results".format(shortname)) |
---|
456 | except Exception as errmsg: |
---|
457 | if GSASIIpath.GetConfigValue('debug'): |
---|
458 | print('\nError reading GPX file:',errmsg) |
---|
459 | import traceback |
---|
460 | print (traceback.format_exc()) |
---|
461 | msg = wx.MessageDialog(G2frame,message="Error reading file "+ |
---|
462 | str(fil)+". This is not a current GSAS-II .gpx file", |
---|
463 | caption="Load Error",style=wx.ICON_ERROR | wx.OK | wx.STAY_ON_TOP) |
---|
464 | msg.ShowModal() |
---|
465 | finally: |
---|
466 | filep.close() |
---|
467 | wx.EndBusyCursor() |
---|
468 | self.GPXtree.Expand(self.root) |
---|
469 | |
---|
470 | def OnDataTreeSelChanged(self,event): |
---|
471 | item = event.GetItem() |
---|
472 | print('selected',item) |
---|
473 | |
---|
474 | #print(self.GPXtree._getTreeItemsList(item)) |
---|
475 | # pltNum = self.G2plotNB.nb.GetSelection() |
---|
476 | # print(pltNum) |
---|
477 | # if pltNum >= 0: #to avoid the startup with no plot! |
---|
478 | # self.G2plotNB.nb.GetPage(pltNum) |
---|
479 | # NewPlot = False |
---|
480 | # else: |
---|
481 | # NewPlot = True |
---|
482 | #if self.getMode() == "Histogram": |
---|
483 | #self.PatternId = self.PickId = item |
---|
484 | #G2plt.PlotPatterns(self,plotType='PWDR',newPlot=NewPlot) |
---|
485 | |
---|
486 | # def OnGPXtreeItemExpanded(self,event): |
---|
487 | # item = event.GetItem() |
---|
488 | # print('expanded',item) |
---|
489 | # print(self.GPXtree._getTreeItemsList(item)) |
---|
490 | # if item == self.root: |
---|
491 | # event.StopPropagation() |
---|
492 | # else: |
---|
493 | # event.Skip(False) |
---|
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 | |
---|
580 | if __name__ == '__main__': |
---|
581 | #if sys.platform == "darwin": |
---|
582 | # application = G2App(0) # create the GUI framework |
---|
583 | #else: |
---|
584 | application = wx.App(0) # create the GUI framework |
---|
585 | try: |
---|
586 | GSASIIpath.SetBinaryPath(True) |
---|
587 | except: |
---|
588 | print('Unable to run with current setup, do you want to update to the') |
---|
589 | try: |
---|
590 | if '2' in platform.python_version_tuple()[0]: |
---|
591 | ans = raw_input("latest GSAS-II version? Update ([Yes]/no): ") |
---|
592 | else: |
---|
593 | ans = input("latest GSAS-II version? Update ([Yes]/no): ") |
---|
594 | except: |
---|
595 | ans = 'no' |
---|
596 | if ans.strip().lower() == "no": |
---|
597 | import sys |
---|
598 | print('Exiting') |
---|
599 | sys.exit() |
---|
600 | print('Updating...') |
---|
601 | GSASIIpath.svnUpdateProcess() |
---|
602 | GSASIIpath.InvokeDebugOpts() |
---|
603 | Frame = main(application) # start the GUI |
---|
604 | argLoadlist = sys.argv[1:] |
---|
605 | for arg in argLoadlist: |
---|
606 | fil = os.path.splitext(arg)[0] + '.gpx' |
---|
607 | if os.path.exists(fil): |
---|
608 | Frame.fileList.append([fil,'GPX']) |
---|
609 | Frame.loadFile(fil) |
---|
610 | else: |
---|
611 | 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 | |
---|
618 | application.MainLoop() |
---|