source: trunk/GSASIIstrMain.py @ 5117

Last change on this file since 5117 was 5117, checked in by vondreele, 9 months ago

complete implementation of EDX data in GSAS-II; will do LeBail? or Pawley refinement
Will not do sample broadening (not even shown at present)
Add column selection for Pawley refinement flag (like similar in Peak List refinement)

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