source: trunk/GSASIIstrMain.py @ 4850

Last change on this file since 4850 was 4850, checked in by vondreele, 8 months ago

fix LeBail? test to skip HKLF histograms in G2strMain
fix DBW access for old gpx files in G2dataGUI
fix Quaternion access for old gpx files in G2plot
add new testSytSym.py routine

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