source: trunk/GSASIIstrMain.py @ 4789

Last change on this file since 4789 was 4789, checked in by toby, 10 months ago

switch to new HessianLSQ; Ouch#4 msg in GUI; misc wx4.1 fixes; docs improvements; switch to new HessianLSQ

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