source: trunk/GSASIIstrMain.py @ 4699

Last change on this file since 4699 was 4699, checked in by toby, 10 months ago

improve refinement warnings

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 51.7 KB
Line 
1# -*- coding: utf-8 -*-
2'''
3*GSASIIstrMain: main structure routine*
4---------------------------------------
5
6'''
7########### SVN repository information ###################
8# $Date: 2021-01-02 17:46:35 +0000 (Sat, 02 Jan 2021) $
9# $Author: toby $
10# $Revision: 4699 $
11# $URL: trunk/GSASIIstrMain.py $
12# $Id: GSASIIstrMain.py 4699 2021-01-02 17:46:35Z toby $
13########### SVN repository information ###################
14from __future__ import division, print_function
15import platform
16import sys
17import os.path as ospath
18import time
19import math
20import copy
21if '2' in platform.python_version_tuple()[0]:
22    import cPickle
23else:
24    import pickle as cPickle
25import numpy as np
26import numpy.linalg as nl
27import scipy.optimize as so
28import GSASIIpath
29GSASIIpath.SetBinaryPath()
30GSASIIpath.SetVersionNumber("$Revision: 4699 $")
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
38import GSASIIfiles as G2fil
39
40sind = lambda x: np.sin(x*np.pi/180.)
41cosd = lambda x: np.cos(x*np.pi/180.)
42tand = lambda x: np.tan(x*np.pi/180.)
43asind = lambda x: 180.*np.arcsin(x)/np.pi
44acosd = lambda x: 180.*np.arccos(x)/np.pi
45atan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
46
47ateln2 = 8.0*math.log(2.0)
48DEBUG = True
49
50def ReportProblems(result,Rvals,varyList):
51    '''Create a message based results from the refinement
52    '''
53    #report on SVD 0's and highly correlated variables
54    msg = ''
55    SVD0 = result[2].get('SVD0')
56    if SVD0 > 0: 
57        msg += 'Warning: There were {} singularities in the Hessian'.format(SVD0)
58    # process singular variables
59    psing = result[2].get('psing',[])
60    if len(psing):
61        if msg: msg += '\n'
62        m = '{} Parameters dropped due to singularities:'.format(len(psing))
63        msg += m
64        G2fil.G2Print(m, mode='warn')
65        m = ''
66        for i,val in enumerate(psing):
67            if i == 0:
68                msg += '\n{}'.format(varyList[val])
69                m = {}'.format(varyList[val])
70            else:
71                if len(m) > 70:
72                    G2fil.G2Print(m, mode='warn')
73                    m = '  '
74                else:
75                    m += ', '
76                m += '{}'.format(varyList[val])
77                if i == 10:
78                    msg += ', {}... see console for full list'.format(varyList[val])
79                elif i > 10:
80                    pass
81                else:
82                    msg += ', {}'.format(varyList[val])
83        if m: G2fil.G2Print(m, mode='warn')
84    #report on highly correlated variables
85    Hcorr = result[2].get('Hcorr',[])
86    if len(Hcorr) > 0: 
87        if msg: msg += '\n'
88        m = 'Note highly correlated parameters:'
89        G2fil.G2Print(m, mode='warn')
90        msg += m
91    elif SVD0 > 0:
92        if msg: msg += '\n'
93        m = 'Check covariance matrix for parameter correlation'
94        G2fil.G2Print(m, mode='warn')
95        msg += m
96    for i,(v1,v2,corr) in enumerate(Hcorr):
97        if corr > .95:
98            stars = '**'
99        else:
100            stars = '   '
101        m = ' {} {} and {} (@{:.2f}%)'.format(
102            stars,varyList[v1],varyList[v2],100.*corr)
103        G2fil.G2Print(m, mode='warn')
104        if i == 5:
105            msg += '\n' + m
106            msg += '\n   ... check console for more'
107        elif i < 5:
108            msg += '\n' + m
109    if msg:
110        if 'msg' not in Rvals: Rvals['msg'] = ''
111        Rvals['msg'] += msg
112
113def RefineCore(Controls,Histograms,Phases,restraintDict,rigidbodyDict,parmDict,varyList,
114    calcControls,pawleyLookup,ifSeq,printFile,dlg,refPlotUpdate=None):
115    '''Core optimization routines, shared between SeqRefine and Refine
116
117    :returns: 5-tuple of ifOk (bool), Rvals (dict), result, covMatrix, sig
118    '''
119    #patch (added Oct 2020) convert variable names for parm limits to G2VarObj
120    import GSASIIscriptable as G2sc
121    G2sc.patchControls(Controls)
122    # end patch
123#    print 'current',varyList
124#    for item in parmDict: print item,parmDict[item] ######### show dict just before refinement
125    ifPrint = True
126    if ifSeq:
127        ifPrint = False
128    Rvals = {}
129    chisq0 = None
130    while True:
131        G2mv.Map2Dict(parmDict,varyList)
132        begin = time.time()
133        values =  np.array(G2stMth.Dict2Values(parmDict, varyList))
134        if np.any(np.isnan(values)):
135            raise G2obj.G2Exception('ERROR - nan found in LS parameters - use Calculate/View LS parms to locate')
136        # test code to compute GOF and save for external repeat
137        #args = ([Histograms,Phases,restraintDict,rigidbodyDict],parmDict,varyList,calcControls,pawleyLookup,dlg)
138        #print '*** before fit chi**2',np.sum(G2stMth.errRefine(values,*args)**2)
139        #fl = open('beforeFit.cpickle','wb')
140        #cPickle.dump(values,fl,1)
141        #cPickle.dump(args[:-1],fl,1)
142        #fl.close()
143        Ftol = Controls['min dM/M']
144        Xtol = Controls['SVDtol']
145        Factor = Controls['shift factor']
146        if 'Jacobian' in Controls['deriv type']:
147            result = so.leastsq(G2stMth.errRefine,values,Dfun=G2stMth.dervRefine,full_output=True,
148                ftol=Ftol,col_deriv=True,factor=Factor,
149                args=([Histograms,Phases,restraintDict,rigidbodyDict],parmDict,varyList,calcControls,pawleyLookup,dlg))
150            ncyc = int(result[2]['nfev']/2)
151            if refPlotUpdate is not None: refPlotUpdate(Histograms)   # update plot after completion
152        elif 'analytic Hessian' in Controls['deriv type']:
153            Lamda = Controls.get('Marquardt',-3)
154            maxCyc = Controls['max cyc']
155            result = G2mth.HessianLSQ(G2stMth.errRefine,values,Hess=G2stMth.HessRefine,ftol=Ftol,xtol=Xtol,maxcyc=maxCyc,Print=ifPrint,lamda=Lamda,
156                args=([Histograms,Phases,restraintDict,rigidbodyDict],parmDict,varyList,calcControls,pawleyLookup,dlg),
157                refPlotUpdate=refPlotUpdate)
158            ncyc = result[2]['num cyc']+1
159            Rvals['lamMax'] = result[2]['lamMax']
160            if 'msg' in result[2]:
161                Rvals['msg'] = result[2]['msg']
162            Controls['Marquardt'] = -3  #reset to default
163            if 'chisq0' in result[2] and chisq0 is None:
164                chisq0 = result[2]['chisq0']
165        elif 'Hessian SVD' in Controls['deriv type']:
166            maxCyc = Controls['max cyc']
167            result = G2mth.HessianSVD(G2stMth.errRefine,values,Hess=G2stMth.HessRefine,ftol=Ftol,xtol=Xtol,maxcyc=maxCyc,Print=ifPrint,
168                args=([Histograms,Phases,restraintDict,rigidbodyDict],parmDict,varyList,calcControls,pawleyLookup,dlg),
169                refPlotUpdate=refPlotUpdate)
170            if result[1] is None:
171                IfOK = False
172                covMatrix = []
173                sig = len(varyList)*[None,]
174                break
175            ncyc = result[2]['num cyc']+1
176            if 'chisq0' in result[2] and chisq0 is None:
177                chisq0 = result[2]['chisq0']
178        else:           #'numeric'
179            result = so.leastsq(G2stMth.errRefine,values,full_output=True,ftol=Ftol,epsfcn=1.e-8,factor=Factor,
180                args=([Histograms,Phases,restraintDict,rigidbodyDict],parmDict,varyList,calcControls,pawleyLookup,dlg))
181            ncyc = 1
182            if len(varyList):
183                ncyc = int(result[2]['nfev']/len(varyList))
184            if refPlotUpdate is not None: refPlotUpdate(Histograms)   # update plot
185#        table = dict(zip(varyList,zip(values,result[0],(result[0]-values))))
186#        for item in table: print item,table[item]               #useful debug - are things shifting?
187        runtime = time.time()-begin
188        Rvals['SVD0'] = result[2].get('SVD0',0)
189        Rvals['converged'] = result[2].get('Converged')
190        Rvals['DelChi2'] = result[2].get('DelChi2',-1.)
191        Rvals['chisq'] = np.sum(result[2]['fvec']**2)
192        G2stMth.Values2Dict(parmDict, varyList, result[0])
193        G2mv.Dict2Map(parmDict,varyList)
194        Rvals['Nobs'] = Histograms['Nobs']
195        Rvals['Rwp'] = np.sqrt(Rvals['chisq']/Histograms['sumwYo'])*100.      #to %
196        Rvals['GOF'] = np.sqrt(Rvals['chisq']/(Histograms['Nobs']-len(varyList)))
197        printFile.write(' Number of function calls: %d No. of observations: %d No. of parameters: %d User rejected: %d Sp. gp. extinct: %d\n'%  \
198            (result[2]['nfev'],Histograms['Nobs'],len(varyList),Histograms['Nrej'],Histograms['Next']))
199        if ncyc:
200            printFile.write(' Refinement time = %8.3fs, %8.3fs/cycle, for %d cycles\n'%(runtime,runtime/ncyc,ncyc))
201        printFile.write(' wR = %7.2f%%, chi**2 = %12.6g, GOF = %6.2f\n'%(Rvals['Rwp'],Rvals['chisq'],Rvals['GOF']))
202        sig = len(varyList)*[None,]
203        if 'None' in str(type(result[1])) and ifSeq:    #this bails out of a sequential refinement on singular matrix
204            IfOK = False
205            covMatrix = []
206            G2fil.G2Print ('Warning: **** Refinement failed - singular matrix ****')
207            if 'Hessian' in Controls['deriv type']:
208                num = len(varyList)-1
209                for i,val in enumerate(np.flipud(result[2]['psing'])):
210                    if val:
211                        G2fil.G2Print('Bad parameter: '+varyList[num-i],mode='warn')
212            else:
213                Ipvt = result[2]['ipvt']
214                for i,ipvt in enumerate(Ipvt):
215                    if not np.sum(result[2]['fjac'],axis=1)[i]:
216                        G2fil.G2Print('Bad parameter: '+varyList[ipvt-1],mode='warn')
217            break
218        IfOK = True
219        try:
220            covMatrix = result[1]*Rvals['GOF']**2
221            sig = np.sqrt(np.diag(covMatrix))
222            Lastshft = result[2].get('Xvec',None)
223            if Lastshft is None:
224                Rvals['Max shft/sig'] = 0.0
225            else:
226                Rvals['Max shft/sig'] = np.max(Lastshft/sig)
227            if np.any(np.isnan(sig)) or not sig.shape:
228                G2fil.G2Print ('*** Least squares aborted - some invalid esds possible ***',mode='error')
229#            table = dict(zip(varyList,zip(values,result[0],(result[0]-values)/sig)))
230#            for item in table: print item,table[item]               #useful debug - are things shifting?
231            # report on refinement issues. Result in Rvals['msg']
232            ReportProblems(result,Rvals,varyList)
233            break                   #refinement succeeded - finish up!
234        except TypeError:
235            # if we get here, no result[1] (covar matrix) was returned or other calc error: refinement failed
236            IfOK = False
237            if not len(varyList):
238                covMatrix = []
239                break
240            if 'Hessian' in Controls['deriv type']:
241                SVD0 = result[2].get('SVD0')
242                if SVD0 == -1:
243                    G2fil.G2Print ('**** Refinement failed - singular matrix ****',mode='error')
244                elif SVD0 == -2:
245                    G2fil.G2Print ('**** Refinement failed - other problem ****',mode='error')
246                elif SVD0 > 0:
247                    G2fil.G2Print ('**** Refinement failed with {} SVD singularities ****'.format(SVD0),mode='error')
248                else:
249                    G2fil.G2Print ('**** Refinement failed ****',mode='error')
250                if result[1] is None:
251                    IfOK = False
252                    covMatrix = []
253                    sig = len(varyList)*[None,]
254                # report on highly correlated variables
255                ReportProblems(result,Rvals,varyList)
256                # process singular variables
257                psing = result[2].get('psing',[])
258                if dlg: break # refining interactively
259                # non-interactive refinement
260                for val in sorted(psing,reverse=True):
261                    G2fil.G2Print ('Removing parameter: '+varyList[val])
262                    del(varyList[val])
263                if not psing: break    # removed variable(s), try again
264            else:
265                G2fil.G2Print ('**** Refinement failed - singular matrix ****',mode='error')
266                Ipvt = result[2]['ipvt']
267                for i,ipvt in enumerate(Ipvt):
268                    if not np.sum(result[2]['fjac'],axis=1)[i]:
269                        G2fil.G2Print ('Removing parameter: '+varyList[ipvt-1])
270                        del(varyList[ipvt-1])
271                        break
272    if IfOK:
273        G2stMth.GetFobsSq(Histograms,Phases,parmDict,calcControls)
274    if chisq0 is not None:
275        Rvals['GOF0'] = np.sqrt(chisq0/(Histograms['Nobs']-len(varyList)))
276    return IfOK,Rvals,result,covMatrix,sig
277
278def Refine(GPXfile,dlg=None,makeBack=True,refPlotUpdate=None):
279    '''Global refinement -- refines to minimize against all histograms.
280    This can be called in one of three ways, from :meth:`GSASIIdataGUI.GSASII.OnRefine` in an
281    interactive refinement, where dlg will be a wx.ProgressDialog, or non-interactively from
282    :meth:`GSASIIscriptable.G2Project.refine` or from :func:`main`, where dlg will be None.
283    '''
284    import GSASIImpsubs as G2mp
285    G2mp.InitMP()
286    import pytexture as ptx
287    ptx.pyqlmninit()            #initialize fortran arrays for spherical harmonics
288
289    printFile = open(ospath.splitext(GPXfile)[0]+'.lst','w')
290    G2stIO.ShowBanner(printFile)
291    varyList = []
292    parmDict = {}
293    G2mv.InitVars()
294    Controls = G2stIO.GetControls(GPXfile)
295    G2stIO.ShowControls(Controls,printFile)
296    calcControls = {}
297    calcControls.update(Controls)
298    constrDict,fixedList = G2stIO.GetConstraints(GPXfile)
299    restraintDict = G2stIO.GetRestraints(GPXfile)
300    Histograms,Phases = G2stIO.GetUsedHistogramsAndPhases(GPXfile)
301    if not Phases:
302        G2fil.G2Print (' *** ERROR - you have no phases to refine! ***')
303        G2fil.G2Print (' *** Refine aborted ***')
304        return False,{'msg':'No phases'}
305    if not Histograms:
306        G2fil.G2Print (' *** ERROR - you have no data to refine with! ***')
307        G2fil.G2Print (' *** Refine aborted ***')
308        return False,{'msg':'No data'}
309    rigidbodyDict = G2stIO.GetRigidBodies(GPXfile)
310    rbIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})
311    rbVary,rbDict = G2stIO.GetRigidBodyModels(rigidbodyDict,pFile=printFile)
312    (Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables,MFtables,
313         maxSSwave) = G2stIO.GetPhaseData(Phases,restraintDict,rbIds,pFile=printFile)
314    calcControls['atomIndx'] = atomIndx
315    calcControls['Natoms'] = Natoms
316    calcControls['FFtables'] = FFtables
317    calcControls['BLtables'] = BLtables
318    calcControls['MFtables'] = MFtables
319    calcControls['maxSSwave'] = maxSSwave
320    hapVary,hapDict,controlDict = G2stIO.GetHistogramPhaseData(Phases,Histograms,pFile=printFile)
321    TwConstr,TwFixed = G2stIO.makeTwinFrConstr(Phases,Histograms,hapVary)
322    constrDict += TwConstr
323    fixedList += TwFixed
324    calcControls.update(controlDict)
325    histVary,histDict,controlDict = G2stIO.GetHistogramData(Histograms,pFile=printFile)
326    calcControls.update(controlDict)
327    varyList = rbVary+phaseVary+hapVary+histVary
328    parmDict.update(rbDict)
329    parmDict.update(phaseDict)
330    parmDict.update(hapDict)
331    parmDict.update(histDict)
332    G2stIO.GetFprime(calcControls,Histograms)
333    # do constraint processing
334    varyListStart = tuple(varyList) # save the original varyList before dependent vars are removed
335    msg = G2mv.EvaluateMultipliers(constrDict,parmDict)
336    if msg:
337        return False,{'msg':'Unable to interpret multiplier(s): '+msg}
338    try:
339        G2mv.GenerateConstraints(varyList,constrDict,fixedList,parmDict)
340        #print(G2mv.VarRemapShow(varyList))
341        #print('DependentVars',G2mv.GetDependentVars())
342        #print('IndependentVars',G2mv.GetIndependentVars())
343    except G2mv.ConstraintException:
344        G2fil.G2Print (' *** ERROR - your constraints are internally inconsistent ***')
345        #errmsg, warnmsg = G2mv.CheckConstraints(varyList,constrDict,fixedList)
346        #print 'Errors',errmsg
347        #if warnmsg: print 'Warnings',warnmsg
348        return False,{'msg':' Constraint error'}
349#    print G2mv.VarRemapShow(varyList)
350
351    # remove frozen vars from refinement
352    if 'parmFrozen' not in Controls:
353        Controls['parmFrozen'] = {}
354    if 'FrozenList' not in Controls['parmFrozen']: 
355        Controls['parmFrozen']['FrozenList'] = []
356    parmFrozenList = Controls['parmFrozen']['FrozenList']
357    frozenList = [i for i in varyList if i in parmFrozenList]
358    if len(frozenList) != 0: 
359        varyList = [i for i in varyList if i not in parmFrozenList]
360        G2fil.G2Print(
361            'Frozen refined variables (due to exceeding limits)\n\t:{}'
362            .format(frozenList))
363       
364    ifSeq = False
365    printFile.write('\n Refinement results:\n')
366    printFile.write(135*'-'+'\n')
367    Rvals = {}
368    try:
369        covData = {}
370        IfOK,Rvals,result,covMatrix,sig = RefineCore(Controls,Histograms,Phases,restraintDict,
371            rigidbodyDict,parmDict,varyList,calcControls,pawleyLookup,ifSeq,printFile,dlg,
372            refPlotUpdate=refPlotUpdate)
373        if IfOK:
374            sigDict = dict(zip(varyList,sig))
375            newCellDict = G2stMth.GetNewCellParms(parmDict,varyList)
376            newAtomDict = G2stMth.ApplyXYZshifts(parmDict,varyList)
377            covData = {'variables':result[0],'varyList':varyList,'sig':sig,'Rvals':Rvals,
378                       'varyListStart':varyListStart,
379                       'covMatrix':covMatrix,'title':GPXfile,'newAtomDict':newAtomDict,
380                       'newCellDict':newCellDict,'freshCOV':True}
381            # add the uncertainties into the esd dictionary (sigDict)
382            sigDict.update(G2mv.ComputeDepESD(covMatrix,varyList,parmDict))
383            # check for variables outside their allowed range, reset and freeze them
384            frozen = dropOOBvars(varyList,parmDict,sigDict,Controls,parmFrozenList)
385            G2mv.PrintIndependentVars(parmDict,varyList,sigDict,pFile=printFile)
386            G2stMth.ApplyRBModels(parmDict,Phases,rigidbodyDict,True)
387            G2stIO.SetRigidBodyModels(parmDict,sigDict,rigidbodyDict,printFile)
388            G2stIO.SetPhaseData(parmDict,sigDict,Phases,rbIds,covData,restraintDict,printFile)
389            G2stIO.PrintISOmodes(printFile,Phases,parmDict,sigDict)
390            G2stIO.SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms,calcControls,pFile=printFile)
391            G2stIO.SetHistogramData(parmDict,sigDict,Histograms,calcControls,pFile=printFile)
392            if len(frozen):
393                if 'msg' in Rvals:
394                    Rvals['msg'] += '\n'
395                else:
396                    Rvals['msg'] = ''
397                msg = ('Warning: {} variable(s) refined outside limits and were frozen ({} total frozen)'
398                    .format(len(frozen),len(parmFrozenList))
399                    )
400                G2fil.G2Print(msg)
401                Rvals['msg'] += msg
402            elif len(parmFrozenList):
403                if 'msg' in Rvals:
404                    Rvals['msg'] += '\n'
405                else:
406                    Rvals['msg'] = ''
407                msg = ('Note: a total of {} variable(s) are frozen due to refining outside limits'
408                    .format(len(parmFrozenList))
409                    )
410                G2fil.G2Print('Note: ',msg)
411                Rvals['msg'] += msg
412            G2stIO.SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases,rigidbodyDict,covData,parmFrozenList,makeBack)
413            printFile.close()
414            G2fil.G2Print (' Refinement results are in file: '+ospath.splitext(GPXfile)[0]+'.lst')
415            G2fil.G2Print (' ***** Refinement successful *****')
416        else:
417            G2fil.G2Print ('****ERROR - Refinement failed',mode='error')
418            if 'msg' in Rvals:
419                G2fil.G2Print ('Note refinement problem:',mode='warn')
420                G2fil.G2Print (Rvals['msg'],mode='warn')
421            raise G2obj.G2Exception('**** ERROR: Refinement failed ****')
422    except G2obj.G2RefineCancel as Msg:
423        printFile.close()
424        G2fil.G2Print (' ***** Refinement stopped *****')
425        if not hasattr(Msg,'msg'): Msg.msg = str(Msg)
426        if 'msg' in Rvals:
427            Rvals['msg'] += '\n'
428            Rvals['msg'] += Msg.msg
429            if not dlg:
430                G2fil.G2Print ('Note refinement problem:',mode='warn')
431                G2fil.G2Print (Rvals['msg'],mode='warn')
432        else:
433            Rvals['msg'] = Msg.msg
434        return False,Rvals
435#    except G2obj.G2Exception as Msg:  # cell metric error, others?
436    except Exception as Msg:  # cell metric error, others?
437        if GSASIIpath.GetConfigValue('debug'):
438            import traceback
439            print(traceback.format_exc())       
440        if not hasattr(Msg,'msg'): Msg.msg = str(Msg)
441        printFile.close()
442        G2fil.G2Print (' ***** Refinement error *****')
443        if 'msg' in Rvals:
444            Rvals['msg'] += '\n\n'
445            Rvals['msg'] += Msg.msg
446            if not dlg:
447                G2fil.G2Print ('Note refinement problem:',mode='warn')
448                G2fil.G2Print (Rvals['msg'],mode='warn')
449        else:
450            Rvals['msg'] = Msg.msg
451        return False,Rvals
452
453#for testing purposes, create a file for testderiv
454    if GSASIIpath.GetConfigValue('debug'):   # and IfOK:
455#needs: values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup
456        fl = open(ospath.splitext(GPXfile)[0]+'.testDeriv','wb')
457        cPickle.dump(result[0],fl,1)
458        cPickle.dump([Histograms,Phases,restraintDict,rigidbodyDict],fl,1)
459        cPickle.dump([constrDict,fixedList,G2mv.GetDependentVars()],fl,1)
460        cPickle.dump(parmDict,fl,1)
461        cPickle.dump(varyList,fl,1)
462        cPickle.dump(calcControls,fl,1)
463        cPickle.dump(pawleyLookup,fl,1)
464        fl.close()
465    if dlg:
466        return True,Rvals
467    elif 'msg' in Rvals:
468        G2fil.G2Print ('Reported from refinement:',mode='warn')
469        G2fil.G2Print (Rvals['msg'],mode='warn')
470
471def phaseCheck(phaseVary,Phases,histogram):
472    '''
473    Removes unused parameters from phase varylist if phase not in histogram
474    #TODO - implement "Fix FXU" for seq refinement here - done?
475    '''
476    NewVary = []
477    for phase in Phases:
478        if histogram not in Phases[phase]['Histograms']: continue
479        if Phases[phase]['Histograms'][histogram]['Use']:
480            pId = Phases[phase]['pId']
481            newVary = [item for item in phaseVary if item.split(':')[0] == str(pId)]
482            FixVals = Phases[phase]['Histograms'][histogram].get('Fix FXU',' ')
483            if 'F' in FixVals:
484                newVary = [item for item in newVary if not 'Afrac' in item]
485            if 'X' in FixVals:
486                newVary = [item for item in newVary if not 'dA' in item]
487            if 'U' in FixVals:
488                newVary = [item for item in newVary if not 'AU' in item]
489            if 'M' in FixVals:
490                newVary = [item for item in newVary if not 'AM' in item]
491            NewVary += newVary
492    return NewVary
493
494def SeqRefine(GPXfile,dlg,refPlotUpdate=None):
495    '''Perform a sequential refinement -- cycles through all selected histgrams,
496    one at a time
497    '''
498    import GSASIImpsubs as G2mp
499    G2mp.InitMP()
500    import pytexture as ptx
501    ptx.pyqlmninit()            #initialize fortran arrays for spherical harmonics
502    msgs = {}
503    printFile = open(ospath.splitext(GPXfile)[0]+'.lst','w')
504    G2fil.G2Print ('Starting Sequential Refinement')
505    G2stIO.ShowBanner(printFile)
506    Controls = G2stIO.GetControls(GPXfile)
507    preFrozenCount = 0
508    for h in Controls['parmFrozen']:
509        if h == 'FrozenList': continue
510        preFrozenCount += len(Controls['parmFrozen'][h])   
511    G2stIO.ShowControls(Controls,printFile,SeqRef=True,
512                            preFrozenCount=preFrozenCount)
513    restraintDict = G2stIO.GetRestraints(GPXfile)
514    Histograms,Phases = G2stIO.GetUsedHistogramsAndPhases(GPXfile)
515    if not Phases:
516        G2fil.G2Print (' *** ERROR - you have no phases to refine! ***')
517        G2fil.G2Print (' *** Refine aborted ***')
518        return False,'No phases'
519    if not Histograms:
520        G2fil.G2Print (' *** ERROR - you have no data to refine with! ***')
521        G2fil.G2Print (' *** Refine aborted ***')
522        return False,'No data'
523    rigidbodyDict = G2stIO.GetRigidBodies(GPXfile)
524    rbIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})
525    rbVary,rbDict = G2stIO.GetRigidBodyModels(rigidbodyDict,pFile=printFile)
526    G2mv.InitVars()
527    (Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables,MFtables,
528         maxSSwave) = G2stIO.GetPhaseData(Phases,restraintDict,rbIds,
529                                    Print=False,pFile=printFile,seqRef=True)
530    for item in phaseVary:
531        if '::A0' in item:
532            G2fil.G2Print ('**** WARNING - lattice parameters should not be refined in a sequential refinement ****')
533            G2fil.G2Print ('****           instead use the Dij parameters for each powder histogram            ****')
534            return False,'Lattice parameter refinement error - see console message'
535        if '::C(' in item:
536            G2fil.G2Print ('**** WARNING - phase texture parameters should not be refined in a sequential refinement ****')
537            G2fil.G2Print ('****           instead use the C(L,N) parameters for each powder histogram               ****')
538            return False,'Phase texture refinement error - see console message'
539    if 'Seq Data' in Controls:
540        histNames = Controls['Seq Data']
541    else: # patch from before Controls['Seq Data'] was implemented?
542        histNames = G2stIO.GetHistogramNames(GPXfile,['PWDR',])
543    if Controls.get('Reverse Seq'):
544        histNames.reverse()
545    SeqResult = G2stIO.GetSeqResult(GPXfile)
546#    SeqResult = {'SeqPseudoVars':{},'SeqParFitEqList':[]}
547    Histo = {}
548    NewparmDict = {}
549    G2stIO.SetupSeqSavePhases(GPXfile)
550    msgs['steepestNum'] = 0
551    msgs['maxshift/sigma'] = []
552    for ihst,histogram in enumerate(histNames):
553        if GSASIIpath.GetConfigValue('Show_timing'): t1 = time.time()
554        G2fil.G2Print('\nRefining with '+str(histogram))
555        G2mv.InitVars()
556        (Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,
557             FFtables,BLtables,MFtables,maxSSwave) = G2stIO.GetPhaseData(
558                 Phases,restraintDict,rbIds,
559                 Print=False,pFile=printFile,seqRef=True)
560        ifPrint = False
561        if dlg:
562            dlg.SetTitle('Residual for histogram '+str(ihst))
563        calcControls = {}
564        calcControls['atomIndx'] = atomIndx
565        calcControls['Natoms'] = Natoms
566        calcControls['FFtables'] = FFtables
567        calcControls['BLtables'] = BLtables
568        calcControls['MFtables'] = MFtables
569        calcControls['maxSSwave'] = maxSSwave
570        if histogram not in Histograms:
571            G2fil.G2Print("Error: not found!")
572            continue
573    #TODO - implement "Fix FXU" for seq refinement here - done?
574        hId = Histograms[histogram]['hId']
575        redphaseVary = phaseCheck(phaseVary,Phases,histogram)
576        Histo = {histogram:Histograms[histogram],}
577        hapVary,hapDict,controlDict = G2stIO.GetHistogramPhaseData(Phases,Histo,Print=False)
578        calcControls.update(controlDict)
579        histVary,histDict,controlDict = G2stIO.GetHistogramData(Histo,False)
580        calcControls.update(controlDict)
581        varyList = rbVary+redphaseVary+hapVary+histVary
582#        if not ihst:
583            # save the initial vary list, but without histogram numbers on parameters
584        saveVaryList = varyList[:]
585        for i,item in enumerate(saveVaryList):
586            items = item.split(':')
587            if items[1]:
588                items[1] = ''
589            item = ':'.join(items)
590            saveVaryList[i] = item
591        if not ihst:
592            SeqResult['varyList'] = saveVaryList
593        else:
594            SeqResult['varyList'] = list(set(SeqResult['varyList']+saveVaryList))
595        parmDict = {}
596        parmDict.update(rbDict)
597        parmDict.update(phaseDict)
598        parmDict.update(hapDict)
599        parmDict.update(histDict)
600        if Controls['Copy2Next']:   # update with parms from last histogram
601            #parmDict.update(NewparmDict) # don't use in case extra entries would cause a problem
602            for parm in NewparmDict:
603                if parm in parmDict:
604                    parmDict[parm] = NewparmDict[parm]
605        elif histogram in SeqResult:  # update phase from last seq ref
606            NewparmDict = SeqResult[histogram].get('parmDict',{})
607            for parm in NewparmDict:
608                if '::' in parm and parm in parmDict:
609                    parmDict[parm] = NewparmDict[parm]
610           
611        G2stIO.GetFprime(calcControls,Histo)
612        # do constraint processing
613        #reload(G2mv) # debug
614        constrDict,fixedList = G2stIO.GetConstraints(GPXfile)
615        varyListStart = tuple(varyList) # save the original varyList before dependent vars are removed
616        msg = G2mv.EvaluateMultipliers(constrDict,parmDict)
617        if msg:
618            return False,'Unable to interpret multiplier(s): '+msg
619        try:
620            groups,parmlist = G2mv.GenerateConstraints(varyList,constrDict,fixedList,parmDict,SeqHist=hId)
621#            if GSASIIpath.GetConfigValue('debug'): print("DBG_"+
622#                G2mv.VarRemapShow(varyList,True))
623#            print('DependentVars',G2mv.GetDependentVars())
624#            print('IndependentVars',G2mv.GetIndependentVars())
625            constraintInfo = (groups,parmlist,constrDict,fixedList,ihst)
626        except G2mv.ConstraintException:
627            G2fil.G2Print (' *** ERROR - your constraints are internally inconsistent ***')
628            #errmsg, warnmsg = G2mv.CheckConstraints(varyList,constrDict,fixedList)
629            #print 'Errors',errmsg
630            #if warnmsg: print 'Warnings',warnmsg
631            return False,' Constraint error'
632        #print G2mv.VarRemapShow(varyList)
633        if not ihst:
634            # first histogram to refine against
635            firstVaryList = []
636            for item in varyList:
637                items = item.split(':')
638                if items[1]:
639                    items[1] = ''
640                item = ':'.join(items)
641                firstVaryList.append(item)
642            newVaryList = firstVaryList
643        else:
644            newVaryList = []
645            for item in varyList:
646                items = item.split(':')
647                if items[1]:
648                    items[1] = ''
649                item = ':'.join(items)
650                newVaryList.append(item)
651        if newVaryList != firstVaryList and Controls['Copy2Next']:
652            # variable lists are expected to match between sequential refinements when Copy2Next is on
653            #print '**** ERROR - variable list for this histogram does not match previous'
654            #print '     Copy of variables is not possible'
655            #print '\ncurrent histogram',histogram,'has',len(newVaryList),'variables'
656            combined = list(set(firstVaryList+newVaryList))
657            c = [var for var in combined if var not in newVaryList]
658            p = [var for var in combined if var not in firstVaryList]
659            G2fil.G2Print('*** Variables change ***')
660            for typ,vars in [('Removed',c),('Added',p)]:
661                line = '  '+typ+': '
662                if vars:
663                    for var in vars:
664                        if len(line) > 70:
665                            G2fil.G2Print(line)
666                            line = '    '
667                        line += var + ', '
668                else:
669                        line += 'none, '
670                G2fil.G2Print(line[:-2])
671            firstVaryList = newVaryList
672
673        ifSeq = True
674        printFile.write('\n Refinement results for histogram id {}: {}\n'
675                            .format(hId,histogram))
676        printFile.write(135*'-'+'\n')
677        # remove frozen vars
678        if 'parmFrozen' not in Controls:
679            Controls['parmFrozen'] = {}
680        if histogram not in Controls['parmFrozen']: 
681            Controls['parmFrozen'][histogram] = []
682        parmFrozenList = Controls['parmFrozen'][histogram]
683        frozenList = [i for i in varyList if i in parmFrozenList]
684        if len(frozenList) != 0: 
685           varyList = [i for i in varyList if i not in parmFrozenList]
686           s = ''
687           for a in frozenList:
688               if s:
689                   s+= ', '
690               s += a
691           printFile.write(
692               ' The following refined variables have previously been frozen due to exceeding limits:\n\t{}\n'
693               .format(s))
694        try:
695            IfOK,Rvals,result,covMatrix,sig = RefineCore(Controls,Histo,Phases,restraintDict,
696                rigidbodyDict,parmDict,varyList,calcControls,pawleyLookup,ifSeq,printFile,dlg,
697                refPlotUpdate=refPlotUpdate)
698            G2fil.G2Print ('  wR = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f, last delta chi = %.4f, last shft/sig = %.4f'%(
699                Rvals['Rwp'],Rvals['chisq'],Rvals['GOF']**2,Rvals['DelChi2'],Rvals.get('Max shft/sig',np.nan)))
700            if Rvals.get('lamMax',0) >= 10.:
701                msgs['steepestNum'] += 1
702            if Rvals.get('Max shft/sig'):
703                msgs['maxshift/sigma'].append(Rvals['Max shft/sig'])
704            # add the uncertainties into the esd dictionary (sigDict)
705            if not IfOK:
706                G2fil.G2Print('***** Sequential refinement failed at histogram '+histogram,mode='warn')
707                break
708            sigDict = dict(zip(varyList,sig))
709            # the uncertainties for dependent constrained parms into the esd dict
710            sigDict.update(G2mv.ComputeDepESD(covMatrix,varyList,parmDict))
711            # check for variables outside their allowed range, reset and freeze them
712            frozen = dropOOBvars(varyList,parmDict,sigDict,Controls,parmFrozenList)
713            msg = None
714            if len(frozen) > 0:
715               msg = ('Hist {}: {} variables were outside limits and were frozen (now {} frozen total)'
716                   .format(ihst,len(frozen),len(parmFrozenList)))
717               G2fil.G2Print(msg)
718               msg = (' {} variables were outside limits and were frozen (now {} frozen total)'
719                   .format(len(frozen),len(parmFrozenList)))
720               for p in frozen:
721                   if p not in varyList:
722                       print('Frozen Warning: {} not in varyList. This should not happen!'.format(p))
723                       continue
724                   i = varyList.index(p)
725                   result[0][i] = parmDict[p]
726                   sig[i] = -0.1
727            # a dict with values & esds for dependent (constrained) parameters - avoid extraneous holds
728            depParmDict = {i:(parmDict[i],sigDict[i]) for i in varyListStart if i in sigDict and i not in varyList}
729            newCellDict = copy.deepcopy(G2stMth.GetNewCellParms(parmDict,varyList))
730            newAtomDict = copy.deepcopy(G2stMth.ApplyXYZshifts(parmDict,varyList))
731            histRefData = {
732                'variables':result[0],'varyList':varyList,'sig':sig,'Rvals':Rvals,
733                'varyListStart':varyListStart,
734                'covMatrix':covMatrix,'title':histogram,'newAtomDict':newAtomDict,
735                'newCellDict':newCellDict,'depParmDict':depParmDict,
736                'constraintInfo':constraintInfo,
737                'parmDict':parmDict,
738                }
739            SeqResult[histogram] = histRefData
740            G2stMth.ApplyRBModels(parmDict,Phases,rigidbodyDict,True)
741#            G2stIO.SetRigidBodyModels(parmDict,sigDict,rigidbodyDict,printFile)
742            G2stIO.SetHistogramPhaseData(parmDict,sigDict,Phases,Histo,None,ifPrint,printFile)
743            G2stIO.SetHistogramData(parmDict,sigDict,Histo,None,ifPrint,printFile,seq=True)
744            G2stIO.SaveUpdatedHistogramsAndPhases(GPXfile,Histo,Phases,rigidbodyDict,histRefData,Controls['parmFrozen'])
745            if msg: 
746                printFile.write(msg+'\n')
747            NewparmDict = {}
748            # make dict of varied parameters in current histogram, renamed to
749            # next histogram, for use in next refinement.
750            if Controls['Copy2Next'] and ihst < len(histNames)-1:
751                hId = Histo[histogram]['hId'] # current histogram
752                nexthId = Histograms[histNames[ihst+1]]['hId']
753                for parm in set(list(varyList)+list(varyListStart)):
754                    items = parm.split(':')
755                    if len(items) < 3: 
756                        continue
757                    if str(hId) in items[1]:
758                        items[1] = str(nexthId)
759                        newparm = ':'.join(items)
760                        NewparmDict[newparm] = parmDict[parm]
761                    else:
762                        if items[2].startswith('dA'): parm = parm.replace(':dA',':A') 
763                        NewparmDict[parm] = parmDict[parm]
764                   
765        except G2obj.G2RefineCancel as Msg:
766            if not hasattr(Msg,'msg'): Msg.msg = str(Msg)
767            printFile.close()
768            G2fil.G2Print (' ***** Refinement stopped *****')
769            return False,Msg.msg
770        except G2obj.G2Exception as Msg:  # cell metric error, others?
771            if not hasattr(Msg,'msg'): Msg.msg = str(Msg)
772            printFile.close()
773            G2fil.G2Print (' ***** Refinement error *****')
774            return False,Msg.msg
775        if GSASIIpath.GetConfigValue('Show_timing'):
776            t2 = time.time()
777            G2fil.G2Print("Fit step time {:.2f} sec.".format(t2-t1))
778            t1 = t2
779    SeqResult['histNames'] = [itm for itm in G2stIO.GetHistogramNames(GPXfile,['PWDR',]) if itm in SeqResult.keys()]
780    try:
781        G2stIO.SetSeqResult(GPXfile,Histograms,SeqResult)
782    except Exception as msg:
783        print('Error reading Sequential results\n',str(msg))
784        if GSASIIpath.GetConfigValue('debug'):
785            import traceback
786            print(traceback.format_exc())       
787    postFrozenCount = 0
788    for h in Controls['parmFrozen']:
789        if h == 'FrozenList': continue
790        postFrozenCount += len(Controls['parmFrozen'][h])
791    if postFrozenCount:
792        msgs['Frozen'] = 'Ending refinement with {} Frozen variables ({} added now)\n'.format(postFrozenCount,postFrozenCount-preFrozenCount)
793        printFile.write('\n'+msgs['Frozen'])
794    printFile.close()
795    G2fil.G2Print (' Sequential refinement results are in file: '+ospath.splitext(GPXfile)[0]+'.lst')
796    G2fil.G2Print (' ***** Sequential refinement successful *****')
797    return True,msgs
798
799def dropOOBvars(varyList,parmDict,sigDict,Controls,parmFrozenList):
800    '''Find variables in the parameters dict that are outside the ranges
801    (in parmMinDict and parmMaxDict) and set them to the limits values.
802    Add any such variables into the list of frozen variable
803    (parmFrozenList). Returns a list of newly frozen variables, if any.
804    '''
805    parmMinDict = Controls.get('parmMinDict',{})
806    parmMaxDict = Controls.get('parmMaxDict',{})
807    freeze = []
808    if parmMinDict or parmMaxDict:
809        for name in varyList:
810            if name not in parmDict: continue
811            n,val = G2obj.prmLookup(name,parmMinDict)
812            if n is not None:
813                if parmDict[name] < parmMinDict[n]:
814                    parmDict[name] = parmMinDict[n]
815                    sigDict[name] = 0.0
816                    freeze.append(name)
817                    continue
818            n,val = G2obj.prmLookup(name,parmMaxDict)
819            if n is not None:
820                if parmDict[name] > parmMaxDict[n]:
821                    parmDict[name] = parmMaxDict[n]
822                    sigDict[name] = 0.0
823                    freeze.append(name)
824                    continue
825        for v in freeze:
826            if v not in parmFrozenList:
827                parmFrozenList.append(v)
828    return freeze
829
830def RetDistAngle(DisAglCtls,DisAglData,dlg=None):
831    '''Compute and return distances and angles
832
833    :param dict DisAglCtls: contains distance/angle radii usually defined using
834       :func:`GSASIIctrlGUI.DisAglDialog`
835    :param dict DisAglData: contains phase data:
836       Items 'OrigAtoms' and 'TargAtoms' contain the atoms to be used
837       for distance/angle origins and atoms to be used as targets.
838       Item 'SGData' has the space group information (see :ref:`Space Group object<SGData_table>`)
839
840    :returns: AtomLabels,DistArray,AngArray where:
841
842      **AtomLabels** is a dict of atom labels, keys are the atom number
843
844      **DistArray** is a dict keyed by the origin atom number where the value is a list
845      of distance entries. The value for each distance is a list containing:
846
847        0) the target atom number (int);
848        1) the unit cell offsets added to x,y & z (tuple of int values)
849        2) the symmetry operator number (which may be modified to indicate centering and center of symmetry)
850        3) an interatomic distance in A (float)
851        4) an uncertainty on the distance in A or 0.0 (float)
852
853      **AngArray** is a dict keyed by the origin (central) atom number where
854      the value is a list of
855      angle entries. The value for each angle entry consists of three values:
856
857        0) a distance item reference for one neighbor (int)
858        1) a distance item reference for a second neighbor (int)
859        2) a angle, uncertainty pair; the s.u. may be zero (tuple of two floats)
860
861      The AngArray distance reference items refer directly to the index of the items in the
862      DistArray item for the list of distances for the central atom.
863    '''
864    import numpy.ma as ma
865
866    SGData = DisAglData['SGData']
867    Cell = DisAglData['Cell']
868
869    Amat,Bmat = G2lat.cell2AB(Cell[:6])
870    covData = {}
871    if len(DisAglData.get('covData',{})):
872        covData = DisAglData['covData']
873        covMatrix = covData['covMatrix']
874        varyList = covData['varyList']
875        pfx = str(DisAglData['pId'])+'::'
876
877    Factor = DisAglCtls['Factors']
878    Radii = dict(zip(DisAglCtls['AtomTypes'],zip(DisAglCtls['BondRadii'],DisAglCtls['AngleRadii'])))
879    indices = (-2,-1,0,1,2)
880    Units = np.array([[h,k,l] for h in indices for k in indices for l in indices])
881    origAtoms = DisAglData['OrigAtoms']
882    targAtoms = DisAglData['TargAtoms']
883    AtomLabels = {}
884    for Oatom in origAtoms:
885        AtomLabels[Oatom[0]] = Oatom[1]
886    for Oatom in targAtoms:
887        AtomLabels[Oatom[0]] = Oatom[1]
888    DistArray = {}
889    AngArray = {}
890    for iO,Oatom in enumerate(origAtoms):
891        DistArray[Oatom[0]] = []
892        AngArray[Oatom[0]] = []
893        OxyzNames = ''
894        IndBlist = []
895        Dist = []
896        Vect = []
897        VectA = []
898        angles = []
899        for Tatom in targAtoms:
900            Xvcov = []
901            TxyzNames = ''
902            if len(DisAglData.get('covData',{})):
903                OxyzNames = [pfx+'dAx:%d'%(Oatom[0]),pfx+'dAy:%d'%(Oatom[0]),pfx+'dAz:%d'%(Oatom[0])]
904                TxyzNames = [pfx+'dAx:%d'%(Tatom[0]),pfx+'dAy:%d'%(Tatom[0]),pfx+'dAz:%d'%(Tatom[0])]
905                Xvcov = G2mth.getVCov(OxyzNames+TxyzNames,varyList,covMatrix)
906            result = G2spc.GenAtom(Tatom[3:6],SGData,False,Move=False)
907            BsumR = (Radii[Oatom[2]][0]+Radii[Tatom[2]][0])*Factor[0]
908            AsumR = (Radii[Oatom[2]][1]+Radii[Tatom[2]][1])*Factor[1]
909            for [Txyz,Top,Tunit,Spn] in result:
910                Dx = (Txyz-np.array(Oatom[3:6]))+Units
911                dx = np.inner(Amat,Dx)
912                dist = ma.masked_less(np.sqrt(np.sum(dx**2,axis=0)),0.5)
913                IndB = ma.nonzero(ma.masked_greater(dist-BsumR,0.))
914                if np.any(IndB):
915                    for indb in IndB:
916                        for i in range(len(indb)):
917                            if str(dx.T[indb][i]) not in IndBlist:
918                                IndBlist.append(str(dx.T[indb][i]))
919                                unit = Units[indb][i]
920                                tunit = (unit[0]+Tunit[0],unit[1]+Tunit[1],unit[2]+Tunit[2])
921                                sig = 0.0
922                                if len(Xvcov):
923                                    pdpx = G2mth.getDistDerv(Oatom[3:6],Tatom[3:6],Amat,unit,Top,SGData)
924                                    sig = np.sqrt(np.inner(pdpx,np.inner(pdpx,Xvcov)))
925                                Dist.append([Oatom[0],Tatom[0],tunit,Top,ma.getdata(dist[indb])[i],sig])
926                                if (Dist[-1][-2]-AsumR) <= 0.:
927                                    Vect.append(dx.T[indb][i]/Dist[-1][-2])
928                                    VectA.append([OxyzNames,np.array(Oatom[3:6]),TxyzNames,np.array(Tatom[3:6]),unit,Top])
929                                else:
930                                    Vect.append([0.,0.,0.])
931                                    VectA.append([])
932        if dlg is not None:
933            dlg.Update(iO,newmsg='Atoms done=%d'%(iO))
934        for D in Dist:
935            DistArray[Oatom[0]].append(D[1:])
936        Vect = np.array(Vect)
937        angles = np.zeros((len(Vect),len(Vect)))
938        angsig = np.zeros((len(Vect),len(Vect)))
939        for i,veca in enumerate(Vect):
940            if np.any(veca):
941                for j,vecb in enumerate(Vect):
942                    if np.any(vecb):
943                        angles[i][j],angsig[i][j] = G2mth.getAngSig(VectA[i],VectA[j],Amat,SGData,covData)
944                        if i <= j: continue
945                        AngArray[Oatom[0]].append((i,j,
946                            G2mth.getAngSig(VectA[i],VectA[j],Amat,SGData,covData)))
947    return AtomLabels,DistArray,AngArray
948
949def PrintDistAngle(DisAglCtls,DisAglData,out=sys.stdout):
950    '''Print distances and angles
951
952    :param dict DisAglCtls: contains distance/angle radii usually defined using
953       :func:`GSASIIctrlGUI.DisAglDialog`
954    :param dict DisAglData: contains phase data:
955       Items 'OrigAtoms' and 'TargAtoms' contain the atoms to be used
956       for distance/angle origins and atoms to be used as targets.
957       Item 'SGData' has the space group information (see :ref:`Space Group object<SGData_table>`)
958    :param file out: file object for output. Defaults to sys.stdout.
959    '''
960    def MyPrint(s):
961        out.write(s+'\n')
962        # print(s,file=out) # use in Python 3
963
964    def ShowBanner(name):
965        MyPrint(80*'*')
966        MyPrint('   Interatomic Distances and Angles for phase '+name)
967        MyPrint((80*'*')+'\n')
968
969    ShowBanner(DisAglCtls['Name'])
970    SGData = DisAglData['SGData']
971    SGtext,SGtable = G2spc.SGPrint(SGData)
972    for line in SGtext: MyPrint(line)
973    if len(SGtable) > 1:
974        for i,item in enumerate(SGtable[::2]):
975            if 2*i+1 == len(SGtable):
976                line = ' %s'%(item.ljust(30))
977            else:
978                line = ' %s %s'%(item.ljust(30),SGtable[2*i+1].ljust(30))
979            MyPrint(line)
980    else:
981        MyPrint(' ( 1)    %s'%(SGtable[0])) #triclinic case
982    Cell = DisAglData['Cell']
983
984    Amat,Bmat = G2lat.cell2AB(Cell[:6])
985    covData = {}
986    if len(DisAglData.get('covData',{})):
987        covData = DisAglData['covData']
988        pfx = str(DisAglData['pId'])+'::'
989        A = G2lat.cell2A(Cell[:6])
990        cellSig = G2stIO.getCellEsd(pfx,SGData,A,covData)
991        names = [' a = ',' b = ',' c = ',' alpha = ',' beta = ',' gamma = ',' Volume = ']
992        valEsd = [G2mth.ValEsd(Cell[i],cellSig[i],True) for i in range(7)]
993        line = '\n Unit cell:'
994        for name,vals in zip(names,valEsd):
995            line += name+vals
996        MyPrint(line)
997    else:
998        MyPrint('\n Unit cell: a = '+('%.5f'%Cell[0])+' b = '+('%.5f'%Cell[1])+' c = '+('%.5f'%Cell[2])+
999            ' alpha = '+('%.3f'%Cell[3])+' beta = '+('%.3f'%Cell[4])+' gamma = '+
1000            ('%.3f'%Cell[5])+' Volume = '+('%.3f'%Cell[6]))
1001
1002    AtomLabels,DistArray,AngArray = RetDistAngle(DisAglCtls,DisAglData)
1003    origAtoms = DisAglData['OrigAtoms']
1004    for Oatom in origAtoms:
1005        i = Oatom[0]
1006        Dist = DistArray[i]
1007        nDist = len(Dist)
1008        angles = np.zeros((nDist,nDist))
1009        angsig = np.zeros((nDist,nDist))
1010        for k,j,tup in AngArray[i]:
1011            angles[k][j],angsig[k][j] = angles[j][k],angsig[j][k] = tup
1012        line = ''
1013        for i,x in enumerate(Oatom[3:6]):
1014            line += ('%12.5f'%x).rstrip('0')
1015        MyPrint('\n Distances & angles for '+Oatom[1]+' at '+line.rstrip())
1016        MyPrint(80*'*')
1017        line = ''
1018        for dist in Dist[:-1]:
1019            line += '%12s'%(AtomLabels[dist[0]].center(12))
1020        MyPrint('  To       cell +(sym. op.)      dist.  '+line.rstrip())
1021        for i,dist in enumerate(Dist):
1022            line = ''
1023            for j,angle in enumerate(angles[i][0:i]):
1024                sig = angsig[i][j]
1025                if angle:
1026                    if sig:
1027                        line += '%12s'%(G2mth.ValEsd(angle,sig,True).center(12))
1028                    else:
1029                        val = '%.3f'%(angle)
1030                        line += '%12s'%(val.center(12))
1031                else:
1032                    line += 12*' '
1033            if dist[4]:            #sig exists!
1034                val = G2mth.ValEsd(dist[3],dist[4])
1035            else:
1036                val = '%8.4f'%(dist[3])
1037            tunit = '[%2d%2d%2d]'% dist[1]
1038            MyPrint((%8s%10s+(%4d) %12s'%(AtomLabels[dist[0]].ljust(8),tunit.ljust(10),dist[2],val.center(12)))+line.rstrip())
1039
1040def DisAglTor(DATData):
1041    'Needs a doc string'
1042    SGData = DATData['SGData']
1043    Cell = DATData['Cell']
1044
1045    Amat,Bmat = G2lat.cell2AB(Cell[:6])
1046    covData = {}
1047    pfx = ''
1048    if 'covData' in DATData:
1049        covData = DATData['covData']
1050        pfx = str(DATData['pId'])+'::'
1051    Datoms = []
1052    Oatoms = []
1053    for i,atom in enumerate(DATData['Datoms']):
1054        symop = atom[-1].split('+')
1055        if len(symop) == 1:
1056            symop.append('0,0,0')
1057        symop[0] = int(symop[0])
1058        symop[1] = eval(symop[1])
1059        atom.append(symop)
1060        Datoms.append(atom)
1061        oatom = DATData['Oatoms'][i]
1062        names = ['','','']
1063        if pfx:
1064            names = [pfx+'dAx:'+str(oatom[0]),pfx+'dAy:'+str(oatom[0]),pfx+'dAz:'+str(oatom[0])]
1065        oatom += [names,]
1066        Oatoms.append(oatom)
1067    atmSeq = [atom[1]+'('+atom[-2]+')' for atom in Datoms]
1068    if DATData['Natoms'] == 4:  #torsion
1069        Tors,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData)
1070        G2fil.G2Print (' Torsion angle for %s atom sequence: %s = %s'%(DATData['Name'],str(atmSeq).replace("'","")[1:-1],G2mth.ValEsd(Tors,sig)))
1071        G2fil.G2Print (' NB: Atom sequence determined by selection order')
1072        return      # done with torsion
1073    elif DATData['Natoms'] == 3:  #angle
1074        Ang,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData)
1075        G2fil.G2Print (' Angle in %s for atom sequence: %s = %s'%(DATData['Name'],str(atmSeq).replace("'","")[1:-1],G2mth.ValEsd(Ang,sig)))
1076        G2fil.G2Print (' NB: Atom sequence determined by selection order')
1077    else:   #2 atoms - distance
1078        Dist,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData)
1079        G2fil.G2Print (' Distance in %s for atom sequence: %s = %s'%(DATData['Name'],str(atmSeq).replace("'","")[1:-1],G2mth.ValEsd(Dist,sig)))
1080
1081def BestPlane(PlaneData):
1082    'Needs a doc string'
1083
1084    def ShowBanner(name):
1085        G2fil.G2Print (80*'*')
1086        G2fil.G2Print ('   Best plane result for phase '+name)
1087        G2fil.G2Print (80*'*','\n')
1088
1089    ShowBanner(PlaneData['Name'])
1090
1091    Cell = PlaneData['Cell']
1092    Amat,Bmat = G2lat.cell2AB(Cell[:6])
1093    Atoms = PlaneData['Atoms']
1094    sumXYZ = np.zeros(3)
1095    XYZ = []
1096    Natoms = len(Atoms)
1097    for atom in Atoms:
1098        xyz = np.array(atom[3:6])
1099        XYZ.append(xyz)
1100        sumXYZ += xyz
1101    sumXYZ /= Natoms
1102    XYZ = np.array(XYZ)-sumXYZ
1103    XYZ = np.inner(Amat,XYZ).T
1104    Zmat = np.zeros((3,3))
1105    for i,xyz in enumerate(XYZ):
1106        Zmat += np.outer(xyz.T,xyz)
1107    G2fil.G2Print (' Selected atoms centered at %10.5f %10.5f %10.5f'%(sumXYZ[0],sumXYZ[1],sumXYZ[2]))
1108    Evec,Emat = nl.eig(Zmat)
1109    Evec = np.sqrt(Evec)/(Natoms-3)
1110    Order = np.argsort(Evec)
1111    XYZ = np.inner(XYZ,Emat.T).T
1112    XYZ = np.array([XYZ[Order[2]],XYZ[Order[1]],XYZ[Order[0]]]).T
1113    G2fil.G2Print (' Atoms in Cartesian best plane coordinates:')
1114    G2fil.G2Print (' Name         X         Y         Z')
1115    for i,xyz in enumerate(XYZ):
1116        G2fil.G2Print (' %6s%10.3f%10.3f%10.3f'%(Atoms[i][1].ljust(6),xyz[0],xyz[1],xyz[2]))
1117    G2fil.G2Print ('\n Best plane RMS X =%8.3f, Y =%8.3f, Z =%8.3f'%(Evec[Order[2]],Evec[Order[1]],Evec[Order[0]]))
1118
1119def main():
1120    'Called to run a refinement when this module is executed '
1121    starttime = time.time()
1122    arg = sys.argv
1123    if len(arg) > 1:
1124        GPXfile = arg[1]
1125        if not ospath.exists(GPXfile):
1126            G2fil.G2Print ('ERROR - '+GPXfile+" doesn't exist!")
1127            exit()
1128    else:
1129        G2fil.G2Print ('ERROR - missing filename')
1130        exit()
1131    # TODO: figure out if this is a sequential refinement and call SeqRefine(GPXfile,None)
1132    Refine(GPXfile,None)
1133    G2fil.G2Print("Done. Execution time {:.2f} sec.".format(time.time()-starttime))
1134
1135if __name__ == '__main__':
1136    GSASIIpath.InvokeDebugOpts()
1137    main()
Note: See TracBrowser for help on using the repository browser.