source: trunk/GSASIIstrMain.py @ 5310

Last change on this file since 5310 was 5310, checked in by vondreele, 15 months ago

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