source: trunk/testDeriv.py @ 4213

Last change on this file since 4213 was 4213, checked in by toby, 22 months ago

multiple small changes to allow docs build under 3.x

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