source: trunk/testDeriv.py @ 4873

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

complete pink Rietveld refinement
fix problem with weights
enhance testDeriv
put new pink parms in parm dictionary
fix problem with Instrument parm display after load
  • 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
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,wxID_SELECTALL,
55] = [wx.NewId() for _init_coll_File_Items in range(5)]
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_SELECTALL,'Select all')
73        self.File.Append(wxID_CLEARSEL,'Clear selections')
74        self.File.Append(wxID_FILEEXIT,'Exit','Exit from testDeriv')
75        self.Bind(wx.EVT_MENU,self.OnTestRead, id=wxID_FILEOPEN)
76        self.Bind(wx.EVT_MENU,self.OnMakePlots,id=wxID_MAKEPLOTS)
77        self.Bind(wx.EVT_MENU,self.ClearSelect,id=wxID_CLEARSEL)
78        self.Bind(wx.EVT_MENU,self.SelectAll,id=wxID_SELECTALL)
79        self.Bind(wx.EVT_MENU,self.OnFileExit, id=wxID_FILEEXIT)
80        self.testDerivMenu.Append(menu=self.File, title='File')
81        self.SetMenuBar(self.testDerivMenu)
82        self.testDerivPanel = wx.ScrolledWindow(self)       
83        self.plotNB = plot.PlotNotebook()
84        self.testFile = ''
85        arg = sys.argv
86        if len(arg) > 1 and arg[1]:
87            try:
88                self.testFile = os.path.splitext(arg[1])[0]+u'.testDeriv'
89            except:
90                self.testFile = os.path.splitext(arg[1])[0]+'.testDeriv'
91            self.TestRead()
92            self.UpdateControls(None)
93       
94    def __init__(self, parent):
95        self._init_ctrls(parent)
96        self.Bind(wx.EVT_CLOSE, self.ExitMain)   
97        self.dirname = ''
98        self.testfile = []
99        self.dataFrame = None
100
101    def ExitMain(self, event):
102        sys.exit()
103       
104    def OnFileExit(self,event):
105        if self.dataFrame:
106            self.dataFrame.Clear() 
107            self.dataFrame.Destroy()
108        self.Close()
109       
110    def SelectAll(self,event):
111        self.use = [True for name in self.names]
112        for i,name in enumerate(self.names):
113            if 'Back' in name:
114                self.use[i] = False
115        self.UpdateControls(event)
116       
117    def ClearSelect(self,event):
118        self.use = [False for i in range(len(self.names))]
119        self.UpdateControls(event)
120
121    def OnTestRead(self,event):
122        dlg = wx.FileDialog(self, 'Open *.testDeriv file',defaultFile='*.testDeriv',
123            wildcard='*.testDeriv')
124        if self.dirname:
125            dlg.SetDirectory(self.dirname)
126        try:
127            if dlg.ShowModal() == wx.ID_OK:
128                self.dirname = dlg.GetDirectory()
129                self.testFile = dlg.GetPath()
130                self.TestRead()
131                self.UpdateControls(event)
132        finally:
133            dlg.Destroy()
134           
135    def TestRead(self):
136        file = open(self.testFile,'rb')
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.