source: trunk/GSASIIstrMain.py @ 3437

Last change on this file since 3437 was 3413, checked in by vondreele, 7 years ago

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