source: trunk/GSASIIstrMain.py @ 4875

Last change on this file since 4875 was 4875, checked in by vondreele, 3 years ago

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