source: trunk/GSASIIstrMain.py @ 5290

Last change on this file since 5290 was 5290, checked in by vondreele, 6 months ago

fix to non 'analytical Hessian'; add 'num cyc' to results[2].

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 59.9 KB
Line 
1# -*- coding: utf-8 -*-
2'''
3*GSASIIstrMain: main structure routine*
4---------------------------------------
5
6'''
7########### SVN repository information ###################
8# $Date: 2022-05-27 13:35:35 +0000 (Fri, 27 May 2022) $
9# $Author: vondreele $
10# $Revision: 5290 $
11# $URL: trunk/GSASIIstrMain.py $
12# $Id: GSASIIstrMain.py 5290 2022-05-27 13:35:35Z 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: 5290 $")
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            result = so.leastsq(G2stMth.errRefine,values,Dfun=G2stMth.dervRefine,full_output=True,
180                ftol=Ftol,col_deriv=True,factor=Factor,
181                args=([Histograms,Phases,restraintDict,rigidbodyDict],parmDict,varyList,calcControls,pawleyLookup,dlg))
182            ncyc = int(result[2]['nfev']/2)
183            result[2]['num cyc'] = ncyc
184            if refPlotUpdate is not None: refPlotUpdate(Histograms)   # update plot after completion
185        elif 'analytic Hessian' in Controls['deriv type']:
186            Lamda = Controls.get('Marquardt',-3)
187            maxCyc = Controls['max cyc']
188            result = G2mth.HessianLSQ(G2stMth.errRefine,values,Hess=G2stMth.HessRefine,ftol=Ftol,xtol=Xtol,maxcyc=maxCyc,Print=ifPrint,lamda=Lamda,
189                args=([Histograms,Phases,restraintDict,rigidbodyDict],parmDict,varyList,calcControls,pawleyLookup,dlg),
190                refPlotUpdate=refPlotUpdate)
191            ncyc = result[2]['num cyc']+1
192            Rvals['lamMax'] = result[2]['lamMax']
193            if 'Ouch#4' in  result[2]:
194                Rvals['Aborted'] = True
195            if 'msg' in result[2]:
196                Rvals['msg'] = result[2]['msg']
197            Controls['Marquardt'] = -3  #reset to default
198            if 'chisq0' in result[2] and chisq0 is None:
199                chisq0 = result[2]['chisq0']
200        elif 'Hessian SVD' in Controls['deriv type']:
201            maxCyc = Controls['max cyc']
202            result = G2mth.HessianSVD(G2stMth.errRefine,values,Hess=G2stMth.HessRefine,ftol=Ftol,xtol=Xtol,maxcyc=maxCyc,Print=ifPrint,
203                args=([Histograms,Phases,restraintDict,rigidbodyDict],parmDict,varyList,calcControls,pawleyLookup,dlg),
204                refPlotUpdate=refPlotUpdate)
205            if result[1] is None:
206                IfOK = False
207                covMatrix = []
208                sig = len(varyList)*[None,]
209                break
210            ncyc = result[2]['num cyc']+1
211            if 'chisq0' in result[2] and chisq0 is None:
212                chisq0 = result[2]['chisq0']
213        else:           #'numeric'
214            result = so.leastsq(G2stMth.errRefine,values,full_output=True,ftol=Ftol,epsfcn=1.e-8,factor=Factor,
215                args=([Histograms,Phases,restraintDict,rigidbodyDict],parmDict,varyList,calcControls,pawleyLookup,dlg))
216            ncyc = 1
217            result[2]['num cyc'] = ncyc
218            if len(varyList):
219                ncyc = int(result[2]['nfev']/len(varyList))
220            if refPlotUpdate is not None: refPlotUpdate(Histograms)   # update plot
221        #table = dict(zip(varyList,zip(values,result[0],(result[0]-values))))
222        #for item in table: print(item,table[item])               #useful debug - are things shifting?
223        runtime = time.time()-begin
224        Rvals['SVD0'] = result[2].get('SVD0',0)
225        Rvals['converged'] = result[2].get('Converged')
226        Rvals['DelChi2'] = result[2].get('DelChi2',-1.)
227        Rvals['chisq'] = np.sum(result[2]['fvec']**2)
228        G2stMth.Values2Dict(parmDict, varyList, result[0])
229        G2mv.Dict2Map(parmDict)
230        Rvals['Nobs'] = Histograms['Nobs']
231        Rvals['Nvars'] = len(varyList)
232        Rvals['Rwp'] = np.sqrt(Rvals['chisq']/Histograms['sumwYo'])*100.      #to %
233        Rvals['GOF'] = np.sqrt(Rvals['chisq']/(Histograms['Nobs']-len(varyList)))
234        printFile.write(' Number of function calls: %d No. of observations: %d No. of parameters: %d User rejected: %d Sp. gp. extinct: %d\n'%  \
235            (result[2]['nfev'],Histograms['Nobs'],len(varyList),Histograms['Nrej'],Histograms['Next']))
236        if ncyc:
237            printFile.write(' Refinement time = %8.3fs, %8.3fs/cycle, for %d cycles\n'%(runtime,runtime/ncyc,ncyc))
238        printFile.write(' wR = %7.2f%%, chi**2 = %12.6g, GOF = %6.2f\n'%(Rvals['Rwp'],Rvals['chisq'],Rvals['GOF']))
239        sig = len(varyList)*[None,]
240        if 'None' in str(type(result[1])) and ifSeq:    #this bails out of a sequential refinement on singular matrix
241            IfOK = False
242            covMatrix = []
243            G2fil.G2Print ('Warning: **** Refinement failed - singular matrix ****')
244            if 'Hessian' in Controls['deriv type']:
245                num = len(varyList)-1
246                # BHT -- I am not sure if this works correctly:
247                for i,val in enumerate(np.flipud(result[2]['psing'])):
248                    if val:
249                        G2fil.G2Print('Bad parameter: '+varyList[num-i],mode='warn')
250            else:
251                Ipvt = result[2]['ipvt']
252                for i,ipvt in enumerate(Ipvt):
253                    if not np.sum(result[2]['fjac'],axis=1)[i]:
254                        G2fil.G2Print('Bad parameter: '+varyList[ipvt-1],mode='warn')
255            break
256        IfOK = True
257        if not len(varyList) or not result[2]['num cyc']:
258            covMatrix = []
259            break
260        try:
261            covMatrix = result[1]*Rvals['GOF']**2
262            sig = np.sqrt(np.diag(covMatrix))
263            Lastshft = result[0]-values     #NOT last shift since values is starting set before current refinement
264            #table = dict(zip(varyList,zip(values,result[0],Lastshft,Lastshft/sig)))
265            #for item in table: print(item,table[item])               #useful debug
266            Rvals['Max shft/sig'] = np.max(np.nan_to_num(Lastshft/sig))
267            if np.any(np.isnan(sig)) or not sig.shape:
268                G2fil.G2Print ('*** Least squares aborted - some invalid esds possible ***',mode='error')
269            else:
270                print('Maximum shift/esd = {:.3f} for all cycles'
271                          .format(Rvals['Max shft/sig']))
272            # report on refinement issues. Result in Rvals['msg']
273            ReportProblems(result,Rvals,varyList)
274            break                   #refinement succeeded - finish up!
275        except TypeError:
276            # if we get here, no result[1] (covar matrix) was returned or other calc error: refinement failed
277            IfOK = False
278            if 'Hessian' in Controls['deriv type']:
279                SVD0 = result[2].get('SVD0')
280                if SVD0 == -1:
281                    G2fil.G2Print ('**** Refinement failed - singular matrix ****',mode='error')
282                elif SVD0 == -2:
283                    G2fil.G2Print ('**** Refinement failed - other problem ****',mode='error')
284                elif SVD0 > 0:
285                    G2fil.G2Print ('**** Refinement failed with {} SVD singularities ****'.format(SVD0),mode='error')
286                else:
287                    G2fil.G2Print ('**** Refinement failed ****',mode='error')
288                if result[1] is None:
289                    IfOK = False
290                    covMatrix = []
291                    sig = len(varyList)*[None,]
292                # report on highly correlated variables
293                ReportProblems(result,Rvals,varyList)
294                # process singular variables
295                if dlg: break # refining interactively
296            else:
297                G2fil.G2Print ('**** Refinement failed - singular matrix ****',mode='error')
298                Ipvt = result[2]['ipvt']
299                for i,ipvt in enumerate(Ipvt):
300                    if not np.sum(result[2]['fjac'],axis=1)[i]:
301                        G2fil.G2Print ('Removing parameter: '+varyList[ipvt-1])
302                        del(varyList[ipvt-1])
303                        break
304    if IfOK:
305        if CheckLeBail(Phases):   # only needed for LeBail extraction
306            G2stMth.errRefine([],[Histograms,Phases,restraintDict,rigidbodyDict],
307                parmDict,[],calcControls,pawleyLookup,dlg)
308        G2stMth.GetFobsSq(Histograms,Phases,parmDict,calcControls)
309    if chisq0 is not None:
310        Rvals['GOF0'] = np.sqrt(chisq0/(Histograms['Nobs']-len(varyList)))
311    return IfOK,Rvals,result,covMatrix,sig,Lastshft
312
313def Refine(GPXfile,dlg=None,makeBack=True,refPlotUpdate=None,newLeBail=False):
314    '''Global refinement -- refines to minimize against all histograms.
315    This can be called in one of three ways, from :meth:`GSASIIdataGUI.GSASII.OnRefine` in an
316    interactive refinement, where dlg will be a wx.ProgressDialog, or non-interactively from
317    :meth:`GSASIIscriptable.G2Project.refine` or from :func:`do_refine`, where dlg will be None.
318    '''
319    import GSASIImpsubs as G2mp
320    G2mp.InitMP()
321    import pytexture as ptx
322    ptx.pyqlmninit()            #initialize fortran arrays for spherical harmonics
323
324    printFile = open(ospath.splitext(GPXfile)[0]+'.lst','w')
325    G2stIO.ShowBanner(printFile)
326    varyList = []
327    parmDict = {}
328    G2mv.InitVars()
329    Controls = G2stIO.GetControls(GPXfile)
330    Controls['newLeBail'] = newLeBail
331    G2stIO.ShowControls(Controls,printFile)
332    calcControls = {}
333    calcControls.update(Controls)
334    constrDict,fixedList = G2stIO.ReadConstraints(GPXfile)
335    restraintDict = G2stIO.GetRestraints(GPXfile)
336    Histograms,Phases = G2stIO.GetUsedHistogramsAndPhases(GPXfile)
337    if not Phases:
338        G2fil.G2Print (' *** ERROR - you have no phases to refine! ***')
339        G2fil.G2Print (' *** Refine aborted ***')
340        return False,{'msg':'No phases'}
341    if not Histograms:
342        G2fil.G2Print (' *** ERROR - you have no data to refine with! ***')
343        G2fil.G2Print (' *** Refine aborted ***')
344        return False,{'msg':'No data'}
345    rigidbodyDict = G2stIO.GetRigidBodies(GPXfile)
346    rbIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})
347    rbVary,rbDict = G2stIO.GetRigidBodyModels(rigidbodyDict,pFile=printFile)
348    (Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables,MFtables,
349         maxSSwave) = G2stIO.GetPhaseData(Phases,restraintDict,rbIds,pFile=printFile)
350    calcControls['atomIndx'] = atomIndx
351    calcControls['Natoms'] = Natoms
352    calcControls['FFtables'] = FFtables
353    calcControls['BLtables'] = BLtables
354    calcControls['MFtables'] = MFtables
355    calcControls['maxSSwave'] = maxSSwave
356    hapVary,hapDict,controlDict = G2stIO.GetHistogramPhaseData(Phases,Histograms,Controls=calcControls,pFile=printFile)
357    TwConstr,TwFixed = G2stIO.makeTwinFrConstr(Phases,Histograms,hapVary)
358    constrDict += TwConstr
359    fixedList += TwFixed
360    calcControls.update(controlDict)
361    histVary,histDict,controlDict = G2stIO.GetHistogramData(Histograms,pFile=printFile)
362    calcControls.update(controlDict)
363    varyList = rbVary+phaseVary+hapVary+histVary
364    parmDict.update(rbDict)
365    parmDict.update(phaseDict)
366    parmDict.update(hapDict)
367    parmDict.update(histDict)
368    G2stIO.GetFprime(calcControls,Histograms)
369    # do constraint processing
370    varyListStart = tuple(varyList) # save the original varyList before dependent vars are removed
371    msg = G2mv.EvaluateMultipliers(constrDict,parmDict)
372    if msg:
373        return False,{'msg':'Unable to interpret multiplier(s): '+msg}
374    try:
375        errmsg,warnmsg,groups,parmlist = G2mv.GenerateConstraints(varyList,constrDict,fixedList,parmDict)
376        G2mv.normParms(parmDict)
377        G2mv.Map2Dict(parmDict,varyList)   # changes varyList
378    except G2mv.ConstraintException:
379        G2fil.G2Print (' *** ERROR - your constraints are internally inconsistent ***')
380        return False,{'msg':' Constraint error'}
381
382    # remove frozen vars from refinement
383    if 'parmFrozen' not in Controls:
384        Controls['parmFrozen'] = {}
385    if 'FrozenList' not in Controls['parmFrozen']: 
386        Controls['parmFrozen']['FrozenList'] = []
387    parmFrozenList = Controls['parmFrozen']['FrozenList']
388    frozenList = [i for i in varyList if i in parmFrozenList]
389    if len(frozenList) != 0: 
390        varyList = [i for i in varyList if i not in parmFrozenList]
391        G2fil.G2Print(
392            'Frozen refined variables (due to exceeding limits)\n\t:{}'
393            .format(frozenList))
394       
395    ifSeq = False
396    printFile.write('\n Refinement results:\n')
397    printFile.write(135*'-'+'\n')
398    Rvals = {}
399    G2mv.Dict2Map(parmDict)  # impose constraints initially
400   
401    try:
402        covData = {}
403        IfOK,Rvals,result,covMatrix,sig,Lastshft = RefineCore(Controls,Histograms,Phases,restraintDict,
404            rigidbodyDict,parmDict,varyList,calcControls,pawleyLookup,ifSeq,printFile,dlg,
405            refPlotUpdate=refPlotUpdate)
406        if IfOK:
407            if len(covMatrix):      #empty for zero cycle refinement
408                sigDict = dict(zip(varyList,sig))
409                newCellDict = G2stMth.GetNewCellParms(parmDict,varyList)
410                newAtomDict = G2stMth.ApplyXYZshifts(parmDict,varyList)
411                covData = {'variables':result[0],'varyList':varyList,'sig':sig,'Rvals':Rvals,
412                           'varyListStart':varyListStart,'Lastshft':Lastshft,
413                           'covMatrix':covMatrix,'title':GPXfile,'newAtomDict':newAtomDict,
414                           'newCellDict':newCellDict,'freshCOV':True}
415                # add indirectly computed uncertainties into the esd dict
416                sigDict.update(G2mv.ComputeDepESD(covMatrix,varyList))
417                G2stIO.PrintIndependentVars(parmDict,varyList,sigDict,pFile=printFile)
418                G2stMth.ApplyRBModels(parmDict,Phases,rigidbodyDict,True)
419                G2stIO.SetRigidBodyModels(parmDict,sigDict,rigidbodyDict,printFile)
420                G2stIO.SetPhaseData(parmDict,sigDict,Phases,rbIds,covData,restraintDict,printFile)
421                G2stIO.SetISOmodes(parmDict,sigDict,Phases,printFile)
422                G2stIO.SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms,calcControls,
423                                             pFile=printFile,covMatrix=covMatrix,varyList=varyList)
424                G2stIO.SetHistogramData(parmDict,sigDict,Histograms,calcControls,pFile=printFile)
425                # check for variables outside their allowed range, reset and freeze them
426                frozen = dropOOBvars(varyList,parmDict,sigDict,Controls,parmFrozenList)
427                # covData['depSig'] = G2stIO.PhFrExtPOSig  # created in G2stIO.SetHistogramData, no longer used?
428                covData['depSigDict'] = {i:(parmDict[i],sigDict[i]) for i in parmDict if i in sigDict}
429                if len(frozen):
430                    if 'msg' in Rvals:
431                        Rvals['msg'] += '\n'
432                    else:
433                        Rvals['msg'] = ''
434                    msg = ('Warning: {} variable(s) refined outside limits and were frozen ({} total frozen)'
435                        .format(len(frozen),len(parmFrozenList))
436                        )
437                    G2fil.G2Print(msg)
438                    Rvals['msg'] += msg
439                elif len(parmFrozenList):
440                    if 'msg' in Rvals:
441                        Rvals['msg'] += '\n'
442                    else:
443                        Rvals['msg'] = ''
444                    msg = ('Note: a total of {} variable(s) are frozen due to refining outside limits'
445                        .format(len(parmFrozenList))
446                        )
447                    G2fil.G2Print('Note: ',msg)
448                    Rvals['msg'] += msg
449            G2stIO.SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases,rigidbodyDict,covData,parmFrozenList,makeBack)
450            printFile.close()
451            G2fil.G2Print (' Refinement results are in file: '+ospath.splitext(GPXfile)[0]+'.lst')
452            G2fil.G2Print (' ***** Refinement successful *****')
453        else:
454            G2fil.G2Print ('****ERROR - Refinement failed',mode='error')
455            if 'msg' in Rvals:
456                G2fil.G2Print ('Note refinement problem:',mode='warn')
457                G2fil.G2Print (Rvals['msg'],mode='warn')
458            raise G2obj.G2Exception('**** ERROR: Refinement failed ****')
459    except G2obj.G2RefineCancel as Msg:
460        printFile.close()
461        G2fil.G2Print (' ***** Refinement stopped *****')
462        if not hasattr(Msg,'msg'): Msg.msg = str(Msg)
463        if 'msg' in Rvals:
464            Rvals['msg'] += '\n'
465            Rvals['msg'] += Msg.msg
466            if not dlg:
467                G2fil.G2Print ('Note refinement problem:',mode='warn')
468                G2fil.G2Print (Rvals['msg'],mode='warn')
469        else:
470            Rvals['msg'] = Msg.msg
471        return False,Rvals
472#    except G2obj.G2Exception as Msg:  # cell metric error, others?
473    except Exception as Msg:  # cell metric error, others?
474        if GSASIIpath.GetConfigValue('debug'):
475            import traceback
476            print(traceback.format_exc())       
477        if not hasattr(Msg,'msg'): Msg.msg = str(Msg)
478        printFile.close()
479        G2fil.G2Print (' ***** Refinement error *****')
480        if 'msg' in Rvals:
481            Rvals['msg'] += '\n\n'
482            Rvals['msg'] += Msg.msg
483            if not dlg:
484                G2fil.G2Print ('Note refinement problem:',mode='warn')
485                G2fil.G2Print (Rvals['msg'],mode='warn')
486        else:
487            Rvals['msg'] = Msg.msg
488        return False,Rvals
489   
490#for testing purposes, create a file for testderiv
491    if GSASIIpath.GetConfigValue('debug'):   # and IfOK:
492#needs: values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup
493        fl = open(ospath.splitext(GPXfile)[0]+'.testDeriv','wb')
494        cPickle.dump(result[0],fl,1)
495        cPickle.dump([Histograms,Phases,restraintDict,rigidbodyDict],fl,1)
496        cPickle.dump([constrDict,fixedList,G2mv.GetDependentVars()],fl,1)
497        cPickle.dump(parmDict,fl,1)
498        cPickle.dump(varyList,fl,1)
499        cPickle.dump(calcControls,fl,1)
500        cPickle.dump(pawleyLookup,fl,1)
501        fl.close()
502    if dlg:
503        return True,Rvals
504    elif 'msg' in Rvals:
505        G2fil.G2Print ('Reported from refinement:',mode='warn')
506        G2fil.G2Print (Rvals['msg'],mode='warn')
507
508def CheckLeBail(Phases):
509    '''Check if there is a LeBail extraction in any histogram
510
511    :returns: True if there is at least one LeBail flag turned on, False otherwise
512    '''
513    for key in Phases:
514        phase = Phases[key]
515        for h in phase['Histograms']:
516            #phase['Histograms'][h]
517            if not phase['Histograms'][h]['Use']: continue
518            try:
519                if phase['Histograms'][h]['LeBail']:
520                     return True
521            except KeyError:    #HKLF & old gpx files
522                pass
523    return False
524       
525def DoLeBail(GPXfile,dlg=None,cycles=10,refPlotUpdate=None):
526    '''Fit LeBail intensities without changes to any other refined parameters.
527    This is a stripped-down version of :func:`Refine` that does not perform
528    any refinement cycles
529    '''
530    import GSASIImpsubs as G2mp
531    G2mp.InitMP()
532    import pytexture as ptx
533    ptx.pyqlmninit()            #initialize fortran arrays for spherical harmonics
534
535    #varyList = []
536    parmDict = {}
537    Controls = G2stIO.GetControls(GPXfile)
538    calcControls = {}
539    calcControls.update(Controls)
540    constrDict,fixedList = G2stIO.ReadConstraints(GPXfile)
541    restraintDict = {}
542    Histograms,Phases = G2stIO.GetUsedHistogramsAndPhases(GPXfile)
543    if not Phases:
544        G2fil.G2Print (' *** ERROR - you have no phases to refine! ***')
545        return False,{'msg':'No phases'}
546    if not Histograms:
547        G2fil.G2Print (' *** ERROR - you have no data to refine with! ***')
548        return False,{'msg':'No data'}
549    if not CheckLeBail(Phases):
550        msg = 'Warning: There are no histograms with LeBail extraction enabled'
551        G2fil.G2Print ('*** '+msg+' ***')
552        return False,{'msg': msg}
553    rigidbodyDict = G2stIO.GetRigidBodies(GPXfile)
554    rbIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})
555    rbVary,rbDict = G2stIO.GetRigidBodyModels(rigidbodyDict,Print=False)
556    (Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables,MFtables,
557         maxSSwave) = G2stIO.GetPhaseData(Phases,restraintDict,rbIds,Print=False)
558    calcControls['atomIndx'] = atomIndx
559    calcControls['Natoms'] = Natoms
560    calcControls['FFtables'] = FFtables
561    calcControls['BLtables'] = BLtables
562    calcControls['MFtables'] = MFtables
563    calcControls['maxSSwave'] = maxSSwave
564    hapVary,hapDict,controlDict = G2stIO.GetHistogramPhaseData(Phases,Histograms,Controls=calcControls,Print=False)
565    calcControls.update(controlDict)
566    histVary,histDict,controlDict = G2stIO.GetHistogramData(Histograms,Print=False)
567    calcControls.update(controlDict)
568    parmDict.update(rbDict)
569    parmDict.update(phaseDict)
570    parmDict.update(hapDict)
571    parmDict.update(histDict)
572    G2stIO.GetFprime(calcControls,Histograms)
573    try:
574        for i in range(cycles):
575            M = G2stMth.errRefine([],[Histograms,Phases,restraintDict,rigidbodyDict],parmDict,[],calcControls,pawleyLookup,dlg)
576            G2stMth.GetFobsSq(Histograms,Phases,parmDict,calcControls)
577            if refPlotUpdate is not None: refPlotUpdate(Histograms,i)
578        Rvals = {}
579        Rvals['chisq'] = np.sum(M**2)
580        Rvals['Nobs'] = Histograms['Nobs']
581        Rvals['Rwp'] = np.sqrt(Rvals['chisq']/Histograms['sumwYo'])*100.      #to %
582        Rvals['GOF'] = np.sqrt(Rvals['chisq']/(Histograms['Nobs'])) # no variables
583
584        covData = {'variables':0,'varyList':[],'sig':[],'Rvals':Rvals,'varyListStart':[],
585            'covMatrix':np.zeros([0,0]),'title':GPXfile,'freshCOV':True}
586          # ??  'newAtomDict':newAtomDict,'newCellDict':newCellDict,
587       
588        G2stIO.SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases,rigidbodyDict,covData,[],True)
589        G2fil.G2Print (' ***** LeBail fit completed *****')
590        return True,Rvals
591    except Exception as Msg:
592        G2fil.G2Print (' ***** LeBail fit error *****')
593        if not hasattr(Msg,'msg'): Msg.msg = str(Msg)
594        if GSASIIpath.GetConfigValue('debug'):
595            import traceback
596            print(traceback.format_exc())       
597        return False,{'msg':Msg.msg}
598
599def phaseCheck(phaseVary,Phases,histogram):
600    '''
601    Removes unused parameters from phase varylist if phase not in histogram
602    for seq refinement removes vars in "Fix FXU" and "FixedSeqVars" here
603    '''
604    NewVary = []
605    for phase in Phases:
606        if histogram not in Phases[phase]['Histograms']: continue
607        if Phases[phase]['Histograms'][histogram]['Use']:
608            pId = Phases[phase]['pId']
609            newVary = [item for item in phaseVary if item.split(':')[0] == str(pId)]
610            FixVals = Phases[phase]['Histograms'][histogram].get('Fix FXU',' ')
611            if 'F' in FixVals:
612                newVary = [item for item in newVary if not 'Afrac' in item]
613            if 'X' in FixVals:
614                newVary = [item for item in newVary if not 'dA' in item]
615            if 'U' in FixVals:
616                newVary = [item for item in newVary if not 'AU' in item]
617            if 'M' in FixVals:
618                newVary = [item for item in newVary if not 'AM' in item]
619            removeVars = Phases[phase]['Histograms'][histogram].get('FixedSeqVars',[])
620            newVary = [item for item in newVary if item not in removeVars]
621            NewVary += newVary
622    return NewVary
623
624def SeqRefine(GPXfile,dlg,refPlotUpdate=None):
625    '''Perform a sequential refinement -- cycles through all selected histgrams,
626    one at a time
627    '''
628    import GSASIImpsubs as G2mp
629    G2mp.InitMP()
630    import pytexture as ptx
631    ptx.pyqlmninit()            #initialize fortran arrays for spherical harmonics
632    msgs = {}
633    printFile = open(ospath.splitext(GPXfile)[0]+'.lst','w')
634    G2fil.G2Print ('Starting Sequential Refinement')
635    G2stIO.ShowBanner(printFile)
636    Controls = G2stIO.GetControls(GPXfile)
637    preFrozenCount = 0
638    for h in Controls['parmFrozen']:
639        if h == 'FrozenList':
640            continue
641        preFrozenCount += len(Controls['parmFrozen'][h])   
642    G2stIO.ShowControls(Controls,printFile,SeqRef=True,preFrozenCount=preFrozenCount)
643    restraintDict = G2stIO.GetRestraints(GPXfile)
644    Histograms,Phases = G2stIO.GetUsedHistogramsAndPhases(GPXfile)
645    if not Phases:
646        G2fil.G2Print (' *** ERROR - you have no phases to refine! ***')
647        G2fil.G2Print (' *** Refine aborted ***')
648        return False,'No phases'
649    if not Histograms:
650        G2fil.G2Print (' *** ERROR - you have no data to refine with! ***')
651        G2fil.G2Print (' *** Refine aborted ***')
652        return False,'No data'
653    rigidbodyDict = G2stIO.GetRigidBodies(GPXfile)
654    rbIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})
655    rbVary,rbDict = G2stIO.GetRigidBodyModels(rigidbodyDict,pFile=printFile)
656    G2mv.InitVars()
657    (Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,
658     FFtables,BLtables,MFtables,maxSSwave) = G2stIO.GetPhaseData(
659             Phases,restraintDict,rbIds,Print=False,pFile=printFile,
660             seqHistName='All')
661    for item in phaseVary:
662        if '::A0' in item:
663            G2fil.G2Print ('**** WARNING - lattice parameters should not be refined in a sequential refinement ****')
664            G2fil.G2Print ('****           instead use the Dij parameters for each powder histogram            ****')
665            return False,'Lattice parameter refinement error - see console message'
666        if '::C(' in item:
667            G2fil.G2Print ('**** WARNING - phase texture parameters should not be refined in a sequential refinement ****')
668            G2fil.G2Print ('****           instead use the C(L,N) parameters for each powder histogram               ****')
669            return False,'Phase texture refinement error - see console message'
670    if 'Seq Data' in Controls:
671        histNames = Controls['Seq Data']
672    else: # patch from before Controls['Seq Data'] was implemented?
673        histNames = G2stIO.GetHistogramNames(GPXfile,['PWDR',])
674    if Controls.get('Reverse Seq'):
675        histNames.reverse()
676    SeqResult = G2stIO.GetSeqResult(GPXfile)
677#    SeqResult = {'SeqPseudoVars':{},'SeqParFitEqList':[]}
678    Histo = {}
679    NewparmDict = {}
680    G2stIO.SetupSeqSavePhases(GPXfile)
681    msgs['steepestNum'] = 0
682    msgs['maxshift/sigma'] = []
683    lasthist = ''
684    for ihst,histogram in enumerate(histNames):
685        if GSASIIpath.GetConfigValue('Show_timing'): t1 = time.time()
686        G2fil.G2Print('\nRefining with '+str(histogram))
687        G2mv.InitVars()
688        (Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,
689         FFtables,BLtables,MFtables,maxSSwave) = G2stIO.GetPhaseData(
690             Phases,restraintDict,rbIds,Print=False,pFile=printFile,
691             seqHistName=histogram)
692        ifPrint = False
693        if dlg:
694            dlg.SetTitle('Residual for histogram '+str(ihst))
695        calcControls = {}
696        calcControls['atomIndx'] = atomIndx
697        calcControls['Natoms'] = Natoms
698        calcControls['FFtables'] = FFtables
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        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.