source: trunk/testDeriv.py @ 5114

Last change on this file since 5114 was 5114, checked in by vondreele, 9 months ago

various fixes for EDX data fitting by LeBail?
various fixes for SAD analysis - plotting updates & importer
removed a couple of unused lambdas from G2sasd

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 11.4 KB
Line 
1# -*- coding: utf-8 -*-
2#testDeriv.py
3'''
4*testDeriv: Check derivative computation*
5=========================================
6
7Use this to check derivatives used in structure least squares
8refinement against numerical values computed in this script.
9
10To use set ``DEBUG=True`` in GSASIIstrMain.py (line 40, as of version
112546); run the least squares - zero cycles is sufficient.  Do the "Save
12Results"; this will write the file testDeriv.dat in the local
13directory.
14
15Then run this program to see plots of derivatives for all
16parameters refined in the last least squares.  Shown will be numerical
17derivatives generated over all observations (including penalty terms)
18and the corresponding analytical ones produced in the least
19squares. They should match. Profiling is also done for function
20calculation & for the 1st selected derivative (rest should be the same).
21'''
22
23import sys
24import os
25import platform
26if '2' in platform.python_version_tuple()[0]:
27    import cPickle
28    import StringIO
29else:
30    import pickle as cPickle
31    import io as StringIO
32import cProfile,pstats
33import wx
34import numpy as np
35import GSASIIpath
36GSASIIpath.SetBinaryPath()
37import GSASIIstrMath as G2stMth
38import GSASIItestplot as plot
39import GSASIImapvars as G2mv
40try:  # fails on doc build
41    import pytexture as ptx
42    ptx.pyqlmninit()            #initialize fortran arrays for spherical harmonics
43except ImportError:
44    pass
45
46try:
47    wx.NewId
48except AttributeError:
49    wx.NewId = wx.NewIdRef
50
51def create(parent):
52    return testDeriv(parent)
53   
54[wxID_FILEEXIT, wxID_FILEOPEN, wxID_MAKEPLOTS, wxID_CLEARSEL,wxID_SELECTALL,
55] = [wx.NewId() for _init_coll_File_Items in range(5)]
56
57def FileDlgFixExt(dlg,file):            #this is needed to fix a problem in linux wx.FileDialog
58    ext = dlg.GetWildcard().split('|')[2*dlg.GetFilterIndex()+1].strip('*')
59    if ext not in file:
60        file += ext
61    return file
62   
63class testDeriv(wx.Frame):
64
65    def _init_ctrls(self, parent):
66        wx.Frame.__init__(self, name='testDeriv', parent=parent,
67            size=wx.Size(800, 250),style=wx.DEFAULT_FRAME_STYLE, title='Test Jacobian Derivatives')
68        self.testDerivMenu = wx.MenuBar()
69        self.File = wx.Menu(title='')
70        self.File.Append(wxID_FILEOPEN,'Open testDeriv file','Open testDeriv')
71        self.File.Append(wxID_MAKEPLOTS,'Make plots','Make derivative plots')
72        self.File.Append(wxID_SELECTALL,'Select all')
73        self.File.Append(wxID_CLEARSEL,'Clear selections')
74        self.File.Append(wxID_FILEEXIT,'Exit','Exit from testDeriv')
75        self.Bind(wx.EVT_MENU,self.OnTestRead, id=wxID_FILEOPEN)
76        self.Bind(wx.EVT_MENU,self.OnMakePlots,id=wxID_MAKEPLOTS)
77        self.Bind(wx.EVT_MENU,self.ClearSelect,id=wxID_CLEARSEL)
78        self.Bind(wx.EVT_MENU,self.SelectAll,id=wxID_SELECTALL)
79        self.Bind(wx.EVT_MENU,self.OnFileExit, id=wxID_FILEEXIT)
80        self.testDerivMenu.Append(menu=self.File, title='File')
81        self.SetMenuBar(self.testDerivMenu)
82        self.testDerivPanel = wx.ScrolledWindow(self)       
83        self.plotNB = plot.PlotNotebook()
84        self.testFile = ''
85        arg = sys.argv
86        if len(arg) > 1 and arg[1]:
87            try:
88                self.testFile = os.path.splitext(arg[1])[0]+u'.testDeriv'
89            except:
90                self.testFile = os.path.splitext(arg[1])[0]+'.testDeriv'
91            self.TestRead()
92            self.UpdateControls(None)
93       
94    def __init__(self, parent):
95        self._init_ctrls(parent)
96        self.Bind(wx.EVT_CLOSE, self.ExitMain)   
97        self.dirname = ''
98        self.testfile = []
99        self.dataFrame = None
100
101    def ExitMain(self, event):
102        sys.exit()
103       
104    def OnFileExit(self,event):
105        if self.dataFrame:
106            self.dataFrame.Clear() 
107            self.dataFrame.Destroy()
108        self.Close()
109       
110    def SelectAll(self,event):
111        self.use = [True for name in self.names]
112        for i,name in enumerate(self.names):
113            if 'Back' in name:
114                self.use[i] = False
115        self.UpdateControls(event)
116       
117    def ClearSelect(self,event):
118        self.use = [False for i in range(len(self.names))]
119        self.UpdateControls(event)
120
121    def OnTestRead(self,event):
122        dlg = wx.FileDialog(self, 'Open *.testDeriv file',defaultFile='*.testDeriv',
123            wildcard='*.testDeriv')
124        if self.dirname:
125            dlg.SetDirectory(self.dirname)
126        try:
127            if dlg.ShowModal() == wx.ID_OK:
128                self.dirname = dlg.GetDirectory()
129                self.testFile = dlg.GetPath()
130                self.TestRead()
131                self.UpdateControls(event)
132        finally:
133            dlg.Destroy()
134           
135    def TestRead(self):
136        file = open(self.testFile,'rb')
137        self.values = cPickle.load(file,encoding='Latin-1')
138        self.HistoPhases = cPickle.load(file,encoding='Latin-1')
139        (self.constrDict,self.fixedList,self.depVarList) = cPickle.load(file,encoding='Latin-1')
140        self.parmDict = cPickle.load(file,encoding='Latin-1')
141        self.varylist = cPickle.load(file,encoding='Latin-1')
142        self.calcControls = cPickle.load(file,encoding='Latin-1')
143        self.pawleyLookup = cPickle.load(file,encoding='Latin-1')
144        self.names = self.varylist+self.depVarList
145        self.use = [False for i in range(len(self.names))]
146        self.delt = [max(abs(self.parmDict[name])*0.0001,1e-6) for name in self.names]
147        for iname,name in enumerate(self.names):
148            if name.split(':')[-1] in ['Shift','DisplaceX','DisplaceY',]:
149                self.delt[iname] = 0.1
150        file.close()
151        msg = G2mv.EvaluateMultipliers(self.constrDict,self.parmDict)
152        if msg:
153            print('Unable to interpret multiplier(s): '+msg)
154            raise Exception
155        G2mv.GenerateConstraints(self.varylist,self.constrDict,self.fixedList,self.parmDict)
156        print(G2mv.VarRemapShow(self.varylist))
157        print('Dependent Vary List:',self.depVarList)
158           
159    def UpdateControls(self,event):
160       
161        def OnItemCk(event):
162            Obj = event.GetEventObject()
163            item = ObjInd[Obj.GetId()]
164            self.use[item] = Obj.GetValue()
165           
166        def OnDelValue(event):
167            event.Skip()
168            Obj = event.GetEventObject()
169            item = ObjInd[Obj.GetId()]
170            try:
171                value = float(Obj.GetValue())
172            except ValueError:
173                value = self.delt[item]
174            self.delt[item] = value
175            Obj.SetValue('%g'%(value))
176       
177        self.testDerivPanel.DestroyChildren()
178        ObjInd = {}
179        names = self.names
180        use = self.use
181        delt = self.delt
182        mainSizer = wx.FlexGridSizer(0,8,5,5)
183        for id,[ck,name,d] in enumerate(zip(use,names,delt)):
184            useVal = wx.CheckBox(self.testDerivPanel,label=name)
185            useVal.SetValue(ck)
186            ObjInd[useVal.GetId()] = id
187            useVal.Bind(wx.EVT_CHECKBOX, OnItemCk)
188            mainSizer.Add(useVal,0)
189            delVal = wx.TextCtrl(self.testDerivPanel,wx.ID_ANY,'%g'%(d),style=wx.TE_PROCESS_ENTER)
190            ObjInd[delVal.GetId()] = id
191            delVal.Bind(wx.EVT_TEXT_ENTER,OnDelValue)
192            delVal.Bind(wx.EVT_KILL_FOCUS,OnDelValue)
193            mainSizer.Add(delVal,0)
194        self.testDerivPanel.SetSizer(mainSizer)   
195        Size = mainSizer.GetMinSize()
196        self.testDerivPanel.SetScrollbars(10,10,Size[0]/10-4,Size[1]/10-1)
197        Size[1] = min(200,Size[1])
198        self.testDerivPanel.SetSize(Size)
199
200    def OnMakePlots(self,event):
201       
202        def test1():
203            fplot = self.plotNB.add('function test').gca()
204            pr = cProfile.Profile()
205            pr.enable()
206            M = G2stMth.errRefine(self.values,self.HistoPhases,
207                self.parmDict,self.varylist,self.calcControls,
208                self.pawleyLookup,None)
209            pr.disable()
210            s = StringIO.StringIO()
211            sortby = 'tottime'
212            ps = pstats.Stats(pr, stream=s).strip_dirs().sort_stats(sortby)
213            print('Profiler of function calculation; top 50% of routines:')
214            ps.print_stats("GSASII",.5)
215            print(s.getvalue())
216            fplot.plot(M,'r',label='M')
217            fplot.legend(loc='best')
218           
219        def test2(name,delt,doProfile):
220            Title = 'derivatives test for '+name
221            names = self.names
222            hplot = self.plotNB.add(Title).gca()
223            if doProfile:
224                pr = cProfile.Profile()
225                pr.enable()
226            #regenerate minimization fxn
227            G2stMth.errRefine(self.values,self.HistoPhases,
228                self.parmDict,self.varylist,self.calcControls,
229                self.pawleyLookup,None)
230            dMdV = G2stMth.dervRefine(self.values,self.HistoPhases,self.parmDict,
231                names,self.calcControls,self.pawleyLookup,None)
232            if doProfile:
233                pr.disable()
234                s = StringIO.StringIO()
235                sortby = 'tottime'
236                ps = pstats.Stats(pr, stream=s).strip_dirs().sort_stats(sortby)
237                ps.print_stats("GSASII",.5)
238                print('Profiler of '+name+' derivative calculation; top 50% of routines:')
239                print(s.getvalue())
240            M2 = dMdV[names.index(name)]
241            hplot.plot(M2,'b',label='analytic deriv')
242            mmin = np.min(dMdV[names.index(name)])
243            mmax = np.max(dMdV[names.index(name)])
244            if name in self.varylist:
245                self.values[self.varylist.index(name)] -= delt
246                M0 = G2stMth.errRefine(self.values,self.HistoPhases,self.parmDict,
247                    names,self.calcControls,self.pawleyLookup,None)
248                self.values[self.varylist.index(name)] += 2.*delt
249                M1 = G2stMth.errRefine(self.values,self.HistoPhases,self.parmDict,
250                    names,self.calcControls,self.pawleyLookup,None)
251                self.values[self.varylist.index(name)] -= delt
252            elif name in self.depVarList:   #in depVarList
253                if 'dA' in name:
254                    name = name.replace('dA','A')
255                    delt *= -1
256                self.parmDict[name] -= delt
257                M0 = G2stMth.errRefine(self.values,self.HistoPhases,self.parmDict,
258                    names,self.calcControls,self.pawleyLookup,None)
259                self.parmDict[name] += 2.*delt
260                M1 = G2stMth.errRefine(self.values,self.HistoPhases,self.parmDict,
261                    names,self.calcControls,self.pawleyLookup,None)
262                self.parmDict[name] -= delt   
263            Mn = (M1-M0)/(2.*abs(delt))
264            print('parameter:',name,self.parmDict[name],delt,mmin,mmax,np.sum(M0),np.sum(M1),np.sum(Mn))
265            hplot.plot(Mn,'r',label='numeric deriv')
266            hplot.legend(loc='best')           
267           
268        while self.plotNB.nb.GetPageCount():
269            self.plotNB.nb.DeletePage(0)
270           
271        test1()
272
273        doProfile = True
274        for use,name,delt in zip(self.use,self.varylist+self.depVarList,self.delt):
275            if use:
276                test2(name,delt,doProfile)
277                doProfile = False
278       
279        self.plotNB.Show()
280       
281class testDerivmain(wx.App):
282    def OnInit(self):
283        self.main = testDeriv(None)
284        self.main.Show()
285        self.SetTopWindow(self.main)
286        return True
287
288def main():
289    'Starts main application to compute and plot derivatives'
290    application = testDerivmain(0)
291    application.MainLoop()
292   
293if __name__ == '__main__':
294    GSASIIpath.InvokeDebugOpts()
295    main()
Note: See TracBrowser for help on using the repository browser.