source: trunk/GSASIIstrMain.py @ 5287

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

improvements to zero cycle refinement & partial calcs.

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