source: trunk/GSASIIstrMain.py @ 5135

Last change on this file since 5135 was 5135, checked in by toby, 9 months ago

show that max(shtf/esd) is from all cycles; fix bug for install path w/underscore

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