source: trunk/testDeriv.py @ 1456

Last change on this file since 1456 was 1456, checked in by vondreele, 7 years ago

consolidate HKLF parts of dervRefine & HessRefine? into new routine dervHKLF

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