source: trunk/testDeriv.py @ 1111

Last change on this file since 1111 was 1111, checked in by vondreele, 8 years ago

small change to FF calc in GSASIIstrMath.py
add a bit of doc test to testDeriv.py

  • Property svn:keywords set to Date Author Revision URL Id
File size: 7.5 KB
Line 
1# -*- coding: utf-8 -*-
2#testDeriv.py
3'''
4Use this to check derivatives used in structure least squares refinement. To use set DEBUG=True
5in GSASIIstrMain.py (line 22 as of version 1110), run the least squares - one cycle is sufficient.
6Do the "Save Results"; this will write the file testDeriv.dat in the local directory7.
7Run this program to see plots of derivatives for all parameters refined in the last least squares.
8Shown will be numerical derivatives generated over all observations (including penalty terms) and the
9corresponding analytical ones produced in the least squares. They should match.
10-----------------------------------------
11'''
12
13import os
14import os.path as ospath
15import sys
16import time
17import cPickle
18import wx
19import numpy as np
20import matplotlib as mpl
21import GSASIIpath
22import GSASIIstrMath as G2stMth
23import GSASIItestplot as plot
24import pytexture as ptx
25ptx.pyqlmninit()            #initialize fortran arrays for spherical harmonics
26
27def create(parent):
28    return testDeriv(parent)
29   
30[wxID_FILEEXIT, wxID_FILEOPEN, wxID_MAKEPLOTS,
31] = [wx.NewId() for _init_coll_File_Items in range(3)]
32
33def FileDlgFixExt(dlg,file):            #this is needed to fix a problem in linux wx.FileDialog
34    ext = dlg.GetWildcard().split('|')[2*dlg.GetFilterIndex()+1].strip('*')
35    if ext not in file:
36        file += ext
37    return file
38   
39class testDeriv(wx.Frame):
40
41    def _init_ctrls(self, parent):
42        wx.Frame.__init__(self, name='testDeriv', parent=parent,
43            size=wx.Size(460, 250),style=wx.DEFAULT_FRAME_STYLE, title='Test Derivatives')
44        self.testDerivMenu = wx.MenuBar()
45        self.File = wx.Menu(title='')
46        self.File.Append(help='Open testDeriv.dat', id=wxID_FILEOPEN,
47             kind=wx.ITEM_NORMAL,text='Open testDeriv.dat file')
48        self.File.Append(help='Make derivative plots',id=wxID_MAKEPLOTS,
49            kind=wx.ITEM_NORMAL,text='Make plots')
50        self.File.Append(help='Exit from testDeriv', id=wxID_FILEEXIT, kind=wx.ITEM_NORMAL,
51            text='Exit')
52        self.Bind(wx.EVT_MENU, self.OnTestRead, id=wxID_FILEOPEN)
53        self.Bind(wx.EVT_MENU,self.OnMakePlots,id=wxID_MAKEPLOTS)
54        self.Bind(wx.EVT_MENU, self.OnFileExit, id=wxID_FILEEXIT)
55        self.testDerivMenu.Append(menu=self.File, title='File')
56        self.SetMenuBar(self.testDerivMenu)
57        self.testDerivPanel = wx.ScrolledWindow(self)       
58        self.plotNB = plot.PlotNotebook()
59       
60    def __init__(self, parent):
61        self._init_ctrls(parent)
62        self.Bind(wx.EVT_CLOSE, self.ExitMain)   
63        self.dirname = ''
64        self.testfile = []
65        self.dataFrame = None
66
67    def ExitMain(self, event):
68        sys.exit()
69       
70    def OnFileExit(self,event):
71        if self.dataFrame:
72            self.dataFrame.Clear() 
73            self.dataFrame.Destroy()
74        self.Close()
75
76    def OnTestRead(self,event):
77        dlg = wx.FileDialog(self, 'Open testDeriv.dat file', '.', 'testDeriv.dat')
78        if self.dirname:
79            dlg.SetDirectory(self.dirname)
80        try:
81            if dlg.ShowModal() == wx.ID_OK:
82                self.dirname = dlg.GetDirectory()
83                testFile = dlg.GetPath()
84                file = open(testFile,'rb')
85                self.values = cPickle.load(file)
86                self.HistoPhases = cPickle.load(file)
87                self.parmDict = cPickle.load(file)
88                self.varylist = cPickle.load(file)
89                self.calcControls = cPickle.load(file)
90                self.pawleyLookup = cPickle.load(file)
91                self.use = [False for i in range(len(self.varylist))]
92                self.delt = [max(abs(self.parmDict[name])*0.001,1e-6) for name in self.varylist]
93                file.close()
94                self.UpdateControls(event)
95        finally:
96            dlg.Destroy()
97           
98    def UpdateControls(self,event):
99       
100        def OnItemCk(event):
101            Obj = event.GetEventObject()
102            item = ObjInd[Obj.GetId()]
103            self.use[item] = Obj.GetValue()
104           
105        def OnDelValue(event):
106            Obj = event.GetEventObject()
107            item = ObjInd[Obj.GetId()]
108            try:
109                value = float(Obj.GetValue())
110            except ValueError:
111                value = self.delt[item]
112            self.delt[item] = value
113            Obj.SetValue('%g'%(value))
114       
115        self.testDerivPanel.DestroyChildren()
116        ObjInd = {}
117        varylist = self.varylist
118        use = self.use
119        delt = self.delt
120        mainSizer = wx.FlexGridSizer(1,8,5,5)
121        for id,[ck,name,d] in enumerate(zip(use,varylist,delt)):
122            useVal = wx.CheckBox(self.testDerivPanel,label=name)
123            useVal.SetValue(ck)
124            ObjInd[useVal.GetId()] = id
125            useVal.Bind(wx.EVT_CHECKBOX, OnItemCk)
126            mainSizer.Add(useVal,0)
127            delVal = wx.TextCtrl(self.testDerivPanel,wx.ID_ANY,'%g'%(d),style=wx.TE_PROCESS_ENTER)
128            ObjInd[delVal.GetId()] = id
129            delVal.Bind(wx.EVT_TEXT_ENTER,OnDelValue)
130            delVal.Bind(wx.EVT_KILL_FOCUS,OnDelValue)
131            mainSizer.Add(delVal,0)
132        mainSizer.Layout()
133        self.testDerivPanel.SetSizer(mainSizer)   
134        Size = mainSizer.Fit(self.testDerivPanel)
135        Size[0] += 40
136        Size[1] = max(Size[1],290) + 35
137        self.testDerivPanel.SetScrollbars(10,10,Size[0]/10-4,Size[1]/10-1)
138        self.testDerivPanel.SetSize(Size)
139
140    def OnMakePlots(self,event):
141       
142        def test1():
143            fplot = self.plotNB.add('function test').gca()
144            M = G2stMth.errRefine(self.values,self.HistoPhases,
145                self.parmDict,self.varylist,self.calcControls,
146                self.pawleyLookup,None)
147            fplot.plot(M,'r',label='M')
148            fplot.legend(loc='best')
149           
150        def test2(name,delt):
151            hplot = self.plotNB.add('derivatives test for '+name).gca()
152            dMdV = G2stMth.dervRefine(self.values,self.HistoPhases,self.parmDict,
153                self.varylist,self.calcControls,self.pawleyLookup,None)
154            hplot.plot(dMdV[self.varylist.index(name)],'b',label='analytic deriv')
155            if name in self.varylist:
156                self.values[self.varylist.index(name)] -= delt
157                M0 = G2stMth.errRefine(self.values,self.HistoPhases,self.parmDict,
158                    self.varylist,self.calcControls,self.pawleyLookup,None)
159                self.values[self.varylist.index(name)] += 2.*delt
160                M1 = G2stMth.errRefine(self.values,self.HistoPhases,self.parmDict,
161                    self.varylist,self.calcControls,self.pawleyLookup,None)
162                self.values[self.varylist.index(name)] -= delt   
163                Mn = (M1-M0)/(2.*delt)
164                hplot.plot(Mn,'r+',label='numeric deriv')
165                hplot.plot(dMdV[self.varylist.index(name)]-Mn,'g',label='diff')
166            hplot.legend(loc='best')
167           
168        while self.plotNB.nb.GetPageCount():
169            self.plotNB.nb.DeletePage(0)
170        test1()
171        for use,name,delt in zip(self.use,self.varylist,self.delt):
172            if use:
173                test2(name,delt)
174       
175        self.plotNB.Show()
176       
177class testDerivmain(wx.App):
178    def OnInit(self):
179        self.main = testDeriv(None)
180        self.main.Show()
181        self.SetTopWindow(self.main)
182        return True
183
184def main():
185    application = testDerivmain(0)
186    application.MainLoop()
187   
188if __name__ == '__main__':
189    main()
Note: See TracBrowser for help on using the repository browser.