source: trunk/GSASIIstrMain.py @ 4693

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

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