source: trunk/testDeriv.py @ 2500

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

protect a dlg in DoIndexPeaks?
fix a lost G2frame.LimitsTable? - need to go back to tree instead
add depVarList to testDeriv file
more work on mag derivs.
modify testDeriv to show depVar derivatives

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