source: trunk/GSASIIstrMain.py @ 4843

Last change on this file since 4843 was 4843, checked in by toby, 7 months ago

fix Le Bail issues: LB before refinement; bug for combined fits; LeBail? -> Le Bail

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 57.8 KB
Line 
1# -*- coding: utf-8 -*-
2'''
3*GSASIIstrMain: main structure routine*
4---------------------------------------
5
6'''
7########### SVN repository information ###################
8# $Date: 2021-03-07 23:46:32 +0000 (Sun, 07 Mar 2021) $
9# $Author: toby $
10# $Revision: 4843 $
11# $URL: trunk/GSASIIstrMain.py $
12# $Id: GSASIIstrMain.py 4843 2021-03-07 23:46:32Z 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: 4843 $")
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            if phase['Histograms'][h]['LeBail']:
512                 return True
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 = 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,
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
995    Amat,Bmat = G2lat.cell2AB(Cell[:6])
996    covData = {}
997    if len(DisAglData.get('covData',{})):
998        covData = DisAglData['covData']
999        covMatrix = covData['covMatrix']
1000        varyList = covData['varyList']
1001        pfx = str(DisAglData['pId'])+'::'
1002
1003    Factor = DisAglCtls['Factors']
1004    Radii = dict(zip(DisAglCtls['AtomTypes'],zip(DisAglCtls['BondRadii'],DisAglCtls['AngleRadii'])))
1005    indices = (-2,-1,0,1,2)
1006    Units = np.array([[h,k,l] for h in indices for k in indices for l in indices])
1007    origAtoms = DisAglData['OrigAtoms']
1008    targAtoms = DisAglData['TargAtoms']
1009    AtomLabels = {}
1010    for Oatom in origAtoms:
1011        AtomLabels[Oatom[0]] = Oatom[1]
1012    for Oatom in targAtoms:
1013        AtomLabels[Oatom[0]] = Oatom[1]
1014    DistArray = {}
1015    AngArray = {}
1016    for iO,Oatom in enumerate(origAtoms):
1017        DistArray[Oatom[0]] = []
1018        AngArray[Oatom[0]] = []
1019        OxyzNames = ''
1020        IndBlist = []
1021        Dist = []
1022        Vect = []
1023        VectA = []
1024        angles = []
1025        for Tatom in targAtoms:
1026            Xvcov = []
1027            TxyzNames = ''
1028            if len(DisAglData.get('covData',{})):
1029                OxyzNames = [pfx+'dAx:%d'%(Oatom[0]),pfx+'dAy:%d'%(Oatom[0]),pfx+'dAz:%d'%(Oatom[0])]
1030                TxyzNames = [pfx+'dAx:%d'%(Tatom[0]),pfx+'dAy:%d'%(Tatom[0]),pfx+'dAz:%d'%(Tatom[0])]
1031                Xvcov = G2mth.getVCov(OxyzNames+TxyzNames,varyList,covMatrix)
1032            result = G2spc.GenAtom(Tatom[3:6],SGData,False,Move=False)
1033            BsumR = (Radii[Oatom[2]][0]+Radii[Tatom[2]][0])*Factor[0]
1034            AsumR = (Radii[Oatom[2]][1]+Radii[Tatom[2]][1])*Factor[1]
1035            for [Txyz,Top,Tunit,Spn] in result:
1036                Dx = (Txyz-np.array(Oatom[3:6]))+Units
1037                dx = np.inner(Amat,Dx)
1038                dist = ma.masked_less(np.sqrt(np.sum(dx**2,axis=0)),0.5)
1039                IndB = ma.nonzero(ma.masked_greater(dist-BsumR,0.))
1040                if np.any(IndB):
1041                    for indb in IndB:
1042                        for i in range(len(indb)):
1043                            if str(dx.T[indb][i]) not in IndBlist:
1044                                IndBlist.append(str(dx.T[indb][i]))
1045                                unit = Units[indb][i]
1046                                tunit = (unit[0]+Tunit[0],unit[1]+Tunit[1],unit[2]+Tunit[2])
1047                                sig = 0.0
1048                                if len(Xvcov):
1049                                    pdpx = G2mth.getDistDerv(Oatom[3:6],Tatom[3:6],Amat,unit,Top,SGData)
1050                                    sig = np.sqrt(np.inner(pdpx,np.inner(pdpx,Xvcov)))
1051                                Dist.append([Oatom[0],Tatom[0],tunit,Top,ma.getdata(dist[indb])[i],sig])
1052                                if (Dist[-1][-2]-AsumR) <= 0.:
1053                                    Vect.append(dx.T[indb][i]/Dist[-1][-2])
1054                                    VectA.append([OxyzNames,np.array(Oatom[3:6]),TxyzNames,np.array(Tatom[3:6]),unit,Top])
1055                                else:
1056                                    Vect.append([0.,0.,0.])
1057                                    VectA.append([])
1058        if dlg is not None:
1059            dlg.Update(iO,newmsg='Atoms done=%d'%(iO))
1060        for D in Dist:
1061            DistArray[Oatom[0]].append(D[1:])
1062        Vect = np.array(Vect)
1063        angles = np.zeros((len(Vect),len(Vect)))
1064        angsig = np.zeros((len(Vect),len(Vect)))
1065        for i,veca in enumerate(Vect):
1066            if np.any(veca):
1067                for j,vecb in enumerate(Vect):
1068                    if np.any(vecb):
1069                        angles[i][j],angsig[i][j] = G2mth.getAngSig(VectA[i],VectA[j],Amat,SGData,covData)
1070                        if i <= j: continue
1071                        AngArray[Oatom[0]].append((i,j,
1072                            G2mth.getAngSig(VectA[i],VectA[j],Amat,SGData,covData)))
1073    return AtomLabels,DistArray,AngArray
1074
1075def PrintDistAngle(DisAglCtls,DisAglData,out=sys.stdout):
1076    '''Print distances and angles
1077
1078    :param dict DisAglCtls: contains distance/angle radii usually defined using
1079       :func:`GSASIIctrlGUI.DisAglDialog`
1080    :param dict DisAglData: contains phase data:
1081       Items 'OrigAtoms' and 'TargAtoms' contain the atoms to be used
1082       for distance/angle origins and atoms to be used as targets.
1083       Item 'SGData' has the space group information (see :ref:`Space Group object<SGData_table>`)
1084    :param file out: file object for output. Defaults to sys.stdout.
1085    '''
1086    def MyPrint(s):
1087        out.write(s+'\n')
1088        # print(s,file=out) # use in Python 3
1089
1090    def ShowBanner(name):
1091        MyPrint(80*'*')
1092        MyPrint('   Interatomic Distances and Angles for phase '+name)
1093        MyPrint((80*'*')+'\n')
1094
1095    ShowBanner(DisAglCtls['Name'])
1096    SGData = DisAglData['SGData']
1097    SGtext,SGtable = G2spc.SGPrint(SGData)
1098    for line in SGtext: MyPrint(line)
1099    if len(SGtable) > 1:
1100        for i,item in enumerate(SGtable[::2]):
1101            if 2*i+1 == len(SGtable):
1102                line = ' %s'%(item.ljust(30))
1103            else:
1104                line = ' %s %s'%(item.ljust(30),SGtable[2*i+1].ljust(30))
1105            MyPrint(line)
1106    else:
1107        MyPrint(' ( 1)    %s'%(SGtable[0])) #triclinic case
1108    Cell = DisAglData['Cell']
1109
1110    Amat,Bmat = G2lat.cell2AB(Cell[:6])
1111    covData = {}
1112    if len(DisAglData.get('covData',{})):
1113        covData = DisAglData['covData']
1114        pfx = str(DisAglData['pId'])+'::'
1115        A = G2lat.cell2A(Cell[:6])
1116        cellSig = G2stIO.getCellEsd(pfx,SGData,A,covData)
1117        names = [' a = ',' b = ',' c = ',' alpha = ',' beta = ',' gamma = ',' Volume = ']
1118        valEsd = [G2mth.ValEsd(Cell[i],cellSig[i],True) for i in range(7)]
1119        line = '\n Unit cell:'
1120        for name,vals in zip(names,valEsd):
1121            line += name+vals
1122        MyPrint(line)
1123    else:
1124        MyPrint('\n Unit cell: a = '+('%.5f'%Cell[0])+' b = '+('%.5f'%Cell[1])+' c = '+('%.5f'%Cell[2])+
1125            ' alpha = '+('%.3f'%Cell[3])+' beta = '+('%.3f'%Cell[4])+' gamma = '+
1126            ('%.3f'%Cell[5])+' Volume = '+('%.3f'%Cell[6]))
1127
1128    AtomLabels,DistArray,AngArray = RetDistAngle(DisAglCtls,DisAglData)
1129    origAtoms = DisAglData['OrigAtoms']
1130    for Oatom in origAtoms:
1131        i = Oatom[0]
1132        Dist = DistArray[i]
1133        nDist = len(Dist)
1134        angles = np.zeros((nDist,nDist))
1135        angsig = np.zeros((nDist,nDist))
1136        for k,j,tup in AngArray[i]:
1137            angles[k][j],angsig[k][j] = angles[j][k],angsig[j][k] = tup
1138        line = ''
1139        for i,x in enumerate(Oatom[3:6]):
1140            line += ('%12.5f'%x).rstrip('0')
1141        MyPrint('\n Distances & angles for '+Oatom[1]+' at '+line.rstrip())
1142        MyPrint(80*'*')
1143        line = ''
1144        for dist in Dist[:-1]:
1145            line += '%12s'%(AtomLabels[dist[0]].center(12))
1146        MyPrint('  To       cell +(sym. op.)      dist.  '+line.rstrip())
1147        BVS = {}
1148        BVdat = {}
1149        Otyp = G2elem.FixValence(Oatom[2]).split('+')[0].split('-')[0]
1150        BVox = [BV for BV in atmdata.BVSoxid[Otyp] if '+' in BV]
1151        if len(BVox):
1152            BVS = {BV:0.0 for BV in BVox}
1153            BVdat = {BV:dict(zip(['O','F','Cl'],atmdata.BVScoeff[BV])) for BV in BVox}
1154            pvline = 'Bond Valence sums for: '
1155        for i,dist in enumerate(Dist):
1156            line = ''
1157            for j,angle in enumerate(angles[i][0:i]):
1158                sig = angsig[i][j]
1159                if angle:
1160                    if sig:
1161                        line += '%12s'%(G2mth.ValEsd(angle,sig,True).center(12))
1162                    else:
1163                        val = '%.3f'%(angle)
1164                        line += '%12s'%(val.center(12))
1165                else:
1166                    line += 12*' '
1167            if dist[4]:            #sig exists!
1168                val = G2mth.ValEsd(dist[3],dist[4])
1169            else:
1170                val = '%8.4f'%(dist[3])
1171            if len(BVox):
1172                Tatm = G2elem.FixValence(DisAglData['TargAtoms'][dist[0]][2]).split('-')[0]
1173                if Tatm in ['O','F','Cl']:
1174                    for BV in BVox:
1175                        BVS[BV] += np.exp((BVdat[BV][Tatm]-dist[3])/0.37)               
1176            tunit = '[%2d%2d%2d]'% dist[1]
1177            MyPrint((%8s%10s+(%4d) %12s'%(AtomLabels[dist[0]].ljust(8),tunit.ljust(10),dist[2],val.center(12)))+line.rstrip())
1178        if len(BVox):
1179            MyPrint(80*'*')
1180            for BV in BVox:
1181                pvline += ' %s: %.2f  '%(BV,BVS[BV])
1182            MyPrint(pvline)
1183
1184def DisAglTor(DATData):
1185    'Needs a doc string'
1186    SGData = DATData['SGData']
1187    Cell = DATData['Cell']
1188
1189    Amat,Bmat = G2lat.cell2AB(Cell[:6])
1190    covData = {}
1191    pfx = ''
1192    if 'covData' in DATData:
1193        covData = DATData['covData']
1194        pfx = str(DATData['pId'])+'::'
1195    Datoms = []
1196    Oatoms = []
1197    for i,atom in enumerate(DATData['Datoms']):
1198        symop = atom[-1].split('+')
1199        if len(symop) == 1:
1200            symop.append('0,0,0')
1201        symop[0] = int(symop[0])
1202        symop[1] = eval(symop[1])
1203        atom.append(symop)
1204        Datoms.append(atom)
1205        oatom = DATData['Oatoms'][i]
1206        names = ['','','']
1207        if pfx:
1208            names = [pfx+'dAx:'+str(oatom[0]),pfx+'dAy:'+str(oatom[0]),pfx+'dAz:'+str(oatom[0])]
1209        oatom += [names,]
1210        Oatoms.append(oatom)
1211    atmSeq = [atom[1]+'('+atom[-2]+')' for atom in Datoms]
1212    if DATData['Natoms'] == 4:  #torsion
1213        Tors,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData)
1214        G2fil.G2Print (' Torsion angle for %s atom sequence: %s = %s'%(DATData['Name'],str(atmSeq).replace("'","")[1:-1],G2mth.ValEsd(Tors,sig)))
1215        G2fil.G2Print (' NB: Atom sequence determined by selection order')
1216        return      # done with torsion
1217    elif DATData['Natoms'] == 3:  #angle
1218        Ang,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData)
1219        G2fil.G2Print (' Angle in %s for atom sequence: %s = %s'%(DATData['Name'],str(atmSeq).replace("'","")[1:-1],G2mth.ValEsd(Ang,sig)))
1220        G2fil.G2Print (' NB: Atom sequence determined by selection order')
1221    else:   #2 atoms - distance
1222        Dist,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData)
1223        G2fil.G2Print (' Distance in %s for atom sequence: %s = %s'%(DATData['Name'],str(atmSeq).replace("'","")[1:-1],G2mth.ValEsd(Dist,sig)))
1224
1225def BestPlane(PlaneData):
1226    'Needs a doc string'
1227
1228    def ShowBanner(name):
1229        G2fil.G2Print (80*'*')
1230        G2fil.G2Print ('   Best plane result for phase '+name)
1231        G2fil.G2Print (80*'*','\n')
1232
1233    ShowBanner(PlaneData['Name'])
1234
1235    Cell = PlaneData['Cell']
1236    Amat,Bmat = G2lat.cell2AB(Cell[:6])
1237    Atoms = PlaneData['Atoms']
1238    sumXYZ = np.zeros(3)
1239    XYZ = []
1240    Natoms = len(Atoms)
1241    for atom in Atoms:
1242        xyz = np.array(atom[3:6])
1243        XYZ.append(xyz)
1244        sumXYZ += xyz
1245    sumXYZ /= Natoms
1246    XYZ = np.array(XYZ)-sumXYZ
1247    XYZ = np.inner(Amat,XYZ).T
1248    Zmat = np.zeros((3,3))
1249    for i,xyz in enumerate(XYZ):
1250        Zmat += np.outer(xyz.T,xyz)
1251    G2fil.G2Print (' Selected atoms centered at %10.5f %10.5f %10.5f'%(sumXYZ[0],sumXYZ[1],sumXYZ[2]))
1252    Evec,Emat = nl.eig(Zmat)
1253    Evec = np.sqrt(Evec)/(Natoms-3)
1254    Order = np.argsort(Evec)
1255    XYZ = np.inner(XYZ,Emat.T).T
1256    XYZ = np.array([XYZ[Order[2]],XYZ[Order[1]],XYZ[Order[0]]]).T
1257    G2fil.G2Print (' Atoms in Cartesian best plane coordinates:')
1258    G2fil.G2Print (' Name         X         Y         Z')
1259    for i,xyz in enumerate(XYZ):
1260        G2fil.G2Print (' %6s%10.3f%10.3f%10.3f'%(Atoms[i][1].ljust(6),xyz[0],xyz[1],xyz[2]))
1261    G2fil.G2Print ('\n Best plane RMS X =%8.3f, Y =%8.3f, Z =%8.3f'%(Evec[Order[2]],Evec[Order[1]],Evec[Order[0]]))
1262
1263def main():
1264    'Called to run a refinement when this module is executed '
1265    starttime = time.time()
1266    arg = sys.argv
1267    if len(arg) > 1:
1268        GPXfile = arg[1]
1269        if not ospath.exists(GPXfile):
1270            G2fil.G2Print ('ERROR - '+GPXfile+" doesn't exist!")
1271            exit()
1272    else:
1273        G2fil.G2Print ('ERROR - missing filename')
1274        exit()
1275    # TODO: figure out if this is a sequential refinement and call SeqRefine(GPXfile,None)
1276    Refine(GPXfile,None)
1277    G2fil.G2Print("Done. Execution time {:.2f} sec.".format(time.time()-starttime))
1278
1279if __name__ == '__main__':
1280    GSASIIpath.InvokeDebugOpts()
1281    main()
Note: See TracBrowser for help on using the repository browser.