source: trunk/testDeriv.py @ 2549

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

fix special position restrictions on Mx,MY,Mz on -1 site.

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