source: trunk/GSASIIstrMain.py @ 4840

Last change on this file since 4840 was 4840, checked in by toby, 9 months ago

LeBail? issues: bug on zero cycles fixed; new LeBail? fit menu item; HessianLSQ: remove extra func eval on 0 cycles; fix noprint options; misc cleanups

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 57.9 KB
Line 
1# -*- coding: utf-8 -*-
2'''
3*GSASIIstrMain: main structure routine*
4---------------------------------------
5
6'''
7########### SVN repository information ###################
8# $Date: 2021-03-05 19:54:34 +0000 (Fri, 05 Mar 2021) $
9# $Author: toby $
10# $Revision: 4840 $
11# $URL: trunk/GSASIIstrMain.py $
12# $Id: GSASIIstrMain.py 4840 2021-03-05 19:54:34Z 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: 4840 $")
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 'Ouch#4' in  result[2]:
191                Rvals['Aborted'] = True
192            if 'msg' in result[2]:
193                Rvals['msg'] = result[2]['msg']
194            Controls['Marquardt'] = -3  #reset to default
195            if 'chisq0' in result[2] and chisq0 is None:
196                chisq0 = result[2]['chisq0']
197        elif 'Hessian SVD' in Controls['deriv type']:
198            maxCyc = Controls['max cyc']
199            result = G2mth.HessianSVD(G2stMth.errRefine,values,Hess=G2stMth.HessRefine,ftol=Ftol,xtol=Xtol,maxcyc=maxCyc,Print=ifPrint,
200                args=([Histograms,Phases,restraintDict,rigidbodyDict],parmDict,varyList,calcControls,pawleyLookup,dlg),
201                refPlotUpdate=refPlotUpdate)
202            if result[1] is None:
203                IfOK = False
204                covMatrix = []
205                sig = len(varyList)*[None,]
206                break
207            ncyc = result[2]['num cyc']+1
208            if 'chisq0' in result[2] and chisq0 is None:
209                chisq0 = result[2]['chisq0']
210        else:           #'numeric'
211            result = so.leastsq(G2stMth.errRefine,values,full_output=True,ftol=Ftol,epsfcn=1.e-8,factor=Factor,
212                args=([Histograms,Phases,restraintDict,rigidbodyDict],parmDict,varyList,calcControls,pawleyLookup,dlg))
213            ncyc = 1
214            if len(varyList):
215                ncyc = int(result[2]['nfev']/len(varyList))
216            if refPlotUpdate is not None: refPlotUpdate(Histograms)   # update plot
217#        table = dict(zip(varyList,zip(values,result[0],(result[0]-values))))
218#        for item in table: print item,table[item]               #useful debug - are things shifting?
219        runtime = time.time()-begin
220        Rvals['SVD0'] = result[2].get('SVD0',0)
221        Rvals['converged'] = result[2].get('Converged')
222        Rvals['DelChi2'] = result[2].get('DelChi2',-1.)
223        Rvals['chisq'] = np.sum(result[2]['fvec']**2)
224        G2stMth.Values2Dict(parmDict, varyList, result[0])
225        G2mv.Dict2Map(parmDict,varyList)
226        Rvals['Nobs'] = Histograms['Nobs']
227        Rvals['Rwp'] = np.sqrt(Rvals['chisq']/Histograms['sumwYo'])*100.      #to %
228        Rvals['GOF'] = np.sqrt(Rvals['chisq']/(Histograms['Nobs']-len(varyList)))
229        printFile.write(' Number of function calls: %d No. of observations: %d No. of parameters: %d User rejected: %d Sp. gp. extinct: %d\n'%  \
230            (result[2]['nfev'],Histograms['Nobs'],len(varyList),Histograms['Nrej'],Histograms['Next']))
231        if ncyc:
232            printFile.write(' Refinement time = %8.3fs, %8.3fs/cycle, for %d cycles\n'%(runtime,runtime/ncyc,ncyc))
233        printFile.write(' wR = %7.2f%%, chi**2 = %12.6g, GOF = %6.2f\n'%(Rvals['Rwp'],Rvals['chisq'],Rvals['GOF']))
234        sig = len(varyList)*[None,]
235        if 'None' in str(type(result[1])) and ifSeq:    #this bails out of a sequential refinement on singular matrix
236            IfOK = False
237            covMatrix = []
238            G2fil.G2Print ('Warning: **** Refinement failed - singular matrix ****')
239            if 'Hessian' in Controls['deriv type']:
240                num = len(varyList)-1
241                # BHT -- I am not sure if this works correctly:
242                for i,val in enumerate(np.flipud(result[2]['psing'])):
243                    if val:
244                        G2fil.G2Print('Bad parameter: '+varyList[num-i],mode='warn')
245            else:
246                Ipvt = result[2]['ipvt']
247                for i,ipvt in enumerate(Ipvt):
248                    if not np.sum(result[2]['fjac'],axis=1)[i]:
249                        G2fil.G2Print('Bad parameter: '+varyList[ipvt-1],mode='warn')
250            break
251        IfOK = True
252        try:
253            covMatrix = result[1]*Rvals['GOF']**2
254            sig = np.sqrt(np.diag(covMatrix))
255            Lastshft = result[2].get('Xvec',None)
256            if Lastshft is None:
257                Rvals['Max shft/sig'] = None
258            else:
259                Rvals['Max shft/sig'] = np.max(np.nan_to_num(Lastshft/sig))
260            if np.any(np.isnan(sig)) or not sig.shape:
261                G2fil.G2Print ('*** Least squares aborted - some invalid esds possible ***',mode='error')
262#            table = dict(zip(varyList,zip(values,result[0],(result[0]-values)/sig)))
263#            for item in table: print item,table[item]               #useful debug - are things shifting?
264            # report on refinement issues. Result in Rvals['msg']
265            ReportProblems(result,Rvals,varyList)
266            break                   #refinement succeeded - finish up!
267        except TypeError:
268            # if we get here, no result[1] (covar matrix) was returned or other calc error: refinement failed
269            IfOK = False
270            if not len(varyList):
271                covMatrix = []
272                break
273            if 'Hessian' in Controls['deriv type']:
274                SVD0 = result[2].get('SVD0')
275                if SVD0 == -1:
276                    G2fil.G2Print ('**** Refinement failed - singular matrix ****',mode='error')
277                elif SVD0 == -2:
278                    G2fil.G2Print ('**** Refinement failed - other problem ****',mode='error')
279                elif SVD0 > 0:
280                    G2fil.G2Print ('**** Refinement failed with {} SVD singularities ****'.format(SVD0),mode='error')
281                else:
282                    G2fil.G2Print ('**** Refinement failed ****',mode='error')
283                if result[1] is None:
284                    IfOK = False
285                    covMatrix = []
286                    sig = len(varyList)*[None,]
287                # report on highly correlated variables
288                ReportProblems(result,Rvals,varyList)
289                # process singular variables
290                if dlg: break # refining interactively
291            else:
292                G2fil.G2Print ('**** Refinement failed - singular matrix ****',mode='error')
293                Ipvt = result[2]['ipvt']
294                for i,ipvt in enumerate(Ipvt):
295                    if not np.sum(result[2]['fjac'],axis=1)[i]:
296                        G2fil.G2Print ('Removing parameter: '+varyList[ipvt-1])
297                        del(varyList[ipvt-1])
298                        break
299    if IfOK:
300        if CheckLeBail(Phases):   # only needed for LeBail extraction
301            G2stMth.errRefine(values,[Histograms,Phases,restraintDict,rigidbodyDict],
302                                parmDict,varyList,calcControls,pawleyLookup,dlg)
303
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 CheckLeBail(Phases):
503    '''Check if there is a LeBail extraction in any histogram
504
505    :returns: True if there is at least one LeBail flag turned on, False otherwise
506    '''
507    for key in Phases:
508        phase = Phases[key]
509        for h in phase['Histograms']:
510            #phase['Histograms'][h]
511            if not phase['Histograms'][h]['Use']: continue
512            if phase['Histograms'][h]['LeBail']:
513                 return True
514    return False
515       
516def DoLeBail(GPXfile,dlg=None,cycles=3,refPlotUpdate=None):
517    '''Fit LeBail intensities without changes to any other refined parameters.
518    This is a stripped-down version of :func:`Refine` that does not perform
519    any refinement cycles
520    '''
521    import GSASIImpsubs as G2mp
522    G2mp.InitMP()
523    import pytexture as ptx
524    ptx.pyqlmninit()            #initialize fortran arrays for spherical harmonics
525
526    #varyList = []
527    parmDict = {}
528    Controls = G2stIO.GetControls(GPXfile)
529    calcControls = {}
530    calcControls.update(Controls)
531    constrDict,fixedList = G2stIO.GetConstraints(GPXfile)
532    restraintDict = {}
533    Histograms,Phases = G2stIO.GetUsedHistogramsAndPhases(GPXfile)
534    if not Phases:
535        G2fil.G2Print (' *** ERROR - you have no phases to refine! ***')
536        return False,{'msg':'No phases'}
537    if not Histograms:
538        G2fil.G2Print (' *** ERROR - you have no data to refine with! ***')
539        return False,{'msg':'No data'}
540    if not CheckLeBail(Phases):
541        msg = 'Warning: There are no histograms with LeBail extraction enabled'
542        G2fil.G2Print ('*** '+msg+' ***')
543        return False,{'msg': msg}
544    rigidbodyDict = G2stIO.GetRigidBodies(GPXfile)
545    rbIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})
546    rbVary,rbDict = G2stIO.GetRigidBodyModels(rigidbodyDict,Print=False)
547    (Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables,MFtables,
548         maxSSwave) = G2stIO.GetPhaseData(Phases,restraintDict,rbIds,Print=False)
549    calcControls['atomIndx'] = atomIndx
550    calcControls['Natoms'] = Natoms
551    calcControls['FFtables'] = FFtables
552    calcControls['BLtables'] = BLtables
553    calcControls['MFtables'] = MFtables
554    calcControls['maxSSwave'] = maxSSwave
555    hapVary,hapDict,controlDict = G2stIO.GetHistogramPhaseData(Phases,Histograms,Print=False)
556    calcControls.update(controlDict)
557    histVary,histDict,controlDict = G2stIO.GetHistogramData(Histograms,Print=False)
558    calcControls.update(controlDict)
559    parmDict.update(rbDict)
560    parmDict.update(phaseDict)
561    parmDict.update(hapDict)
562    parmDict.update(histDict)
563    G2stIO.GetFprime(calcControls,Histograms)
564    parmFrozenList = Controls['parmFrozen']['FrozenList']
565    try:
566        for i in range(cycles):
567            M = G2stMth.errRefine([],[Histograms,Phases,restraintDict,rigidbodyDict],parmDict,[],calcControls,pawleyLookup,dlg)
568            G2stMth.GetFobsSq(Histograms,Phases,parmDict,calcControls)
569            if refPlotUpdate is not None: refPlotUpdate(Histograms,i)
570        Rvals = {}
571        Rvals['chisq'] = np.sum(M**2)
572        Rvals['Nobs'] = Histograms['Nobs']
573        Rvals['Rwp'] = np.sqrt(Rvals['chisq']/Histograms['sumwYo'])*100.      #to %
574        Rvals['GOF'] = np.sqrt(Rvals['chisq']/(Histograms['Nobs'])) # no variables
575
576        covData = {'variables':0,'varyList':[],'sig':[],'Rvals':Rvals,
577                       'varyListStart':[],
578                       'covMatrix':np.zeros([0,0]),'title':GPXfile,
579                       #'newAtomDict':newAtomDict,'newCellDict':newCellDict,
580                       'freshCOV':True}
581       
582        G2stIO.SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases,rigidbodyDict,covData,parmFrozenList,True)
583        G2fil.G2Print (' ***** LeBail fit completed *****')
584        return True,Rvals
585    except Exception as Msg:
586        G2fil.G2Print (' ***** LeBail fit error *****')
587        if not hasattr(Msg,'msg'): Msg.msg = str(Msg)
588        if GSASIIpath.GetConfigValue('debug'):
589            import traceback
590            print(traceback.format_exc())       
591        return False,{'msg':Msg.msg}
592
593def phaseCheck(phaseVary,Phases,histogram):
594    '''
595    Removes unused parameters from phase varylist if phase not in histogram
596    for seq refinement removes vars in "Fix FXU" and "FixedSeqVars" here
597    '''
598    NewVary = []
599    for phase in Phases:
600        if histogram not in Phases[phase]['Histograms']: continue
601        if Phases[phase]['Histograms'][histogram]['Use']:
602            pId = Phases[phase]['pId']
603            newVary = [item for item in phaseVary if item.split(':')[0] == str(pId)]
604            FixVals = Phases[phase]['Histograms'][histogram].get('Fix FXU',' ')
605            if 'F' in FixVals:
606                newVary = [item for item in newVary if not 'Afrac' in item]
607            if 'X' in FixVals:
608                newVary = [item for item in newVary if not 'dA' in item]
609            if 'U' in FixVals:
610                newVary = [item for item in newVary if not 'AU' in item]
611            if 'M' in FixVals:
612                newVary = [item for item in newVary if not 'AM' in item]
613            removeVars = Phases[phase]['Histograms'][histogram].get('FixedSeqVars',[])
614            newVary = [item for item in newVary if item not in removeVars]
615            NewVary += newVary
616    return NewVary
617
618def SeqRefine(GPXfile,dlg,refPlotUpdate=None):
619    '''Perform a sequential refinement -- cycles through all selected histgrams,
620    one at a time
621    '''
622    import GSASIImpsubs as G2mp
623    G2mp.InitMP()
624    import pytexture as ptx
625    ptx.pyqlmninit()            #initialize fortran arrays for spherical harmonics
626    msgs = {}
627    printFile = open(ospath.splitext(GPXfile)[0]+'.lst','w')
628    G2fil.G2Print ('Starting Sequential Refinement')
629    G2stIO.ShowBanner(printFile)
630    Controls = G2stIO.GetControls(GPXfile)
631    preFrozenCount = 0
632    for h in Controls['parmFrozen']:
633        if h == 'FrozenList': continue
634        preFrozenCount += len(Controls['parmFrozen'][h])   
635    G2stIO.ShowControls(Controls,printFile,SeqRef=True,
636                            preFrozenCount=preFrozenCount)
637    restraintDict = G2stIO.GetRestraints(GPXfile)
638    Histograms,Phases = G2stIO.GetUsedHistogramsAndPhases(GPXfile)
639    if not Phases:
640        G2fil.G2Print (' *** ERROR - you have no phases to refine! ***')
641        G2fil.G2Print (' *** Refine aborted ***')
642        return False,'No phases'
643    if not Histograms:
644        G2fil.G2Print (' *** ERROR - you have no data to refine with! ***')
645        G2fil.G2Print (' *** Refine aborted ***')
646        return False,'No data'
647    rigidbodyDict = G2stIO.GetRigidBodies(GPXfile)
648    rbIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})
649    rbVary,rbDict = G2stIO.GetRigidBodyModels(rigidbodyDict,pFile=printFile)
650    G2mv.InitVars()
651    (Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables,MFtables,
652         maxSSwave) = G2stIO.GetPhaseData(Phases,restraintDict,rbIds,
653                                    Print=False,pFile=printFile,seqRef=True)
654    for item in phaseVary:
655        if '::A0' in item:
656            G2fil.G2Print ('**** WARNING - lattice parameters should not be refined in a sequential refinement ****')
657            G2fil.G2Print ('****           instead use the Dij parameters for each powder histogram            ****')
658            return False,'Lattice parameter refinement error - see console message'
659        if '::C(' in item:
660            G2fil.G2Print ('**** WARNING - phase texture parameters should not be refined in a sequential refinement ****')
661            G2fil.G2Print ('****           instead use the C(L,N) parameters for each powder histogram               ****')
662            return False,'Phase texture refinement error - see console message'
663    if 'Seq Data' in Controls:
664        histNames = Controls['Seq Data']
665    else: # patch from before Controls['Seq Data'] was implemented?
666        histNames = G2stIO.GetHistogramNames(GPXfile,['PWDR',])
667    if Controls.get('Reverse Seq'):
668        histNames.reverse()
669    SeqResult = G2stIO.GetSeqResult(GPXfile)
670#    SeqResult = {'SeqPseudoVars':{},'SeqParFitEqList':[]}
671    Histo = {}
672    NewparmDict = {}
673    G2stIO.SetupSeqSavePhases(GPXfile)
674    msgs['steepestNum'] = 0
675    msgs['maxshift/sigma'] = []
676    for ihst,histogram in enumerate(histNames):
677        if GSASIIpath.GetConfigValue('Show_timing'): t1 = time.time()
678        G2fil.G2Print('\nRefining with '+str(histogram))
679        G2mv.InitVars()
680        (Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,
681             FFtables,BLtables,MFtables,maxSSwave) = G2stIO.GetPhaseData(
682                 Phases,restraintDict,rbIds,
683                 Print=False,pFile=printFile,seqRef=True)
684        ifPrint = False
685        if dlg:
686            dlg.SetTitle('Residual for histogram '+str(ihst))
687        calcControls = {}
688        calcControls['atomIndx'] = atomIndx
689        calcControls['Natoms'] = Natoms
690        calcControls['FFtables'] = FFtables
691        calcControls['BLtables'] = BLtables
692        calcControls['MFtables'] = MFtables
693        calcControls['maxSSwave'] = maxSSwave
694        if histogram not in Histograms:
695            G2fil.G2Print("Error: not found!")
696            raise G2obj.G2Exception("refining with invalid histogram {}".
697                                        format(histogram))
698        hId = Histograms[histogram]['hId']
699        redphaseVary = phaseCheck(phaseVary,Phases,histogram)
700        Histo = {histogram:Histograms[histogram],}
701        hapVary,hapDict,controlDict = G2stIO.GetHistogramPhaseData(Phases,Histo,Print=False)
702        calcControls.update(controlDict)
703        histVary,histDict,controlDict = G2stIO.GetHistogramData(Histo,False)
704        calcControls.update(controlDict)
705        varyList = rbVary+redphaseVary+hapVary+histVary
706#        if not ihst:
707            # save the initial vary list, but without histogram numbers on parameters
708        saveVaryList = varyList[:]
709        for i,item in enumerate(saveVaryList):
710            items = item.split(':')
711            if items[1]:
712                items[1] = ''
713            item = ':'.join(items)
714            saveVaryList[i] = item
715        if not ihst:
716            SeqResult['varyList'] = saveVaryList
717        else:
718            SeqResult['varyList'] = list(set(SeqResult['varyList']+saveVaryList))
719        parmDict = {}
720        parmDict.update(rbDict)
721        parmDict.update(phaseDict)
722        parmDict.update(hapDict)
723        parmDict.update(histDict)
724        if Controls['Copy2Next']:   # update with parms from last histogram
725            #parmDict.update(NewparmDict) # don't use in case extra entries would cause a problem
726            for parm in NewparmDict:
727                if parm in parmDict:
728                    parmDict[parm] = NewparmDict[parm]
729        elif histogram in SeqResult:  # update phase from last seq ref
730            NewparmDict = SeqResult[histogram].get('parmDict',{})
731            for parm in NewparmDict:
732                if '::' in parm and parm in parmDict:
733                    parmDict[parm] = NewparmDict[parm]
734           
735        G2stIO.GetFprime(calcControls,Histo)
736        # do constraint processing
737        #reload(G2mv) # debug
738        constrDict,fixedList = G2stIO.GetConstraints(GPXfile)
739        varyListStart = tuple(varyList) # save the original varyList before dependent vars are removed
740        msg = G2mv.EvaluateMultipliers(constrDict,parmDict)
741        if msg:
742            return False,'Unable to interpret multiplier(s): '+msg
743        try:
744            groups,parmlist = G2mv.GenerateConstraints(varyList,constrDict,fixedList,parmDict,SeqHist=hId)
745#            if GSASIIpath.GetConfigValue('debug'): print("DBG_"+
746#                G2mv.VarRemapShow(varyList,True))
747#            print('DependentVars',G2mv.GetDependentVars())
748#            print('IndependentVars',G2mv.GetIndependentVars())
749            constraintInfo = (groups,parmlist,constrDict,fixedList,ihst)
750        except G2mv.ConstraintException:
751            G2fil.G2Print (' *** ERROR - your constraints are internally inconsistent ***')
752            #errmsg, warnmsg = G2mv.CheckConstraints(varyList,constrDict,fixedList)
753            #print 'Errors',errmsg
754            #if warnmsg: print 'Warnings',warnmsg
755            return False,' Constraint error'
756        #print G2mv.VarRemapShow(varyList)
757        if not ihst:
758            # first histogram to refine against
759            firstVaryList = []
760            for item in varyList:
761                items = item.split(':')
762                if items[1]:
763                    items[1] = ''
764                item = ':'.join(items)
765                firstVaryList.append(item)
766            newVaryList = firstVaryList
767        else:
768            newVaryList = []
769            for item in varyList:
770                items = item.split(':')
771                if items[1]:
772                    items[1] = ''
773                item = ':'.join(items)
774                newVaryList.append(item)
775        if newVaryList != firstVaryList and Controls['Copy2Next']:
776            # variable lists are expected to match between sequential refinements when Copy2Next is on
777            #print '**** ERROR - variable list for this histogram does not match previous'
778            #print '     Copy of variables is not possible'
779            #print '\ncurrent histogram',histogram,'has',len(newVaryList),'variables'
780            combined = list(set(firstVaryList+newVaryList))
781            c = [var for var in combined if var not in newVaryList]
782            p = [var for var in combined if var not in firstVaryList]
783            G2fil.G2Print('*** Variables change ***')
784            for typ,vars in [('Removed',c),('Added',p)]:
785                line = '  '+typ+': '
786                if vars:
787                    for var in vars:
788                        if len(line) > 70:
789                            G2fil.G2Print(line)
790                            line = '    '
791                        line += var + ', '
792                else:
793                        line += 'none, '
794                G2fil.G2Print(line[:-2])
795            firstVaryList = newVaryList
796
797        ifSeq = True
798        printFile.write('\n Refinement results for histogram id {}: {}\n'
799                            .format(hId,histogram))
800        printFile.write(135*'-'+'\n')
801        # remove frozen vars
802        if 'parmFrozen' not in Controls:
803            Controls['parmFrozen'] = {}
804        if histogram not in Controls['parmFrozen']: 
805            Controls['parmFrozen'][histogram] = []
806        parmFrozenList = Controls['parmFrozen'][histogram]
807        frozenList = [i for i in varyList if i in parmFrozenList]
808        if len(frozenList) != 0: 
809           varyList = [i for i in varyList if i not in parmFrozenList]
810           s = ''
811           for a in frozenList:
812               if s:
813                   s+= ', '
814               s += a
815           printFile.write(
816               ' The following refined variables have previously been frozen due to exceeding limits:\n\t{}\n'
817               .format(s))
818        try:
819            IfOK,Rvals,result,covMatrix,sig = RefineCore(Controls,Histo,Phases,restraintDict,
820                rigidbodyDict,parmDict,varyList,calcControls,pawleyLookup,ifSeq,printFile,dlg,
821                refPlotUpdate=refPlotUpdate)
822            try:
823                shft = '%.4f'% Rvals['Max shft/sig']
824            except:
825                shft = '?'
826            G2fil.G2Print ('  wR = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f, last delta chi = %.4f, last shft/sig = %s'%(
827                Rvals['Rwp'],Rvals['chisq'],Rvals['GOF']**2,Rvals['DelChi2'],shft))
828            if Rvals.get('lamMax',0) >= 10.:
829                msgs['steepestNum'] += 1
830            if Rvals.get('Max shft/sig'):
831                msgs['maxshift/sigma'].append(Rvals['Max shft/sig'])
832            # add the uncertainties into the esd dictionary (sigDict)
833            if not IfOK:
834                G2fil.G2Print('***** Sequential refinement failed at histogram '+histogram,mode='warn')
835                break
836            sigDict = dict(zip(varyList,sig))
837            # the uncertainties for dependent constrained parms into the esd dict
838            sigDict.update(G2mv.ComputeDepESD(covMatrix,varyList,parmDict))
839            # check for variables outside their allowed range, reset and freeze them
840            frozen = dropOOBvars(varyList,parmDict,sigDict,Controls,parmFrozenList)
841            msg = None
842            if len(frozen) > 0:
843               msg = ('Hist {}: {} variables were outside limits and were frozen (now {} frozen total)'
844                   .format(ihst,len(frozen),len(parmFrozenList)))
845               G2fil.G2Print(msg)
846               msg = (' {} variables were outside limits and were frozen (now {} frozen total)'
847                   .format(len(frozen),len(parmFrozenList)))
848               for p in frozen:
849                   if p not in varyList:
850                       print('Frozen Warning: {} not in varyList. This should not happen!'.format(p))
851                       continue
852                   i = varyList.index(p)
853                   result[0][i] = parmDict[p]
854                   sig[i] = -0.1
855            # a dict with values & esds for dependent (constrained) parameters - avoid extraneous holds
856            depParmDict = {i:(parmDict[i],sigDict[i]) for i in varyListStart if i in sigDict and i not in varyList}
857            newCellDict = copy.deepcopy(G2stMth.GetNewCellParms(parmDict,varyList))
858            newAtomDict = copy.deepcopy(G2stMth.ApplyXYZshifts(parmDict,varyList))
859            histRefData = {
860                'variables':result[0],'varyList':varyList,'sig':sig,'Rvals':Rvals,
861                'varyListStart':varyListStart,
862                'covMatrix':covMatrix,'title':histogram,'newAtomDict':newAtomDict,
863                'newCellDict':newCellDict,'depParmDict':depParmDict,
864                'constraintInfo':constraintInfo,
865                'parmDict':parmDict,
866                }
867            SeqResult[histogram] = histRefData
868            G2stMth.ApplyRBModels(parmDict,Phases,rigidbodyDict,True)
869#            G2stIO.SetRigidBodyModels(parmDict,sigDict,rigidbodyDict,printFile)
870            G2stIO.SetHistogramPhaseData(parmDict,sigDict,Phases,Histo,None,ifPrint,printFile)
871            G2stIO.SetHistogramData(parmDict,sigDict,Histo,None,ifPrint,printFile,seq=True)
872            G2stIO.SaveUpdatedHistogramsAndPhases(GPXfile,Histo,Phases,rigidbodyDict,histRefData,Controls['parmFrozen'])
873            if msg: 
874                printFile.write(msg+'\n')
875            NewparmDict = {}
876            # make dict of varied parameters in current histogram, renamed to
877            # next histogram, for use in next refinement.
878            if Controls['Copy2Next'] and ihst < len(histNames)-1:
879                hId = Histo[histogram]['hId'] # current histogram
880                nexthId = Histograms[histNames[ihst+1]]['hId']
881                for parm in set(list(varyList)+list(varyListStart)):
882                    items = parm.split(':')
883                    if len(items) < 3: 
884                        continue
885                    if str(hId) in items[1]:
886                        items[1] = str(nexthId)
887                        newparm = ':'.join(items)
888                        NewparmDict[newparm] = parmDict[parm]
889                    else:
890                        if items[2].startswith('dA'): parm = parm.replace(':dA',':A') 
891                        NewparmDict[parm] = parmDict[parm]
892                   
893        except G2obj.G2RefineCancel as Msg:
894            if not hasattr(Msg,'msg'): Msg.msg = str(Msg)
895            printFile.close()
896            G2fil.G2Print (' ***** Refinement stopped *****')
897            return False,Msg.msg
898        except G2obj.G2Exception as Msg:  # cell metric error, others?
899            if not hasattr(Msg,'msg'): Msg.msg = str(Msg)
900            printFile.close()
901            G2fil.G2Print (' ***** Refinement error *****')
902            return False,Msg.msg
903        if GSASIIpath.GetConfigValue('Show_timing'):
904            t2 = time.time()
905            G2fil.G2Print("Fit step time {:.2f} sec.".format(t2-t1))
906            t1 = t2
907    SeqResult['histNames'] = [itm for itm in G2stIO.GetHistogramNames(GPXfile,['PWDR',]) if itm in SeqResult.keys()]
908    try:
909        G2stIO.SetSeqResult(GPXfile,Histograms,SeqResult)
910    except Exception as msg:
911        print('Error reading Sequential results\n',str(msg))
912        if GSASIIpath.GetConfigValue('debug'):
913            import traceback
914            print(traceback.format_exc())       
915    postFrozenCount = 0
916    for h in Controls['parmFrozen']:
917        if h == 'FrozenList': continue
918        postFrozenCount += len(Controls['parmFrozen'][h])
919    if postFrozenCount:
920        msgs['Frozen'] = 'Ending refinement with {} Frozen variables ({} added now)\n'.format(postFrozenCount,postFrozenCount-preFrozenCount)
921        printFile.write('\n'+msgs['Frozen'])
922    printFile.close()
923    G2fil.G2Print (' Sequential refinement results are in file: '+ospath.splitext(GPXfile)[0]+'.lst')
924    G2fil.G2Print (' ***** Sequential refinement successful *****')
925    return True,msgs
926
927def dropOOBvars(varyList,parmDict,sigDict,Controls,parmFrozenList):
928    '''Find variables in the parameters dict that are outside the ranges
929    (in parmMinDict and parmMaxDict) and set them to the limits values.
930    Add any such variables into the list of frozen variable
931    (parmFrozenList). Returns a list of newly frozen variables, if any.
932    '''
933    parmMinDict = Controls.get('parmMinDict',{})
934    parmMaxDict = Controls.get('parmMaxDict',{})
935    freeze = []
936    if parmMinDict or parmMaxDict:
937        for name in varyList:
938            if name not in parmDict: continue
939            n,val = G2obj.prmLookup(name,parmMinDict)
940            if n is not None:
941                if parmDict[name] < parmMinDict[n]:
942                    parmDict[name] = parmMinDict[n]
943                    sigDict[name] = 0.0
944                    freeze.append(name)
945                    continue
946            n,val = G2obj.prmLookup(name,parmMaxDict)
947            if n is not None:
948                if parmDict[name] > parmMaxDict[n]:
949                    parmDict[name] = parmMaxDict[n]
950                    sigDict[name] = 0.0
951                    freeze.append(name)
952                    continue
953        for v in freeze:
954            if v not in parmFrozenList:
955                parmFrozenList.append(v)
956    return freeze
957
958def RetDistAngle(DisAglCtls,DisAglData,dlg=None):
959    '''Compute and return distances and angles
960
961    :param dict DisAglCtls: contains distance/angle radii usually defined using
962       :func:`GSASIIctrlGUI.DisAglDialog`
963    :param dict DisAglData: contains phase data:
964       Items 'OrigAtoms' and 'TargAtoms' contain the atoms to be used
965       for distance/angle origins and atoms to be used as targets.
966       Item 'SGData' has the space group information (see :ref:`Space Group object<SGData_table>`)
967
968    :returns: AtomLabels,DistArray,AngArray where:
969
970      **AtomLabels** is a dict of atom labels, keys are the atom number
971
972      **DistArray** is a dict keyed by the origin atom number where the value is a list
973      of distance entries. The value for each distance is a list containing:
974
975        0) the target atom number (int);
976        1) the unit cell offsets added to x,y & z (tuple of int values)
977        2) the symmetry operator number (which may be modified to indicate centering and center of symmetry)
978        3) an interatomic distance in A (float)
979        4) an uncertainty on the distance in A or 0.0 (float)
980
981      **AngArray** is a dict keyed by the origin (central) atom number where
982      the value is a list of
983      angle entries. The value for each angle entry consists of three values:
984
985        0) a distance item reference for one neighbor (int)
986        1) a distance item reference for a second neighbor (int)
987        2) a angle, uncertainty pair; the s.u. may be zero (tuple of two floats)
988
989      The AngArray distance reference items refer directly to the index of the items in the
990      DistArray item for the list of distances for the central atom.
991    '''
992    import numpy.ma as ma
993
994    SGData = DisAglData['SGData']
995    Cell = DisAglData['Cell']
996
997    Amat,Bmat = G2lat.cell2AB(Cell[:6])
998    covData = {}
999    if len(DisAglData.get('covData',{})):
1000        covData = DisAglData['covData']
1001        covMatrix = covData['covMatrix']
1002        varyList = covData['varyList']
1003        pfx = str(DisAglData['pId'])+'::'
1004
1005    Factor = DisAglCtls['Factors']
1006    Radii = dict(zip(DisAglCtls['AtomTypes'],zip(DisAglCtls['BondRadii'],DisAglCtls['AngleRadii'])))
1007    indices = (-2,-1,0,1,2)
1008    Units = np.array([[h,k,l] for h in indices for k in indices for l in indices])
1009    origAtoms = DisAglData['OrigAtoms']
1010    targAtoms = DisAglData['TargAtoms']
1011    AtomLabels = {}
1012    for Oatom in origAtoms:
1013        AtomLabels[Oatom[0]] = Oatom[1]
1014    for Oatom in targAtoms:
1015        AtomLabels[Oatom[0]] = Oatom[1]
1016    DistArray = {}
1017    AngArray = {}
1018    for iO,Oatom in enumerate(origAtoms):
1019        DistArray[Oatom[0]] = []
1020        AngArray[Oatom[0]] = []
1021        OxyzNames = ''
1022        IndBlist = []
1023        Dist = []
1024        Vect = []
1025        VectA = []
1026        angles = []
1027        for Tatom in targAtoms:
1028            Xvcov = []
1029            TxyzNames = ''
1030            if len(DisAglData.get('covData',{})):
1031                OxyzNames = [pfx+'dAx:%d'%(Oatom[0]),pfx+'dAy:%d'%(Oatom[0]),pfx+'dAz:%d'%(Oatom[0])]
1032                TxyzNames = [pfx+'dAx:%d'%(Tatom[0]),pfx+'dAy:%d'%(Tatom[0]),pfx+'dAz:%d'%(Tatom[0])]
1033                Xvcov = G2mth.getVCov(OxyzNames+TxyzNames,varyList,covMatrix)
1034            result = G2spc.GenAtom(Tatom[3:6],SGData,False,Move=False)
1035            BsumR = (Radii[Oatom[2]][0]+Radii[Tatom[2]][0])*Factor[0]
1036            AsumR = (Radii[Oatom[2]][1]+Radii[Tatom[2]][1])*Factor[1]
1037            for [Txyz,Top,Tunit,Spn] in result:
1038                Dx = (Txyz-np.array(Oatom[3:6]))+Units
1039                dx = np.inner(Amat,Dx)
1040                dist = ma.masked_less(np.sqrt(np.sum(dx**2,axis=0)),0.5)
1041                IndB = ma.nonzero(ma.masked_greater(dist-BsumR,0.))
1042                if np.any(IndB):
1043                    for indb in IndB:
1044                        for i in range(len(indb)):
1045                            if str(dx.T[indb][i]) not in IndBlist:
1046                                IndBlist.append(str(dx.T[indb][i]))
1047                                unit = Units[indb][i]
1048                                tunit = (unit[0]+Tunit[0],unit[1]+Tunit[1],unit[2]+Tunit[2])
1049                                sig = 0.0
1050                                if len(Xvcov):
1051                                    pdpx = G2mth.getDistDerv(Oatom[3:6],Tatom[3:6],Amat,unit,Top,SGData)
1052                                    sig = np.sqrt(np.inner(pdpx,np.inner(pdpx,Xvcov)))
1053                                Dist.append([Oatom[0],Tatom[0],tunit,Top,ma.getdata(dist[indb])[i],sig])
1054                                if (Dist[-1][-2]-AsumR) <= 0.:
1055                                    Vect.append(dx.T[indb][i]/Dist[-1][-2])
1056                                    VectA.append([OxyzNames,np.array(Oatom[3:6]),TxyzNames,np.array(Tatom[3:6]),unit,Top])
1057                                else:
1058                                    Vect.append([0.,0.,0.])
1059                                    VectA.append([])
1060        if dlg is not None:
1061            dlg.Update(iO,newmsg='Atoms done=%d'%(iO))
1062        for D in Dist:
1063            DistArray[Oatom[0]].append(D[1:])
1064        Vect = np.array(Vect)
1065        angles = np.zeros((len(Vect),len(Vect)))
1066        angsig = np.zeros((len(Vect),len(Vect)))
1067        for i,veca in enumerate(Vect):
1068            if np.any(veca):
1069                for j,vecb in enumerate(Vect):
1070                    if np.any(vecb):
1071                        angles[i][j],angsig[i][j] = G2mth.getAngSig(VectA[i],VectA[j],Amat,SGData,covData)
1072                        if i <= j: continue
1073                        AngArray[Oatom[0]].append((i,j,
1074                            G2mth.getAngSig(VectA[i],VectA[j],Amat,SGData,covData)))
1075    return AtomLabels,DistArray,AngArray
1076
1077def PrintDistAngle(DisAglCtls,DisAglData,out=sys.stdout):
1078    '''Print distances and angles
1079
1080    :param dict DisAglCtls: contains distance/angle radii usually defined using
1081       :func:`GSASIIctrlGUI.DisAglDialog`
1082    :param dict DisAglData: contains phase data:
1083       Items 'OrigAtoms' and 'TargAtoms' contain the atoms to be used
1084       for distance/angle origins and atoms to be used as targets.
1085       Item 'SGData' has the space group information (see :ref:`Space Group object<SGData_table>`)
1086    :param file out: file object for output. Defaults to sys.stdout.
1087    '''
1088    def MyPrint(s):
1089        out.write(s+'\n')
1090        # print(s,file=out) # use in Python 3
1091
1092    def ShowBanner(name):
1093        MyPrint(80*'*')
1094        MyPrint('   Interatomic Distances and Angles for phase '+name)
1095        MyPrint((80*'*')+'\n')
1096
1097    ShowBanner(DisAglCtls['Name'])
1098    SGData = DisAglData['SGData']
1099    SGtext,SGtable = G2spc.SGPrint(SGData)
1100    for line in SGtext: MyPrint(line)
1101    if len(SGtable) > 1:
1102        for i,item in enumerate(SGtable[::2]):
1103            if 2*i+1 == len(SGtable):
1104                line = ' %s'%(item.ljust(30))
1105            else:
1106                line = ' %s %s'%(item.ljust(30),SGtable[2*i+1].ljust(30))
1107            MyPrint(line)
1108    else:
1109        MyPrint(' ( 1)    %s'%(SGtable[0])) #triclinic case
1110    Cell = DisAglData['Cell']
1111
1112    Amat,Bmat = G2lat.cell2AB(Cell[:6])
1113    covData = {}
1114    if len(DisAglData.get('covData',{})):
1115        covData = DisAglData['covData']
1116        pfx = str(DisAglData['pId'])+'::'
1117        A = G2lat.cell2A(Cell[:6])
1118        cellSig = G2stIO.getCellEsd(pfx,SGData,A,covData)
1119        names = [' a = ',' b = ',' c = ',' alpha = ',' beta = ',' gamma = ',' Volume = ']
1120        valEsd = [G2mth.ValEsd(Cell[i],cellSig[i],True) for i in range(7)]
1121        line = '\n Unit cell:'
1122        for name,vals in zip(names,valEsd):
1123            line += name+vals
1124        MyPrint(line)
1125    else:
1126        MyPrint('\n Unit cell: a = '+('%.5f'%Cell[0])+' b = '+('%.5f'%Cell[1])+' c = '+('%.5f'%Cell[2])+
1127            ' alpha = '+('%.3f'%Cell[3])+' beta = '+('%.3f'%Cell[4])+' gamma = '+
1128            ('%.3f'%Cell[5])+' Volume = '+('%.3f'%Cell[6]))
1129
1130    AtomLabels,DistArray,AngArray = RetDistAngle(DisAglCtls,DisAglData)
1131    origAtoms = DisAglData['OrigAtoms']
1132    for Oatom in origAtoms:
1133        i = Oatom[0]
1134        Dist = DistArray[i]
1135        nDist = len(Dist)
1136        angles = np.zeros((nDist,nDist))
1137        angsig = np.zeros((nDist,nDist))
1138        for k,j,tup in AngArray[i]:
1139            angles[k][j],angsig[k][j] = angles[j][k],angsig[j][k] = tup
1140        line = ''
1141        for i,x in enumerate(Oatom[3:6]):
1142            line += ('%12.5f'%x).rstrip('0')
1143        MyPrint('\n Distances & angles for '+Oatom[1]+' at '+line.rstrip())
1144        MyPrint(80*'*')
1145        line = ''
1146        for dist in Dist[:-1]:
1147            line += '%12s'%(AtomLabels[dist[0]].center(12))
1148        MyPrint('  To       cell +(sym. op.)      dist.  '+line.rstrip())
1149        BVS = {}
1150        BVdat = {}
1151        Otyp = G2elem.FixValence(Oatom[2]).split('+')[0].split('-')[0]
1152        BVox = [BV for BV in atmdata.BVSoxid[Otyp] if '+' in BV]
1153        if len(BVox):
1154            BVS = {BV:0.0 for BV in BVox}
1155            BVdat = {BV:dict(zip(['O','F','Cl'],atmdata.BVScoeff[BV])) for BV in BVox}
1156            pvline = 'Bond Valence sums for: '
1157        for i,dist in enumerate(Dist):
1158            line = ''
1159            for j,angle in enumerate(angles[i][0:i]):
1160                sig = angsig[i][j]
1161                if angle:
1162                    if sig:
1163                        line += '%12s'%(G2mth.ValEsd(angle,sig,True).center(12))
1164                    else:
1165                        val = '%.3f'%(angle)
1166                        line += '%12s'%(val.center(12))
1167                else:
1168                    line += 12*' '
1169            if dist[4]:            #sig exists!
1170                val = G2mth.ValEsd(dist[3],dist[4])
1171            else:
1172                val = '%8.4f'%(dist[3])
1173            if len(BVox):
1174                Tatm = G2elem.FixValence(DisAglData['TargAtoms'][dist[0]][2]).split('-')[0]
1175                if Tatm in ['O','F','Cl']:
1176                    for BV in BVox:
1177                        BVS[BV] += np.exp((BVdat[BV][Tatm]-dist[3])/0.37)               
1178            tunit = '[%2d%2d%2d]'% dist[1]
1179            MyPrint((%8s%10s+(%4d) %12s'%(AtomLabels[dist[0]].ljust(8),tunit.ljust(10),dist[2],val.center(12)))+line.rstrip())
1180        if len(BVox):
1181            MyPrint(80*'*')
1182            for BV in BVox:
1183                pvline += ' %s: %.2f  '%(BV,BVS[BV])
1184            MyPrint(pvline)
1185
1186def DisAglTor(DATData):
1187    'Needs a doc string'
1188    SGData = DATData['SGData']
1189    Cell = DATData['Cell']
1190
1191    Amat,Bmat = G2lat.cell2AB(Cell[:6])
1192    covData = {}
1193    pfx = ''
1194    if 'covData' in DATData:
1195        covData = DATData['covData']
1196        pfx = str(DATData['pId'])+'::'
1197    Datoms = []
1198    Oatoms = []
1199    for i,atom in enumerate(DATData['Datoms']):
1200        symop = atom[-1].split('+')
1201        if len(symop) == 1:
1202            symop.append('0,0,0')
1203        symop[0] = int(symop[0])
1204        symop[1] = eval(symop[1])
1205        atom.append(symop)
1206        Datoms.append(atom)
1207        oatom = DATData['Oatoms'][i]
1208        names = ['','','']
1209        if pfx:
1210            names = [pfx+'dAx:'+str(oatom[0]),pfx+'dAy:'+str(oatom[0]),pfx+'dAz:'+str(oatom[0])]
1211        oatom += [names,]
1212        Oatoms.append(oatom)
1213    atmSeq = [atom[1]+'('+atom[-2]+')' for atom in Datoms]
1214    if DATData['Natoms'] == 4:  #torsion
1215        Tors,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData)
1216        G2fil.G2Print (' Torsion angle for %s atom sequence: %s = %s'%(DATData['Name'],str(atmSeq).replace("'","")[1:-1],G2mth.ValEsd(Tors,sig)))
1217        G2fil.G2Print (' NB: Atom sequence determined by selection order')
1218        return      # done with torsion
1219    elif DATData['Natoms'] == 3:  #angle
1220        Ang,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData)
1221        G2fil.G2Print (' Angle in %s for atom sequence: %s = %s'%(DATData['Name'],str(atmSeq).replace("'","")[1:-1],G2mth.ValEsd(Ang,sig)))
1222        G2fil.G2Print (' NB: Atom sequence determined by selection order')
1223    else:   #2 atoms - distance
1224        Dist,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData)
1225        G2fil.G2Print (' Distance in %s for atom sequence: %s = %s'%(DATData['Name'],str(atmSeq).replace("'","")[1:-1],G2mth.ValEsd(Dist,sig)))
1226
1227def BestPlane(PlaneData):
1228    'Needs a doc string'
1229
1230    def ShowBanner(name):
1231        G2fil.G2Print (80*'*')
1232        G2fil.G2Print ('   Best plane result for phase '+name)
1233        G2fil.G2Print (80*'*','\n')
1234
1235    ShowBanner(PlaneData['Name'])
1236
1237    Cell = PlaneData['Cell']
1238    Amat,Bmat = G2lat.cell2AB(Cell[:6])
1239    Atoms = PlaneData['Atoms']
1240    sumXYZ = np.zeros(3)
1241    XYZ = []
1242    Natoms = len(Atoms)
1243    for atom in Atoms:
1244        xyz = np.array(atom[3:6])
1245        XYZ.append(xyz)
1246        sumXYZ += xyz
1247    sumXYZ /= Natoms
1248    XYZ = np.array(XYZ)-sumXYZ
1249    XYZ = np.inner(Amat,XYZ).T
1250    Zmat = np.zeros((3,3))
1251    for i,xyz in enumerate(XYZ):
1252        Zmat += np.outer(xyz.T,xyz)
1253    G2fil.G2Print (' Selected atoms centered at %10.5f %10.5f %10.5f'%(sumXYZ[0],sumXYZ[1],sumXYZ[2]))
1254    Evec,Emat = nl.eig(Zmat)
1255    Evec = np.sqrt(Evec)/(Natoms-3)
1256    Order = np.argsort(Evec)
1257    XYZ = np.inner(XYZ,Emat.T).T
1258    XYZ = np.array([XYZ[Order[2]],XYZ[Order[1]],XYZ[Order[0]]]).T
1259    G2fil.G2Print (' Atoms in Cartesian best plane coordinates:')
1260    G2fil.G2Print (' Name         X         Y         Z')
1261    for i,xyz in enumerate(XYZ):
1262        G2fil.G2Print (' %6s%10.3f%10.3f%10.3f'%(Atoms[i][1].ljust(6),xyz[0],xyz[1],xyz[2]))
1263    G2fil.G2Print ('\n Best plane RMS X =%8.3f, Y =%8.3f, Z =%8.3f'%(Evec[Order[2]],Evec[Order[1]],Evec[Order[0]]))
1264
1265def main():
1266    'Called to run a refinement when this module is executed '
1267    starttime = time.time()
1268    arg = sys.argv
1269    if len(arg) > 1:
1270        GPXfile = arg[1]
1271        if not ospath.exists(GPXfile):
1272            G2fil.G2Print ('ERROR - '+GPXfile+" doesn't exist!")
1273            exit()
1274    else:
1275        G2fil.G2Print ('ERROR - missing filename')
1276        exit()
1277    # TODO: figure out if this is a sequential refinement and call SeqRefine(GPXfile,None)
1278    Refine(GPXfile,None)
1279    G2fil.G2Print("Done. Execution time {:.2f} sec.".format(time.time()-starttime))
1280
1281if __name__ == '__main__':
1282    GSASIIpath.InvokeDebugOpts()
1283    main()
Note: See TracBrowser for help on using the repository browser.