source: trunk/GSASIIstrMain.py @ 3913

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

Update Rietveld cycle-by-cycle

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