source: trunk/GSASIIstrMain.py @ 3845

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

more sequential timing; doc fixes for scriptable

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