source: trunk/GSASIIstrMain.py @ 3136

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

make GSAS-II python 3.6 compliant & preserve python 2.7 use;changes:
do from future import division, print_function for all GSAS-II py sources
all menu items revised to be py 2.7/3.6 compliant
all wx.OPEN --> wx.FD_OPEN in file dialogs
all integer divides (typically for image pixel math) made explicit with ; ambiguous ones made floats as appropriate
all print "stuff" --> print (stuff)
all print >> pFile,'stuff' --> pFile.writeCIFtemplate('stuff')
all read file opens made explicit 'r' or 'rb'
all cPickle imports made for py2.7 or 3.6 as cPickle or _pickle; test for '2' platform.version_tuple[0] for py 2.7
define cPickleload to select load(fp) or load(fp,encoding='latin-1') for loading gpx files; provides cross compatibility between py 2.7/3.6 gpx files
make dict.keys() as explicit list(dict.keys()) as needed (NB: possible source of remaining py3.6 bugs)
make zip(a,b) as explicit list(zip(a,b)) as needed (NB: possible source of remaining py3.6 bugs)
select unichr/chr according test for '2' platform.version_tuple[0] for py 2.7 (G2pwdGUI * G2plot) for special characters
select wg.EVT_GRID_CELL_CHANGE (classic) or wg.EVT_GRID_CELL_CHANGED (phoenix) in grid Bind
maxint --> maxsize; used in random number stuff
raise Exception,"stuff" --> raise Exception("stuff")
wx 'classic' sizer.DeleteWindows?() or 'phoenix' sizer.Clear(True)
wx 'classic' SetToolTipString?(text) or 'phoenix' SetToolTip?(wx.ToolTip?(text)); define SetToolTipString?(self,text) to handle the choice in plots
status.SetFields? --> status.SetStatusText?
'classic' AddSimpleTool? or 'phoenix' self.AddTool? for plot toolbar; Bind different as well
define GetItemPydata? as it doesn't exist in wx 'phoenix'
allow python versions 2.7 & 3.6 to run GSAS-II
Bind override commented out - no logging capability (NB: remove all logging code?)
all import ContentsValidator? open filename & test if valid then close; filepointer removed from Reader
binary importers (mostly images) test for 'byte' type & convert as needed to satisfy py 3.6 str/byte rules

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 36.2 KB
Line 
1# -*- coding: utf-8 -*-
2'''
3*GSASIIstrMain: main structure routine*
4---------------------------------------
5
6'''
7########### SVN repository information ###################
8# $Date: 2017-10-23 16:39:16 +0000 (Mon, 23 Oct 2017) $
9# $Author: vondreele $
10# $Revision: 3136 $
11# $URL: trunk/GSASIIstrMain.py $
12# $Id: GSASIIstrMain.py 3136 2017-10-23 16:39:16Z vondreele $
13########### SVN repository information ###################
14from __future__ import division, print_function
15import platform
16import sys
17import os.path as ospath
18import time
19import math
20import copy
21if '2' in platform.python_version_tuple()[0]:
22    import cPickle
23else:
24    import pickle as cPickle
25import numpy as np
26import numpy.linalg as nl
27import scipy.optimize as so
28import GSASIIpath
29GSASIIpath.SetBinaryPath()
30GSASIIpath.SetVersionNumber("$Revision: 3136 $")
31import GSASIIlattice as G2lat
32import GSASIIspc as G2spc
33import GSASIImapvars as G2mv
34import GSASIImath as G2mth
35import GSASIIstrIO as G2stIO
36import GSASIIstrMath as G2stMth
37import GSASIIobj as G2obj
38
39sind = lambda x: np.sin(x*np.pi/180.)
40cosd = lambda x: np.cos(x*np.pi/180.)
41tand = lambda x: np.tan(x*np.pi/180.)
42asind = lambda x: 180.*np.arcsin(x)/np.pi
43acosd = lambda x: 180.*np.arccos(x)/np.pi
44atan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
45
46ateln2 = 8.0*math.log(2.0)
47DEBUG = True
48
49def RefineCore(Controls,Histograms,Phases,restraintDict,rigidbodyDict,parmDict,varyList,
50    calcControls,pawleyLookup,ifPrint,printFile,dlg):
51    '''Core optimization routines, shared between SeqRefine and Refine
52
53    :returns: 5-tuple of ifOk (bool), Rvals (dict), result, covMatrix, sig
54    '''
55#    print 'current',varyList
56#    for item in parmDict: print item,parmDict[item] ######### show dict just before refinement
57    G2mv.Map2Dict(parmDict,varyList)
58    Rvals = {}
59    while True:
60        begin = time.time()
61        values =  np.array(G2stMth.Dict2Values(parmDict, varyList))
62        if np.any(np.isnan(values)):
63            raise G2obj.G2Exception('ERROR - nan found in LS parameters - use Calculate/View LS parms to locate')
64        # test code to compute GOF and save for external repeat
65        #args = ([Histograms,Phases,restraintDict,rigidbodyDict],parmDict,varyList,calcControls,pawleyLookup,dlg)
66        #print '*** before fit chi**2',np.sum(G2stMth.errRefine(values,*args)**2)
67        #fl = open('beforeFit.cpickle','wb')
68        #cPickle.dump(values,fl,1)
69        #cPickle.dump(args[:-1],fl,1)
70        #fl.close()
71        Ftol = Controls['min dM/M']
72        Xtol = Controls['SVDtol']
73        Factor = Controls['shift factor']
74        if 'Jacobian' in Controls['deriv type']:
75            result = so.leastsq(G2stMth.errRefine,values,Dfun=G2stMth.dervRefine,full_output=True,
76                ftol=Ftol,col_deriv=True,factor=Factor,
77                args=([Histograms,Phases,restraintDict,rigidbodyDict],parmDict,varyList,calcControls,pawleyLookup,dlg))
78            ncyc = int(result[2]['nfev']/2)
79        elif 'analytic Hessian' in Controls['deriv type']:
80            Lamda = Controls.get('Marquardt',-3)
81            maxCyc = Controls['max cyc']
82            result = G2mth.HessianLSQ(G2stMth.errRefine,values,Hess=G2stMth.HessRefine,ftol=Ftol,xtol=Xtol,maxcyc=maxCyc,Print=ifPrint,lamda=Lamda,
83                args=([Histograms,Phases,restraintDict,rigidbodyDict],parmDict,varyList,calcControls,pawleyLookup,dlg))
84            ncyc = result[2]['num cyc']+1
85            Rvals['lamMax'] = result[2]['lamMax']
86            Controls['Marquardt'] = -3  #reset to default
87        elif 'Hessian SVD' in Controls['deriv type']:
88            maxCyc = Controls['max cyc']
89            result = G2mth.HessianSVD(G2stMth.errRefine,values,Hess=G2stMth.HessRefine,ftol=Ftol,xtol=Xtol,maxcyc=maxCyc,Print=ifPrint,
90                args=([Histograms,Phases,restraintDict,rigidbodyDict],parmDict,varyList,calcControls,pawleyLookup,dlg))
91            ncyc = result[2]['num cyc']+1
92        else:           #'numeric'
93            result = so.leastsq(G2stMth.errRefine,values,full_output=True,ftol=Ftol,epsfcn=1.e-8,factor=Factor,
94                args=([Histograms,Phases,restraintDict,rigidbodyDict],parmDict,varyList,calcControls,pawleyLookup,dlg))
95            ncyc = 1
96            if len(varyList):
97                ncyc = int(result[2]['nfev']/len(varyList))
98#        table = dict(zip(varyList,zip(values,result[0],(result[0]-values))))
99#        for item in table: print item,table[item]               #useful debug - are things shifting?
100        runtime = time.time()-begin
101        Rvals['SVD0'] = result[2].get('SVD0',0)
102        Rvals['converged'] = result[2].get('Converged')
103        Rvals['DelChi2'] = result[2].get('DelChi2',-1.)
104        Rvals['chisq'] = np.sum(result[2]['fvec']**2)
105        G2stMth.Values2Dict(parmDict, varyList, result[0])
106        G2mv.Dict2Map(parmDict,varyList)
107        Rvals['Nobs'] = Histograms['Nobs']
108        Rvals['Rwp'] = np.sqrt(Rvals['chisq']/Histograms['sumwYo'])*100.      #to %
109        Rvals['GOF'] = np.sqrt(Rvals['chisq']/(Histograms['Nobs']-len(varyList)))
110        printFile.write(' Number of function calls: %d No. of observations: %d No. of parameters: %d User rejected: %d Sp. gp. extinct: %d\n'%  \
111            (result[2]['nfev'],Histograms['Nobs'],len(varyList),Histograms['Nrej'],Histograms['Next']))
112        printFile.write(' Refinement time = %8.3fs, %8.3fs/cycle, for %d cycles\n'%(runtime,runtime/ncyc,ncyc))
113        printFile.write(' wR = %7.2f%%, chi**2 = %12.6g, GOF = %6.2f\n'%(Rvals['Rwp'],Rvals['chisq'],Rvals['GOF']))
114        IfOK = True
115        try:
116            covMatrix = result[1]*Rvals['GOF']**2
117            sig = np.sqrt(np.diag(covMatrix))
118            if np.any(np.isnan(sig)) or not sig.shape:
119                print ('*** Least squares aborted - some invalid esds possible ***')
120#            table = dict(zip(varyList,zip(values,result[0],(result[0]-values)/sig)))
121#            for item in table: print item,table[item]               #useful debug - are things shifting?
122            break                   #refinement succeeded - finish up!
123        except TypeError:          #result[1] is None on singular matrix or LinAlgError
124            IfOK = False
125            if not len(varyList):
126                covMatrix = []
127                sig = []
128                break
129            print ('**** Refinement failed - singular matrix ****')
130            if 'Hessian' in Controls['deriv type']:
131                num = len(varyList)-1
132                for i,val in enumerate(np.flipud(result[2]['psing'])):
133                    if val:
134                        print ('Removing parameter: '+varyList[num-i])
135                        del(varyList[num-i])
136            else:
137                Ipvt = result[2]['ipvt']
138                for i,ipvt in enumerate(Ipvt):
139                    if not np.sum(result[2]['fjac'],axis=1)[i]:
140                        print ('Removing parameter: '+varyList[ipvt-1])
141                        del(varyList[ipvt-1])
142                        break
143    if IfOK:
144        G2stMth.GetFobsSq(Histograms,Phases,parmDict,calcControls)
145    return IfOK,Rvals,result,covMatrix,sig
146
147def Refine(GPXfile,dlg=None,makeBack=True):
148    'Global refinement -- refines to minimize against all histograms'
149    import GSASIImpsubs as G2mp
150    G2mp.InitMP()
151    import pytexture as ptx
152    ptx.pyqlmninit()            #initialize fortran arrays for spherical harmonics
153
154    printFile = open(ospath.splitext(GPXfile)[0]+'.lst','w')
155    G2stIO.ShowBanner(printFile)
156    varyList = []
157    parmDict = {}
158    G2mv.InitVars()
159    Controls = G2stIO.GetControls(GPXfile)
160    G2stIO.ShowControls(Controls,printFile)
161    calcControls = {}
162    calcControls.update(Controls)
163    constrDict,fixedList = G2stIO.GetConstraints(GPXfile)
164    restraintDict = G2stIO.GetRestraints(GPXfile)
165    Histograms,Phases = G2stIO.GetUsedHistogramsAndPhases(GPXfile)
166    if not Phases:
167        print (' *** ERROR - you have no phases to refine! ***')
168        print (' *** Refine aborted ***')
169        return False,'No phases'
170    if not Histograms:
171        print (' *** ERROR - you have no data to refine with! ***')
172        print (' *** Refine aborted ***')
173        return False,'No data'
174    rigidbodyDict = G2stIO.GetRigidBodies(GPXfile)
175    rbIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})
176    rbVary,rbDict = G2stIO.GetRigidBodyModels(rigidbodyDict,pFile=printFile)
177    Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables,MFtables,maxSSwave = \
178        G2stIO.GetPhaseData(Phases,restraintDict,rbIds,pFile=printFile)
179    calcControls['atomIndx'] = atomIndx
180    calcControls['Natoms'] = Natoms
181    calcControls['FFtables'] = FFtables
182    calcControls['BLtables'] = BLtables
183    calcControls['MFtables'] = MFtables
184    calcControls['maxSSwave'] = maxSSwave
185    hapVary,hapDict,controlDict = G2stIO.GetHistogramPhaseData(Phases,Histograms,pFile=printFile)
186    TwConstr,TwFixed = G2stIO.makeTwinFrConstr(Phases,Histograms,hapVary)
187    constrDict += TwConstr
188    fixedList += TwFixed
189    calcControls.update(controlDict)
190    histVary,histDict,controlDict = G2stIO.GetHistogramData(Histograms,pFile=printFile)
191    calcControls.update(controlDict)
192    varyList = rbVary+phaseVary+hapVary+histVary
193    parmDict.update(rbDict)
194    parmDict.update(phaseDict)
195    parmDict.update(hapDict)
196    parmDict.update(histDict)
197    G2stIO.GetFprime(calcControls,Histograms)
198    # do constraint processing
199    varyListStart = tuple(varyList) # save the original varyList before dependent vars are removed
200    try:
201        groups,parmlist = G2mv.GroupConstraints(constrDict)
202        G2mv.GenerateConstraints(groups,parmlist,varyList,constrDict,fixedList,parmDict)
203        #print G2mv.VarRemapShow(varyList)
204        #print 'DependentVars',G2mv.GetDependentVars()
205        #print 'IndependentVars',G2mv.GetIndependentVars()
206    except:
207        print (' *** ERROR - your constraints are internally inconsistent ***')
208        #errmsg, warnmsg = G2mv.CheckConstraints(varyList,constrDict,fixedList)
209        #print 'Errors',errmsg
210        #if warnmsg: print 'Warnings',warnmsg
211        return False,' Constraint error'
212#    print G2mv.VarRemapShow(varyList)
213
214    ifPrint = True
215    printFile.write('\n Refinement results:\n')
216    printFile.write(135*'-'+'\n')
217    if True:
218#    try:
219        covData = {}
220        IfOK,Rvals,result,covMatrix,sig = RefineCore(Controls,Histograms,Phases,restraintDict,
221            rigidbodyDict,parmDict,varyList,calcControls,pawleyLookup,ifPrint,printFile,dlg)
222        if IfOK:
223            sigDict = dict(zip(varyList,sig))
224            newCellDict = G2stMth.GetNewCellParms(parmDict,varyList)
225            newAtomDict = G2stMth.ApplyXYZshifts(parmDict,varyList)
226            covData = {'variables':result[0],'varyList':varyList,'sig':sig,'Rvals':Rvals,
227                       'varyListStart':varyListStart,
228                       'covMatrix':covMatrix,'title':GPXfile,'newAtomDict':newAtomDict,
229                       'newCellDict':newCellDict,'freshCOV':True}
230            # add the uncertainties into the esd dictionary (sigDict)
231            sigDict.update(G2mv.ComputeDepESD(covMatrix,varyList,parmDict))
232            G2mv.PrintIndependentVars(parmDict,varyList,sigDict,pFile=printFile)
233            G2stMth.ApplyRBModels(parmDict,Phases,rigidbodyDict,True)
234            G2stIO.SetRigidBodyModels(parmDict,sigDict,rigidbodyDict,printFile)
235            G2stIO.SetPhaseData(parmDict,sigDict,Phases,rbIds,covData,restraintDict,printFile)
236            G2stIO.SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms,calcControls['FFtables'],pFile=printFile)
237            G2stIO.SetHistogramData(parmDict,sigDict,Histograms,calcControls['FFtables'],pFile=printFile)
238            G2stIO.SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases,rigidbodyDict,covData,makeBack)
239            printFile.close()
240            print (' Refinement results are in file: '+ospath.splitext(GPXfile)[0]+'.lst')
241            print (' ***** Refinement successful *****')
242        else:
243            print ('****ERROR - Refinement failed')
244            raise G2obj.G2Exception('****ERROR - Refinement failed')
245#    except G2obj.G2Exception(Msg):
246#        printFile.close()
247#        return False,Msg.msg
248
249#for testing purposes!!!
250    if DEBUG and IfOK:
251#needs: values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup
252        fl = open('testDeriv.dat','wb')
253        cPickle.dump(result[0],fl,1)
254        cPickle.dump([Histograms,Phases,restraintDict,rigidbodyDict],fl,1)
255        cPickle.dump([constrDict,fixedList,G2mv.GetDependentVars()],fl,1)
256        cPickle.dump(parmDict,fl,1)
257        cPickle.dump(varyList,fl,1)
258        cPickle.dump(calcControls,fl,1)
259        cPickle.dump(pawleyLookup,fl,1)
260        fl.close()
261    if dlg:
262        return True,Rvals
263
264def phaseCheck(phaseVary,Phases,histogram):
265    '''
266    Removes unused parameters from phase varylist if phase not in histogram
267    #TODO - implement "Fix FXU" for seq refinement here - done?
268    '''
269    NewVary = []
270    for phase in Phases:
271        if histogram not in Phases[phase]['Histograms']: continue
272        if Phases[phase]['Histograms'][histogram]['Use']:
273            pId = Phases[phase]['pId']
274            newVary = [item for item in phaseVary if item.split(':')[0] == str(pId)]
275            FixVals = Phases[phase]['Histograms'][histogram].get('Fix FXU',' ')
276            if 'F' in FixVals:
277                newVary = [item for item in newVary if not 'AF' in item]
278            if 'X' in FixVals:
279                newVary = [item for item in newVary if not 'dA' in item]
280            if 'U' in FixVals:
281                newVary = [item for item in newVary if not 'AU' in item]
282            if 'M' in FixVals:
283                newVary = [item for item in newVary if not 'AM' in item]
284            NewVary += newVary
285    return NewVary
286
287def SeqRefine(GPXfile,dlg,PlotFunction=None,G2frame=None):
288    '''Perform a sequential refinement -- cycles through all selected histgrams,
289    one at a time
290    '''
291    import GSASIImpsubs as G2mp
292    G2mp.InitMP()
293    import pytexture as ptx
294    ptx.pyqlmninit()            #initialize fortran arrays for spherical harmonics
295
296    printFile = open(ospath.splitext(GPXfile)[0]+'.lst','w')
297    print ('Starting Sequential Refinement')
298    G2stIO.ShowBanner(printFile)
299    Controls = G2stIO.GetControls(GPXfile)
300    G2stIO.ShowControls(Controls,printFile,SeqRef=True)
301    restraintDict = G2stIO.GetRestraints(GPXfile)
302    Histograms,Phases = G2stIO.GetUsedHistogramsAndPhases(GPXfile)
303    if not Phases:
304        print (' *** ERROR - you have no phases to refine! ***')
305        print (' *** Refine aborted ***')
306        return False,'No phases'
307    if not Histograms:
308        print (' *** ERROR - you have no data to refine with! ***')
309        print (' *** Refine aborted ***')
310        return False,'No data'
311    rigidbodyDict = G2stIO.GetRigidBodies(GPXfile)
312    rbIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})
313    rbVary,rbDict = G2stIO.GetRigidBodyModels(rigidbodyDict,pFile=printFile)
314    Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables,MFtables,maxSSwave = \
315        G2stIO.GetPhaseData(Phases,restraintDict,rbIds,False,printFile,seqRef=True)
316    for item in phaseVary:
317        if '::A0' in item:
318            print ('**** WARNING - lattice parameters should not be refined in a sequential refinement ****')
319            print ('****           instead use the Dij parameters for each powder histogram            ****')
320            return False,'Lattice parameter refinement error - see console message'
321        if '::C(' in item:
322            print ('**** WARNING - phase texture parameters should not be refined in a sequential refinement ****')
323            print ('****           instead use the C(L,N) parameters for each powder histogram               ****')
324            return False,'Phase texture refinement error - see console message'
325    if 'Seq Data' in Controls:
326        histNames = Controls['Seq Data']
327    else:
328        histNames = G2stIO.GetHistogramNames(GPXfile,['PWDR',])
329    if Controls.get('Reverse Seq'):
330        histNames.reverse()
331    SeqResult = G2stIO.GetSeqResult(GPXfile)
332#    SeqResult = {'SeqPseudoVars':{},'SeqParFitEqList':[]}
333    makeBack = True
334    Histo = {}
335    NewparmDict = {}
336    for ihst,histogram in enumerate(histNames):
337        print('\nRefining with '+str(histogram))
338        ifPrint = False
339        if dlg:
340            dlg.SetTitle('Residual for histogram '+str(ihst))
341        calcControls = {}
342        calcControls['atomIndx'] = atomIndx
343        calcControls['Natoms'] = Natoms
344        calcControls['FFtables'] = FFtables
345        calcControls['BLtables'] = BLtables
346        calcControls['MFtables'] = MFtables
347        calcControls['maxSSwave'] = maxSSwave
348        if histogram not in Histograms:
349            print("Error: not found!")
350            continue
351    #TODO - implement "Fix FXU" for seq refinement here - done?
352        redphaseVary = phaseCheck(phaseVary,Phases,histogram)
353        Histo = {histogram:Histograms[histogram],}
354        hapVary,hapDict,controlDict = G2stIO.GetHistogramPhaseData(Phases,Histo,Print=False)
355        calcControls.update(controlDict)
356        histVary,histDict,controlDict = G2stIO.GetHistogramData(Histo,False)
357        calcControls.update(controlDict)
358        varyList = rbVary+redphaseVary+hapVary+histVary
359#        if not ihst:
360            # save the initial vary list, but without histogram numbers on parameters
361        saveVaryList = varyList[:]
362        for i,item in enumerate(saveVaryList):
363            items = item.split(':')
364            if items[1]:
365                items[1] = ''
366            item = ':'.join(items)
367            saveVaryList[i] = item
368        if not ihst:
369            SeqResult['varyList'] = saveVaryList
370        else:
371            SeqResult['varyList'] = list(set(SeqResult['varyList']+saveVaryList))
372        parmDict = {}
373        parmDict.update(phaseDict)
374        parmDict.update(hapDict)
375        parmDict.update(histDict)
376        if Controls['Copy2Next']:
377            #parmDict.update(NewparmDict) # don't use in case extra entries would cause a problem
378            for parm in NewparmDict:
379                if parm in parmDict:
380                    parmDict[parm] = NewparmDict[parm]
381        G2stIO.GetFprime(calcControls,Histo)
382        # do constraint processing
383        #reload(G2mv) # debug
384        G2mv.InitVars()
385        constrDict,fixedList = G2stIO.GetConstraints(GPXfile)
386        varyListStart = tuple(varyList) # save the original varyList before dependent vars are removed
387        try:
388            groups,parmlist = G2mv.GroupConstraints(constrDict)
389            G2mv.GenerateConstraints(groups,parmlist,varyList,constrDict,fixedList,parmDict,SeqHist=ihst)
390            constraintInfo = (groups,parmlist,constrDict,fixedList,ihst)
391        except:
392            print (' *** ERROR - your constraints are internally inconsistent ***')
393            #errmsg, warnmsg = G2mv.CheckConstraints(varyList,constrDict,fixedList)
394            #print 'Errors',errmsg
395            #if warnmsg: print 'Warnings',warnmsg
396            return False,' Constraint error'
397        #print G2mv.VarRemapShow(varyList)
398        if not ihst:
399            # first histogram to refine against
400            firstVaryList = []
401            for item in varyList:
402                items = item.split(':')
403                if items[1]:
404                    items[1] = ''
405                item = ':'.join(items)
406                firstVaryList.append(item)
407            newVaryList = firstVaryList
408        else:
409            newVaryList = []
410            for item in varyList:
411                items = item.split(':')
412                if items[1]:
413                    items[1] = ''
414                item = ':'.join(items)
415                newVaryList.append(item)
416        if newVaryList != firstVaryList and Controls['Copy2Next']:
417            # variable lists are expected to match between sequential refinements when Copy2Next is on
418            #print '**** ERROR - variable list for this histogram does not match previous'
419            #print '     Copy of variables is not possible'
420            #print '\ncurrent histogram',histogram,'has',len(newVaryList),'variables'
421            combined = list(set(firstVaryList+newVaryList))
422            c = [var for var in combined if var not in newVaryList]
423            p = [var for var in combined if var not in firstVaryList]
424            print('*** Variables change ***')
425            for typ,vars in [('Removed',c),('Added',p)]:
426                line = '  '+typ+': '
427                if vars:
428                    for var in vars:
429                        if len(line) > 70:
430                            print(line)
431                            line = '    '
432                        line += var + ', '
433                else:
434                        line += 'none, '
435                print(line[:-2])
436            firstVaryList = newVaryList
437
438        ifPrint = False
439        printFile.write('\n Refinement results for histogram: %s\n'%histogram)
440        printFile.write(135*'-'+'\n')
441        if True:
442#        try:
443            IfOK,Rvals,result,covMatrix,sig = RefineCore(Controls,Histo,Phases,restraintDict,
444                rigidbodyDict,parmDict,varyList,calcControls,pawleyLookup,ifPrint,printFile,dlg)
445            if PlotFunction:
446                PlotFunction(G2frame,Histo[histogram]['Data'],histogram)
447
448            print ('  wR = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f, last delta chi = %.4f'%(
449                Rvals['Rwp'],Rvals['chisq'],Rvals['GOF']**2,Rvals['DelChi2']))
450            # add the uncertainties into the esd dictionary (sigDict)
451            sigDict = dict(zip(varyList,sig))
452            # the uncertainties for dependent constrained parms into the esd dict
453            sigDict.update(G2mv.ComputeDepESD(covMatrix,varyList,parmDict))
454
455            # a dict with values & esds for dependent (constrained) parameters - avoid extraneous holds
456            depParmDict = {i:(parmDict[i],sigDict[i]) for i in varyListStart if i in sigDict and i not in varyList}
457            newCellDict = copy.deepcopy(G2stMth.GetNewCellParms(parmDict,varyList))
458            newAtomDict = copy.deepcopy(G2stMth.ApplyXYZshifts(parmDict,varyList))
459            histRefData = {
460                'variables':result[0],'varyList':varyList,'sig':sig,'Rvals':Rvals,
461                'varyListStart':varyListStart,
462                'covMatrix':covMatrix,'title':histogram,'newAtomDict':newAtomDict,
463                'newCellDict':newCellDict,'depParmDict':depParmDict,
464                'constraintInfo':constraintInfo,
465                'parmDict':parmDict}
466            SeqResult[histogram] = histRefData
467            G2stMth.ApplyRBModels(parmDict,Phases,rigidbodyDict,True)
468    #        G2stIO.SetRigidBodyModels(parmDict,sigDict,rigidbodyDict,printFile)
469            G2stIO.SetHistogramPhaseData(parmDict,sigDict,Phases,Histo,None,ifPrint,printFile)
470            G2stIO.SetHistogramData(parmDict,sigDict,Histo,None,ifPrint,printFile)
471            G2stIO.SetUsedHistogramsAndPhases(GPXfile,Histo,Phases,rigidbodyDict,histRefData,makeBack)
472            makeBack = False
473            NewparmDict = {}
474            # make dict of varied parameters in current histogram, renamed to
475            # next histogram, for use in next refinement.
476            if Controls['Copy2Next'] and ihst < len(histNames)-1:
477                hId = Histo[histogram]['hId'] # current histogram
478                nexthId = Histograms[histNames[ihst+1]]['hId']
479                for parm in set(list(varyList)+list(varyListStart)):
480                    items = parm.split(':')
481                    if len(items) < 3: continue
482                    if str(hId) in items[1]:
483                        items[1] = str(nexthId)
484                        newparm = ':'.join(items)
485                        NewparmDict[newparm] = parmDict[parm]
486#        except G2obj.G2Exception(Msg):
487#            printFile.close()
488#            print (' ***** Refinement aborted *****')
489#            return False,Msg.msg
490    if Controls.get('Reverse Seq'):
491        histNames.reverse()
492    SeqResult['histNames'] = histNames
493    G2stIO.SetSeqResult(GPXfile,Histograms,SeqResult)
494    printFile.close()
495    print (' Sequential refinement results are in file: '+ospath.splitext(GPXfile)[0]+'.lst')
496    print (' ***** Sequential refinement successful *****')
497    return True,'Success'
498
499def RetDistAngle(DisAglCtls,DisAglData):
500    '''Compute and return distances and angles
501
502    :param dict DisAglCtls: contains distance/angle radii usually defined using
503       :func:`GSASIIctrlGUI.DisAglDialog`
504    :param dict DisAglData: contains phase data:
505       Items 'OrigAtoms' and 'TargAtoms' contain the atoms to be used
506       for distance/angle origins and atoms to be used as targets.
507       Item 'SGData' has the space group information (see :ref:`Space Group object<SGData_table>`)
508
509    :returns: AtomLabels,DistArray,AngArray where:
510
511      **AtomLabels** is a dict of atom labels, keys are the atom number
512
513      **DistArray** is a dict keyed by the origin atom number where the value is a list
514      of distance entries. The value for each distance is a list containing:
515
516        0) the target atom number (int);
517        1) the unit cell offsets added to x,y & z (tuple of int values)
518        2) the symmetry operator number (which may be modified to indicate centering and center of symmetry)
519        3) an interatomic distance in A (float)
520        4) an uncertainty on the distance in A or 0.0 (float)
521
522      **AngArray** is a dict keyed by the origin (central) atom number where
523      the value is a list of
524      angle entries. The value for each angle entry consists of three values:
525
526        0) a distance item reference for one neighbor (int)
527        1) a distance item reference for a second neighbor (int)
528        2) a angle, uncertainty pair; the s.u. may be zero (tuple of two floats)
529
530      The AngArray distance reference items refer directly to the index of the items in the
531      DistArray item for the list of distances for the central atom.
532    '''
533    import numpy.ma as ma
534
535    SGData = DisAglData['SGData']
536    Cell = DisAglData['Cell']
537
538    Amat,Bmat = G2lat.cell2AB(Cell[:6])
539    covData = {}
540    if 'covData' in DisAglData:
541        covData = DisAglData['covData']
542        covMatrix = covData['covMatrix']
543        varyList = covData['varyList']
544        pfx = str(DisAglData['pId'])+'::'
545
546    Factor = DisAglCtls['Factors']
547    Radii = dict(zip(DisAglCtls['AtomTypes'],zip(DisAglCtls['BondRadii'],DisAglCtls['AngleRadii'])))
548    indices = (-2,-1,0,1,2)
549    Units = np.array([[h,k,l] for h in indices for k in indices for l in indices])
550    origAtoms = DisAglData['OrigAtoms']
551    targAtoms = DisAglData['TargAtoms']
552    AtomLabels = {}
553    for Oatom in origAtoms:
554        AtomLabels[Oatom[0]] = Oatom[1]
555    for Oatom in targAtoms:
556        AtomLabels[Oatom[0]] = Oatom[1]
557    DistArray = {}
558    AngArray = {}
559    for Oatom in origAtoms:
560        DistArray[Oatom[0]] = []
561        AngArray[Oatom[0]] = []
562        OxyzNames = ''
563        IndBlist = []
564        Dist = []
565        Vect = []
566        VectA = []
567        angles = []
568        for Tatom in targAtoms:
569            Xvcov = []
570            TxyzNames = ''
571            if 'covData' in DisAglData:
572                OxyzNames = [pfx+'dAx:%d'%(Oatom[0]),pfx+'dAy:%d'%(Oatom[0]),pfx+'dAz:%d'%(Oatom[0])]
573                TxyzNames = [pfx+'dAx:%d'%(Tatom[0]),pfx+'dAy:%d'%(Tatom[0]),pfx+'dAz:%d'%(Tatom[0])]
574                Xvcov = G2mth.getVCov(OxyzNames+TxyzNames,varyList,covMatrix)
575            result = G2spc.GenAtom(Tatom[3:6],SGData,False,Move=False)
576            BsumR = (Radii[Oatom[2]][0]+Radii[Tatom[2]][0])*Factor[0]
577            AsumR = (Radii[Oatom[2]][1]+Radii[Tatom[2]][1])*Factor[1]
578            for Txyz,Top,Tunit in result:
579                Dx = (Txyz-np.array(Oatom[3:6]))+Units
580                dx = np.inner(Amat,Dx)
581                dist = ma.masked_less(np.sqrt(np.sum(dx**2,axis=0)),0.5)
582                IndB = ma.nonzero(ma.masked_greater(dist-BsumR,0.))
583                if np.any(IndB):
584                    for indb in IndB:
585                        for i in range(len(indb)):
586                            if str(dx.T[indb][i]) not in IndBlist:
587                                IndBlist.append(str(dx.T[indb][i]))
588                                unit = Units[indb][i]
589                                tunit = (unit[0]+Tunit[0],unit[1]+Tunit[1],unit[2]+Tunit[2])
590                                pdpx = G2mth.getDistDerv(Oatom[3:6],Tatom[3:6],Amat,unit,Top,SGData)
591                                sig = 0.0
592                                if len(Xvcov):
593                                    sig = np.sqrt(np.inner(pdpx,np.inner(pdpx,Xvcov)))
594                                Dist.append([Oatom[0],Tatom[0],tunit,Top,ma.getdata(dist[indb])[i],sig])
595                                if (Dist[-1][-2]-AsumR) <= 0.:
596                                    Vect.append(dx.T[indb][i]/Dist[-1][-2])
597                                    VectA.append([OxyzNames,np.array(Oatom[3:6]),TxyzNames,np.array(Tatom[3:6]),unit,Top])
598                                else:
599                                    Vect.append([0.,0.,0.])
600                                    VectA.append([])
601        for D in Dist:
602            DistArray[Oatom[0]].append(D[1:])
603        Vect = np.array(Vect)
604        angles = np.zeros((len(Vect),len(Vect)))
605        angsig = np.zeros((len(Vect),len(Vect)))
606        for i,veca in enumerate(Vect):
607            if np.any(veca):
608                for j,vecb in enumerate(Vect):
609                    if np.any(vecb):
610                        angles[i][j],angsig[i][j] = G2mth.getAngSig(VectA[i],VectA[j],Amat,SGData,covData)
611                        if i <= j: continue
612                        AngArray[Oatom[0]].append((i,j,
613                            G2mth.getAngSig(VectA[i],VectA[j],Amat,SGData,covData)))
614    return AtomLabels,DistArray,AngArray
615
616def PrintDistAngle(DisAglCtls,DisAglData,out=sys.stdout):
617    '''Print distances and angles
618
619    :param dict DisAglCtls: contains distance/angle radii usually defined using
620       :func:`GSASIIctrlGUI.DisAglDialog`
621    :param dict DisAglData: contains phase data:
622       Items 'OrigAtoms' and 'TargAtoms' contain the atoms to be used
623       for distance/angle origins and atoms to be used as targets.
624       Item 'SGData' has the space group information (see :ref:`Space Group object<SGData_table>`)
625    :param file out: file object for output. Defaults to sys.stdout.
626    '''
627    def MyPrint(s):
628        out.write(s+'\n')
629        # print(s,file=out) # use in Python 3
630
631    def ShowBanner(name):
632        MyPrint(80*'*')
633        MyPrint('   Interatomic Distances and Angles for phase '+name)
634        MyPrint((80*'*')+'\n')
635
636    ShowBanner(DisAglCtls['Name'])
637    SGData = DisAglData['SGData']
638    SGtext,SGtable = G2spc.SGPrint(SGData)
639    for line in SGtext: MyPrint(line)
640    if len(SGtable) > 1:
641        for i,item in enumerate(SGtable[::2]):
642            if 2*i+1 == len(SGtable):
643                line = ' %s'%(item.ljust(30))
644            else:
645                line = ' %s %s'%(item.ljust(30),SGtable[2*i+1].ljust(30))
646            MyPrint(line)
647    else:
648        MyPrint(' ( 1)    %s'%(SGtable[0])) #triclinic case
649    Cell = DisAglData['Cell']
650
651    Amat,Bmat = G2lat.cell2AB(Cell[:6])
652    covData = {}
653    if 'covData' in DisAglData:
654        covData = DisAglData['covData']
655        pfx = str(DisAglData['pId'])+'::'
656        A = G2lat.cell2A(Cell[:6])
657        cellSig = G2stIO.getCellEsd(pfx,SGData,A,covData)
658        names = [' a = ',' b = ',' c = ',' alpha = ',' beta = ',' gamma = ',' Volume = ']
659        valEsd = [G2mth.ValEsd(Cell[i],cellSig[i],True) for i in range(7)]
660        line = '\n Unit cell:'
661        for name,vals in zip(names,valEsd):
662            line += name+vals
663        MyPrint(line)
664    else:
665        MyPrint('\n Unit cell: a = '+('%.5f'%Cell[0])+' b = '+('%.5f'%Cell[1])+' c = '+('%.5f'%Cell[2])+
666            ' alpha = '+('%.3f'%Cell[3])+' beta = '+('%.3f'%Cell[4])+' gamma = '+
667            ('%.3f'%Cell[5])+' Volume = '+('%.3f'%Cell[6]))
668
669    AtomLabels,DistArray,AngArray = RetDistAngle(DisAglCtls,DisAglData)
670    origAtoms = DisAglData['OrigAtoms']
671    for Oatom in origAtoms:
672        i = Oatom[0]
673        Dist = DistArray[i]
674        nDist = len(Dist)
675        angles = np.zeros((nDist,nDist))
676        angsig = np.zeros((nDist,nDist))
677        for k,j,tup in AngArray[i]:
678            angles[k][j],angsig[k][j] = angles[j][k],angsig[j][k] = tup
679        line = ''
680        for i,x in enumerate(Oatom[3:6]):
681            line += ('%12.5f'%x).rstrip('0')
682        MyPrint('\n Distances & angles for '+Oatom[1]+' at '+line.rstrip())
683        MyPrint(80*'*')
684        line = ''
685        for dist in Dist[:-1]:
686            line += '%12s'%(AtomLabels[dist[0]].center(12))
687        MyPrint('  To       cell +(sym. op.)      dist.  '+line.rstrip())
688        for i,dist in enumerate(Dist):
689            line = ''
690            for j,angle in enumerate(angles[i][0:i]):
691                sig = angsig[i][j]
692                if angle:
693                    if sig:
694                        line += '%12s'%(G2mth.ValEsd(angle,sig,True).center(12))
695                    else:
696                        val = '%.3f'%(angle)
697                        line += '%12s'%(val.center(12))
698                else:
699                    line += 12*' '
700            if dist[4]:            #sig exists!
701                val = G2mth.ValEsd(dist[3],dist[4])
702            else:
703                val = '%8.4f'%(dist[3])
704            tunit = '[%2d%2d%2d]'% dist[1]
705            MyPrint((%8s%10s+(%4d) %12s'%(AtomLabels[dist[0]].ljust(8),tunit.ljust(10),dist[2],val.center(12)))+line.rstrip())
706
707def DisAglTor(DATData):
708    'Needs a doc string'
709    SGData = DATData['SGData']
710    Cell = DATData['Cell']
711
712    Amat,Bmat = G2lat.cell2AB(Cell[:6])
713    covData = {}
714    pfx = ''
715    if 'covData' in DATData:
716        covData = DATData['covData']
717        pfx = str(DATData['pId'])+'::'
718    Datoms = []
719    Oatoms = []
720    for i,atom in enumerate(DATData['Datoms']):
721        symop = atom[-1].split('+')
722        if len(symop) == 1:
723            symop.append('0,0,0')
724        symop[0] = int(symop[0])
725        symop[1] = eval(symop[1])
726        atom.append(symop)
727        Datoms.append(atom)
728        oatom = DATData['Oatoms'][i]
729        names = ['','','']
730        if pfx:
731            names = [pfx+'dAx:'+str(oatom[0]),pfx+'dAy:'+str(oatom[0]),pfx+'dAz:'+str(oatom[0])]
732        oatom += [names,]
733        Oatoms.append(oatom)
734    atmSeq = [atom[1]+'('+atom[-2]+')' for atom in Datoms]
735    if DATData['Natoms'] == 4:  #torsion
736        Tors,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData)
737        print (' Torsion angle for %s atom sequence: %s = %s'%(DATData['Name'],str(atmSeq).replace("'","")[1:-1],G2mth.ValEsd(Tors,sig)))
738        print (' NB: Atom sequence determined by selection order')
739        return      # done with torsion
740    elif DATData['Natoms'] == 3:  #angle
741        Ang,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData)
742        print (' Angle in %s for atom sequence: %s = %s'%(DATData['Name'],str(atmSeq).replace("'","")[1:-1],G2mth.ValEsd(Ang,sig)))
743        print (' NB: Atom sequence determined by selection order')
744    else:   #2 atoms - distance
745        Dist,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData)
746        print (' Distance in %s for atom sequence: %s = %s'%(DATData['Name'],str(atmSeq).replace("'","")[1:-1],G2mth.ValEsd(Dist,sig)))
747
748def BestPlane(PlaneData):
749    'Needs a doc string'
750
751    def ShowBanner(name):
752        print (80*'*')
753        print ('   Best plane result for phase '+name)
754        print (80*'*','\n')
755
756    ShowBanner(PlaneData['Name'])
757
758    Cell = PlaneData['Cell']
759    Amat,Bmat = G2lat.cell2AB(Cell[:6])
760    Atoms = PlaneData['Atoms']
761    sumXYZ = np.zeros(3)
762    XYZ = []
763    Natoms = len(Atoms)
764    for atom in Atoms:
765        xyz = np.array(atom[3:6])
766        XYZ.append(xyz)
767        sumXYZ += xyz
768    sumXYZ /= Natoms
769    XYZ = np.array(XYZ)-sumXYZ
770    XYZ = np.inner(Amat,XYZ).T
771    Zmat = np.zeros((3,3))
772    for i,xyz in enumerate(XYZ):
773        Zmat += np.outer(xyz.T,xyz)
774    print (' Selected atoms centered at %10.5f %10.5f %10.5f'%(sumXYZ[0],sumXYZ[1],sumXYZ[2]))
775    Evec,Emat = nl.eig(Zmat)
776    Evec = np.sqrt(Evec)/(Natoms-3)
777    Order = np.argsort(Evec)
778    XYZ = np.inner(XYZ,Emat.T).T
779    XYZ = np.array([XYZ[Order[2]],XYZ[Order[1]],XYZ[Order[0]]]).T
780    print (' Atoms in Cartesian best plane coordinates:')
781    print (' Name         X         Y         Z')
782    for i,xyz in enumerate(XYZ):
783        print (' %6s%10.3f%10.3f%10.3f'%(Atoms[i][1].ljust(6),xyz[0],xyz[1],xyz[2]))
784    print ('\n Best plane RMS X =%8.3f, Y =%8.3f, Z =%8.3f'%(Evec[Order[2]],Evec[Order[1]],Evec[Order[0]]))
785
786def main():
787    'Called to run a refinement when this module is executed '
788    starttime = time.time()
789    arg = sys.argv
790    if len(arg) > 1:
791        GPXfile = arg[1]
792        if not ospath.exists(GPXfile):
793            print ('ERROR - '+GPXfile+" doesn't exist!")
794            exit()
795    else:
796        print ('ERROR - missing filename')
797        exit()
798    # TODO: figure out if this is a sequential refinement and call SeqRefine(GPXfile,None)
799    Refine(GPXfile,None)
800    print("Done. Execution time {:.2f} sec.".format(time.time()-starttime))
801
802if __name__ == '__main__':
803    GSASIIpath.InvokeDebugOpts()
804    main()
Note: See TracBrowser for help on using the repository browser.