source: trunk/testDeriv.py @ 3369

Last change on this file since 3369 was 3369, checked in by vondreele, 5 years ago

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