source: trunk/GSASIIstrMain.py @ 4783

Last change on this file since 4783 was 4783, checked in by toby, 11 months ago

wx4.1 updates; new CIF data item; scripting: peak fit details returned; improve multistring gui

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