source: trunk/GSASIIstrMain.py @ 4189

Last change on this file since 4189 was 4189, checked in by vondreele, 2 years ago

revert back on progress bar change

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 40.3 KB
Line 
1# -*- coding: utf-8 -*-
2'''
3*GSASIIstrMain: main structure routine*
4---------------------------------------
5
6'''
7########### SVN repository information ###################
8# $Date: 2019-11-01 10:08:20 +0000 (Fri, 01 Nov 2019) $
9# $Author: vondreele $
10# $Revision: 4189 $
11# $URL: trunk/GSASIIstrMain.py $
12# $Id: GSASIIstrMain.py 4189 2019-11-01 10:08:20Z 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: 4189 $")
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
39
40sind = lambda x: np.sin(x*np.pi/180.)
41cosd = lambda x: np.cos(x*np.pi/180.)
42tand = lambda x: np.tan(x*np.pi/180.)
43asind = lambda x: 180.*np.arcsin(x)/np.pi
44acosd = lambda x: 180.*np.arccos(x)/np.pi
45atan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
46
47ateln2 = 8.0*math.log(2.0)
48DEBUG = True
49
50def RefineCore(Controls,Histograms,Phases,restraintDict,rigidbodyDict,parmDict,varyList,
51    calcControls,pawleyLookup,ifSeq,printFile,dlg,refPlotUpdate=None):
52    '''Core optimization routines, shared between SeqRefine and Refine
53
54    :returns: 5-tuple of ifOk (bool), Rvals (dict), result, covMatrix, sig
55    '''
56#    print 'current',varyList
57#    for item in parmDict: print item,parmDict[item] ######### show dict just before refinement
58    G2mv.Map2Dict(parmDict,varyList)
59    ifPrint = True
60    if ifSeq:
61        ifPrint = False
62    Rvals = {}
63    while True:
64        begin = time.time()
65        values =  np.array(G2stMth.Dict2Values(parmDict, varyList))
66        if np.any(np.isnan(values)):
67            raise G2obj.G2Exception('ERROR - nan found in LS parameters - use Calculate/View LS parms to locate')
68        # test code to compute GOF and save for external repeat
69        #args = ([Histograms,Phases,restraintDict,rigidbodyDict],parmDict,varyList,calcControls,pawleyLookup,dlg)
70        #print '*** before fit chi**2',np.sum(G2stMth.errRefine(values,*args)**2)
71        #fl = open('beforeFit.cpickle','wb')
72        #cPickle.dump(values,fl,1)
73        #cPickle.dump(args[:-1],fl,1)
74        #fl.close()
75        Ftol = Controls['min dM/M']
76        Xtol = Controls['SVDtol']
77        Factor = Controls['shift factor']
78        if 'Jacobian' in Controls['deriv type']:
79            result = so.leastsq(G2stMth.errRefine,values,Dfun=G2stMth.dervRefine,full_output=True,
80                ftol=Ftol,col_deriv=True,factor=Factor,
81                args=([Histograms,Phases,restraintDict,rigidbodyDict],parmDict,varyList,calcControls,pawleyLookup,dlg))
82            ncyc = int(result[2]['nfev']/2)
83            if refPlotUpdate is not None: refPlotUpdate(Histograms)   # update plot after completion
84        elif 'analytic Hessian' in Controls['deriv type']:
85            Lamda = Controls.get('Marquardt',-3)
86            maxCyc = Controls['max cyc']
87            result = G2mth.HessianLSQ(G2stMth.errRefine,values,Hess=G2stMth.HessRefine,ftol=Ftol,xtol=Xtol,maxcyc=maxCyc,Print=ifPrint,lamda=Lamda,
88                args=([Histograms,Phases,restraintDict,rigidbodyDict],parmDict,varyList,calcControls,pawleyLookup,dlg),
89                refPlotUpdate=refPlotUpdate)
90            ncyc = result[2]['num cyc']+1
91            Rvals['lamMax'] = result[2]['lamMax']
92            Controls['Marquardt'] = -3  #reset to default
93        elif 'Hessian SVD' in Controls['deriv type']:
94            maxCyc = Controls['max cyc']
95            result = G2mth.HessianSVD(G2stMth.errRefine,values,Hess=G2stMth.HessRefine,ftol=Ftol,xtol=Xtol,maxcyc=maxCyc,Print=ifPrint,
96                args=([Histograms,Phases,restraintDict,rigidbodyDict],parmDict,varyList,calcControls,pawleyLookup,dlg),
97                refPlotUpdate=refPlotUpdate)
98            if result[1] is None:
99                IfOK = False
100                covMatrix = []
101                sig = len(varyList)*[None,]
102                break
103            ncyc = result[2]['num cyc']+1
104        else:           #'numeric'
105            result = so.leastsq(G2stMth.errRefine,values,full_output=True,ftol=Ftol,epsfcn=1.e-8,factor=Factor,
106                args=([Histograms,Phases,restraintDict,rigidbodyDict],parmDict,varyList,calcControls,pawleyLookup,dlg))
107            ncyc = 1
108            if len(varyList):
109                ncyc = int(result[2]['nfev']/len(varyList))
110            if refPlotUpdate is not None: refPlotUpdate(Histograms)   # update plot
111#        table = dict(zip(varyList,zip(values,result[0],(result[0]-values))))
112#        for item in table: print item,table[item]               #useful debug - are things shifting?
113        runtime = time.time()-begin
114        Rvals['SVD0'] = result[2].get('SVD0',0)
115        Rvals['converged'] = result[2].get('Converged')
116        Rvals['DelChi2'] = result[2].get('DelChi2',-1.)
117        Rvals['chisq'] = np.sum(result[2]['fvec']**2)
118        G2stMth.Values2Dict(parmDict, varyList, result[0])
119        G2mv.Dict2Map(parmDict,varyList)
120        Rvals['Nobs'] = Histograms['Nobs']
121        Rvals['Rwp'] = np.sqrt(Rvals['chisq']/Histograms['sumwYo'])*100.      #to %
122        Rvals['GOF'] = np.sqrt(Rvals['chisq']/(Histograms['Nobs']-len(varyList)))
123        printFile.write(' Number of function calls: %d No. of observations: %d No. of parameters: %d User rejected: %d Sp. gp. extinct: %d\n'%  \
124            (result[2]['nfev'],Histograms['Nobs'],len(varyList),Histograms['Nrej'],Histograms['Next']))
125        if ncyc:
126            printFile.write(' Refinement time = %8.3fs, %8.3fs/cycle, for %d cycles\n'%(runtime,runtime/ncyc,ncyc))
127        printFile.write(' wR = %7.2f%%, chi**2 = %12.6g, GOF = %6.2f\n'%(Rvals['Rwp'],Rvals['chisq'],Rvals['GOF']))
128        sig = len(varyList)*[None,]
129        if 'None' in str(type(result[1])) and ifSeq:    #this bails out of a sequential refinement on singular matrix
130            IfOK = False
131            covMatrix = []
132            G2fil.G2Print ('Warning: **** Refinement failed - singular matrix ****')
133            if 'Hessian' in Controls['deriv type']:
134                num = len(varyList)-1
135                for i,val in enumerate(np.flipud(result[2]['psing'])):
136                    if val:
137                        G2fil.G2Print('Bad parameter: '+varyList[num-i],mode='warn')
138            else:
139                Ipvt = result[2]['ipvt']
140                for i,ipvt in enumerate(Ipvt):
141                    if not np.sum(result[2]['fjac'],axis=1)[i]:
142                        G2fil.G2Print('Bad parameter: '+varyList[ipvt-1],mode='warn')
143            break
144        IfOK = True
145        try:
146            covMatrix = result[1]*Rvals['GOF']**2
147            sig = np.sqrt(np.diag(covMatrix))
148            if np.any(np.isnan(sig)) or not sig.shape:
149                G2fil.G2Print ('*** Least squares aborted - some invalid esds possible ***',mode='error')
150#            table = dict(zip(varyList,zip(values,result[0],(result[0]-values)/sig)))
151#            for item in table: print item,table[item]               #useful debug - are things shifting?
152            break                   #refinement succeeded - finish up!
153        except TypeError:          #result[1] is None on singular matrix or LinAlgError
154            IfOK = False
155            if not len(varyList):
156                covMatrix = []
157                break
158            G2fil.G2Print ('**** Refinement failed - singular matrix ****',mode='error')
159            if 'Hessian' in Controls['deriv type']:
160                if result[1] is None:
161                    IfOK = False
162                    covMatrix = []
163                    sig = len(varyList)*[None,]
164                    break
165                num = len(varyList)-1
166                for i,val in enumerate(np.flipud(result[2]['psing'])):
167                    if val:
168                        G2fil.G2Print ('Removing parameter: '+varyList[num-i])
169                        del(varyList[num-i])
170            else:
171                Ipvt = result[2]['ipvt']
172                for i,ipvt in enumerate(Ipvt):
173                    if not np.sum(result[2]['fjac'],axis=1)[i]:
174                        G2fil.G2Print ('Removing parameter: '+varyList[ipvt-1])
175                        del(varyList[ipvt-1])
176                        break
177    if IfOK:
178        G2stMth.GetFobsSq(Histograms,Phases,parmDict,calcControls)
179    return IfOK,Rvals,result,covMatrix,sig
180
181def Refine(GPXfile,dlg=None,makeBack=True,refPlotUpdate=None):
182    'Global refinement -- refines to minimize against all histograms'
183    import GSASIImpsubs as G2mp
184    G2mp.InitMP()
185    import pytexture as ptx
186    ptx.pyqlmninit()            #initialize fortran arrays for spherical harmonics
187
188    printFile = open(ospath.splitext(GPXfile)[0]+'.lst','w')
189    G2stIO.ShowBanner(printFile)
190    varyList = []
191    parmDict = {}
192    G2mv.InitVars()
193    Controls = G2stIO.GetControls(GPXfile)
194    G2stIO.ShowControls(Controls,printFile)
195    calcControls = {}
196    calcControls.update(Controls)
197    constrDict,fixedList = G2stIO.GetConstraints(GPXfile)
198    restraintDict = G2stIO.GetRestraints(GPXfile)
199    Histograms,Phases = G2stIO.GetUsedHistogramsAndPhases(GPXfile)
200    if not Phases:
201        G2fil.G2Print (' *** ERROR - you have no phases to refine! ***')
202        G2fil.G2Print (' *** Refine aborted ***')
203        return False,'No phases'
204    if not Histograms:
205        G2fil.G2Print (' *** ERROR - you have no data to refine with! ***')
206        G2fil.G2Print (' *** Refine aborted ***')
207        return False,'No data'
208    rigidbodyDict = G2stIO.GetRigidBodies(GPXfile)
209    rbIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})
210    rbVary,rbDict = G2stIO.GetRigidBodyModels(rigidbodyDict,pFile=printFile)
211    (Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables,MFtables,
212         maxSSwave) = G2stIO.GetPhaseData(Phases,restraintDict,rbIds,pFile=printFile)
213    calcControls['atomIndx'] = atomIndx
214    calcControls['Natoms'] = Natoms
215    calcControls['FFtables'] = FFtables
216    calcControls['BLtables'] = BLtables
217    calcControls['MFtables'] = MFtables
218    calcControls['maxSSwave'] = maxSSwave
219    hapVary,hapDict,controlDict = G2stIO.GetHistogramPhaseData(Phases,Histograms,pFile=printFile)
220    TwConstr,TwFixed = G2stIO.makeTwinFrConstr(Phases,Histograms,hapVary)
221    constrDict += TwConstr
222    fixedList += TwFixed
223    calcControls.update(controlDict)
224    histVary,histDict,controlDict = G2stIO.GetHistogramData(Histograms,pFile=printFile)
225    calcControls.update(controlDict)
226    varyList = rbVary+phaseVary+hapVary+histVary
227    parmDict.update(rbDict)
228    parmDict.update(phaseDict)
229    parmDict.update(hapDict)
230    parmDict.update(histDict)
231    G2stIO.GetFprime(calcControls,Histograms)
232    # do constraint processing
233    varyListStart = tuple(varyList) # save the original varyList before dependent vars are removed
234    msg = G2mv.EvaluateMultipliers(constrDict,parmDict)
235    if msg:
236        return False,'Unable to interpret multiplier(s): '+msg
237    try:
238        G2mv.GenerateConstraints(varyList,constrDict,fixedList,parmDict)
239        #print(G2mv.VarRemapShow(varyList))
240        #print('DependentVars',G2mv.GetDependentVars())
241        #print('IndependentVars',G2mv.GetIndependentVars())
242    except G2mv.ConstraintException:
243        G2fil.G2Print (' *** ERROR - your constraints are internally inconsistent ***')
244        #errmsg, warnmsg = G2mv.CheckConstraints(varyList,constrDict,fixedList)
245        #print 'Errors',errmsg
246        #if warnmsg: print 'Warnings',warnmsg
247        return False,' Constraint error'
248#    print G2mv.VarRemapShow(varyList)
249
250    ifSeq = False
251    printFile.write('\n Refinement results:\n')
252    printFile.write(135*'-'+'\n')
253    try:
254        covData = {}
255        IfOK,Rvals,result,covMatrix,sig = RefineCore(Controls,Histograms,Phases,restraintDict,
256            rigidbodyDict,parmDict,varyList,calcControls,pawleyLookup,ifSeq,printFile,dlg,
257            refPlotUpdate=refPlotUpdate)
258        if IfOK:
259            sigDict = dict(zip(varyList,sig))
260            newCellDict = G2stMth.GetNewCellParms(parmDict,varyList)
261            newAtomDict = G2stMth.ApplyXYZshifts(parmDict,varyList)
262            covData = {'variables':result[0],'varyList':varyList,'sig':sig,'Rvals':Rvals,
263                       'varyListStart':varyListStart,
264                       'covMatrix':covMatrix,'title':GPXfile,'newAtomDict':newAtomDict,
265                       'newCellDict':newCellDict,'freshCOV':True}
266            # add the uncertainties into the esd dictionary (sigDict)
267            sigDict.update(G2mv.ComputeDepESD(covMatrix,varyList,parmDict))
268            G2mv.PrintIndependentVars(parmDict,varyList,sigDict,pFile=printFile)
269            G2stMth.ApplyRBModels(parmDict,Phases,rigidbodyDict,True)
270            G2stIO.SetRigidBodyModels(parmDict,sigDict,rigidbodyDict,printFile)
271            G2stIO.SetPhaseData(parmDict,sigDict,Phases,rbIds,covData,restraintDict,printFile)
272            G2stIO.SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms,calcControls['FFtables'],pFile=printFile)
273            G2stIO.SetHistogramData(parmDict,sigDict,Histograms,calcControls['FFtables'],pFile=printFile)
274            G2stIO.SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases,rigidbodyDict,covData,makeBack)
275            printFile.close()
276            G2fil.G2Print (' Refinement results are in file: '+ospath.splitext(GPXfile)[0]+'.lst')
277            G2fil.G2Print (' ***** Refinement successful *****')
278        else:
279            G2fil.G2Print ('****ERROR - Refinement failed')
280            raise G2obj.G2Exception('****ERROR - Refinement failed')
281    except G2obj.G2RefineCancel as Msg:
282        printFile.close()
283        G2fil.G2Print (' ***** Refinement stopped *****')
284        return False,Msg.msg
285    except G2obj.G2Exception as Msg:  # cell metric error, others?
286        printFile.close()
287        G2fil.G2Print (' ***** Refinement error *****')
288        return False,Msg.msg
289
290#for testing purposes, create a file for testderiv
291    if GSASIIpath.GetConfigValue('debug'):   # and IfOK:
292#needs: values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup
293        fl = open(ospath.splitext(GPXfile)[0]+'.testDeriv','wb')
294        cPickle.dump(result[0],fl,1)
295        cPickle.dump([Histograms,Phases,restraintDict,rigidbodyDict],fl,1)
296        cPickle.dump([constrDict,fixedList,G2mv.GetDependentVars()],fl,1)
297        cPickle.dump(parmDict,fl,1)
298        cPickle.dump(varyList,fl,1)
299        cPickle.dump(calcControls,fl,1)
300        cPickle.dump(pawleyLookup,fl,1)
301        fl.close()
302    if dlg:
303        return True,Rvals
304
305def phaseCheck(phaseVary,Phases,histogram):
306    '''
307    Removes unused parameters from phase varylist if phase not in histogram
308    #TODO - implement "Fix FXU" for seq refinement here - done?
309    '''
310    NewVary = []
311    for phase in Phases:
312        if histogram not in Phases[phase]['Histograms']: continue
313        if Phases[phase]['Histograms'][histogram]['Use']:
314            pId = Phases[phase]['pId']
315            newVary = [item for item in phaseVary if item.split(':')[0] == str(pId)]
316            FixVals = Phases[phase]['Histograms'][histogram].get('Fix FXU',' ')
317            if 'F' in FixVals:
318                newVary = [item for item in newVary if not 'Afrac' in item]
319            if 'X' in FixVals:
320                newVary = [item for item in newVary if not 'dA' in item]
321            if 'U' in FixVals:
322                newVary = [item for item in newVary if not 'AU' in item]
323            if 'M' in FixVals:
324                newVary = [item for item in newVary if not 'AM' in item]
325            NewVary += newVary
326    return NewVary
327
328def SeqRefine(GPXfile,dlg,refPlotUpdate=None):
329    '''Perform a sequential refinement -- cycles through all selected histgrams,
330    one at a time
331    '''
332    import GSASIImpsubs as G2mp
333    G2mp.InitMP()
334    import pytexture as ptx
335    ptx.pyqlmninit()            #initialize fortran arrays for spherical harmonics
336
337    printFile = open(ospath.splitext(GPXfile)[0]+'.lst','w')
338    G2fil.G2Print ('Starting Sequential Refinement')
339    G2stIO.ShowBanner(printFile)
340    Controls = G2stIO.GetControls(GPXfile)
341    G2stIO.ShowControls(Controls,printFile,SeqRef=True)
342    restraintDict = G2stIO.GetRestraints(GPXfile)
343    Histograms,Phases = G2stIO.GetUsedHistogramsAndPhases(GPXfile)
344    if not Phases:
345        G2fil.G2Print (' *** ERROR - you have no phases to refine! ***')
346        G2fil.G2Print (' *** Refine aborted ***')
347        return False,'No phases'
348    if not Histograms:
349        G2fil.G2Print (' *** ERROR - you have no data to refine with! ***')
350        G2fil.G2Print (' *** Refine aborted ***')
351        return False,'No data'
352    rigidbodyDict = G2stIO.GetRigidBodies(GPXfile)
353    rbIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})
354    rbVary,rbDict = G2stIO.GetRigidBodyModels(rigidbodyDict,pFile=printFile)
355    G2mv.InitVars()
356    (Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables,MFtables,
357         maxSSwave) = G2stIO.GetPhaseData(Phases,restraintDict,rbIds,
358                                    Print=False,pFile=printFile,seqRef=True)
359    for item in phaseVary:
360        if '::A0' in item:
361            G2fil.G2Print ('**** WARNING - lattice parameters should not be refined in a sequential refinement ****')
362            G2fil.G2Print ('****           instead use the Dij parameters for each powder histogram            ****')
363            return False,'Lattice parameter refinement error - see console message'
364        if '::C(' in item:
365            G2fil.G2Print ('**** WARNING - phase texture parameters should not be refined in a sequential refinement ****')
366            G2fil.G2Print ('****           instead use the C(L,N) parameters for each powder histogram               ****')
367            return False,'Phase texture refinement error - see console message'
368    if 'Seq Data' in Controls:
369        histNames = Controls['Seq Data']
370    else: # patch from before Controls['Seq Data'] was implemented?
371        histNames = G2stIO.GetHistogramNames(GPXfile,['PWDR',])
372    if Controls.get('Reverse Seq'):
373        histNames.reverse()
374    SeqResult = G2stIO.GetSeqResult(GPXfile)
375#    SeqResult = {'SeqPseudoVars':{},'SeqParFitEqList':[]}
376    Histo = {}
377    NewparmDict = {}
378    G2stIO.SetupSeqSavePhases(GPXfile)
379    for ihst,histogram in enumerate(histNames):
380        if GSASIIpath.GetConfigValue('Show_timing'): t1 = time.time()
381        G2fil.G2Print('\nRefining with '+str(histogram))
382        G2mv.InitVars()
383        (Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,
384             FFtables,BLtables,MFtables,maxSSwave) = G2stIO.GetPhaseData(
385                 Phases,restraintDict,rbIds,
386                 Print=False,pFile=printFile,seqRef=True)
387        ifPrint = False
388        if dlg:
389            dlg.SetTitle('Residual for histogram '+str(ihst))
390        calcControls = {}
391        calcControls['atomIndx'] = atomIndx
392        calcControls['Natoms'] = Natoms
393        calcControls['FFtables'] = FFtables
394        calcControls['BLtables'] = BLtables
395        calcControls['MFtables'] = MFtables
396        calcControls['maxSSwave'] = maxSSwave
397        if histogram not in Histograms:
398            G2fil.G2Print("Error: not found!")
399            continue
400    #TODO - implement "Fix FXU" for seq refinement here - done?
401        hId = Histograms[histogram]['hId']
402        redphaseVary = phaseCheck(phaseVary,Phases,histogram)
403        Histo = {histogram:Histograms[histogram],}
404        hapVary,hapDict,controlDict = G2stIO.GetHistogramPhaseData(Phases,Histo,Print=False)
405        calcControls.update(controlDict)
406        histVary,histDict,controlDict = G2stIO.GetHistogramData(Histo,False)
407        calcControls.update(controlDict)
408        varyList = rbVary+redphaseVary+hapVary+histVary
409#        if not ihst:
410            # save the initial vary list, but without histogram numbers on parameters
411        saveVaryList = varyList[:]
412        for i,item in enumerate(saveVaryList):
413            items = item.split(':')
414            if items[1]:
415                items[1] = ''
416            item = ':'.join(items)
417            saveVaryList[i] = item
418        if not ihst:
419            SeqResult['varyList'] = saveVaryList
420        else:
421            SeqResult['varyList'] = list(set(SeqResult['varyList']+saveVaryList))
422        parmDict = {}
423        parmDict.update(rbDict)
424        parmDict.update(phaseDict)
425        parmDict.update(hapDict)
426        parmDict.update(histDict)
427        if Controls['Copy2Next']:   # update with parms from last histogram
428            #parmDict.update(NewparmDict) # don't use in case extra entries would cause a problem
429            for parm in NewparmDict:
430                if parm in parmDict:
431                    parmDict[parm] = NewparmDict[parm]
432        elif histogram in SeqResult:  # update phase from last seq ref
433            NewparmDict = SeqResult[histogram].get('parmDict',{})
434            for parm in NewparmDict:
435                if '::' in parm and parm in parmDict:
436                    parmDict[parm] = NewparmDict[parm]
437           
438        G2stIO.GetFprime(calcControls,Histo)
439        # do constraint processing
440        #reload(G2mv) # debug
441        constrDict,fixedList = G2stIO.GetConstraints(GPXfile)
442        varyListStart = tuple(varyList) # save the original varyList before dependent vars are removed
443        msg = G2mv.EvaluateMultipliers(constrDict,parmDict)
444        if msg:
445            return False,'Unable to interpret multiplier(s): '+msg
446        try:
447            groups,parmlist = G2mv.GenerateConstraints(varyList,constrDict,fixedList,parmDict,SeqHist=hId)
448#            if GSASIIpath.GetConfigValue('debug'): print("DBG_"+
449#                G2mv.VarRemapShow(varyList,True))
450#            print('DependentVars',G2mv.GetDependentVars())
451#            print('IndependentVars',G2mv.GetIndependentVars())
452            constraintInfo = (groups,parmlist,constrDict,fixedList,ihst)
453        except G2mv.ConstraintException:
454            G2fil.G2Print (' *** ERROR - your constraints are internally inconsistent ***')
455            #errmsg, warnmsg = G2mv.CheckConstraints(varyList,constrDict,fixedList)
456            #print 'Errors',errmsg
457            #if warnmsg: print 'Warnings',warnmsg
458            return False,' Constraint error'
459        #print G2mv.VarRemapShow(varyList)
460        if not ihst:
461            # first histogram to refine against
462            firstVaryList = []
463            for item in varyList:
464                items = item.split(':')
465                if items[1]:
466                    items[1] = ''
467                item = ':'.join(items)
468                firstVaryList.append(item)
469            newVaryList = firstVaryList
470        else:
471            newVaryList = []
472            for item in varyList:
473                items = item.split(':')
474                if items[1]:
475                    items[1] = ''
476                item = ':'.join(items)
477                newVaryList.append(item)
478        if newVaryList != firstVaryList and Controls['Copy2Next']:
479            # variable lists are expected to match between sequential refinements when Copy2Next is on
480            #print '**** ERROR - variable list for this histogram does not match previous'
481            #print '     Copy of variables is not possible'
482            #print '\ncurrent histogram',histogram,'has',len(newVaryList),'variables'
483            combined = list(set(firstVaryList+newVaryList))
484            c = [var for var in combined if var not in newVaryList]
485            p = [var for var in combined if var not in firstVaryList]
486            G2fil.G2Print('*** Variables change ***')
487            for typ,vars in [('Removed',c),('Added',p)]:
488                line = '  '+typ+': '
489                if vars:
490                    for var in vars:
491                        if len(line) > 70:
492                            G2fil.G2Print(line)
493                            line = '    '
494                        line += var + ', '
495                else:
496                        line += 'none, '
497                G2fil.G2Print(line[:-2])
498            firstVaryList = newVaryList
499
500        ifSeq = True
501        printFile.write('\n Refinement results for histogram: %s\n'%histogram)
502        printFile.write(135*'-'+'\n')
503        try:
504            IfOK,Rvals,result,covMatrix,sig = RefineCore(Controls,Histo,Phases,restraintDict,
505                rigidbodyDict,parmDict,varyList,calcControls,pawleyLookup,ifSeq,printFile,dlg,
506                refPlotUpdate=refPlotUpdate)
507            G2fil.G2Print ('  wR = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f, last delta chi = %.4f'%(
508                Rvals['Rwp'],Rvals['chisq'],Rvals['GOF']**2,Rvals['DelChi2']))
509            # add the uncertainties into the esd dictionary (sigDict)
510            if not IfOK:
511                G2fil.G2Print('***** Sequential refinement failed at histogram '+histogram,mode='warn')
512                break
513            sigDict = dict(zip(varyList,sig))
514            # the uncertainties for dependent constrained parms into the esd dict
515            sigDict.update(G2mv.ComputeDepESD(covMatrix,varyList,parmDict))
516
517            # a dict with values & esds for dependent (constrained) parameters - avoid extraneous holds
518            depParmDict = {i:(parmDict[i],sigDict[i]) for i in varyListStart if i in sigDict and i not in varyList}
519            newCellDict = copy.deepcopy(G2stMth.GetNewCellParms(parmDict,varyList))
520            newAtomDict = copy.deepcopy(G2stMth.ApplyXYZshifts(parmDict,varyList))
521            histRefData = {
522                'variables':result[0],'varyList':varyList,'sig':sig,'Rvals':Rvals,
523                'varyListStart':varyListStart,
524                'covMatrix':covMatrix,'title':histogram,'newAtomDict':newAtomDict,
525                'newCellDict':newCellDict,'depParmDict':depParmDict,
526                'constraintInfo':constraintInfo,
527                'parmDict':parmDict}
528            SeqResult[histogram] = histRefData
529            G2stMth.ApplyRBModels(parmDict,Phases,rigidbodyDict,True)
530#            G2stIO.SetRigidBodyModels(parmDict,sigDict,rigidbodyDict,printFile)
531            G2stIO.SetHistogramPhaseData(parmDict,sigDict,Phases,Histo,None,ifPrint,printFile)
532            G2stIO.SetHistogramData(parmDict,sigDict,Histo,None,ifPrint,printFile)
533            G2stIO.SaveUpdatedHistogramsAndPhases(GPXfile,Histo,Phases,rigidbodyDict,histRefData)
534            NewparmDict = {}
535            # make dict of varied parameters in current histogram, renamed to
536            # next histogram, for use in next refinement.
537            if Controls['Copy2Next'] and ihst < len(histNames)-1:
538                hId = Histo[histogram]['hId'] # current histogram
539                nexthId = Histograms[histNames[ihst+1]]['hId']
540                for parm in set(list(varyList)+list(varyListStart)):
541                    items = parm.split(':')
542                    if len(items) < 3: 
543                        continue
544                    if str(hId) in items[1]:
545                        items[1] = str(nexthId)
546                        newparm = ':'.join(items)
547                        NewparmDict[newparm] = parmDict[parm]
548                    else:
549                        if items[2].startswith('dA'): parm = parm.replace(':dA',':A') 
550                        NewparmDict[parm] = parmDict[parm]
551                   
552        except G2obj.G2RefineCancel as Msg:
553            printFile.close()
554            G2fil.G2Print (' ***** Refinement stopped *****')
555            return False,Msg.msg
556        except G2obj.G2Exception as Msg:  # cell metric error, others?
557            printFile.close()
558            G2fil.G2Print (' ***** Refinement error *****')
559            return False,Msg.msg
560        if GSASIIpath.GetConfigValue('Show_timing'):
561            t2 = time.time()
562            G2fil.G2Print("Fit step time {:.2f} sec.".format(t2-t1))
563            t1 = t2
564    SeqResult['histNames'] = [itm for itm in G2stIO.GetHistogramNames(GPXfile,['PWDR',]) if itm in SeqResult.keys()]
565    G2stIO.SetSeqResult(GPXfile,Histograms,SeqResult)
566    printFile.close()
567    G2fil.G2Print (' Sequential refinement results are in file: '+ospath.splitext(GPXfile)[0]+'.lst')
568    G2fil.G2Print (' ***** Sequential refinement successful *****')
569    return True,'Success'
570
571def RetDistAngle(DisAglCtls,DisAglData):
572    '''Compute and return distances and angles
573
574    :param dict DisAglCtls: contains distance/angle radii usually defined using
575       :func:`GSASIIctrlGUI.DisAglDialog`
576    :param dict DisAglData: contains phase data:
577       Items 'OrigAtoms' and 'TargAtoms' contain the atoms to be used
578       for distance/angle origins and atoms to be used as targets.
579       Item 'SGData' has the space group information (see :ref:`Space Group object<SGData_table>`)
580
581    :returns: AtomLabels,DistArray,AngArray where:
582
583      **AtomLabels** is a dict of atom labels, keys are the atom number
584
585      **DistArray** is a dict keyed by the origin atom number where the value is a list
586      of distance entries. The value for each distance is a list containing:
587
588        0) the target atom number (int);
589        1) the unit cell offsets added to x,y & z (tuple of int values)
590        2) the symmetry operator number (which may be modified to indicate centering and center of symmetry)
591        3) an interatomic distance in A (float)
592        4) an uncertainty on the distance in A or 0.0 (float)
593
594      **AngArray** is a dict keyed by the origin (central) atom number where
595      the value is a list of
596      angle entries. The value for each angle entry consists of three values:
597
598        0) a distance item reference for one neighbor (int)
599        1) a distance item reference for a second neighbor (int)
600        2) a angle, uncertainty pair; the s.u. may be zero (tuple of two floats)
601
602      The AngArray distance reference items refer directly to the index of the items in the
603      DistArray item for the list of distances for the central atom.
604    '''
605    import numpy.ma as ma
606
607    SGData = DisAglData['SGData']
608    Cell = DisAglData['Cell']
609
610    Amat,Bmat = G2lat.cell2AB(Cell[:6])
611    covData = {}
612    if 'covData' in DisAglData:
613        covData = DisAglData['covData']
614        covMatrix = covData['covMatrix']
615        varyList = covData['varyList']
616        pfx = str(DisAglData['pId'])+'::'
617
618    Factor = DisAglCtls['Factors']
619    Radii = dict(zip(DisAglCtls['AtomTypes'],zip(DisAglCtls['BondRadii'],DisAglCtls['AngleRadii'])))
620    indices = (-2,-1,0,1,2)
621    Units = np.array([[h,k,l] for h in indices for k in indices for l in indices])
622    origAtoms = DisAglData['OrigAtoms']
623    targAtoms = DisAglData['TargAtoms']
624    AtomLabels = {}
625    for Oatom in origAtoms:
626        AtomLabels[Oatom[0]] = Oatom[1]
627    for Oatom in targAtoms:
628        AtomLabels[Oatom[0]] = Oatom[1]
629    DistArray = {}
630    AngArray = {}
631    for Oatom in origAtoms:
632        DistArray[Oatom[0]] = []
633        AngArray[Oatom[0]] = []
634        OxyzNames = ''
635        IndBlist = []
636        Dist = []
637        Vect = []
638        VectA = []
639        angles = []
640        for Tatom in targAtoms:
641            Xvcov = []
642            TxyzNames = ''
643            if 'covData' in DisAglData:
644                OxyzNames = [pfx+'dAx:%d'%(Oatom[0]),pfx+'dAy:%d'%(Oatom[0]),pfx+'dAz:%d'%(Oatom[0])]
645                TxyzNames = [pfx+'dAx:%d'%(Tatom[0]),pfx+'dAy:%d'%(Tatom[0]),pfx+'dAz:%d'%(Tatom[0])]
646                Xvcov = G2mth.getVCov(OxyzNames+TxyzNames,varyList,covMatrix)
647            result = G2spc.GenAtom(Tatom[3:6],SGData,False,Move=False)
648            BsumR = (Radii[Oatom[2]][0]+Radii[Tatom[2]][0])*Factor[0]
649            AsumR = (Radii[Oatom[2]][1]+Radii[Tatom[2]][1])*Factor[1]
650            for [Txyz,Top,Tunit,Spn] in result:
651                Dx = (Txyz-np.array(Oatom[3:6]))+Units
652                dx = np.inner(Amat,Dx)
653                dist = ma.masked_less(np.sqrt(np.sum(dx**2,axis=0)),0.5)
654                IndB = ma.nonzero(ma.masked_greater(dist-BsumR,0.))
655                if np.any(IndB):
656                    for indb in IndB:
657                        for i in range(len(indb)):
658                            if str(dx.T[indb][i]) not in IndBlist:
659                                IndBlist.append(str(dx.T[indb][i]))
660                                unit = Units[indb][i]
661                                tunit = (unit[0]+Tunit[0],unit[1]+Tunit[1],unit[2]+Tunit[2])
662                                pdpx = G2mth.getDistDerv(Oatom[3:6],Tatom[3:6],Amat,unit,Top,SGData)
663                                sig = 0.0
664                                if len(Xvcov):
665                                    sig = np.sqrt(np.inner(pdpx,np.inner(pdpx,Xvcov)))
666                                Dist.append([Oatom[0],Tatom[0],tunit,Top,ma.getdata(dist[indb])[i],sig])
667                                if (Dist[-1][-2]-AsumR) <= 0.:
668                                    Vect.append(dx.T[indb][i]/Dist[-1][-2])
669                                    VectA.append([OxyzNames,np.array(Oatom[3:6]),TxyzNames,np.array(Tatom[3:6]),unit,Top])
670                                else:
671                                    Vect.append([0.,0.,0.])
672                                    VectA.append([])
673        for D in Dist:
674            DistArray[Oatom[0]].append(D[1:])
675        Vect = np.array(Vect)
676        angles = np.zeros((len(Vect),len(Vect)))
677        angsig = np.zeros((len(Vect),len(Vect)))
678        for i,veca in enumerate(Vect):
679            if np.any(veca):
680                for j,vecb in enumerate(Vect):
681                    if np.any(vecb):
682                        angles[i][j],angsig[i][j] = G2mth.getAngSig(VectA[i],VectA[j],Amat,SGData,covData)
683                        if i <= j: continue
684                        AngArray[Oatom[0]].append((i,j,
685                            G2mth.getAngSig(VectA[i],VectA[j],Amat,SGData,covData)))
686    return AtomLabels,DistArray,AngArray
687
688def PrintDistAngle(DisAglCtls,DisAglData,out=sys.stdout):
689    '''Print distances and angles
690
691    :param dict DisAglCtls: contains distance/angle radii usually defined using
692       :func:`GSASIIctrlGUI.DisAglDialog`
693    :param dict DisAglData: contains phase data:
694       Items 'OrigAtoms' and 'TargAtoms' contain the atoms to be used
695       for distance/angle origins and atoms to be used as targets.
696       Item 'SGData' has the space group information (see :ref:`Space Group object<SGData_table>`)
697    :param file out: file object for output. Defaults to sys.stdout.
698    '''
699    def MyPrint(s):
700        out.write(s+'\n')
701        # print(s,file=out) # use in Python 3
702
703    def ShowBanner(name):
704        MyPrint(80*'*')
705        MyPrint('   Interatomic Distances and Angles for phase '+name)
706        MyPrint((80*'*')+'\n')
707
708    ShowBanner(DisAglCtls['Name'])
709    SGData = DisAglData['SGData']
710    SGtext,SGtable = G2spc.SGPrint(SGData)
711    for line in SGtext: MyPrint(line)
712    if len(SGtable) > 1:
713        for i,item in enumerate(SGtable[::2]):
714            if 2*i+1 == len(SGtable):
715                line = ' %s'%(item.ljust(30))
716            else:
717                line = ' %s %s'%(item.ljust(30),SGtable[2*i+1].ljust(30))
718            MyPrint(line)
719    else:
720        MyPrint(' ( 1)    %s'%(SGtable[0])) #triclinic case
721    Cell = DisAglData['Cell']
722
723    Amat,Bmat = G2lat.cell2AB(Cell[:6])
724    covData = {}
725    if 'covData' in DisAglData:
726        covData = DisAglData['covData']
727        pfx = str(DisAglData['pId'])+'::'
728        A = G2lat.cell2A(Cell[:6])
729        cellSig = G2stIO.getCellEsd(pfx,SGData,A,covData)
730        names = [' a = ',' b = ',' c = ',' alpha = ',' beta = ',' gamma = ',' Volume = ']
731        valEsd = [G2mth.ValEsd(Cell[i],cellSig[i],True) for i in range(7)]
732        line = '\n Unit cell:'
733        for name,vals in zip(names,valEsd):
734            line += name+vals
735        MyPrint(line)
736    else:
737        MyPrint('\n Unit cell: a = '+('%.5f'%Cell[0])+' b = '+('%.5f'%Cell[1])+' c = '+('%.5f'%Cell[2])+
738            ' alpha = '+('%.3f'%Cell[3])+' beta = '+('%.3f'%Cell[4])+' gamma = '+
739            ('%.3f'%Cell[5])+' Volume = '+('%.3f'%Cell[6]))
740
741    AtomLabels,DistArray,AngArray = RetDistAngle(DisAglCtls,DisAglData)
742    origAtoms = DisAglData['OrigAtoms']
743    for Oatom in origAtoms:
744        i = Oatom[0]
745        Dist = DistArray[i]
746        nDist = len(Dist)
747        angles = np.zeros((nDist,nDist))
748        angsig = np.zeros((nDist,nDist))
749        for k,j,tup in AngArray[i]:
750            angles[k][j],angsig[k][j] = angles[j][k],angsig[j][k] = tup
751        line = ''
752        for i,x in enumerate(Oatom[3:6]):
753            line += ('%12.5f'%x).rstrip('0')
754        MyPrint('\n Distances & angles for '+Oatom[1]+' at '+line.rstrip())
755        MyPrint(80*'*')
756        line = ''
757        for dist in Dist[:-1]:
758            line += '%12s'%(AtomLabels[dist[0]].center(12))
759        MyPrint('  To       cell +(sym. op.)      dist.  '+line.rstrip())
760        for i,dist in enumerate(Dist):
761            line = ''
762            for j,angle in enumerate(angles[i][0:i]):
763                sig = angsig[i][j]
764                if angle:
765                    if sig:
766                        line += '%12s'%(G2mth.ValEsd(angle,sig,True).center(12))
767                    else:
768                        val = '%.3f'%(angle)
769                        line += '%12s'%(val.center(12))
770                else:
771                    line += 12*' '
772            if dist[4]:            #sig exists!
773                val = G2mth.ValEsd(dist[3],dist[4])
774            else:
775                val = '%8.4f'%(dist[3])
776            tunit = '[%2d%2d%2d]'% dist[1]
777            MyPrint((%8s%10s+(%4d) %12s'%(AtomLabels[dist[0]].ljust(8),tunit.ljust(10),dist[2],val.center(12)))+line.rstrip())
778
779def DisAglTor(DATData):
780    'Needs a doc string'
781    SGData = DATData['SGData']
782    Cell = DATData['Cell']
783
784    Amat,Bmat = G2lat.cell2AB(Cell[:6])
785    covData = {}
786    pfx = ''
787    if 'covData' in DATData:
788        covData = DATData['covData']
789        pfx = str(DATData['pId'])+'::'
790    Datoms = []
791    Oatoms = []
792    for i,atom in enumerate(DATData['Datoms']):
793        symop = atom[-1].split('+')
794        if len(symop) == 1:
795            symop.append('0,0,0')
796        symop[0] = int(symop[0])
797        symop[1] = eval(symop[1])
798        atom.append(symop)
799        Datoms.append(atom)
800        oatom = DATData['Oatoms'][i]
801        names = ['','','']
802        if pfx:
803            names = [pfx+'dAx:'+str(oatom[0]),pfx+'dAy:'+str(oatom[0]),pfx+'dAz:'+str(oatom[0])]
804        oatom += [names,]
805        Oatoms.append(oatom)
806    atmSeq = [atom[1]+'('+atom[-2]+')' for atom in Datoms]
807    if DATData['Natoms'] == 4:  #torsion
808        Tors,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData)
809        G2fil.G2Print (' Torsion angle for %s atom sequence: %s = %s'%(DATData['Name'],str(atmSeq).replace("'","")[1:-1],G2mth.ValEsd(Tors,sig)))
810        G2fil.G2Print (' NB: Atom sequence determined by selection order')
811        return      # done with torsion
812    elif DATData['Natoms'] == 3:  #angle
813        Ang,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData)
814        G2fil.G2Print (' Angle in %s for atom sequence: %s = %s'%(DATData['Name'],str(atmSeq).replace("'","")[1:-1],G2mth.ValEsd(Ang,sig)))
815        G2fil.G2Print (' NB: Atom sequence determined by selection order')
816    else:   #2 atoms - distance
817        Dist,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData)
818        G2fil.G2Print (' Distance in %s for atom sequence: %s = %s'%(DATData['Name'],str(atmSeq).replace("'","")[1:-1],G2mth.ValEsd(Dist,sig)))
819
820def BestPlane(PlaneData):
821    'Needs a doc string'
822
823    def ShowBanner(name):
824        G2fil.G2Print (80*'*')
825        G2fil.G2Print ('   Best plane result for phase '+name)
826        G2fil.G2Print (80*'*','\n')
827
828    ShowBanner(PlaneData['Name'])
829
830    Cell = PlaneData['Cell']
831    Amat,Bmat = G2lat.cell2AB(Cell[:6])
832    Atoms = PlaneData['Atoms']
833    sumXYZ = np.zeros(3)
834    XYZ = []
835    Natoms = len(Atoms)
836    for atom in Atoms:
837        xyz = np.array(atom[3:6])
838        XYZ.append(xyz)
839        sumXYZ += xyz
840    sumXYZ /= Natoms
841    XYZ = np.array(XYZ)-sumXYZ
842    XYZ = np.inner(Amat,XYZ).T
843    Zmat = np.zeros((3,3))
844    for i,xyz in enumerate(XYZ):
845        Zmat += np.outer(xyz.T,xyz)
846    G2fil.G2Print (' Selected atoms centered at %10.5f %10.5f %10.5f'%(sumXYZ[0],sumXYZ[1],sumXYZ[2]))
847    Evec,Emat = nl.eig(Zmat)
848    Evec = np.sqrt(Evec)/(Natoms-3)
849    Order = np.argsort(Evec)
850    XYZ = np.inner(XYZ,Emat.T).T
851    XYZ = np.array([XYZ[Order[2]],XYZ[Order[1]],XYZ[Order[0]]]).T
852    G2fil.G2Print (' Atoms in Cartesian best plane coordinates:')
853    G2fil.G2Print (' Name         X         Y         Z')
854    for i,xyz in enumerate(XYZ):
855        G2fil.G2Print (' %6s%10.3f%10.3f%10.3f'%(Atoms[i][1].ljust(6),xyz[0],xyz[1],xyz[2]))
856    G2fil.G2Print ('\n Best plane RMS X =%8.3f, Y =%8.3f, Z =%8.3f'%(Evec[Order[2]],Evec[Order[1]],Evec[Order[0]]))
857
858def main():
859    'Called to run a refinement when this module is executed '
860    starttime = time.time()
861    arg = sys.argv
862    if len(arg) > 1:
863        GPXfile = arg[1]
864        if not ospath.exists(GPXfile):
865            G2fil.G2Print ('ERROR - '+GPXfile+" doesn't exist!")
866            exit()
867    else:
868        G2fil.G2Print ('ERROR - missing filename')
869        exit()
870    # TODO: figure out if this is a sequential refinement and call SeqRefine(GPXfile,None)
871    Refine(GPXfile,None)
872    G2fil.G2Print("Done. Execution time {:.2f} sec.".format(time.time()-starttime))
873
874if __name__ == '__main__':
875    GSASIIpath.InvokeDebugOpts()
876    main()
Note: See TracBrowser for help on using the repository browser.