source: trunk/GSASIIstrMain.py @ 5484

Last change on this file since 5484 was 5484, checked in by toby, 8 months ago

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