source: trunk/testDeriv.py @ 3415

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

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