source: trunk/testDeriv.py @ 3307

Last change on this file since 3307 was 3232, checked in by vondreele, 7 years ago

fix magnetic math in G2strMath
  • 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
29GSASIIpath.SetBinaryPath()
30import GSASIIstrMath as G2stMth
31import GSASIItestplot as plot
32import GSASIImapvars as G2mv
33import pytexture as ptx
34ptx.pyqlmninit()            #initialize fortran arrays for spherical harmonics
35
36def create(parent):
37    return testDeriv(parent)
38   
39[wxID_FILEEXIT, wxID_FILEOPEN, wxID_MAKEPLOTS,
40] = [wx.NewId() for _init_coll_File_Items in range(3)]
41
42def FileDlgFixExt(dlg,file):            #this is needed to fix a problem in linux wx.FileDialog
43    ext = dlg.GetWildcard().split('|')[2*dlg.GetFilterIndex()+1].strip('*')
44    if ext not in file:
45        file += ext
46    return file
47   
48class testDeriv(wx.Frame):
49
50    def _init_ctrls(self, parent):
51        wx.Frame.__init__(self, name='testDeriv', parent=parent,
52            size=wx.Size(800, 250),style=wx.DEFAULT_FRAME_STYLE, title='Test Jacobian Derivatives')
53        self.testDerivMenu = wx.MenuBar()
54        self.File = wx.Menu(title='')
55        self.File.Append(help='Open testDeriv.dat', id=wxID_FILEOPEN,
56             kind=wx.ITEM_NORMAL,text='Open testDeriv.dat file')
57        self.File.Append(help='Make derivative plots',id=wxID_MAKEPLOTS,
58            kind=wx.ITEM_NORMAL,text='Make plots')
59        self.File.Append(help='Exit from testDeriv', id=wxID_FILEEXIT, kind=wx.ITEM_NORMAL,
60            text='Exit')
61        self.Bind(wx.EVT_MENU, self.OnTestRead, id=wxID_FILEOPEN)
62        self.Bind(wx.EVT_MENU,self.OnMakePlots,id=wxID_MAKEPLOTS)
63        self.Bind(wx.EVT_MENU, self.OnFileExit, id=wxID_FILEEXIT)
64        self.testDerivMenu.Append(menu=self.File, title='File')
65        self.SetMenuBar(self.testDerivMenu)
66        self.testDerivPanel = wx.ScrolledWindow(self)       
67        self.plotNB = plot.PlotNotebook()
68       
69    def __init__(self, parent):
70        self._init_ctrls(parent)
71        self.Bind(wx.EVT_CLOSE, self.ExitMain)   
72        self.dirname = ''
73        self.testfile = []
74        self.dataFrame = None
75
76    def ExitMain(self, event):
77        sys.exit()
78       
79    def OnFileExit(self,event):
80        if self.dataFrame:
81            self.dataFrame.Clear() 
82            self.dataFrame.Destroy()
83        self.Close()
84
85    def OnTestRead(self,event):
86        dlg = wx.FileDialog(self, 'Open testDeriv.dat file',defaultFile='testDeriv.dat',
87            wildcard='testDeriv.dat')
88        if self.dirname:
89            dlg.SetDirectory(self.dirname)
90        try:
91            if dlg.ShowModal() == wx.ID_OK:
92                self.dirname = dlg.GetDirectory()
93                testFile = dlg.GetPath()
94                file = open(testFile,'rb')
95                self.values = cPickle.load(file)
96                self.HistoPhases = cPickle.load(file)
97                (self.constrDict,self.fixedList,self.depVarList) = cPickle.load(file)
98                self.parmDict = cPickle.load(file)
99                self.varylist = cPickle.load(file)
100                self.calcControls = cPickle.load(file)
101                self.pawleyLookup = cPickle.load(file)
102                self.use = [False for i in range(len(self.varylist+self.depVarList))]
103                self.delt = [max(abs(self.parmDict[name])*0.0001,1e-6) for name in self.varylist+self.depVarList]
104                file.close()
105                groups,parmlist = G2mv.GroupConstraints(self.constrDict)
106                G2mv.GenerateConstraints(groups,parmlist,self.varylist,self.constrDict,self.fixedList,self.parmDict)
107                self.UpdateControls(event)
108                print(G2mv.VarRemapShow(self.varylist))
109                print('Dependent Vary List:',self.depVarList)
110        finally:
111            dlg.Destroy()
112           
113    def UpdateControls(self,event):
114       
115        def OnItemCk(event):
116            Obj = event.GetEventObject()
117            item = ObjInd[Obj.GetId()]
118            self.use[item] = Obj.GetValue()
119           
120        def OnDelValue(event):
121            event.Skip()
122            Obj = event.GetEventObject()
123            item = ObjInd[Obj.GetId()]
124            try:
125                value = float(Obj.GetValue())
126            except ValueError:
127                value = self.delt[item]
128            self.delt[item] = value
129            Obj.SetValue('%g'%(value))
130       
131        self.testDerivPanel.DestroyChildren()
132        ObjInd = {}
133        varylist = self.varylist
134        depVarList = self.depVarList
135        use = self.use
136        delt = self.delt
137        mainSizer = wx.FlexGridSizer(0,8,5,5)
138        for id,[ck,name,d] in enumerate(zip(use,varylist+depVarList,delt)):
139            useVal = wx.CheckBox(self.testDerivPanel,label=name)
140            useVal.SetValue(ck)
141            ObjInd[useVal.GetId()] = id
142            useVal.Bind(wx.EVT_CHECKBOX, OnItemCk)
143            mainSizer.Add(useVal,0)
144            delVal = wx.TextCtrl(self.testDerivPanel,wx.ID_ANY,'%g'%(d),style=wx.TE_PROCESS_ENTER)
145            ObjInd[delVal.GetId()] = id
146            delVal.Bind(wx.EVT_TEXT_ENTER,OnDelValue)
147            delVal.Bind(wx.EVT_KILL_FOCUS,OnDelValue)
148            mainSizer.Add(delVal,0)
149        self.testDerivPanel.SetSizer(mainSizer)   
150        Size = mainSizer.GetMinSize()
151        self.testDerivPanel.SetScrollbars(10,10,Size[0]/10-4,Size[1]/10-1)
152        Size[1] = min(200,Size[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.