source: trunk/GSASIIstrMain.py @ 3838

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

fix previous

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