source: trunk/testDeriv.py @ 2548

Last change on this file since 2548 was 2548, checked in by vondreele, 6 years ago

put profiler into testDeriv.py
  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 9.8 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 22, as of version
111110); run the least squares - one cycle 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.
20'''
21
22import sys
23import time
24import cPickle
25import cProfile,pstats,StringIO
26import wx
27import numpy as np
28import GSASIIpath
29import GSASIIstrMath as G2stMth
30import GSASIItestplot as plot
31import GSASIImapvars as G2mv
32import pytexture as ptx
33ptx.pyqlmninit()            #initialize fortran arrays for spherical harmonics
34
35def create(parent):
36    return testDeriv(parent)
37   
38[wxID_FILEEXIT, wxID_FILEOPEN, wxID_MAKEPLOTS,
39] = [wx.NewId() for _init_coll_File_Items in range(3)]
40
41def FileDlgFixExt(dlg,file):            #this is needed to fix a problem in linux wx.FileDialog
42    ext = dlg.GetWildcard().split('|')[2*dlg.GetFilterIndex()+1].strip('*')
43    if ext not in file:
44        file += ext
45    return file
46   
47class testDeriv(wx.Frame):
48
49    def _init_ctrls(self, parent):
50        wx.Frame.__init__(self, name='testDeriv', parent=parent,
51            size=wx.Size(800, 250),style=wx.DEFAULT_FRAME_STYLE, title='Test Jacobian Derivatives')
52        self.testDerivMenu = wx.MenuBar()
53        self.File = wx.Menu(title='')
54        self.File.Append(help='Open testDeriv.dat', id=wxID_FILEOPEN,
55             kind=wx.ITEM_NORMAL,text='Open testDeriv.dat file')
56        self.File.Append(help='Make derivative plots',id=wxID_MAKEPLOTS,
57            kind=wx.ITEM_NORMAL,text='Make plots')
58        self.File.Append(help='Exit from testDeriv', id=wxID_FILEEXIT, kind=wx.ITEM_NORMAL,
59            text='Exit')
60        self.Bind(wx.EVT_MENU, self.OnTestRead, id=wxID_FILEOPEN)
61        self.Bind(wx.EVT_MENU,self.OnMakePlots,id=wxID_MAKEPLOTS)
62        self.Bind(wx.EVT_MENU, self.OnFileExit, id=wxID_FILEEXIT)
63        self.testDerivMenu.Append(menu=self.File, title='File')
64        self.SetMenuBar(self.testDerivMenu)
65        self.testDerivPanel = wx.ScrolledWindow(self)       
66        self.plotNB = plot.PlotNotebook()
67       
68    def __init__(self, parent):
69        self._init_ctrls(parent)
70        self.Bind(wx.EVT_CLOSE, self.ExitMain)   
71        self.dirname = ''
72        self.testfile = []
73        self.dataFrame = None
74
75    def ExitMain(self, event):
76        sys.exit()
77       
78    def OnFileExit(self,event):
79        if self.dataFrame:
80            self.dataFrame.Clear() 
81            self.dataFrame.Destroy()
82        self.Close()
83
84    def OnTestRead(self,event):
85        dlg = wx.FileDialog(self, 'Open testDeriv.dat file',defaultFile='testDeriv.dat',
86            wildcard='testDeriv.dat')
87        if self.dirname:
88            dlg.SetDirectory(self.dirname)
89        try:
90            if dlg.ShowModal() == wx.ID_OK:
91                self.dirname = dlg.GetDirectory()
92                testFile = dlg.GetPath()
93                file = open(testFile,'rb')
94                self.values = cPickle.load(file)
95                self.HistoPhases = cPickle.load(file)
96                (self.constrDict,self.fixedList,self.depVarList) = cPickle.load(file)
97                self.parmDict = cPickle.load(file)
98                self.varylist = cPickle.load(file)
99                self.calcControls = cPickle.load(file)
100                self.pawleyLookup = cPickle.load(file)
101                self.use = [False for i in range(len(self.varylist+self.depVarList))]
102                self.delt = [max(abs(self.parmDict[name])*0.0001,1e-6) for name in self.varylist+self.depVarList]
103                file.close()
104                groups,parmlist = G2mv.GroupConstraints(self.constrDict)
105                G2mv.GenerateConstraints(groups,parmlist,self.varylist,self.constrDict,self.fixedList,self.parmDict)
106                self.UpdateControls(event)
107                print G2mv.VarRemapShow(self.varylist)
108                print 'Dependent Vary List:',self.depVarList
109        finally:
110            dlg.Destroy()
111           
112    def UpdateControls(self,event):
113       
114        def OnItemCk(event):
115            Obj = event.GetEventObject()
116            item = ObjInd[Obj.GetId()]
117            self.use[item] = Obj.GetValue()
118           
119        def OnDelValue(event):
120            event.Skip()
121            Obj = event.GetEventObject()
122            item = ObjInd[Obj.GetId()]
123            try:
124                value = float(Obj.GetValue())
125            except ValueError:
126                value = self.delt[item]
127            self.delt[item] = value
128            Obj.SetValue('%g'%(value))
129       
130        self.testDerivPanel.DestroyChildren()
131        ObjInd = {}
132        varylist = self.varylist
133        depVarList = self.depVarList
134        use = self.use
135        delt = self.delt
136        mainSizer = wx.FlexGridSizer(0,8,5,5)
137        for id,[ck,name,d] in enumerate(zip(use,varylist+depVarList,delt)):
138            useVal = wx.CheckBox(self.testDerivPanel,label=name)
139            useVal.SetValue(ck)
140            ObjInd[useVal.GetId()] = id
141            useVal.Bind(wx.EVT_CHECKBOX, OnItemCk)
142            mainSizer.Add(useVal,0)
143            delVal = wx.TextCtrl(self.testDerivPanel,wx.ID_ANY,'%g'%(d),style=wx.TE_PROCESS_ENTER)
144            ObjInd[delVal.GetId()] = id
145            delVal.Bind(wx.EVT_TEXT_ENTER,OnDelValue)
146            delVal.Bind(wx.EVT_KILL_FOCUS,OnDelValue)
147            mainSizer.Add(delVal,0)
148        self.testDerivPanel.SetSizer(mainSizer)   
149        Size = mainSizer.Fit(self.testDerivPanel)
150        Size[0] = 800
151        Size[1] = max(Size[1],290) + 35
152        self.testDerivPanel.SetScrollbars(10,10,Size[0]/10-4,Size[1]/10-1)
153        self.testDerivPanel.SetSize(Size)
154
155    def OnMakePlots(self,event):
156       
157        def test1():
158            fplot = self.plotNB.add('function test').gca()
159            pr = cProfile.Profile()
160            pr.enable()
161            M = G2stMth.errRefine(self.values,self.HistoPhases,
162                self.parmDict,self.varylist,self.calcControls,
163                self.pawleyLookup,None)
164            pr.disable()
165            s = StringIO.StringIO()
166            sortby = 'tottime'
167            ps = pstats.Stats(pr, stream=s).strip_dirs().sort_stats(sortby)
168            print 'Profiler of function calculation; top 50% of routines:'
169            ps.print_stats("GSASII",.5)
170            print s.getvalue()
171            fplot.plot(M,'r',label='M')
172            fplot.legend(loc='best')
173           
174        def test2(name,delt,doProfile):
175            Title = 'derivatives test for '+name
176            varyList = self.varylist+self.depVarList
177            hplot = self.plotNB.add(Title).gca()
178            if doProfile:
179                pr = cProfile.Profile()
180                pr.enable()
181            dMdV = G2stMth.dervRefine(self.values,self.HistoPhases,self.parmDict,
182                varyList,self.calcControls,self.pawleyLookup,None)
183            if doProfile:
184                pr.disable()
185                s = StringIO.StringIO()
186                sortby = 'tottime'
187                ps = pstats.Stats(pr, stream=s).strip_dirs().sort_stats(sortby)
188                ps.print_stats("GSASII",.5)
189                print 'Profiler of '+name+' derivative calculation; top 50% of routines:'
190                print s.getvalue()
191            M2 = dMdV[varyList.index(name)]
192            hplot.plot(M2,'b',label='analytic deriv')
193            mmin = np.min(dMdV[varyList.index(name)])
194            mmax = np.max(dMdV[varyList.index(name)])
195            print 'parameter:',name,self.parmDict[name],delt,mmin,mmax
196            if name in self.varylist:
197                self.values[self.varylist.index(name)] -= delt
198                M0 = G2stMth.errRefine(self.values,self.HistoPhases,self.parmDict,
199                    varyList,self.calcControls,self.pawleyLookup,None)
200                self.values[self.varylist.index(name)] += 2.*delt
201                M1 = G2stMth.errRefine(self.values,self.HistoPhases,self.parmDict,
202                    varyList,self.calcControls,self.pawleyLookup,None)
203                self.values[self.varylist.index(name)] -= delt
204            elif name in self.depVarList:   #in depVarList
205                if 'dA' in name:
206                    name = name.replace('dA','A')
207                    delt *= -1
208                self.parmDict[name] -= delt
209                M0 = G2stMth.errRefine(self.values,self.HistoPhases,self.parmDict,
210                    varyList,self.calcControls,self.pawleyLookup,None)
211                self.parmDict[name] += 2.*delt
212                M1 = G2stMth.errRefine(self.values,self.HistoPhases,self.parmDict,
213                    varyList,self.calcControls,self.pawleyLookup,None)
214                self.parmDict[name] -= delt   
215            Mn = (M1-M0)/(2.*abs(delt))
216            hplot.plot(Mn,'r',label='numeric deriv')
217            hplot.plot(M2-Mn,'g',label='diff')
218#            GSASIIpath.IPyBreak()
219            hplot.legend(loc='best')           
220           
221        while self.plotNB.nb.GetPageCount():
222            self.plotNB.nb.DeletePage(0)
223           
224        test1()
225
226        doProfile = True
227        for use,name,delt in zip(self.use,self.varylist+self.depVarList,self.delt):
228            if use:
229                test2(name,delt,doProfile)
230                doProfile = False
231       
232        self.plotNB.Show()
233       
234class testDerivmain(wx.App):
235    def OnInit(self):
236        self.main = testDeriv(None)
237        self.main.Show()
238        self.SetTopWindow(self.main)
239        return True
240
241def main():
242    'Starts main application to compute and plot derivatives'
243    application = testDerivmain(0)
244    application.MainLoop()
245   
246if __name__ == '__main__':
247    GSASIIpath.InvokeDebugOpts()
248    main()
Note: See TracBrowser for help on using the repository browser.