source: trunk/testDeriv.py @ 1835

Last change on this file since 1835 was 1835, checked in by vondreele, 6 years ago

fix problem with deleting/inserting data into phases
fix plot problem from Ddata
fix problem with testDeriv & constraints

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 7.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 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#        mainSizer.Layout()
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            hplot = self.plotNB.add('derivatives test for '+name).gca()
166            dMdV = G2stMth.dervRefine(self.values,self.HistoPhases,self.parmDict,
167                self.varylist,self.calcControls,self.pawleyLookup,None)
168            hplot.plot(dMdV[self.varylist.index(name)],'b',label='analytic deriv')
169            if name in self.varylist:
170                print 'parameter:',name,self.values[self.varylist.index(name)],delt
171                self.values[self.varylist.index(name)] -= delt
172                M0 = G2stMth.errRefine(self.values,self.HistoPhases,self.parmDict,
173                    self.varylist,self.calcControls,self.pawleyLookup,None)
174                self.values[self.varylist.index(name)] += 2.*delt
175                M1 = G2stMth.errRefine(self.values,self.HistoPhases,self.parmDict,
176                    self.varylist,self.calcControls,self.pawleyLookup,None)
177                self.values[self.varylist.index(name)] -= delt   
178                Mn = (M1-M0)/(2.*delt)
179                hplot.plot(Mn,'r',label='numeric deriv')
180                hplot.plot(dMdV[self.varylist.index(name)]-Mn,'g',label='diff')
181            hplot.legend(loc='best')
182           
183        while self.plotNB.nb.GetPageCount():
184            self.plotNB.nb.DeletePage(0)
185        test1()
186        for use,name,delt in zip(self.use,self.varylist,self.delt):
187            if use:
188                test2(name,delt)
189       
190        self.plotNB.Show()
191       
192class testDerivmain(wx.App):
193    def OnInit(self):
194        self.main = testDeriv(None)
195        self.main.Show()
196        self.SetTopWindow(self.main)
197        return True
198
199def main():
200    'Starts main application to compute and plot derivatives'
201    application = testDerivmain(0)
202    application.MainLoop()
203   
204if __name__ == '__main__':
205    main()
Note: See TracBrowser for help on using the repository browser.