source: trunk/GSASIIstrMain.py @ 4009

Last change on this file since 4009 was 4009, checked in by toby, 4 years ago

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