source: trunk/testDeriv.py @ 1860

Last change on this file since 1860 was 1860, checked in by vondreele, 6 years ago
  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 8.0 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) = 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))]
103                self.delt = [max(abs(self.parmDict[name])*0.001,1e-6) for name in self.varylist]
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        finally:
110            dlg.Destroy()
111           
112    def UpdateControls(self,event):
113       
114        def OnItemCk(event):
115            Obj = event.GetEventObject()
116            item = ObjInd[Obj.GetId()]
117            self.use[item] = Obj.GetValue()
118           
119        def OnDelValue(event):
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        use = self.use
133        delt = self.delt
134        mainSizer = wx.FlexGridSizer(0,8,5,5)
135        for id,[ck,name,d] in enumerate(zip(use,varylist,delt)):
136            useVal = wx.CheckBox(self.testDerivPanel,label=name)
137            useVal.SetValue(ck)
138            ObjInd[useVal.GetId()] = id
139            useVal.Bind(wx.EVT_CHECKBOX, OnItemCk)
140            mainSizer.Add(useVal,0)
141            delVal = wx.TextCtrl(self.testDerivPanel,wx.ID_ANY,'%g'%(d),style=wx.TE_PROCESS_ENTER)
142            ObjInd[delVal.GetId()] = id
143            delVal.Bind(wx.EVT_TEXT_ENTER,OnDelValue)
144            delVal.Bind(wx.EVT_KILL_FOCUS,OnDelValue)
145            mainSizer.Add(delVal,0)
146        self.testDerivPanel.SetSizer(mainSizer)   
147        Size = mainSizer.Fit(self.testDerivPanel)
148        Size[0] = 800
149        Size[1] = max(Size[1],290) + 35
150        self.testDerivPanel.SetScrollbars(10,10,Size[0]/10-4,Size[1]/10-1)
151        self.testDerivPanel.SetSize(Size)
152
153    def OnMakePlots(self,event):
154       
155        def test1():
156            fplot = self.plotNB.add('function test').gca()
157            M = G2stMth.errRefine(self.values,self.HistoPhases,
158                self.parmDict,self.varylist,self.calcControls,
159                self.pawleyLookup,None)
160            fplot.plot(M,'r',label='M')
161            fplot.legend(loc='best')
162           
163        def test2(name,delt):
164           
165            Title = 'derivatives test for '+name
166            hplot = self.plotNB.add(Title).gca()
167            dMdV = G2stMth.dervRefine(self.values,self.HistoPhases,self.parmDict,
168                self.varylist,self.calcControls,self.pawleyLookup,None)
169            hplot.plot(dMdV[self.varylist.index(name)],'b',label='analytic deriv')
170            if name in self.varylist:
171                print 'parameter:',name,self.values[self.varylist.index(name)],delt
172                self.values[self.varylist.index(name)] -= delt
173                M0 = G2stMth.errRefine(self.values,self.HistoPhases,self.parmDict,
174                    self.varylist,self.calcControls,self.pawleyLookup,None)
175                self.values[self.varylist.index(name)] += 2.*delt
176                M1 = G2stMth.errRefine(self.values,self.HistoPhases,self.parmDict,
177                    self.varylist,self.calcControls,self.pawleyLookup,None)
178                self.values[self.varylist.index(name)] -= delt   
179                Mn = (M1-M0)/(2.*delt)
180                hplot.plot(Mn,'r',label='numeric deriv')
181                hplot.plot(dMdV[self.varylist.index(name)]-Mn,'g',label='diff')
182            hplot.legend(loc='best')           
183           
184        while self.plotNB.nb.GetPageCount():
185            self.plotNB.nb.DeletePage(0)
186        test1()
187        for use,name,delt in zip(self.use,self.varylist,self.delt):
188            if use:
189                test2(name,delt)
190       
191        self.plotNB.Show()
192       
193class testDerivmain(wx.App):
194    def OnInit(self):
195        self.main = testDeriv(None)
196        self.main.Show()
197        self.SetTopWindow(self.main)
198        return True
199
200def main():
201    'Starts main application to compute and plot derivatives'
202    application = testDerivmain(0)
203    application.MainLoop()
204   
205if __name__ == '__main__':
206    main()
Note: See TracBrowser for help on using the repository browser.