source: trunk/GSASIIstrMain.py @ 3686

Last change on this file since 3686 was 3579, checked in by vondreele, 7 years ago

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