source: trunk/GSASIIstrMain.py @ 3068

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

fix parameter display to show nan (crashed before)
trap with error if any parameter is nan
skip postprocessing for failed refinements (IfOK=False)
fix Jacobian to use GetProfileDervMP

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 35.2 KB
Line 
1# -*- coding: utf-8 -*-
2'''
3*GSASIIstrMain: main structure routine*
4---------------------------------------
5
6'''
7########### SVN repository information ###################
8# $Date: 2017-09-12 15:54:45 +0000 (Tue, 12 Sep 2017) $
9# $Author: vondreele $
10# $Revision: 3068 $
11# $URL: trunk/GSASIIstrMain.py $
12# $Id: GSASIIstrMain.py 3068 2017-09-12 15:54:45Z 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: 3068 $")
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    '''
264    pIds = []
265    for phase in Phases:
266        if histogram not in Phases[phase]['Histograms']: continue
267        if Phases[phase]['Histograms'][histogram]['Use']:
268            pIds.append(str(Phases[phase]['pId']))
269    return [item for item in phaseVary if item.split(':')[0] in pIds]
270
271def SeqRefine(GPXfile,dlg,PlotFunction=None,G2frame=None):
272    '''Perform a sequential refinement -- cycles through all selected histgrams,
273    one at a time
274    '''
275    import GSASIImpsubs as G2mp
276    G2mp.InitMP()
277    import pytexture as ptx
278    ptx.pyqlmninit()            #initialize fortran arrays for spherical harmonics
279
280    printFile = open(ospath.splitext(GPXfile)[0]+'.lst','w')
281    print 'Starting Sequential Refinement'
282    G2stIO.ShowBanner(printFile)
283    Controls = G2stIO.GetControls(GPXfile)
284    G2stIO.ShowControls(Controls,printFile,SeqRef=True)
285    restraintDict = G2stIO.GetRestraints(GPXfile)
286    Histograms,Phases = G2stIO.GetUsedHistogramsAndPhases(GPXfile)
287    if not Phases:
288        print ' *** ERROR - you have no phases to refine! ***'
289        print ' *** Refine aborted ***'
290        return False,'No phases'
291    if not Histograms:
292        print ' *** ERROR - you have no data to refine with! ***'
293        print ' *** Refine aborted ***'
294        return False,'No data'
295    rigidbodyDict = G2stIO.GetRigidBodies(GPXfile)
296    rbIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})
297    rbVary,rbDict = G2stIO.GetRigidBodyModels(rigidbodyDict,pFile=printFile)
298    Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables,MFtables,maxSSwave = \
299        G2stIO.GetPhaseData(Phases,restraintDict,rbIds,False,printFile,seqRef=True)
300    for item in phaseVary:
301        if '::A0' in item:
302            print '**** WARNING - lattice parameters should not be refined in a sequential refinement ****'
303            print '****           instead use the Dij parameters for each powder histogram            ****'
304            return False,'Lattice parameter refinement error - see console message'
305        if '::C(' in item:
306            print '**** WARNING - phase texture parameters should not be refined in a sequential refinement ****'
307            print '****           instead use the C(L,N) parameters for each powder histogram               ****'
308            return False,'Phase texture refinement error - see console message'
309    if 'Seq Data' in Controls:
310        histNames = Controls['Seq Data']
311    else:
312        histNames = G2stIO.GetHistogramNames(GPXfile,['PWDR',])
313    if Controls.get('Reverse Seq'):
314        histNames.reverse()
315    SeqResult = G2stIO.GetSeqResult(GPXfile)
316#    SeqResult = {'SeqPseudoVars':{},'SeqParFitEqList':[]}
317    makeBack = True
318    Histo = {}
319    NewparmDict = {}
320    for ihst,histogram in enumerate(histNames):
321        print('\nRefining with '+str(histogram))
322        ifPrint = False
323        if dlg:
324            dlg.SetTitle('Residual for histogram '+str(ihst))
325        calcControls = {}
326        calcControls['atomIndx'] = atomIndx
327        calcControls['Natoms'] = Natoms
328        calcControls['FFtables'] = FFtables
329        calcControls['BLtables'] = BLtables
330        calcControls['MFtables'] = MFtables
331        calcControls['maxSSwave'] = maxSSwave
332        if histogram not in Histograms:
333            print("Error: not found!")
334            continue
335        redphaseVary = phaseCheck(phaseVary,Phases,histogram)
336        Histo = {histogram:Histograms[histogram],}
337        hapVary,hapDict,controlDict = G2stIO.GetHistogramPhaseData(Phases,Histo,Print=False)
338        calcControls.update(controlDict)
339        histVary,histDict,controlDict = G2stIO.GetHistogramData(Histo,False)
340        calcControls.update(controlDict)
341        varyList = rbVary+redphaseVary+hapVary+histVary
342#        if not ihst:
343            # save the initial vary list, but without histogram numbers on parameters
344        saveVaryList = varyList[:]
345        for i,item in enumerate(saveVaryList):
346            items = item.split(':')
347            if items[1]:
348                items[1] = ''
349            item = ':'.join(items)
350            saveVaryList[i] = item
351        if not ihst:
352            SeqResult['varyList'] = saveVaryList
353        else:
354            SeqResult['varyList'] = list(set(SeqResult['varyList']+saveVaryList))
355        parmDict = {}
356        parmDict.update(phaseDict)
357        parmDict.update(hapDict)
358        parmDict.update(histDict)
359        if Controls['Copy2Next']:
360            #parmDict.update(NewparmDict) # don't use in case extra entries would cause a problem
361            for parm in NewparmDict:
362                if parm in parmDict:
363                    parmDict[parm] = NewparmDict[parm]
364        G2stIO.GetFprime(calcControls,Histo)
365        # do constraint processing
366        #reload(G2mv) # debug
367        G2mv.InitVars()
368        constrDict,fixedList = G2stIO.GetConstraints(GPXfile)
369        varyListStart = tuple(varyList) # save the original varyList before dependent vars are removed
370        try:
371            groups,parmlist = G2mv.GroupConstraints(constrDict)
372            G2mv.GenerateConstraints(groups,parmlist,varyList,constrDict,fixedList,parmDict,SeqHist=ihst)
373            constraintInfo = (groups,parmlist,constrDict,fixedList,ihst)
374        except:
375            print ' *** ERROR - your constraints are internally inconsistent ***'
376            #errmsg, warnmsg = G2mv.CheckConstraints(varyList,constrDict,fixedList)
377            #print 'Errors',errmsg
378            #if warnmsg: print 'Warnings',warnmsg
379            return False,' Constraint error'
380        #print G2mv.VarRemapShow(varyList)
381        if not ihst:
382            # first histogram to refine against
383            firstVaryList = []
384            for item in varyList:
385                items = item.split(':')
386                if items[1]:
387                    items[1] = ''
388                item = ':'.join(items)
389                firstVaryList.append(item)
390            newVaryList = firstVaryList
391        else:
392            newVaryList = []
393            for item in varyList:
394                items = item.split(':')
395                if items[1]:
396                    items[1] = ''
397                item = ':'.join(items)
398                newVaryList.append(item)
399        if newVaryList != firstVaryList and Controls['Copy2Next']:
400            # variable lists are expected to match between sequential refinements when Copy2Next is on
401            #print '**** ERROR - variable list for this histogram does not match previous'
402            #print '     Copy of variables is not possible'
403            #print '\ncurrent histogram',histogram,'has',len(newVaryList),'variables'
404            combined = list(set(firstVaryList+newVaryList))
405            c = [var for var in combined if var not in newVaryList]
406            p = [var for var in combined if var not in firstVaryList]
407            print('*** Variables change ***')
408            for typ,vars in [('Removed',c),('Added',p)]:
409                line = '  '+typ+': '
410                if vars:
411                    for var in vars:
412                        if len(line) > 70:
413                            print(line)
414                            line = '    '
415                        line += var + ', '
416                else:
417                        line += 'none, '
418                print(line[:-2])
419            firstVaryList = newVaryList
420
421        ifPrint = False
422        print >>printFile,'\n Refinement results for histogram: v'+histogram
423        print >>printFile,135*'-'
424        try:
425            IfOK,Rvals,result,covMatrix,sig = RefineCore(Controls,Histo,Phases,restraintDict,
426                rigidbodyDict,parmDict,varyList,calcControls,pawleyLookup,ifPrint,printFile,dlg)
427            if PlotFunction:
428                PlotFunction(G2frame,Histo[histogram]['Data'],histogram)
429
430            print '  wR = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f, last delta chi = %.4f'%(
431                Rvals['Rwp'],Rvals['chisq'],Rvals['GOF']**2,Rvals['DelChi2'])
432            # add the uncertainties into the esd dictionary (sigDict)
433            sigDict = dict(zip(varyList,sig))
434            # the uncertainties for dependent constrained parms into the esd dict
435            sigDict.update(G2mv.ComputeDepESD(covMatrix,varyList,parmDict))
436
437            # a dict with values & esds for dependent (constrained) parameters - avoid extraneous holds
438            depParmDict = {i:(parmDict[i],sigDict[i]) for i in varyListStart if i in sigDict and i not in varyList}
439            newCellDict = copy.deepcopy(G2stMth.GetNewCellParms(parmDict,varyList))
440            newAtomDict = copy.deepcopy(G2stMth.ApplyXYZshifts(parmDict,varyList))
441            histRefData = {
442                'variables':result[0],'varyList':varyList,'sig':sig,'Rvals':Rvals,
443                'varyListStart':varyListStart,
444                'covMatrix':covMatrix,'title':histogram,'newAtomDict':newAtomDict,
445                'newCellDict':newCellDict,'depParmDict':depParmDict,
446                'constraintInfo':constraintInfo,
447                'parmDict':parmDict}
448            SeqResult[histogram] = histRefData
449            G2stMth.ApplyRBModels(parmDict,Phases,rigidbodyDict,True)
450    #        G2stIO.SetRigidBodyModels(parmDict,sigDict,rigidbodyDict,printFile)
451            G2stIO.SetHistogramPhaseData(parmDict,sigDict,Phases,Histo,None,ifPrint,printFile)
452            G2stIO.SetHistogramData(parmDict,sigDict,Histo,None,ifPrint,printFile)
453            G2stIO.SetUsedHistogramsAndPhases(GPXfile,Histo,Phases,rigidbodyDict,histRefData,makeBack)
454            makeBack = False
455            NewparmDict = {}
456            # make dict of varied parameters in current histogram, renamed to
457            # next histogram, for use in next refinement.
458            if Controls['Copy2Next'] and ihst < len(histNames)-1:
459                hId = Histo[histogram]['hId'] # current histogram
460                nexthId = Histograms[histNames[ihst+1]]['hId']
461                for parm in set(list(varyList)+list(varyListStart)):
462                    items = parm.split(':')
463                    if len(items) < 3: continue
464                    if str(hId) in items[1]:
465                        items[1] = str(nexthId)
466                        newparm = ':'.join(items)
467                        NewparmDict[newparm] = parmDict[parm]
468        except G2obj.G2Exception,Msg:
469            printFile.close()
470            print ' ***** Refinement aborted *****'
471            return False,Msg.msg
472    if Controls.get('Reverse Seq'):
473        histNames.reverse()
474    SeqResult['histNames'] = histNames
475    G2stIO.SetSeqResult(GPXfile,Histograms,SeqResult)
476    printFile.close()
477    print ' Sequential refinement results are in file: '+ospath.splitext(GPXfile)[0]+'.lst'
478    print ' ***** Sequential refinement successful *****'
479    return True,'Success'
480
481def RetDistAngle(DisAglCtls,DisAglData):
482    '''Compute and return distances and angles
483
484    :param dict DisAglCtls: contains distance/angle radii usually defined using
485       :func:`GSASIIctrlGUI.DisAglDialog`
486    :param dict DisAglData: contains phase data:
487       Items 'OrigAtoms' and 'TargAtoms' contain the atoms to be used
488       for distance/angle origins and atoms to be used as targets.
489       Item 'SGData' has the space group information (see :ref:`Space Group object<SGData_table>`)
490
491    :returns: AtomLabels,DistArray,AngArray where:
492
493      **AtomLabels** is a dict of atom labels, keys are the atom number
494
495      **DistArray** is a dict keyed by the origin atom number where the value is a list
496      of distance entries. The value for each distance is a list containing:
497
498        0) the target atom number (int);
499        1) the unit cell offsets added to x,y & z (tuple of int values)
500        2) the symmetry operator number (which may be modified to indicate centering and center of symmetry)
501        3) an interatomic distance in A (float)
502        4) an uncertainty on the distance in A or 0.0 (float)
503
504      **AngArray** is a dict keyed by the origin (central) atom number where
505      the value is a list of
506      angle entries. The value for each angle entry consists of three values:
507
508        0) a distance item reference for one neighbor (int)
509        1) a distance item reference for a second neighbor (int)
510        2) a angle, uncertainty pair; the s.u. may be zero (tuple of two floats)
511
512      The AngArray distance reference items refer directly to the index of the items in the
513      DistArray item for the list of distances for the central atom.
514    '''
515    import numpy.ma as ma
516
517    SGData = DisAglData['SGData']
518    Cell = DisAglData['Cell']
519
520    Amat,Bmat = G2lat.cell2AB(Cell[:6])
521    covData = {}
522    if 'covData' in DisAglData:
523        covData = DisAglData['covData']
524        covMatrix = covData['covMatrix']
525        varyList = covData['varyList']
526        pfx = str(DisAglData['pId'])+'::'
527
528    Factor = DisAglCtls['Factors']
529    Radii = dict(zip(DisAglCtls['AtomTypes'],zip(DisAglCtls['BondRadii'],DisAglCtls['AngleRadii'])))
530    indices = (-2,-1,0,1,2)
531    Units = np.array([[h,k,l] for h in indices for k in indices for l in indices])
532    origAtoms = DisAglData['OrigAtoms']
533    targAtoms = DisAglData['TargAtoms']
534    AtomLabels = {}
535    for Oatom in origAtoms:
536        AtomLabels[Oatom[0]] = Oatom[1]
537    for Oatom in targAtoms:
538        AtomLabels[Oatom[0]] = Oatom[1]
539    DistArray = {}
540    AngArray = {}
541    for Oatom in origAtoms:
542        DistArray[Oatom[0]] = []
543        AngArray[Oatom[0]] = []
544        OxyzNames = ''
545        IndBlist = []
546        Dist = []
547        Vect = []
548        VectA = []
549        angles = []
550        for Tatom in targAtoms:
551            Xvcov = []
552            TxyzNames = ''
553            if 'covData' in DisAglData:
554                OxyzNames = [pfx+'dAx:%d'%(Oatom[0]),pfx+'dAy:%d'%(Oatom[0]),pfx+'dAz:%d'%(Oatom[0])]
555                TxyzNames = [pfx+'dAx:%d'%(Tatom[0]),pfx+'dAy:%d'%(Tatom[0]),pfx+'dAz:%d'%(Tatom[0])]
556                Xvcov = G2mth.getVCov(OxyzNames+TxyzNames,varyList,covMatrix)
557            result = G2spc.GenAtom(Tatom[3:6],SGData,False,Move=False)
558            BsumR = (Radii[Oatom[2]][0]+Radii[Tatom[2]][0])*Factor[0]
559            AsumR = (Radii[Oatom[2]][1]+Radii[Tatom[2]][1])*Factor[1]
560            for Txyz,Top,Tunit in result:
561                Dx = (Txyz-np.array(Oatom[3:6]))+Units
562                dx = np.inner(Amat,Dx)
563                dist = ma.masked_less(np.sqrt(np.sum(dx**2,axis=0)),0.5)
564                IndB = ma.nonzero(ma.masked_greater(dist-BsumR,0.))
565                if np.any(IndB):
566                    for indb in IndB:
567                        for i in range(len(indb)):
568                            if str(dx.T[indb][i]) not in IndBlist:
569                                IndBlist.append(str(dx.T[indb][i]))
570                                unit = Units[indb][i]
571                                tunit = (unit[0]+Tunit[0],unit[1]+Tunit[1],unit[2]+Tunit[2])
572                                pdpx = G2mth.getDistDerv(Oatom[3:6],Tatom[3:6],Amat,unit,Top,SGData)
573                                sig = 0.0
574                                if len(Xvcov):
575                                    sig = np.sqrt(np.inner(pdpx,np.inner(pdpx,Xvcov)))
576                                Dist.append([Oatom[0],Tatom[0],tunit,Top,ma.getdata(dist[indb])[i],sig])
577                                if (Dist[-1][-2]-AsumR) <= 0.:
578                                    Vect.append(dx.T[indb][i]/Dist[-1][-2])
579                                    VectA.append([OxyzNames,np.array(Oatom[3:6]),TxyzNames,np.array(Tatom[3:6]),unit,Top])
580                                else:
581                                    Vect.append([0.,0.,0.])
582                                    VectA.append([])
583        for D in Dist:
584            DistArray[Oatom[0]].append(D[1:])
585        Vect = np.array(Vect)
586        angles = np.zeros((len(Vect),len(Vect)))
587        angsig = np.zeros((len(Vect),len(Vect)))
588        for i,veca in enumerate(Vect):
589            if np.any(veca):
590                for j,vecb in enumerate(Vect):
591                    if np.any(vecb):
592                        angles[i][j],angsig[i][j] = G2mth.getAngSig(VectA[i],VectA[j],Amat,SGData,covData)
593                        if i <= j: continue
594                        AngArray[Oatom[0]].append((i,j,
595                            G2mth.getAngSig(VectA[i],VectA[j],Amat,SGData,covData)))
596    return AtomLabels,DistArray,AngArray
597
598def PrintDistAngle(DisAglCtls,DisAglData,out=sys.stdout):
599    '''Print distances and angles
600
601    :param dict DisAglCtls: contains distance/angle radii usually defined using
602       :func:`GSASIIctrlGUI.DisAglDialog`
603    :param dict DisAglData: contains phase data:
604       Items 'OrigAtoms' and 'TargAtoms' contain the atoms to be used
605       for distance/angle origins and atoms to be used as targets.
606       Item 'SGData' has the space group information (see :ref:`Space Group object<SGData_table>`)
607    :param file out: file object for output. Defaults to sys.stdout.
608    '''
609    def MyPrint(s):
610        out.write(s+'\n')
611        # print(s,file=out) # use in Python 3
612
613    def ShowBanner(name):
614        MyPrint(80*'*')
615        MyPrint('   Interatomic Distances and Angles for phase '+name)
616        MyPrint((80*'*')+'\n')
617
618    ShowBanner(DisAglCtls['Name'])
619    SGData = DisAglData['SGData']
620    SGtext,SGtable = G2spc.SGPrint(SGData)
621    for line in SGtext: MyPrint(line)
622    if len(SGtable) > 1:
623        for i,item in enumerate(SGtable[::2]):
624            if 2*i+1 == len(SGtable):
625                line = ' %s'%(item.ljust(30))
626            else:
627                line = ' %s %s'%(item.ljust(30),SGtable[2*i+1].ljust(30))
628            MyPrint(line)
629    else:
630        MyPrint(' ( 1)    %s'%(SGtable[0])) #triclinic case
631    Cell = DisAglData['Cell']
632
633    Amat,Bmat = G2lat.cell2AB(Cell[:6])
634    covData = {}
635    if 'covData' in DisAglData:
636        covData = DisAglData['covData']
637        pfx = str(DisAglData['pId'])+'::'
638        A = G2lat.cell2A(Cell[:6])
639        cellSig = G2stIO.getCellEsd(pfx,SGData,A,covData)
640        names = [' a = ',' b = ',' c = ',' alpha = ',' beta = ',' gamma = ',' Volume = ']
641        valEsd = [G2mth.ValEsd(Cell[i],cellSig[i],True) for i in range(7)]
642        line = '\n Unit cell:'
643        for name,vals in zip(names,valEsd):
644            line += name+vals
645        MyPrint(line)
646    else:
647        MyPrint('\n Unit cell: a = '+('%.5f'%Cell[0])+' b = '+('%.5f'%Cell[1])+' c = '+('%.5f'%Cell[2])+
648            ' alpha = '+('%.3f'%Cell[3])+' beta = '+('%.3f'%Cell[4])+' gamma = '+
649            ('%.3f'%Cell[5])+' Volume = '+('%.3f'%Cell[6]))
650
651    AtomLabels,DistArray,AngArray = RetDistAngle(DisAglCtls,DisAglData)
652    origAtoms = DisAglData['OrigAtoms']
653    for Oatom in origAtoms:
654        i = Oatom[0]
655        Dist = DistArray[i]
656        nDist = len(Dist)
657        angles = np.zeros((nDist,nDist))
658        angsig = np.zeros((nDist,nDist))
659        for k,j,tup in AngArray[i]:
660            angles[k][j],angsig[k][j] = angles[j][k],angsig[j][k] = tup
661        line = ''
662        for i,x in enumerate(Oatom[3:6]):
663            line += ('%12.5f'%x).rstrip('0')
664        MyPrint('\n Distances & angles for '+Oatom[1]+' at '+line.rstrip())
665        MyPrint(80*'*')
666        line = ''
667        for dist in Dist[:-1]:
668            line += '%12s'%(AtomLabels[dist[0]].center(12))
669        MyPrint('  To       cell +(sym. op.)      dist.  '+line.rstrip())
670        for i,dist in enumerate(Dist):
671            line = ''
672            for j,angle in enumerate(angles[i][0:i]):
673                sig = angsig[i][j]
674                if angle:
675                    if sig:
676                        line += '%12s'%(G2mth.ValEsd(angle,sig,True).center(12))
677                    else:
678                        val = '%.3f'%(angle)
679                        line += '%12s'%(val.center(12))
680                else:
681                    line += 12*' '
682            if dist[4]:            #sig exists!
683                val = G2mth.ValEsd(dist[3],dist[4])
684            else:
685                val = '%8.4f'%(dist[3])
686            tunit = '[%2d%2d%2d]'% dist[1]
687            MyPrint((%8s%10s+(%4d) %12s'%(AtomLabels[dist[0]].ljust(8),tunit.ljust(10),dist[2],val.center(12)))+line.rstrip())
688
689def DisAglTor(DATData):
690    'Needs a doc string'
691    SGData = DATData['SGData']
692    Cell = DATData['Cell']
693
694    Amat,Bmat = G2lat.cell2AB(Cell[:6])
695    covData = {}
696    pfx = ''
697    if 'covData' in DATData:
698        covData = DATData['covData']
699        pfx = str(DATData['pId'])+'::'
700    Datoms = []
701    Oatoms = []
702    for i,atom in enumerate(DATData['Datoms']):
703        symop = atom[-1].split('+')
704        if len(symop) == 1:
705            symop.append('0,0,0')
706        symop[0] = int(symop[0])
707        symop[1] = eval(symop[1])
708        atom.append(symop)
709        Datoms.append(atom)
710        oatom = DATData['Oatoms'][i]
711        names = ['','','']
712        if pfx:
713            names = [pfx+'dAx:'+str(oatom[0]),pfx+'dAy:'+str(oatom[0]),pfx+'dAz:'+str(oatom[0])]
714        oatom += [names,]
715        Oatoms.append(oatom)
716    atmSeq = [atom[1]+'('+atom[-2]+')' for atom in Datoms]
717    if DATData['Natoms'] == 4:  #torsion
718        Tors,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData)
719        print ' Torsion angle for '+DATData['Name']+' atom sequence: ',atmSeq,'=',G2mth.ValEsd(Tors,sig)
720        print ' NB: Atom sequence determined by selection order'
721        return      # done with torsion
722    elif DATData['Natoms'] == 3:  #angle
723        Ang,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData)
724        print ' Angle in '+DATData['Name']+' for atom sequence: ',atmSeq,'=',G2mth.ValEsd(Ang,sig)
725        print ' NB: Atom sequence determined by selection order'
726    else:   #2 atoms - distance
727        Dist,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData)
728        print ' Distance in '+DATData['Name']+' for atom sequence: ',atmSeq,'=',G2mth.ValEsd(Dist,sig)
729
730def BestPlane(PlaneData):
731    'Needs a doc string'
732
733    def ShowBanner(name):
734        print 80*'*'
735        print '   Best plane result for phase '+name
736        print 80*'*','\n'
737
738    ShowBanner(PlaneData['Name'])
739
740    Cell = PlaneData['Cell']
741    Amat,Bmat = G2lat.cell2AB(Cell[:6])
742    Atoms = PlaneData['Atoms']
743    sumXYZ = np.zeros(3)
744    XYZ = []
745    Natoms = len(Atoms)
746    for atom in Atoms:
747        xyz = np.array(atom[3:6])
748        XYZ.append(xyz)
749        sumXYZ += xyz
750    sumXYZ /= Natoms
751    XYZ = np.array(XYZ)-sumXYZ
752    XYZ = np.inner(Amat,XYZ).T
753    Zmat = np.zeros((3,3))
754    for i,xyz in enumerate(XYZ):
755        Zmat += np.outer(xyz.T,xyz)
756    print ' Selected atoms centered at %10.5f %10.5f %10.5f'%(sumXYZ[0],sumXYZ[1],sumXYZ[2])
757    Evec,Emat = nl.eig(Zmat)
758    Evec = np.sqrt(Evec)/(Natoms-3)
759    Order = np.argsort(Evec)
760    XYZ = np.inner(XYZ,Emat.T).T
761    XYZ = np.array([XYZ[Order[2]],XYZ[Order[1]],XYZ[Order[0]]]).T
762    print ' Atoms in Cartesian best plane coordinates:'
763    print ' Name         X         Y         Z'
764    for i,xyz in enumerate(XYZ):
765        print ' %6s%10.3f%10.3f%10.3f'%(Atoms[i][1].ljust(6),xyz[0],xyz[1],xyz[2])
766    print '\n Best plane RMS X =%8.3f, Y =%8.3f, Z =%8.3f'%(Evec[Order[2]],Evec[Order[1]],Evec[Order[0]])
767
768def main():
769    'Called to run a refinement when this module is executed '
770    starttime = time.time()
771    arg = sys.argv
772    if len(arg) > 1:
773        GPXfile = arg[1]
774        if not ospath.exists(GPXfile):
775            print 'ERROR - ',GPXfile," doesn't exist!"
776            exit()
777    else:
778        print 'ERROR - missing filename'
779        exit()
780    # TODO: figure out if this is a sequential refinement and call SeqRefine(GPXfile,None)
781    Refine(GPXfile,None)
782    print("Done. Execution time {:.2f} sec.".format(time.time()-starttime))
783
784if __name__ == '__main__':
785    GSASIIpath.InvokeDebugOpts()
786    main()
Note: See TracBrowser for help on using the repository browser.