source: trunk/testDeriv.py @ 3897

Last change on this file since 3897 was 3828, checked in by toby, 6 years ago

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