source: trunk/testDeriv.py @ 3166

Last change on this file since 3166 was 3166, checked in by toby, 8 years ago

Misc Py3 fixes
  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 9.9 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 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.GetMinSize()
150        self.testDerivPanel.SetScrollbars(10,10,Size[0]/10-4,Size[1]/10-1)
151        Size[1] = min(200,Size[1])
152        self.testDerivPanel.SetSize(Size)
153
154    def OnMakePlots(self,event):
155       
156        def test1():
157            fplot = self.plotNB.add('function test').gca()
158            pr = cProfile.Profile()
159            pr.enable()
160            M = G2stMth.errRefine(self.values,self.HistoPhases,
161                self.parmDict,self.varylist,self.calcControls,
162                self.pawleyLookup,None)
163            pr.disable()
164            s = StringIO.StringIO()
165            sortby = 'tottime'
166            ps = pstats.Stats(pr, stream=s).strip_dirs().sort_stats(sortby)
167            print('Profiler of function calculation; top 50% of routines:')
168            ps.print_stats("GSASII",.5)
169            print(s.getvalue())
170            fplot.plot(M,'r',label='M')
171            fplot.legend(loc='best')
172           
173        def test2(name,delt,doProfile):
174            Title = 'derivatives test for '+name
175            varyList = self.varylist+self.depVarList
176            hplot = self.plotNB.add(Title).gca()
177            if doProfile:
178                pr = cProfile.Profile()
179                pr.enable()
180            dMdV = G2stMth.dervRefine(self.values,self.HistoPhases,self.parmDict,
181                varyList,self.calcControls,self.pawleyLookup,None)
182            if doProfile:
183                pr.disable()
184                s = StringIO.StringIO()
185                sortby = 'tottime'
186                ps = pstats.Stats(pr, stream=s).strip_dirs().sort_stats(sortby)
187                ps.print_stats("GSASII",.5)
188                print('Profiler of '+name+' derivative calculation; top 50% of routines:')
189                print(s.getvalue())
190            M2 = dMdV[varyList.index(name)]
191            hplot.plot(M2,'b',label='analytic deriv')
192            mmin = np.min(dMdV[varyList.index(name)])
193            mmax = np.max(dMdV[varyList.index(name)])
194            print('parameter:',name,self.parmDict[name],delt,mmin,mmax)
195            if name in self.varylist:
196                self.values[self.varylist.index(name)] -= delt
197                M0 = G2stMth.errRefine(self.values,self.HistoPhases,self.parmDict,
198                    varyList,self.calcControls,self.pawleyLookup,None)
199                self.values[self.varylist.index(name)] += 2.*delt
200                M1 = G2stMth.errRefine(self.values,self.HistoPhases,self.parmDict,
201                    varyList,self.calcControls,self.pawleyLookup,None)
202                self.values[self.varylist.index(name)] -= delt
203            elif name in self.depVarList:   #in depVarList
204                if 'dA' in name:
205                    name = name.replace('dA','A')
206                    delt *= -1
207                self.parmDict[name] -= delt
208                M0 = G2stMth.errRefine(self.values,self.HistoPhases,self.parmDict,
209                    varyList,self.calcControls,self.pawleyLookup,None)
210                self.parmDict[name] += 2.*delt
211                M1 = G2stMth.errRefine(self.values,self.HistoPhases,self.parmDict,
212                    varyList,self.calcControls,self.pawleyLookup,None)
213                self.parmDict[name] -= delt   
214            Mn = (M1-M0)/(2.*abs(delt))
215            hplot.plot(Mn,'r',label='numeric deriv')
216            hplot.plot(M2-Mn,'g',label='diff')
217#            GSASIIpath.IPyBreak()
218            hplot.legend(loc='best')           
219           
220        while self.plotNB.nb.GetPageCount():
221            self.plotNB.nb.DeletePage(0)
222           
223        test1()
224
225        doProfile = True
226        for use,name,delt in zip(self.use,self.varylist+self.depVarList,self.delt):
227            if use:
228                test2(name,delt,doProfile)
229                doProfile = False
230       
231        self.plotNB.Show()
232       
233class testDerivmain(wx.App):
234    def OnInit(self):
235        self.main = testDeriv(None)
236        self.main.Show()
237        self.SetTopWindow(self.main)
238        return True
239
240def main():
241    'Starts main application to compute and plot derivatives'
242    application = testDerivmain(0)
243    application.MainLoop()
244   
245if __name__ == '__main__':
246    GSASIIpath.InvokeDebugOpts()
247    main()
Note: See TracBrowser for help on using the repository browser.