source: trunk/testDeriv.py

Last change on this file was 5706, checked in by vondreele, 18 months ago

change deformation model to have only one kappa for valence orbitals & one for deformation orbitals.
A few fixes for testDeriv
  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 13.0 KB
RevLine 
[1299]1# -*- coding: utf-8 -*-
2#testDeriv.py
[5577]3########### SVN repository information ###################
4# $Date: 2023-12-12 18:48:03 +0000 (Tue, 12 Dec 2023) $
5# $Author: toby $
6# $Revision: 5706 $
7# $URL: trunk/testDeriv.py $
8# $Id: testDeriv.py 5706 2023-12-12 18:48:03Z toby $
9########### SVN repository information ###################
[1299]10'''
[2550]11To use set ``DEBUG=True`` in GSASIIstrMain.py (line 40, as of version
122546); run the least squares - zero cycles is sufficient.  Do the "Save
[1299]13Results"; this will write the file testDeriv.dat in the local
14directory.
15
16Then run this program to see plots of derivatives for all
17parameters refined in the last least squares.  Shown will be numerical
18derivatives generated over all observations (including penalty terms)
19and the corresponding analytical ones produced in the least
[2550]20squares. They should match. Profiling is also done for function
21calculation & for the 1st selected derivative (rest should be the same).
[1299]22'''
23
24import sys
[3415]25import os
[3340]26import platform
[5145]27import copy
[3340]28if '2' in platform.python_version_tuple()[0]:
29    import cPickle
30    import StringIO
31else:
[3828]32    import pickle as cPickle
[3340]33    import io as StringIO
34import cProfile,pstats
[1299]35import wx
36import numpy as np
37import GSASIIpath
[3232]38GSASIIpath.SetBinaryPath()
[1299]39import GSASIIstrMath as G2stMth
40import GSASIItestplot as plot
[1493]41import GSASIImapvars as G2mv
[4213]42try:  # fails on doc build
43    import pytexture as ptx
44    ptx.pyqlmninit()            #initialize fortran arrays for spherical harmonics
45except ImportError:
46    pass
[1299]47
[3826]48try:
[5145]49    NewId = wx.NewIdRef
[3826]50except AttributeError:
[5145]51    NewId = wx.NewId
[4521]52[wxID_FILEEXIT, wxID_FILEOPEN, wxID_MAKEPLOTS, wxID_CLEARSEL,wxID_SELECTALL,
[5145]53] = [NewId() for _init_coll_File_Items in range(5)]
[1299]54
55def FileDlgFixExt(dlg,file):            #this is needed to fix a problem in linux wx.FileDialog
56    ext = dlg.GetWildcard().split('|')[2*dlg.GetFilterIndex()+1].strip('*')
57    if ext not in file:
58        file += ext
59    return file
60   
61class testDeriv(wx.Frame):
62
63    def _init_ctrls(self, parent):
64        wx.Frame.__init__(self, name='testDeriv', parent=parent,
[1483]65            size=wx.Size(800, 250),style=wx.DEFAULT_FRAME_STYLE, title='Test Jacobian Derivatives')
[1299]66        self.testDerivMenu = wx.MenuBar()
67        self.File = wx.Menu(title='')
[5145]68        self.File.Append(wxID_FILEOPEN,'Open testDeriv file\tCtrl+O','Open testDeriv')
69        self.File.Append(wxID_MAKEPLOTS,'Make plots\tCtrl+P','Make derivative plots')
70        self.File.Append(wxID_SELECTALL,'Select all\tCtrl+S')
71        self.File.Append(wxID_CLEARSEL,'Clear selections\tCtrl+C')
72        self.File.Append(wxID_FILEEXIT,'Exit\tALT+F4','Exit from testDeriv')
[3774]73        self.Bind(wx.EVT_MENU,self.OnTestRead, id=wxID_FILEOPEN)
[1299]74        self.Bind(wx.EVT_MENU,self.OnMakePlots,id=wxID_MAKEPLOTS)
[3774]75        self.Bind(wx.EVT_MENU,self.ClearSelect,id=wxID_CLEARSEL)
[4521]76        self.Bind(wx.EVT_MENU,self.SelectAll,id=wxID_SELECTALL)
[3774]77        self.Bind(wx.EVT_MENU,self.OnFileExit, id=wxID_FILEEXIT)
[1299]78        self.testDerivMenu.Append(menu=self.File, title='File')
79        self.SetMenuBar(self.testDerivMenu)
[5145]80        self.testDerivPanel = wx.ScrolledWindow(self)
[1299]81        self.plotNB = plot.PlotNotebook()
[3415]82        self.testFile = ''
83        arg = sys.argv
84        if len(arg) > 1 and arg[1]:
85            try:
86                self.testFile = os.path.splitext(arg[1])[0]+u'.testDeriv'
87            except:
88                self.testFile = os.path.splitext(arg[1])[0]+'.testDeriv'
89            self.TestRead()
90            self.UpdateControls(None)
[1299]91       
92    def __init__(self, parent):
93        self._init_ctrls(parent)
94        self.Bind(wx.EVT_CLOSE, self.ExitMain)   
95        self.dirname = ''
96        self.testfile = []
97        self.dataFrame = None
[5145]98        self.timingOn = False
[1299]99
100    def ExitMain(self, event):
101        sys.exit()
102       
103    def OnFileExit(self,event):
104        if self.dataFrame:
105            self.dataFrame.Clear() 
106            self.dataFrame.Destroy()
107        self.Close()
[3774]108       
[4521]109    def SelectAll(self,event):
110        self.use = [True for name in self.names]
111        for i,name in enumerate(self.names):
112            if 'Back' in name:
113                self.use[i] = False
114        self.UpdateControls(event)
115       
[3774]116    def ClearSelect(self,event):
117        self.use = [False for i in range(len(self.names))]
118        self.UpdateControls(event)
[1299]119
120    def OnTestRead(self,event):
[3369]121        dlg = wx.FileDialog(self, 'Open *.testDeriv file',defaultFile='*.testDeriv',
122            wildcard='*.testDeriv')
[1299]123        if self.dirname:
124            dlg.SetDirectory(self.dirname)
125        try:
126            if dlg.ShowModal() == wx.ID_OK:
127                self.dirname = dlg.GetDirectory()
[3415]128                self.testFile = dlg.GetPath()
129                self.TestRead()
[1299]130                self.UpdateControls(event)
131        finally:
132            dlg.Destroy()
133           
[3415]134    def TestRead(self):
135        file = open(self.testFile,'rb')
[4521]136        self.values = cPickle.load(file,encoding='Latin-1')
137        self.HistoPhases = cPickle.load(file,encoding='Latin-1')
138        (self.constrDict,self.fixedList,self.depVarList) = cPickle.load(file,encoding='Latin-1')
139        self.parmDict = cPickle.load(file,encoding='Latin-1')
140        self.varylist = cPickle.load(file,encoding='Latin-1')
141        self.calcControls = cPickle.load(file,encoding='Latin-1')
142        self.pawleyLookup = cPickle.load(file,encoding='Latin-1')
[3774]143        self.names = self.varylist+self.depVarList
144        self.use = [False for i in range(len(self.names))]
145        self.delt = [max(abs(self.parmDict[name])*0.0001,1e-6) for name in self.names]
146        for iname,name in enumerate(self.names):
147            if name.split(':')[-1] in ['Shift','DisplaceX','DisplaceY',]:
148                self.delt[iname] = 0.1
[3415]149        file.close()
[5145]150        G2mv.InitVars()
[3711]151        msg = G2mv.EvaluateMultipliers(self.constrDict,self.parmDict)
152        if msg:
153            print('Unable to interpret multiplier(s): '+msg)
154            raise Exception
155        G2mv.GenerateConstraints(self.varylist,self.constrDict,self.fixedList,self.parmDict)
[3415]156        print(G2mv.VarRemapShow(self.varylist))
157        print('Dependent Vary List:',self.depVarList)
[5145]158        G2mv.Map2Dict(self.parmDict,copy.copy(self.varylist))   # compute independent params, N.B. changes varylist
159        G2mv.Dict2Map(self.parmDict) # imposes constraints on dependent values
160
[1299]161    def UpdateControls(self,event):
162        def OnItemCk(event):
163            Obj = event.GetEventObject()
164            item = ObjInd[Obj.GetId()]
165            self.use[item] = Obj.GetValue()
166           
167        def OnDelValue(event):
[2500]168            event.Skip()
[1299]169            Obj = event.GetEventObject()
170            item = ObjInd[Obj.GetId()]
171            try:
172                value = float(Obj.GetValue())
173            except ValueError:
174                value = self.delt[item]
175            self.delt[item] = value
176            Obj.SetValue('%g'%(value))
177       
[5145]178        if self.testDerivPanel.GetSizer():
179            self.testDerivPanel.GetSizer().Clear(True)
[1299]180        ObjInd = {}
181        use = self.use
182        delt = self.delt
[5145]183        topSizer = wx.BoxSizer(wx.VERTICAL)
184        self.timingVal = wx.CheckBox(self.testDerivPanel,label='Show Execution Profiling')
185        topSizer.Add(self.timingVal,0)
186        topSizer.Add((-1,10))
[1378]187        mainSizer = wx.FlexGridSizer(0,8,5,5)
[5145]188        for id,[ck,name,d] in enumerate(zip(use,self.names,delt)):
[1299]189            useVal = wx.CheckBox(self.testDerivPanel,label=name)
190            useVal.SetValue(ck)
191            ObjInd[useVal.GetId()] = id
192            useVal.Bind(wx.EVT_CHECKBOX, OnItemCk)
193            mainSizer.Add(useVal,0)
194            delVal = wx.TextCtrl(self.testDerivPanel,wx.ID_ANY,'%g'%(d),style=wx.TE_PROCESS_ENTER)
195            ObjInd[delVal.GetId()] = id
196            delVal.Bind(wx.EVT_TEXT_ENTER,OnDelValue)
197            delVal.Bind(wx.EVT_KILL_FOCUS,OnDelValue)
198            mainSizer.Add(delVal,0)
[5145]199        topSizer.Add(mainSizer,0)
200        self.testDerivPanel.SetSizer(topSizer)   
201        Size = topSizer.GetMinSize()
202        self.testDerivPanel.SetScrollbars(10,10,int(Size[0]/10-4),int(Size[1]/10-1))
[5509]203        Size[1] = max(200,Size[1])
[5145]204        Size[0] += 20
205        self.SetSize(Size)
[1299]206
207    def OnMakePlots(self,event):
208       
209        def test1():
210            fplot = self.plotNB.add('function test').gca()
[5145]211            if self.timingOn:
212                pr = cProfile.Profile()
213                pr.enable()
[1299]214            M = G2stMth.errRefine(self.values,self.HistoPhases,
215                self.parmDict,self.varylist,self.calcControls,
216                self.pawleyLookup,None)
[5145]217            if self.timingOn:
218                pr.disable()
219                s = StringIO.StringIO()
220                sortby = 'tottime'
221                ps = pstats.Stats(pr, stream=s).strip_dirs().sort_stats(sortby)
222                print('Profiler of function calculation; top 50% of routines:')
223                ps.print_stats("GSASII",.5)
224                print(s.getvalue())
[1299]225            fplot.plot(M,'r',label='M')
226            fplot.legend(loc='best')
227           
[2548]228        def test2(name,delt,doProfile):
[1860]229            Title = 'derivatives test for '+name
[5145]230            ind = self.names.index(name)
[1860]231            hplot = self.plotNB.add(Title).gca()
[5145]232            if doProfile and self.timingOn:
[2548]233                pr = cProfile.Profile()
234                pr.enable()
[3774]235            #regenerate minimization fxn
236            G2stMth.errRefine(self.values,self.HistoPhases,
237                self.parmDict,self.varylist,self.calcControls,
238                self.pawleyLookup,None)
[1299]239            dMdV = G2stMth.dervRefine(self.values,self.HistoPhases,self.parmDict,
[5145]240                self.names,self.calcControls,self.pawleyLookup,None)
241            if doProfile and self.timingOn:
[2548]242                pr.disable()
243                s = StringIO.StringIO()
244                sortby = 'tottime'
245                ps = pstats.Stats(pr, stream=s).strip_dirs().sort_stats(sortby)
246                ps.print_stats("GSASII",.5)
[3166]247                print('Profiler of '+name+' derivative calculation; top 50% of routines:')
248                print(s.getvalue())
[5145]249            M2 = dMdV[ind]
[2508]250            hplot.plot(M2,'b',label='analytic deriv')
[5145]251            mmin = np.min(dMdV[ind])
252            mmax = np.max(dMdV[ind])
[2539]253            if name in self.varylist:
[5145]254                ind = self.varylist.index(name)
255                orig = copy.copy(self.parmDict)  # save parmDict before changes
256                self.parmDict[name] = self.values[ind] =  self.values[ind] - delt
257                G2mv.Dict2Map(self.parmDict)
258                first = True
259                for i in self.parmDict:
[5706]260                    if 'UVmat' in i:
261                        continue
[5145]262                    if orig[i] != self.parmDict[i] and i != name:
263                        if first:
[5506]264                            print('Propagated changes from this shift')
[5145]265                            print(name,orig[name],self.parmDict[name],orig[name]-self.parmDict[name])
266                            print('are:')
267                            first = False
268                        print(i,orig[i],self.parmDict[i],orig[i]-self.parmDict[i])
[2539]269                M0 = G2stMth.errRefine(self.values,self.HistoPhases,self.parmDict,
[5145]270                    self.names,self.calcControls,self.pawleyLookup,None)
271                self.parmDict[name] = self.values[ind] =  self.values[ind] + 2.*delt
272                G2mv.Dict2Map(self.parmDict)
[2539]273                M1 = G2stMth.errRefine(self.values,self.HistoPhases,self.parmDict,
[5145]274                    self.names,self.calcControls,self.pawleyLookup,None)
275                self.parmDict[name] = self.values[ind] =  self.values[ind] - delt
276                G2mv.Dict2Map(self.parmDict)
[2539]277            elif name in self.depVarList:   #in depVarList
278                if 'dA' in name:
279                    name = name.replace('dA','A')
[5145]280                    #delt *= -1  # why???
[2539]281                self.parmDict[name] -= delt
[5145]282                G2mv.Dict2Map(self.parmDict)
[2539]283                M0 = G2stMth.errRefine(self.values,self.HistoPhases,self.parmDict,
[5145]284                        self.names,self.calcControls,self.pawleyLookup,None)
[2539]285                self.parmDict[name] += 2.*delt
[5145]286                G2mv.Dict2Map(self.parmDict)
[2539]287                M1 = G2stMth.errRefine(self.values,self.HistoPhases,self.parmDict,
[5145]288                        self.names,self.calcControls,self.pawleyLookup,None)
[2539]289                self.parmDict[name] -= delt   
[5145]290                G2mv.Dict2Map(self.parmDict)
[2539]291            Mn = (M1-M0)/(2.*abs(delt))
[5114]292            print('parameter:',name,self.parmDict[name],delt,mmin,mmax,np.sum(M0),np.sum(M1),np.sum(Mn))
[2539]293            hplot.plot(Mn,'r',label='numeric deriv')
[1860]294            hplot.legend(loc='best')           
[1299]295           
296        while self.plotNB.nb.GetPageCount():
297            self.plotNB.nb.DeletePage(0)
[2548]298           
[2501]299        test1()
[5145]300        self.timingOn = self.timingVal.GetValue()
[2548]301
302        doProfile = True
[5145]303        for use,name,delt in zip(self.use,self.names,self.delt):
[1299]304            if use:
[2548]305                test2(name,delt,doProfile)
306                doProfile = False
[1299]307       
308        self.plotNB.Show()
309       
310def main():
311    'Starts main application to compute and plot derivatives'
[5145]312    application = wx.App(0)
313    application.main = testDeriv(None)
314    application.main.Show()
315    application.SetTopWindow(application.main)
[1299]316    application.MainLoop()
317   
318if __name__ == '__main__':
[2500]319    GSASIIpath.InvokeDebugOpts()
[1299]320    main()
Note: See TracBrowser for help on using the repository browser.