source: trunk/GSASIIstrMain.py @ 3095

Last change on this file since 3095 was 3095, checked in by vondreele, 4 years ago

add new confit parm: Show_timing
enable editing of the Use column in seq results table; doesn't seem to affect plot tho.
always print Found SVD zeros in seq refinement as well
add local directory to path for finding fit equation functions
add a pretranslate vector to atom transformations - still needs work tho.
couple of TODOs in G2strMain for fix parameters in seq refinements
change 'debug' to 'Show_timing' in G2strMath

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