source: trunk/GSASIIstrMain.py @ 4691

Last change on this file since 4691 was 4691, checked in by toby, 3 years ago

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