source: trunk/GSASIIstrMain.py

Last change on this file was 5661, checked in by toby, 2 days ago

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