source: trunk/GSASIIstrMain.py @ 4819

Last change on this file since 4819 was 4819, checked in by toby, 2 years ago

seq ref: fix & edit phase vars for 1 hist; spyder conveniences; improve filter; fix bkg copy flags bug

fail

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