source: trunk/testDeriv.py @ 3774

Last change on this file since 3774 was 3774, checked in by vondreele, 3 years ago

fix super indexing problem in transposeHKLF
fix reflection generation for incommensurate mag case in G2lattice & G2pwd
clean up non Fourier modulation calcs & remove analytic derivative stuff (now numeric)
fix uij derivative bug
work on incommensurate magnetic sturcture factors - not working yet
clean up of testDeriv - better choices for delt & force reflection regeneration before each derivative test

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