source: trunk/testDeriv.py @ 2539

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

fix to unicode character in directory name problem
numerous unused parameters removed - more to come

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 8.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 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 time
24import cPickle
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            M = G2stMth.errRefine(self.values,self.HistoPhases,
159                self.parmDict,self.varylist,self.calcControls,
160                self.pawleyLookup,None)
161            fplot.plot(M,'r',label='M')
162            fplot.legend(loc='best')
163           
164        def test2(name,delt):
165            Title = 'derivatives test for '+name
166            varyList = self.varylist+self.depVarList
167            hplot = self.plotNB.add(Title).gca()
168            dMdV = G2stMth.dervRefine(self.values,self.HistoPhases,self.parmDict,
169                varyList,self.calcControls,self.pawleyLookup,None)
170            M2 = dMdV[varyList.index(name)]
171            hplot.plot(M2,'b',label='analytic deriv')
172            mmin = np.min(dMdV[varyList.index(name)])
173            mmax = np.max(dMdV[varyList.index(name)])
174            print 'parameter:',name,self.parmDict[name],delt,mmin,mmax
175            if name in self.varylist:
176                self.values[self.varylist.index(name)] -= delt
177                M0 = G2stMth.errRefine(self.values,self.HistoPhases,self.parmDict,
178                    varyList,self.calcControls,self.pawleyLookup,None)
179                self.values[self.varylist.index(name)] += 2.*delt
180                M1 = G2stMth.errRefine(self.values,self.HistoPhases,self.parmDict,
181                    varyList,self.calcControls,self.pawleyLookup,None)
182                self.values[self.varylist.index(name)] -= delt
183            elif name in self.depVarList:   #in depVarList
184                if 'dA' in name:
185                    name = name.replace('dA','A')
186                    delt *= -1
187                self.parmDict[name] -= delt
188                M0 = G2stMth.errRefine(self.values,self.HistoPhases,self.parmDict,
189                    varyList,self.calcControls,self.pawleyLookup,None)
190                self.parmDict[name] += 2.*delt
191                M1 = G2stMth.errRefine(self.values,self.HistoPhases,self.parmDict,
192                    varyList,self.calcControls,self.pawleyLookup,None)
193                self.parmDict[name] -= delt   
194            Mn = (M1-M0)/(2.*abs(delt))
195            hplot.plot(Mn,'r',label='numeric deriv')
196            hplot.plot(M2-Mn,'g',label='diff')
197#            GSASIIpath.IPyBreak()
198            hplot.legend(loc='best')           
199           
200        while self.plotNB.nb.GetPageCount():
201            self.plotNB.nb.DeletePage(0)
202        test1()
203        for use,name,delt in zip(self.use,self.varylist+self.depVarList,self.delt):
204            if use:
205                test2(name,delt)
206       
207        self.plotNB.Show()
208       
209class testDerivmain(wx.App):
210    def OnInit(self):
211        self.main = testDeriv(None)
212        self.main.Show()
213        self.SetTopWindow(self.main)
214        return True
215
216def main():
217    'Starts main application to compute and plot derivatives'
218    application = testDerivmain(0)
219    application.MainLoop()
220   
221if __name__ == '__main__':
222    GSASIIpath.InvokeDebugOpts()
223    main()
Note: See TracBrowser for help on using the repository browser.