source: trunk/testDeriv.py @ 2508

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

fix move difference curve problem
mag derivs are better but still not right - refinements work (sort of)
fix to testDeriv to show atom pos derivs for depVarList parameters

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 9.1 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',defaultFile='testDeriv.dat',
88            wildcard='testDeriv.dat')
89        if self.dirname:
90            dlg.SetDirectory(self.dirname)
91        try:
92            if dlg.ShowModal() == wx.ID_OK:
93                self.dirname = dlg.GetDirectory()
94                testFile = dlg.GetPath()
95                file = open(testFile,'rb')
96                self.values = cPickle.load(file)
97                self.HistoPhases = cPickle.load(file)
98                (self.constrDict,self.fixedList,self.depVarList) = cPickle.load(file)
99                self.parmDict = cPickle.load(file)
100                self.varylist = cPickle.load(file)
101                self.calcControls = cPickle.load(file)
102                self.pawleyLookup = cPickle.load(file)
103                self.use = [False for i in range(len(self.varylist+self.depVarList))]
104                self.delt = [max(abs(self.parmDict[name])*0.001,1e-6) for name in self.varylist+self.depVarList]
105                file.close()
106                groups,parmlist = G2mv.GroupConstraints(self.constrDict)
107                G2mv.GenerateConstraints(groups,parmlist,self.varylist,self.constrDict,self.fixedList,self.parmDict)
108                self.UpdateControls(event)
109                print G2mv.VarRemapShow(self.varylist)
110                print 'Dependent Vary List:',self.depVarList
111        finally:
112            dlg.Destroy()
113           
114    def UpdateControls(self,event):
115       
116        def OnItemCk(event):
117            Obj = event.GetEventObject()
118            item = ObjInd[Obj.GetId()]
119            self.use[item] = Obj.GetValue()
120           
121        def OnDelValue(event):
122            event.Skip()
123            Obj = event.GetEventObject()
124            item = ObjInd[Obj.GetId()]
125            try:
126                value = float(Obj.GetValue())
127            except ValueError:
128                value = self.delt[item]
129            self.delt[item] = value
130            Obj.SetValue('%g'%(value))
131       
132        self.testDerivPanel.DestroyChildren()
133        ObjInd = {}
134        varylist = self.varylist
135        depVarList = self.depVarList
136        use = self.use
137        delt = self.delt
138        mainSizer = wx.FlexGridSizer(0,8,5,5)
139        for id,[ck,name,d] in enumerate(zip(use,varylist+depVarList,delt)):
140            useVal = wx.CheckBox(self.testDerivPanel,label=name)
141            useVal.SetValue(ck)
142            ObjInd[useVal.GetId()] = id
143            useVal.Bind(wx.EVT_CHECKBOX, OnItemCk)
144            mainSizer.Add(useVal,0)
145            delVal = wx.TextCtrl(self.testDerivPanel,wx.ID_ANY,'%g'%(d),style=wx.TE_PROCESS_ENTER)
146            ObjInd[delVal.GetId()] = id
147            delVal.Bind(wx.EVT_TEXT_ENTER,OnDelValue)
148            delVal.Bind(wx.EVT_KILL_FOCUS,OnDelValue)
149            mainSizer.Add(delVal,0)
150        self.testDerivPanel.SetSizer(mainSizer)   
151        Size = mainSizer.Fit(self.testDerivPanel)
152        Size[0] = 800
153        Size[1] = max(Size[1],290) + 35
154        self.testDerivPanel.SetScrollbars(10,10,Size[0]/10-4,Size[1]/10-1)
155        self.testDerivPanel.SetSize(Size)
156
157    def OnMakePlots(self,event):
158       
159        def test1():
160            fplot = self.plotNB.add('function test').gca()
161            M = G2stMth.errRefine(self.values,self.HistoPhases,
162                self.parmDict,self.varylist,self.calcControls,
163                self.pawleyLookup,None)
164            fplot.plot(M,'r',label='M')
165            fplot.legend(loc='best')
166           
167        def test2(name,delt):
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            M2 = dMdV[varyList.index(name)]
174            hplot.plot(M2,'b',label='analytic deriv')
175            if name in varyList:
176                mmin = np.min(dMdV[varyList.index(name)])
177                mmax = np.max(dMdV[varyList.index(name)])
178                print 'parameter:',name,self.parmDict[name],delt,mmin,mmax
179                if name in self.varylist:
180                    self.values[self.varylist.index(name)] -= delt
181                    M0 = G2stMth.errRefine(self.values,self.HistoPhases,self.parmDict,
182                        varyList,self.calcControls,self.pawleyLookup,None)
183                    self.values[self.varylist.index(name)] += 2.*delt
184                    M1 = G2stMth.errRefine(self.values,self.HistoPhases,self.parmDict,
185                        varyList,self.calcControls,self.pawleyLookup,None)
186                    self.values[self.varylist.index(name)] -= delt
187                elif name in self.depVarList:   #in depVarList
188                    if 'dA' in name:
189                        name = name.replace('dA','A')
190                        delt *= -1
191                    self.parmDict[name] -= delt
192                    M0 = G2stMth.errRefine(self.values,self.HistoPhases,self.parmDict,
193                        varyList,self.calcControls,self.pawleyLookup,None)
194                    self.parmDict[name] += 2.*delt
195                    M1 = G2stMth.errRefine(self.values,self.HistoPhases,self.parmDict,
196                        varyList,self.calcControls,self.pawleyLookup,None)
197                    self.parmDict[name] -= delt   
198                Mn = (M1-M0)/(2.*abs(delt))
199                hplot.plot(Mn,'r',label='numeric deriv')
200                hplot.plot(M2-Mn,'g',label='diff')
201#            GSASIIpath.IPyBreak()
202            hplot.legend(loc='best')           
203           
204        while self.plotNB.nb.GetPageCount():
205            self.plotNB.nb.DeletePage(0)
206        test1()
207        for use,name,delt in zip(self.use,self.varylist+self.depVarList,self.delt):
208            if use:
209                test2(name,delt)
210       
211        self.plotNB.Show()
212       
213class testDerivmain(wx.App):
214    def OnInit(self):
215        self.main = testDeriv(None)
216        self.main.Show()
217        self.SetTopWindow(self.main)
218        return True
219
220def main():
221    'Starts main application to compute and plot derivatives'
222    application = testDerivmain(0)
223    application.MainLoop()
224   
225if __name__ == '__main__':
226    GSASIIpath.InvokeDebugOpts()
227    main()
Note: See TracBrowser for help on using the repository browser.