source: trunk/GSASIIstrMain.py @ 4021

Last change on this file since 4021 was 4021, checked in by toby, 4 years ago

implement filter for screen messages; start to replace print() with G2Print(); scripting changes for PDF
  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 39.0 KB
Line 
1# -*- coding: utf-8 -*-
2'''
3*GSASIIstrMain: main structure routine*
4---------------------------------------
5
6'''
7########### SVN repository information ###################
8# $Date: 2019-06-12 03:05:07 +0000 (Wed, 12 Jun 2019) $
9# $Author: toby $
10# $Revision: 4021 $
11# $URL: trunk/GSASIIstrMain.py $
12# $Id: GSASIIstrMain.py 4021 2019-06-12 03:05:07Z toby $
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: 4021 $")
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,maxSSwave = \
212        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    if True:
254#    try:
255        covData = {}
256        IfOK,Rvals,result,covMatrix,sig = RefineCore(Controls,Histograms,Phases,restraintDict,
257            rigidbodyDict,parmDict,varyList,calcControls,pawleyLookup,ifSeq,printFile,dlg,
258            refPlotUpdate=refPlotUpdate)
259        if IfOK:
260            sigDict = dict(zip(varyList,sig))
261            newCellDict = G2stMth.GetNewCellParms(parmDict,varyList)
262            newAtomDict = G2stMth.ApplyXYZshifts(parmDict,varyList)
263            covData = {'variables':result[0],'varyList':varyList,'sig':sig,'Rvals':Rvals,
264                       'varyListStart':varyListStart,
265                       'covMatrix':covMatrix,'title':GPXfile,'newAtomDict':newAtomDict,
266                       'newCellDict':newCellDict,'freshCOV':True}
267            # add the uncertainties into the esd dictionary (sigDict)
268            sigDict.update(G2mv.ComputeDepESD(covMatrix,varyList,parmDict))
269            G2mv.PrintIndependentVars(parmDict,varyList,sigDict,pFile=printFile)
270            G2stMth.ApplyRBModels(parmDict,Phases,rigidbodyDict,True)
271            G2stIO.SetRigidBodyModels(parmDict,sigDict,rigidbodyDict,printFile)
272            G2stIO.SetPhaseData(parmDict,sigDict,Phases,rbIds,covData,restraintDict,printFile)
273            G2stIO.SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms,calcControls['FFtables'],pFile=printFile)
274            G2stIO.SetHistogramData(parmDict,sigDict,Histograms,calcControls['FFtables'],pFile=printFile)
275            G2stIO.SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases,rigidbodyDict,covData,makeBack)
276            printFile.close()
277            G2fil.G2Print (' Refinement results are in file: '+ospath.splitext(GPXfile)[0]+'.lst')
278            G2fil.G2Print (' ***** Refinement successful *****')
279        else:
280            G2fil.G2Print ('****ERROR - Refinement failed')
281            raise G2obj.G2Exception('****ERROR - Refinement failed')
282#    except G2obj.G2Exception(Msg):
283#        printFile.close()
284#        return False,Msg.msg
285
286#for testing purposes, create a file for testderiv
287    if GSASIIpath.GetConfigValue('debug'):   # and IfOK:
288#needs: values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup
289        fl = open(ospath.splitext(GPXfile)[0]+'.testDeriv','wb')
290        cPickle.dump(result[0],fl,1)
291        cPickle.dump([Histograms,Phases,restraintDict,rigidbodyDict],fl,1)
292        cPickle.dump([constrDict,fixedList,G2mv.GetDependentVars()],fl,1)
293        cPickle.dump(parmDict,fl,1)
294        cPickle.dump(varyList,fl,1)
295        cPickle.dump(calcControls,fl,1)
296        cPickle.dump(pawleyLookup,fl,1)
297        fl.close()
298    if dlg:
299        return True,Rvals
300
301def phaseCheck(phaseVary,Phases,histogram):
302    '''
303    Removes unused parameters from phase varylist if phase not in histogram
304    #TODO - implement "Fix FXU" for seq refinement here - done?
305    '''
306    NewVary = []
307    for phase in Phases:
308        if histogram not in Phases[phase]['Histograms']: continue
309        if Phases[phase]['Histograms'][histogram]['Use']:
310            pId = Phases[phase]['pId']
311            newVary = [item for item in phaseVary if item.split(':')[0] == str(pId)]
312            FixVals = Phases[phase]['Histograms'][histogram].get('Fix FXU',' ')
313            if 'F' in FixVals:
314                newVary = [item for item in newVary if not 'Afrac' in item]
315            if 'X' in FixVals:
316                newVary = [item for item in newVary if not 'dA' in item]
317            if 'U' in FixVals:
318                newVary = [item for item in newVary if not 'AU' in item]
319            if 'M' in FixVals:
320                newVary = [item for item in newVary if not 'AM' in item]
321            NewVary += newVary
322    return NewVary
323
324def SeqRefine(GPXfile,dlg,refPlotUpdate=None):
325    '''Perform a sequential refinement -- cycles through all selected histgrams,
326    one at a time
327    '''
328    import GSASIImpsubs as G2mp
329    G2mp.InitMP()
330    import pytexture as ptx
331    ptx.pyqlmninit()            #initialize fortran arrays for spherical harmonics
332
333    printFile = open(ospath.splitext(GPXfile)[0]+'.lst','w')
334    G2fil.G2Print ('Starting Sequential Refinement')
335    G2stIO.ShowBanner(printFile)
336    Controls = G2stIO.GetControls(GPXfile)
337    G2stIO.ShowControls(Controls,printFile,SeqRef=True)
338    restraintDict = G2stIO.GetRestraints(GPXfile)
339    Histograms,Phases = G2stIO.GetUsedHistogramsAndPhases(GPXfile)
340    if not Phases:
341        G2fil.G2Print (' *** ERROR - you have no phases to refine! ***')
342        G2fil.G2Print (' *** Refine aborted ***')
343        return False,'No phases'
344    if not Histograms:
345        G2fil.G2Print (' *** ERROR - you have no data to refine with! ***')
346        G2fil.G2Print (' *** Refine aborted ***')
347        return False,'No data'
348    rigidbodyDict = G2stIO.GetRigidBodies(GPXfile)
349    rbIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})
350    rbVary,rbDict = G2stIO.GetRigidBodyModels(rigidbodyDict,pFile=printFile)
351    Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables,MFtables,maxSSwave = \
352        G2stIO.GetPhaseData(Phases,restraintDict,rbIds,False,printFile,seqRef=True)
353    for item in phaseVary:
354        if '::A0' in item:
355            G2fil.G2Print ('**** WARNING - lattice parameters should not be refined in a sequential refinement ****')
356            G2fil.G2Print ('****           instead use the Dij parameters for each powder histogram            ****')
357            return False,'Lattice parameter refinement error - see console message'
358        if '::C(' in item:
359            G2fil.G2Print ('**** WARNING - phase texture parameters should not be refined in a sequential refinement ****')
360            G2fil.G2Print ('****           instead use the C(L,N) parameters for each powder histogram               ****')
361            return False,'Phase texture refinement error - see console message'
362    if 'Seq Data' in Controls:
363        histNames = Controls['Seq Data']
364    else: # patch from before Controls['Seq Data'] was implemented?
365        histNames = G2stIO.GetHistogramNames(GPXfile,['PWDR',])
366    if Controls.get('Reverse Seq'):
367        histNames.reverse()
368    SeqResult = G2stIO.GetSeqResult(GPXfile)
369#    SeqResult = {'SeqPseudoVars':{},'SeqParFitEqList':[]}
370    makeBack = True
371    Histo = {}
372    NewparmDict = {}
373    G2stIO.SetupSeqSavePhases(GPXfile)
374    for ihst,histogram in enumerate(histNames):
375        if GSASIIpath.GetConfigValue('Show_timing'): t1 = time.time()
376        G2fil.G2Print('\nRefining with '+str(histogram))
377        ifPrint = False
378        if dlg:
379            dlg.SetTitle('Residual for histogram '+str(ihst))
380        calcControls = {}
381        calcControls['atomIndx'] = atomIndx
382        calcControls['Natoms'] = Natoms
383        calcControls['FFtables'] = FFtables
384        calcControls['BLtables'] = BLtables
385        calcControls['MFtables'] = MFtables
386        calcControls['maxSSwave'] = maxSSwave
387        if histogram not in Histograms:
388            G2fil.G2Print("Error: not found!")
389            continue
390    #TODO - implement "Fix FXU" for seq refinement here - done?
391        hId = Histograms[histogram]['hId']
392        redphaseVary = phaseCheck(phaseVary,Phases,histogram)
393        Histo = {histogram:Histograms[histogram],}
394        hapVary,hapDict,controlDict = G2stIO.GetHistogramPhaseData(Phases,Histo,Print=False)
395        calcControls.update(controlDict)
396        histVary,histDict,controlDict = G2stIO.GetHistogramData(Histo,False)
397        calcControls.update(controlDict)
398        varyList = rbVary+redphaseVary+hapVary+histVary
399#        if not ihst:
400            # save the initial vary list, but without histogram numbers on parameters
401        saveVaryList = varyList[:]
402        for i,item in enumerate(saveVaryList):
403            items = item.split(':')
404            if items[1]:
405                items[1] = ''
406            item = ':'.join(items)
407            saveVaryList[i] = item
408        if not ihst:
409            SeqResult['varyList'] = saveVaryList
410        else:
411            SeqResult['varyList'] = list(set(SeqResult['varyList']+saveVaryList))
412        parmDict = {}
413        parmDict.update(rbDict)
414        parmDict.update(phaseDict)
415        parmDict.update(hapDict)
416        parmDict.update(histDict)
417        if Controls['Copy2Next']:
418            #parmDict.update(NewparmDict) # don't use in case extra entries would cause a problem
419            for parm in NewparmDict:
420                if parm in parmDict:
421                    parmDict[parm] = NewparmDict[parm]
422        G2stIO.GetFprime(calcControls,Histo)
423        # do constraint processing
424        #reload(G2mv) # debug
425        G2mv.InitVars()
426        constrDict,fixedList = G2stIO.GetConstraints(GPXfile)
427        varyListStart = tuple(varyList) # save the original varyList before dependent vars are removed
428        msg = G2mv.EvaluateMultipliers(constrDict,parmDict)
429        if msg:
430            return False,'Unable to interpret multiplier(s): '+msg
431        try:
432            groups,parmlist = G2mv.GenerateConstraints(varyList,constrDict,fixedList,parmDict,SeqHist=hId)
433#            if GSASIIpath.GetConfigValue('debug'): print("DBG_"+
434#                G2mv.VarRemapShow(varyList,True))
435            constraintInfo = (groups,parmlist,constrDict,fixedList,ihst)
436        except G2mv.ConstraintException:
437            G2fil.G2Print (' *** ERROR - your constraints are internally inconsistent ***')
438            #errmsg, warnmsg = G2mv.CheckConstraints(varyList,constrDict,fixedList)
439            #print 'Errors',errmsg
440            #if warnmsg: print 'Warnings',warnmsg
441            return False,' Constraint error'
442        #print G2mv.VarRemapShow(varyList)
443        if not ihst:
444            # first histogram to refine against
445            firstVaryList = []
446            for item in varyList:
447                items = item.split(':')
448                if items[1]:
449                    items[1] = ''
450                item = ':'.join(items)
451                firstVaryList.append(item)
452            newVaryList = firstVaryList
453        else:
454            newVaryList = []
455            for item in varyList:
456                items = item.split(':')
457                if items[1]:
458                    items[1] = ''
459                item = ':'.join(items)
460                newVaryList.append(item)
461        if newVaryList != firstVaryList and Controls['Copy2Next']:
462            # variable lists are expected to match between sequential refinements when Copy2Next is on
463            #print '**** ERROR - variable list for this histogram does not match previous'
464            #print '     Copy of variables is not possible'
465            #print '\ncurrent histogram',histogram,'has',len(newVaryList),'variables'
466            combined = list(set(firstVaryList+newVaryList))
467            c = [var for var in combined if var not in newVaryList]
468            p = [var for var in combined if var not in firstVaryList]
469            G2fil.G2Print('*** Variables change ***')
470            for typ,vars in [('Removed',c),('Added',p)]:
471                line = '  '+typ+': '
472                if vars:
473                    for var in vars:
474                        if len(line) > 70:
475                            G2fil.G2Print(line)
476                            line = '    '
477                        line += var + ', '
478                else:
479                        line += 'none, '
480                G2fil.G2Print(line[:-2])
481            firstVaryList = newVaryList
482
483        ifSeq = True
484        printFile.write('\n Refinement results for histogram: %s\n'%histogram)
485        printFile.write(135*'-'+'\n')
486        if True:
487#        try:
488            IfOK,Rvals,result,covMatrix,sig = RefineCore(Controls,Histo,Phases,restraintDict,
489                rigidbodyDict,parmDict,varyList,calcControls,pawleyLookup,ifSeq,printFile,dlg,
490                refPlotUpdate=refPlotUpdate)
491            G2fil.G2Print ('  wR = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f, last delta chi = %.4f'%(
492                Rvals['Rwp'],Rvals['chisq'],Rvals['GOF']**2,Rvals['DelChi2']))
493            # add the uncertainties into the esd dictionary (sigDict)
494            if not IfOK:
495                G2fil.G2Print('***** Sequential refinement failed at histogram '+histogram,mode='warn')
496                break
497            sigDict = dict(zip(varyList,sig))
498            # the uncertainties for dependent constrained parms into the esd dict
499            sigDict.update(G2mv.ComputeDepESD(covMatrix,varyList,parmDict))
500
501            # a dict with values & esds for dependent (constrained) parameters - avoid extraneous holds
502            depParmDict = {i:(parmDict[i],sigDict[i]) for i in varyListStart if i in sigDict and i not in varyList}
503            newCellDict = copy.deepcopy(G2stMth.GetNewCellParms(parmDict,varyList))
504            newAtomDict = copy.deepcopy(G2stMth.ApplyXYZshifts(parmDict,varyList))
505            histRefData = {
506                'variables':result[0],'varyList':varyList,'sig':sig,'Rvals':Rvals,
507                'varyListStart':varyListStart,
508                'covMatrix':covMatrix,'title':histogram,'newAtomDict':newAtomDict,
509                'newCellDict':newCellDict,'depParmDict':depParmDict,
510                'constraintInfo':constraintInfo,
511                'parmDict':parmDict}
512            SeqResult[histogram] = histRefData
513            G2stMth.ApplyRBModels(parmDict,Phases,rigidbodyDict,True)
514#            G2stIO.SetRigidBodyModels(parmDict,sigDict,rigidbodyDict,printFile)
515            G2stIO.SetHistogramPhaseData(parmDict,sigDict,Phases,Histo,None,ifPrint,printFile)
516            G2stIO.SetHistogramData(parmDict,sigDict,Histo,None,ifPrint,printFile)
517            G2stIO.SaveUpdatedHistogramsAndPhases(GPXfile,Histo,Phases,rigidbodyDict,histRefData)
518            makeBack = False
519            NewparmDict = {}
520            # make dict of varied parameters in current histogram, renamed to
521            # next histogram, for use in next refinement.
522            if Controls['Copy2Next'] and ihst < len(histNames)-1:
523                hId = Histo[histogram]['hId'] # current histogram
524                nexthId = Histograms[histNames[ihst+1]]['hId']
525                for parm in set(list(varyList)+list(varyListStart)):
526                    items = parm.split(':')
527                    if len(items) < 3: continue
528                    if str(hId) in items[1]:
529                        items[1] = str(nexthId)
530                        newparm = ':'.join(items)
531                        NewparmDict[newparm] = parmDict[parm]
532#        except G2obj.G2Exception(Msg):
533#            printFile.close()
534#            G2fil.G2Print (' ***** Refinement aborted *****')
535#            return False,Msg.msg
536        if GSASIIpath.GetConfigValue('Show_timing'):
537            t2 = time.time()
538            G2fil.G2Print("Fit step time {:.2f} sec.".format(t2-t1))
539            t1 = t2
540    SeqResult['histNames'] = [itm for itm in G2stIO.GetHistogramNames(GPXfile,['PWDR',]) if itm in SeqResult.keys()]
541    G2stIO.SetSeqResult(GPXfile,Histograms,SeqResult)
542    printFile.close()
543    G2fil.G2Print (' Sequential refinement results are in file: '+ospath.splitext(GPXfile)[0]+'.lst')
544    G2fil.G2Print (' ***** Sequential refinement successful *****')
545    return True,'Success'
546
547def RetDistAngle(DisAglCtls,DisAglData):
548    '''Compute and return distances and angles
549
550    :param dict DisAglCtls: contains distance/angle radii usually defined using
551       :func:`GSASIIctrlGUI.DisAglDialog`
552    :param dict DisAglData: contains phase data:
553       Items 'OrigAtoms' and 'TargAtoms' contain the atoms to be used
554       for distance/angle origins and atoms to be used as targets.
555       Item 'SGData' has the space group information (see :ref:`Space Group object<SGData_table>`)
556
557    :returns: AtomLabels,DistArray,AngArray where:
558
559      **AtomLabels** is a dict of atom labels, keys are the atom number
560
561      **DistArray** is a dict keyed by the origin atom number where the value is a list
562      of distance entries. The value for each distance is a list containing:
563
564        0) the target atom number (int);
565        1) the unit cell offsets added to x,y & z (tuple of int values)
566        2) the symmetry operator number (which may be modified to indicate centering and center of symmetry)
567        3) an interatomic distance in A (float)
568        4) an uncertainty on the distance in A or 0.0 (float)
569
570      **AngArray** is a dict keyed by the origin (central) atom number where
571      the value is a list of
572      angle entries. The value for each angle entry consists of three values:
573
574        0) a distance item reference for one neighbor (int)
575        1) a distance item reference for a second neighbor (int)
576        2) a angle, uncertainty pair; the s.u. may be zero (tuple of two floats)
577
578      The AngArray distance reference items refer directly to the index of the items in the
579      DistArray item for the list of distances for the central atom.
580    '''
581    import numpy.ma as ma
582
583    SGData = DisAglData['SGData']
584    Cell = DisAglData['Cell']
585
586    Amat,Bmat = G2lat.cell2AB(Cell[:6])
587    covData = {}
588    if 'covData' in DisAglData:
589        covData = DisAglData['covData']
590        covMatrix = covData['covMatrix']
591        varyList = covData['varyList']
592        pfx = str(DisAglData['pId'])+'::'
593
594    Factor = DisAglCtls['Factors']
595    Radii = dict(zip(DisAglCtls['AtomTypes'],zip(DisAglCtls['BondRadii'],DisAglCtls['AngleRadii'])))
596    indices = (-2,-1,0,1,2)
597    Units = np.array([[h,k,l] for h in indices for k in indices for l in indices])
598    origAtoms = DisAglData['OrigAtoms']
599    targAtoms = DisAglData['TargAtoms']
600    AtomLabels = {}
601    for Oatom in origAtoms:
602        AtomLabels[Oatom[0]] = Oatom[1]
603    for Oatom in targAtoms:
604        AtomLabels[Oatom[0]] = Oatom[1]
605    DistArray = {}
606    AngArray = {}
607    for Oatom in origAtoms:
608        DistArray[Oatom[0]] = []
609        AngArray[Oatom[0]] = []
610        OxyzNames = ''
611        IndBlist = []
612        Dist = []
613        Vect = []
614        VectA = []
615        angles = []
616        for Tatom in targAtoms:
617            Xvcov = []
618            TxyzNames = ''
619            if 'covData' in DisAglData:
620                OxyzNames = [pfx+'dAx:%d'%(Oatom[0]),pfx+'dAy:%d'%(Oatom[0]),pfx+'dAz:%d'%(Oatom[0])]
621                TxyzNames = [pfx+'dAx:%d'%(Tatom[0]),pfx+'dAy:%d'%(Tatom[0]),pfx+'dAz:%d'%(Tatom[0])]
622                Xvcov = G2mth.getVCov(OxyzNames+TxyzNames,varyList,covMatrix)
623            result = G2spc.GenAtom(Tatom[3:6],SGData,False,Move=False)
624            BsumR = (Radii[Oatom[2]][0]+Radii[Tatom[2]][0])*Factor[0]
625            AsumR = (Radii[Oatom[2]][1]+Radii[Tatom[2]][1])*Factor[1]
626            for [Txyz,Top,Tunit,Spn] in result:
627                Dx = (Txyz-np.array(Oatom[3:6]))+Units
628                dx = np.inner(Amat,Dx)
629                dist = ma.masked_less(np.sqrt(np.sum(dx**2,axis=0)),0.5)
630                IndB = ma.nonzero(ma.masked_greater(dist-BsumR,0.))
631                if np.any(IndB):
632                    for indb in IndB:
633                        for i in range(len(indb)):
634                            if str(dx.T[indb][i]) not in IndBlist:
635                                IndBlist.append(str(dx.T[indb][i]))
636                                unit = Units[indb][i]
637                                tunit = (unit[0]+Tunit[0],unit[1]+Tunit[1],unit[2]+Tunit[2])
638                                pdpx = G2mth.getDistDerv(Oatom[3:6],Tatom[3:6],Amat,unit,Top,SGData)
639                                sig = 0.0
640                                if len(Xvcov):
641                                    sig = np.sqrt(np.inner(pdpx,np.inner(pdpx,Xvcov)))
642                                Dist.append([Oatom[0],Tatom[0],tunit,Top,ma.getdata(dist[indb])[i],sig])
643                                if (Dist[-1][-2]-AsumR) <= 0.:
644                                    Vect.append(dx.T[indb][i]/Dist[-1][-2])
645                                    VectA.append([OxyzNames,np.array(Oatom[3:6]),TxyzNames,np.array(Tatom[3:6]),unit,Top])
646                                else:
647                                    Vect.append([0.,0.,0.])
648                                    VectA.append([])
649        for D in Dist:
650            DistArray[Oatom[0]].append(D[1:])
651        Vect = np.array(Vect)
652        angles = np.zeros((len(Vect),len(Vect)))
653        angsig = np.zeros((len(Vect),len(Vect)))
654        for i,veca in enumerate(Vect):
655            if np.any(veca):
656                for j,vecb in enumerate(Vect):
657                    if np.any(vecb):
658                        angles[i][j],angsig[i][j] = G2mth.getAngSig(VectA[i],VectA[j],Amat,SGData,covData)
659                        if i <= j: continue
660                        AngArray[Oatom[0]].append((i,j,
661                            G2mth.getAngSig(VectA[i],VectA[j],Amat,SGData,covData)))
662    return AtomLabels,DistArray,AngArray
663
664def PrintDistAngle(DisAglCtls,DisAglData,out=sys.stdout):
665    '''Print distances and angles
666
667    :param dict DisAglCtls: contains distance/angle radii usually defined using
668       :func:`GSASIIctrlGUI.DisAglDialog`
669    :param dict DisAglData: contains phase data:
670       Items 'OrigAtoms' and 'TargAtoms' contain the atoms to be used
671       for distance/angle origins and atoms to be used as targets.
672       Item 'SGData' has the space group information (see :ref:`Space Group object<SGData_table>`)
673    :param file out: file object for output. Defaults to sys.stdout.
674    '''
675    def MyPrint(s):
676        out.write(s+'\n')
677        # print(s,file=out) # use in Python 3
678
679    def ShowBanner(name):
680        MyPrint(80*'*')
681        MyPrint('   Interatomic Distances and Angles for phase '+name)
682        MyPrint((80*'*')+'\n')
683
684    ShowBanner(DisAglCtls['Name'])
685    SGData = DisAglData['SGData']
686    SGtext,SGtable = G2spc.SGPrint(SGData)
687    for line in SGtext: MyPrint(line)
688    if len(SGtable) > 1:
689        for i,item in enumerate(SGtable[::2]):
690            if 2*i+1 == len(SGtable):
691                line = ' %s'%(item.ljust(30))
692            else:
693                line = ' %s %s'%(item.ljust(30),SGtable[2*i+1].ljust(30))
694            MyPrint(line)
695    else:
696        MyPrint(' ( 1)    %s'%(SGtable[0])) #triclinic case
697    Cell = DisAglData['Cell']
698
699    Amat,Bmat = G2lat.cell2AB(Cell[:6])
700    covData = {}
701    if 'covData' in DisAglData:
702        covData = DisAglData['covData']
703        pfx = str(DisAglData['pId'])+'::'
704        A = G2lat.cell2A(Cell[:6])
705        cellSig = G2stIO.getCellEsd(pfx,SGData,A,covData)
706        names = [' a = ',' b = ',' c = ',' alpha = ',' beta = ',' gamma = ',' Volume = ']
707        valEsd = [G2mth.ValEsd(Cell[i],cellSig[i],True) for i in range(7)]
708        line = '\n Unit cell:'
709        for name,vals in zip(names,valEsd):
710            line += name+vals
711        MyPrint(line)
712    else:
713        MyPrint('\n Unit cell: a = '+('%.5f'%Cell[0])+' b = '+('%.5f'%Cell[1])+' c = '+('%.5f'%Cell[2])+
714            ' alpha = '+('%.3f'%Cell[3])+' beta = '+('%.3f'%Cell[4])+' gamma = '+
715            ('%.3f'%Cell[5])+' Volume = '+('%.3f'%Cell[6]))
716
717    AtomLabels,DistArray,AngArray = RetDistAngle(DisAglCtls,DisAglData)
718    origAtoms = DisAglData['OrigAtoms']
719    for Oatom in origAtoms:
720        i = Oatom[0]
721        Dist = DistArray[i]
722        nDist = len(Dist)
723        angles = np.zeros((nDist,nDist))
724        angsig = np.zeros((nDist,nDist))
725        for k,j,tup in AngArray[i]:
726            angles[k][j],angsig[k][j] = angles[j][k],angsig[j][k] = tup
727        line = ''
728        for i,x in enumerate(Oatom[3:6]):
729            line += ('%12.5f'%x).rstrip('0')
730        MyPrint('\n Distances & angles for '+Oatom[1]+' at '+line.rstrip())
731        MyPrint(80*'*')
732        line = ''
733        for dist in Dist[:-1]:
734            line += '%12s'%(AtomLabels[dist[0]].center(12))
735        MyPrint('  To       cell +(sym. op.)      dist.  '+line.rstrip())
736        for i,dist in enumerate(Dist):
737            line = ''
738            for j,angle in enumerate(angles[i][0:i]):
739                sig = angsig[i][j]
740                if angle:
741                    if sig:
742                        line += '%12s'%(G2mth.ValEsd(angle,sig,True).center(12))
743                    else:
744                        val = '%.3f'%(angle)
745                        line += '%12s'%(val.center(12))
746                else:
747                    line += 12*' '
748            if dist[4]:            #sig exists!
749                val = G2mth.ValEsd(dist[3],dist[4])
750            else:
751                val = '%8.4f'%(dist[3])
752            tunit = '[%2d%2d%2d]'% dist[1]
753            MyPrint(('  %8s%10s+(%4d) %12s'%(AtomLabels[dist[0]].ljust(8),tunit.ljust(10),dist[2],val.center(12)))+line.rstrip())
754
755def DisAglTor(DATData):
756    'Needs a doc string'
757    SGData = DATData['SGData']
758    Cell = DATData['Cell']
759
760    Amat,Bmat = G2lat.cell2AB(Cell[:6])
761    covData = {}
762    pfx = ''
763    if 'covData' in DATData:
764        covData = DATData['covData']
765        pfx = str(DATData['pId'])+'::'
766    Datoms = []
767    Oatoms = []
768    for i,atom in enumerate(DATData['Datoms']):
769        symop = atom[-1].split('+')
770        if len(symop) == 1:
771            symop.append('0,0,0')
772        symop[0] = int(symop[0])
773        symop[1] = eval(symop[1])
774        atom.append(symop)
775        Datoms.append(atom)
776        oatom = DATData['Oatoms'][i]
777        names = ['','','']
778        if pfx:
779            names = [pfx+'dAx:'+str(oatom[0]),pfx+'dAy:'+str(oatom[0]),pfx+'dAz:'+str(oatom[0])]
780        oatom += [names,]
781        Oatoms.append(oatom)
782    atmSeq = [atom[1]+'('+atom[-2]+')' for atom in Datoms]
783    if DATData['Natoms'] == 4:  #torsion
784        Tors,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData)
785        G2fil.G2Print (' Torsion angle for %s atom sequence: %s = %s'%(DATData['Name'],str(atmSeq).replace("'","")[1:-1],G2mth.ValEsd(Tors,sig)))
786        G2fil.G2Print (' NB: Atom sequence determined by selection order')
787        return      # done with torsion
788    elif DATData['Natoms'] == 3:  #angle
789        Ang,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData)
790        G2fil.G2Print (' Angle in %s for atom sequence: %s = %s'%(DATData['Name'],str(atmSeq).replace("'","")[1:-1],G2mth.ValEsd(Ang,sig)))
791        G2fil.G2Print (' NB: Atom sequence determined by selection order')
792    else:   #2 atoms - distance
793        Dist,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData)
794        G2fil.G2Print (' Distance in %s for atom sequence: %s = %s'%(DATData['Name'],str(atmSeq).replace("'","")[1:-1],G2mth.ValEsd(Dist,sig)))
795
796def BestPlane(PlaneData):
797    'Needs a doc string'
798
799    def ShowBanner(name):
800        G2fil.G2Print (80*'*')
801        G2fil.G2Print ('   Best plane result for phase '+name)
802        G2fil.G2Print (80*'*','\n')
803
804    ShowBanner(PlaneData['Name'])
805
806    Cell = PlaneData['Cell']
807    Amat,Bmat = G2lat.cell2AB(Cell[:6])
808    Atoms = PlaneData['Atoms']
809    sumXYZ = np.zeros(3)
810    XYZ = []
811    Natoms = len(Atoms)
812    for atom in Atoms:
813        xyz = np.array(atom[3:6])
814        XYZ.append(xyz)
815        sumXYZ += xyz
816    sumXYZ /= Natoms
817    XYZ = np.array(XYZ)-sumXYZ
818    XYZ = np.inner(Amat,XYZ).T
819    Zmat = np.zeros((3,3))
820    for i,xyz in enumerate(XYZ):
821        Zmat += np.outer(xyz.T,xyz)
822    G2fil.G2Print (' Selected atoms centered at %10.5f %10.5f %10.5f'%(sumXYZ[0],sumXYZ[1],sumXYZ[2]))
823    Evec,Emat = nl.eig(Zmat)
824    Evec = np.sqrt(Evec)/(Natoms-3)
825    Order = np.argsort(Evec)
826    XYZ = np.inner(XYZ,Emat.T).T
827    XYZ = np.array([XYZ[Order[2]],XYZ[Order[1]],XYZ[Order[0]]]).T
828    G2fil.G2Print (' Atoms in Cartesian best plane coordinates:')
829    G2fil.G2Print (' Name         X         Y         Z')
830    for i,xyz in enumerate(XYZ):
831        G2fil.G2Print (' %6s%10.3f%10.3f%10.3f'%(Atoms[i][1].ljust(6),xyz[0],xyz[1],xyz[2]))
832    G2fil.G2Print ('\n Best plane RMS X =%8.3f, Y =%8.3f, Z =%8.3f'%(Evec[Order[2]],Evec[Order[1]],Evec[Order[0]]))
833
834def main():
835    'Called to run a refinement when this module is executed '
836    starttime = time.time()
837    arg = sys.argv
838    if len(arg) > 1:
839        GPXfile = arg[1]
840        if not ospath.exists(GPXfile):
841            G2fil.G2Print ('ERROR - '+GPXfile+" doesn't exist!")
842            exit()
843    else:
844        G2fil.G2Print ('ERROR - missing filename')
845        exit()
846    # TODO: figure out if this is a sequential refinement and call SeqRefine(GPXfile,None)
847    Refine(GPXfile,None)
848    G2fil.G2Print("Done. Execution time {:.2f} sec.".format(time.time()-starttime))
849
850if __name__ == '__main__':
851    GSASIIpath.InvokeDebugOpts()
852    main()
Note: See TracBrowser for help on using the repository browser.