source: trunk/GSASIIstrMain.py @ 3447

Last change on this file since 3447 was 3447, checked in by toby, 3 years ago

fix seq ref constraint numbering bug; fix use if HAP vars in PSvars & equations; complete RepaintHistogramInfo? w/Py3 + misc Py3 fixes; improve initialization after project read

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 37.4 KB
Line 
1# -*- coding: utf-8 -*-
2'''
3*GSASIIstrMain: main structure routine*
4---------------------------------------
5
6'''
7########### SVN repository information ###################
8# $Date: 2018-06-24 00:46:31 +0000 (Sun, 24 Jun 2018) $
9# $Author: toby $
10# $Revision: 3447 $
11# $URL: trunk/GSASIIstrMain.py $
12# $Id: GSASIIstrMain.py 3447 2018-06-24 00:46:31Z 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: 3447 $")
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        hId = Histograms[histogram]['hId']
371        redphaseVary = phaseCheck(phaseVary,Phases,histogram)
372        Histo = {histogram:Histograms[histogram],}
373        hapVary,hapDict,controlDict = G2stIO.GetHistogramPhaseData(Phases,Histo,Print=False)
374        calcControls.update(controlDict)
375        histVary,histDict,controlDict = G2stIO.GetHistogramData(Histo,False)
376        calcControls.update(controlDict)
377        varyList = rbVary+redphaseVary+hapVary+histVary
378#        if not ihst:
379            # save the initial vary list, but without histogram numbers on parameters
380        saveVaryList = varyList[:]
381        for i,item in enumerate(saveVaryList):
382            items = item.split(':')
383            if items[1]:
384                items[1] = ''
385            item = ':'.join(items)
386            saveVaryList[i] = item
387        if not ihst:
388            SeqResult['varyList'] = saveVaryList
389        else:
390            SeqResult['varyList'] = list(set(SeqResult['varyList']+saveVaryList))
391        parmDict = {}
392        parmDict.update(phaseDict)
393        parmDict.update(hapDict)
394        parmDict.update(histDict)
395        if Controls['Copy2Next']:
396            #parmDict.update(NewparmDict) # don't use in case extra entries would cause a problem
397            for parm in NewparmDict:
398                if parm in parmDict:
399                    parmDict[parm] = NewparmDict[parm]
400        G2stIO.GetFprime(calcControls,Histo)
401        # do constraint processing
402        #reload(G2mv) # debug
403        G2mv.InitVars()
404        constrDict,fixedList = G2stIO.GetConstraints(GPXfile)
405        varyListStart = tuple(varyList) # save the original varyList before dependent vars are removed
406        try:
407            groups,parmlist = G2mv.GroupConstraints(constrDict)
408            G2mv.GenerateConstraints(groups,parmlist,varyList,constrDict,fixedList,parmDict,SeqHist=hId)
409#            if GSASIIpath.GetConfigValue('debug'): print("DBG_"+
410#                G2mv.VarRemapShow(varyList,True))
411            constraintInfo = (groups,parmlist,constrDict,fixedList,ihst)
412        except G2mv.ConstraintException:
413            print (' *** ERROR - your constraints are internally inconsistent ***')
414            #errmsg, warnmsg = G2mv.CheckConstraints(varyList,constrDict,fixedList)
415            #print 'Errors',errmsg
416            #if warnmsg: print 'Warnings',warnmsg
417            return False,' Constraint error'
418        #print G2mv.VarRemapShow(varyList)
419        if not ihst:
420            # first histogram to refine against
421            firstVaryList = []
422            for item in varyList:
423                items = item.split(':')
424                if items[1]:
425                    items[1] = ''
426                item = ':'.join(items)
427                firstVaryList.append(item)
428            newVaryList = firstVaryList
429        else:
430            newVaryList = []
431            for item in varyList:
432                items = item.split(':')
433                if items[1]:
434                    items[1] = ''
435                item = ':'.join(items)
436                newVaryList.append(item)
437        if newVaryList != firstVaryList and Controls['Copy2Next']:
438            # variable lists are expected to match between sequential refinements when Copy2Next is on
439            #print '**** ERROR - variable list for this histogram does not match previous'
440            #print '     Copy of variables is not possible'
441            #print '\ncurrent histogram',histogram,'has',len(newVaryList),'variables'
442            combined = list(set(firstVaryList+newVaryList))
443            c = [var for var in combined if var not in newVaryList]
444            p = [var for var in combined if var not in firstVaryList]
445            print('*** Variables change ***')
446            for typ,vars in [('Removed',c),('Added',p)]:
447                line = '  '+typ+': '
448                if vars:
449                    for var in vars:
450                        if len(line) > 70:
451                            print(line)
452                            line = '    '
453                        line += var + ', '
454                else:
455                        line += 'none, '
456                print(line[:-2])
457            firstVaryList = newVaryList
458
459        ifSeq = True
460        printFile.write('\n Refinement results for histogram: %s\n'%histogram)
461        printFile.write(135*'-'+'\n')
462        if True:
463#        try:
464            IfOK,Rvals,result,covMatrix,sig = RefineCore(Controls,Histo,Phases,restraintDict,
465                rigidbodyDict,parmDict,varyList,calcControls,pawleyLookup,ifSeq,printFile,dlg)
466            if PlotFunction:
467                PlotFunction(G2frame,Histo[histogram]['Data'],histogram)
468
469            print ('  wR = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f, last delta chi = %.4f'%(
470                Rvals['Rwp'],Rvals['chisq'],Rvals['GOF']**2,Rvals['DelChi2']))
471            # add the uncertainties into the esd dictionary (sigDict)
472            if not IfOK:
473                print('***** Sequential refinement failed at histogram '+histogram)
474                break
475            sigDict = dict(zip(varyList,sig))
476            # the uncertainties for dependent constrained parms into the esd dict
477            sigDict.update(G2mv.ComputeDepESD(covMatrix,varyList,parmDict))
478
479            # a dict with values & esds for dependent (constrained) parameters - avoid extraneous holds
480            depParmDict = {i:(parmDict[i],sigDict[i]) for i in varyListStart if i in sigDict and i not in varyList}
481            newCellDict = copy.deepcopy(G2stMth.GetNewCellParms(parmDict,varyList))
482            newAtomDict = copy.deepcopy(G2stMth.ApplyXYZshifts(parmDict,varyList))
483            histRefData = {
484                'variables':result[0],'varyList':varyList,'sig':sig,'Rvals':Rvals,
485                'varyListStart':varyListStart,
486                'covMatrix':covMatrix,'title':histogram,'newAtomDict':newAtomDict,
487                'newCellDict':newCellDict,'depParmDict':depParmDict,
488                'constraintInfo':constraintInfo,
489                'parmDict':parmDict}
490            SeqResult[histogram] = histRefData
491            G2stMth.ApplyRBModels(parmDict,Phases,rigidbodyDict,True)
492    #        G2stIO.SetRigidBodyModels(parmDict,sigDict,rigidbodyDict,printFile)
493            G2stIO.SetHistogramPhaseData(parmDict,sigDict,Phases,Histo,None,ifPrint,printFile)
494            G2stIO.SetHistogramData(parmDict,sigDict,Histo,None,ifPrint,printFile)
495            G2stIO.SetUsedHistogramsAndPhases(GPXfile,Histo,Phases,rigidbodyDict,histRefData,makeBack)
496            makeBack = False
497            NewparmDict = {}
498            # make dict of varied parameters in current histogram, renamed to
499            # next histogram, for use in next refinement.
500            if Controls['Copy2Next'] and ihst < len(histNames)-1:
501                hId = Histo[histogram]['hId'] # current histogram
502                nexthId = Histograms[histNames[ihst+1]]['hId']
503                for parm in set(list(varyList)+list(varyListStart)):
504                    items = parm.split(':')
505                    if len(items) < 3: continue
506                    if str(hId) in items[1]:
507                        items[1] = str(nexthId)
508                        newparm = ':'.join(items)
509                        NewparmDict[newparm] = parmDict[parm]
510#        except G2obj.G2Exception(Msg):
511#            printFile.close()
512#            print (' ***** Refinement aborted *****')
513#            return False,Msg.msg
514    SeqResult['histNames'] = [itm for itm in G2stIO.GetHistogramNames(GPXfile,['PWDR',]) if itm in SeqResult.keys()]
515    G2stIO.SetSeqResult(GPXfile,Histograms,SeqResult)
516    printFile.close()
517    print (' Sequential refinement results are in file: '+ospath.splitext(GPXfile)[0]+'.lst')
518    print (' ***** Sequential refinement successful *****')
519    return True,'Success'
520
521def RetDistAngle(DisAglCtls,DisAglData):
522    '''Compute and return distances and angles
523
524    :param dict DisAglCtls: contains distance/angle radii usually defined using
525       :func:`GSASIIctrlGUI.DisAglDialog`
526    :param dict DisAglData: contains phase data:
527       Items 'OrigAtoms' and 'TargAtoms' contain the atoms to be used
528       for distance/angle origins and atoms to be used as targets.
529       Item 'SGData' has the space group information (see :ref:`Space Group object<SGData_table>`)
530
531    :returns: AtomLabels,DistArray,AngArray where:
532
533      **AtomLabels** is a dict of atom labels, keys are the atom number
534
535      **DistArray** is a dict keyed by the origin atom number where the value is a list
536      of distance entries. The value for each distance is a list containing:
537
538        0) the target atom number (int);
539        1) the unit cell offsets added to x,y & z (tuple of int values)
540        2) the symmetry operator number (which may be modified to indicate centering and center of symmetry)
541        3) an interatomic distance in A (float)
542        4) an uncertainty on the distance in A or 0.0 (float)
543
544      **AngArray** is a dict keyed by the origin (central) atom number where
545      the value is a list of
546      angle entries. The value for each angle entry consists of three values:
547
548        0) a distance item reference for one neighbor (int)
549        1) a distance item reference for a second neighbor (int)
550        2) a angle, uncertainty pair; the s.u. may be zero (tuple of two floats)
551
552      The AngArray distance reference items refer directly to the index of the items in the
553      DistArray item for the list of distances for the central atom.
554    '''
555    import numpy.ma as ma
556
557    SGData = DisAglData['SGData']
558    Cell = DisAglData['Cell']
559
560    Amat,Bmat = G2lat.cell2AB(Cell[:6])
561    covData = {}
562    if 'covData' in DisAglData:
563        covData = DisAglData['covData']
564        covMatrix = covData['covMatrix']
565        varyList = covData['varyList']
566        pfx = str(DisAglData['pId'])+'::'
567
568    Factor = DisAglCtls['Factors']
569    Radii = dict(zip(DisAglCtls['AtomTypes'],zip(DisAglCtls['BondRadii'],DisAglCtls['AngleRadii'])))
570    indices = (-2,-1,0,1,2)
571    Units = np.array([[h,k,l] for h in indices for k in indices for l in indices])
572    origAtoms = DisAglData['OrigAtoms']
573    targAtoms = DisAglData['TargAtoms']
574    AtomLabels = {}
575    for Oatom in origAtoms:
576        AtomLabels[Oatom[0]] = Oatom[1]
577    for Oatom in targAtoms:
578        AtomLabels[Oatom[0]] = Oatom[1]
579    DistArray = {}
580    AngArray = {}
581    for Oatom in origAtoms:
582        DistArray[Oatom[0]] = []
583        AngArray[Oatom[0]] = []
584        OxyzNames = ''
585        IndBlist = []
586        Dist = []
587        Vect = []
588        VectA = []
589        angles = []
590        for Tatom in targAtoms:
591            Xvcov = []
592            TxyzNames = ''
593            if 'covData' in DisAglData:
594                OxyzNames = [pfx+'dAx:%d'%(Oatom[0]),pfx+'dAy:%d'%(Oatom[0]),pfx+'dAz:%d'%(Oatom[0])]
595                TxyzNames = [pfx+'dAx:%d'%(Tatom[0]),pfx+'dAy:%d'%(Tatom[0]),pfx+'dAz:%d'%(Tatom[0])]
596                Xvcov = G2mth.getVCov(OxyzNames+TxyzNames,varyList,covMatrix)
597            result = G2spc.GenAtom(Tatom[3:6],SGData,False,Move=False)
598            BsumR = (Radii[Oatom[2]][0]+Radii[Tatom[2]][0])*Factor[0]
599            AsumR = (Radii[Oatom[2]][1]+Radii[Tatom[2]][1])*Factor[1]
600            for [Txyz,Top,Tunit,Spn] in result:
601                Dx = (Txyz-np.array(Oatom[3:6]))+Units
602                dx = np.inner(Amat,Dx)
603                dist = ma.masked_less(np.sqrt(np.sum(dx**2,axis=0)),0.5)
604                IndB = ma.nonzero(ma.masked_greater(dist-BsumR,0.))
605                if np.any(IndB):
606                    for indb in IndB:
607                        for i in range(len(indb)):
608                            if str(dx.T[indb][i]) not in IndBlist:
609                                IndBlist.append(str(dx.T[indb][i]))
610                                unit = Units[indb][i]
611                                tunit = (unit[0]+Tunit[0],unit[1]+Tunit[1],unit[2]+Tunit[2])
612                                pdpx = G2mth.getDistDerv(Oatom[3:6],Tatom[3:6],Amat,unit,Top,SGData)
613                                sig = 0.0
614                                if len(Xvcov):
615                                    sig = np.sqrt(np.inner(pdpx,np.inner(pdpx,Xvcov)))
616                                Dist.append([Oatom[0],Tatom[0],tunit,Top,ma.getdata(dist[indb])[i],sig])
617                                if (Dist[-1][-2]-AsumR) <= 0.:
618                                    Vect.append(dx.T[indb][i]/Dist[-1][-2])
619                                    VectA.append([OxyzNames,np.array(Oatom[3:6]),TxyzNames,np.array(Tatom[3:6]),unit,Top])
620                                else:
621                                    Vect.append([0.,0.,0.])
622                                    VectA.append([])
623        for D in Dist:
624            DistArray[Oatom[0]].append(D[1:])
625        Vect = np.array(Vect)
626        angles = np.zeros((len(Vect),len(Vect)))
627        angsig = np.zeros((len(Vect),len(Vect)))
628        for i,veca in enumerate(Vect):
629            if np.any(veca):
630                for j,vecb in enumerate(Vect):
631                    if np.any(vecb):
632                        angles[i][j],angsig[i][j] = G2mth.getAngSig(VectA[i],VectA[j],Amat,SGData,covData)
633                        if i <= j: continue
634                        AngArray[Oatom[0]].append((i,j,
635                            G2mth.getAngSig(VectA[i],VectA[j],Amat,SGData,covData)))
636    return AtomLabels,DistArray,AngArray
637
638def PrintDistAngle(DisAglCtls,DisAglData,out=sys.stdout):
639    '''Print distances and angles
640
641    :param dict DisAglCtls: contains distance/angle radii usually defined using
642       :func:`GSASIIctrlGUI.DisAglDialog`
643    :param dict DisAglData: contains phase data:
644       Items 'OrigAtoms' and 'TargAtoms' contain the atoms to be used
645       for distance/angle origins and atoms to be used as targets.
646       Item 'SGData' has the space group information (see :ref:`Space Group object<SGData_table>`)
647    :param file out: file object for output. Defaults to sys.stdout.
648    '''
649    def MyPrint(s):
650        out.write(s+'\n')
651        # print(s,file=out) # use in Python 3
652
653    def ShowBanner(name):
654        MyPrint(80*'*')
655        MyPrint('   Interatomic Distances and Angles for phase '+name)
656        MyPrint((80*'*')+'\n')
657
658    ShowBanner(DisAglCtls['Name'])
659    SGData = DisAglData['SGData']
660    SGtext,SGtable = G2spc.SGPrint(SGData)
661    for line in SGtext: MyPrint(line)
662    if len(SGtable) > 1:
663        for i,item in enumerate(SGtable[::2]):
664            if 2*i+1 == len(SGtable):
665                line = ' %s'%(item.ljust(30))
666            else:
667                line = ' %s %s'%(item.ljust(30),SGtable[2*i+1].ljust(30))
668            MyPrint(line)
669    else:
670        MyPrint(' ( 1)    %s'%(SGtable[0])) #triclinic case
671    Cell = DisAglData['Cell']
672
673    Amat,Bmat = G2lat.cell2AB(Cell[:6])
674    covData = {}
675    if 'covData' in DisAglData:
676        covData = DisAglData['covData']
677        pfx = str(DisAglData['pId'])+'::'
678        A = G2lat.cell2A(Cell[:6])
679        cellSig = G2stIO.getCellEsd(pfx,SGData,A,covData)
680        names = [' a = ',' b = ',' c = ',' alpha = ',' beta = ',' gamma = ',' Volume = ']
681        valEsd = [G2mth.ValEsd(Cell[i],cellSig[i],True) for i in range(7)]
682        line = '\n Unit cell:'
683        for name,vals in zip(names,valEsd):
684            line += name+vals
685        MyPrint(line)
686    else:
687        MyPrint('\n Unit cell: a = '+('%.5f'%Cell[0])+' b = '+('%.5f'%Cell[1])+' c = '+('%.5f'%Cell[2])+
688            ' alpha = '+('%.3f'%Cell[3])+' beta = '+('%.3f'%Cell[4])+' gamma = '+
689            ('%.3f'%Cell[5])+' Volume = '+('%.3f'%Cell[6]))
690
691    AtomLabels,DistArray,AngArray = RetDistAngle(DisAglCtls,DisAglData)
692    origAtoms = DisAglData['OrigAtoms']
693    for Oatom in origAtoms:
694        i = Oatom[0]
695        Dist = DistArray[i]
696        nDist = len(Dist)
697        angles = np.zeros((nDist,nDist))
698        angsig = np.zeros((nDist,nDist))
699        for k,j,tup in AngArray[i]:
700            angles[k][j],angsig[k][j] = angles[j][k],angsig[j][k] = tup
701        line = ''
702        for i,x in enumerate(Oatom[3:6]):
703            line += ('%12.5f'%x).rstrip('0')
704        MyPrint('\n Distances & angles for '+Oatom[1]+' at '+line.rstrip())
705        MyPrint(80*'*')
706        line = ''
707        for dist in Dist[:-1]:
708            line += '%12s'%(AtomLabels[dist[0]].center(12))
709        MyPrint('  To       cell +(sym. op.)      dist.  '+line.rstrip())
710        for i,dist in enumerate(Dist):
711            line = ''
712            for j,angle in enumerate(angles[i][0:i]):
713                sig = angsig[i][j]
714                if angle:
715                    if sig:
716                        line += '%12s'%(G2mth.ValEsd(angle,sig,True).center(12))
717                    else:
718                        val = '%.3f'%(angle)
719                        line += '%12s'%(val.center(12))
720                else:
721                    line += 12*' '
722            if dist[4]:            #sig exists!
723                val = G2mth.ValEsd(dist[3],dist[4])
724            else:
725                val = '%8.4f'%(dist[3])
726            tunit = '[%2d%2d%2d]'% dist[1]
727            MyPrint((%8s%10s+(%4d) %12s'%(AtomLabels[dist[0]].ljust(8),tunit.ljust(10),dist[2],val.center(12)))+line.rstrip())
728
729def DisAglTor(DATData):
730    'Needs a doc string'
731    SGData = DATData['SGData']
732    Cell = DATData['Cell']
733
734    Amat,Bmat = G2lat.cell2AB(Cell[:6])
735    covData = {}
736    pfx = ''
737    if 'covData' in DATData:
738        covData = DATData['covData']
739        pfx = str(DATData['pId'])+'::'
740    Datoms = []
741    Oatoms = []
742    for i,atom in enumerate(DATData['Datoms']):
743        symop = atom[-1].split('+')
744        if len(symop) == 1:
745            symop.append('0,0,0')
746        symop[0] = int(symop[0])
747        symop[1] = eval(symop[1])
748        atom.append(symop)
749        Datoms.append(atom)
750        oatom = DATData['Oatoms'][i]
751        names = ['','','']
752        if pfx:
753            names = [pfx+'dAx:'+str(oatom[0]),pfx+'dAy:'+str(oatom[0]),pfx+'dAz:'+str(oatom[0])]
754        oatom += [names,]
755        Oatoms.append(oatom)
756    atmSeq = [atom[1]+'('+atom[-2]+')' for atom in Datoms]
757    if DATData['Natoms'] == 4:  #torsion
758        Tors,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData)
759        print (' Torsion angle for %s atom sequence: %s = %s'%(DATData['Name'],str(atmSeq).replace("'","")[1:-1],G2mth.ValEsd(Tors,sig)))
760        print (' NB: Atom sequence determined by selection order')
761        return      # done with torsion
762    elif DATData['Natoms'] == 3:  #angle
763        Ang,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData)
764        print (' Angle in %s for atom sequence: %s = %s'%(DATData['Name'],str(atmSeq).replace("'","")[1:-1],G2mth.ValEsd(Ang,sig)))
765        print (' NB: Atom sequence determined by selection order')
766    else:   #2 atoms - distance
767        Dist,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData)
768        print (' Distance in %s for atom sequence: %s = %s'%(DATData['Name'],str(atmSeq).replace("'","")[1:-1],G2mth.ValEsd(Dist,sig)))
769
770def BestPlane(PlaneData):
771    'Needs a doc string'
772
773    def ShowBanner(name):
774        print (80*'*')
775        print ('   Best plane result for phase '+name)
776        print (80*'*','\n')
777
778    ShowBanner(PlaneData['Name'])
779
780    Cell = PlaneData['Cell']
781    Amat,Bmat = G2lat.cell2AB(Cell[:6])
782    Atoms = PlaneData['Atoms']
783    sumXYZ = np.zeros(3)
784    XYZ = []
785    Natoms = len(Atoms)
786    for atom in Atoms:
787        xyz = np.array(atom[3:6])
788        XYZ.append(xyz)
789        sumXYZ += xyz
790    sumXYZ /= Natoms
791    XYZ = np.array(XYZ)-sumXYZ
792    XYZ = np.inner(Amat,XYZ).T
793    Zmat = np.zeros((3,3))
794    for i,xyz in enumerate(XYZ):
795        Zmat += np.outer(xyz.T,xyz)
796    print (' Selected atoms centered at %10.5f %10.5f %10.5f'%(sumXYZ[0],sumXYZ[1],sumXYZ[2]))
797    Evec,Emat = nl.eig(Zmat)
798    Evec = np.sqrt(Evec)/(Natoms-3)
799    Order = np.argsort(Evec)
800    XYZ = np.inner(XYZ,Emat.T).T
801    XYZ = np.array([XYZ[Order[2]],XYZ[Order[1]],XYZ[Order[0]]]).T
802    print (' Atoms in Cartesian best plane coordinates:')
803    print (' Name         X         Y         Z')
804    for i,xyz in enumerate(XYZ):
805        print (' %6s%10.3f%10.3f%10.3f'%(Atoms[i][1].ljust(6),xyz[0],xyz[1],xyz[2]))
806    print ('\n Best plane RMS X =%8.3f, Y =%8.3f, Z =%8.3f'%(Evec[Order[2]],Evec[Order[1]],Evec[Order[0]]))
807
808def main():
809    'Called to run a refinement when this module is executed '
810    starttime = time.time()
811    arg = sys.argv
812    if len(arg) > 1:
813        GPXfile = arg[1]
814        if not ospath.exists(GPXfile):
815            print ('ERROR - '+GPXfile+" doesn't exist!")
816            exit()
817    else:
818        print ('ERROR - missing filename')
819        exit()
820    # TODO: figure out if this is a sequential refinement and call SeqRefine(GPXfile,None)
821    Refine(GPXfile,None)
822    print("Done. Execution time {:.2f} sec.".format(time.time()-starttime))
823
824if __name__ == '__main__':
825    GSASIIpath.InvokeDebugOpts()
826    main()
Note: See TracBrowser for help on using the repository browser.