source: trunk/testDeriv.py @ 1113

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

fix Babinet derivatives for SC & powder

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