source: trunk/testDeriv.py @ 3711

Last change on this file since 3711 was 3711, checked in by toby, 3 years ago

major constraint update: move conflicting equivs to constraints; allow a formula for multiplier; update docs extensively. New routine EvaluateMultipliers? needed to evaluate formulae

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