source: trunk/GSASIIstrMain.py @ 4894

Last change on this file since 4894 was 4894, checked in by toby, 6 months ago

bug on aborted refinement; more on FPA

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