source: trunk/GSASIIstrMain.py @ 3839

Last change on this file since 3839 was 3839, checked in by vondreele, 4 years ago

add some timing for plot patterns - controlled by showTiming
more inn. mag. st.fctr.stuff
fix timing prints

  • 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 21:56:01 +0000 (Wed, 06 Mar 2019) $
9# $Author: vondreele $
10# $Revision: 3839 $
11# $URL: trunk/GSASIIstrMain.py $
12# $Id: GSASIIstrMain.py 3839 2019-03-06 21:56:01Z vondreele $
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: 3839 $")
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('debug'): 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('debug'):
418            t2 = time.time()
419            print("#1 debug time {:.2f} sec.".format(t2-t1))
420            t1 = t2
421        G2mv.InitVars()
422        constrDict,fixedList = G2stIO.GetConstraints(GPXfile)
423        if GSASIIpath.GetConfigValue('debug'):
424            t2 = time.time()
425            print("#2 debug 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('debug'):
487            t2 = time.time()
488            print("#3 debug 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 PlotFunction:
495                PlotFunction(G2frame,Histo[histogram]['Data'],histogram)
496
497            print ('  wR = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f, last delta chi = %.4f'%(
498                Rvals['Rwp'],Rvals['chisq'],Rvals['GOF']**2,Rvals['DelChi2']))
499            # add the uncertainties into the esd dictionary (sigDict)
500            if not IfOK:
501                print('***** Sequential refinement failed at histogram '+histogram)
502                break
503            sigDict = dict(zip(varyList,sig))
504            # the uncertainties for dependent constrained parms into the esd dict
505            sigDict.update(G2mv.ComputeDepESD(covMatrix,varyList,parmDict))
506
507            # a dict with values & esds for dependent (constrained) parameters - avoid extraneous holds
508            depParmDict = {i:(parmDict[i],sigDict[i]) for i in varyListStart if i in sigDict and i not in varyList}
509            newCellDict = copy.deepcopy(G2stMth.GetNewCellParms(parmDict,varyList))
510            newAtomDict = copy.deepcopy(G2stMth.ApplyXYZshifts(parmDict,varyList))
511            histRefData = {
512                'variables':result[0],'varyList':varyList,'sig':sig,'Rvals':Rvals,
513                'varyListStart':varyListStart,
514                'covMatrix':covMatrix,'title':histogram,'newAtomDict':newAtomDict,
515                'newCellDict':newCellDict,'depParmDict':depParmDict,
516                'constraintInfo':constraintInfo,
517                'parmDict':parmDict}
518            SeqResult[histogram] = histRefData
519            if GSASIIpath.GetConfigValue('debug'):
520                t2 = time.time()
521                print("#4 debug time {:.2f} sec.".format(t2-t1))
522                t1 = t2
523            G2stMth.ApplyRBModels(parmDict,Phases,rigidbodyDict,True)
524    #        G2stIO.SetRigidBodyModels(parmDict,sigDict,rigidbodyDict,printFile)
525            G2stIO.SetHistogramPhaseData(parmDict,sigDict,Phases,Histo,None,ifPrint,printFile)
526            G2stIO.SetHistogramData(parmDict,sigDict,Histo,None,ifPrint,printFile)
527            G2stIO.SetUsedHistogramsAndPhases(GPXfile,Histo,Phases,rigidbodyDict,histRefData,makeBack)
528            if GSASIIpath.GetConfigValue('debug'):
529                t2 = time.time()
530                print("#5 debug time {:.2f} sec.".format(t2-t1))
531                t1 = t2
532            makeBack = False
533            NewparmDict = {}
534            # make dict of varied parameters in current histogram, renamed to
535            # next histogram, for use in next refinement.
536            if Controls['Copy2Next'] and ihst < len(histNames)-1:
537                hId = Histo[histogram]['hId'] # current histogram
538                nexthId = Histograms[histNames[ihst+1]]['hId']
539                for parm in set(list(varyList)+list(varyListStart)):
540                    items = parm.split(':')
541                    if len(items) < 3: continue
542                    if str(hId) in items[1]:
543                        items[1] = str(nexthId)
544                        newparm = ':'.join(items)
545                        NewparmDict[newparm] = parmDict[parm]
546#        except G2obj.G2Exception(Msg):
547#            printFile.close()
548#            print (' ***** Refinement aborted *****')
549#            return False,Msg.msg
550        if GSASIIpath.GetConfigValue('debug'):
551            t2 = time.time()
552            print("#6 debug time {:.2f} sec.".format(t2-t1))
553            t1 = t2
554    SeqResult['histNames'] = [itm for itm in G2stIO.GetHistogramNames(GPXfile,['PWDR',]) if itm in SeqResult.keys()]
555    G2stIO.SetSeqResult(GPXfile,Histograms,SeqResult)
556    printFile.close()
557    print (' Sequential refinement results are in file: '+ospath.splitext(GPXfile)[0]+'.lst')
558    print (' ***** Sequential refinement successful *****')
559    return True,'Success'
560
561def RetDistAngle(DisAglCtls,DisAglData):
562    '''Compute and return distances and angles
563
564    :param dict DisAglCtls: contains distance/angle radii usually defined using
565       :func:`GSASIIctrlGUI.DisAglDialog`
566    :param dict DisAglData: contains phase data:
567       Items 'OrigAtoms' and 'TargAtoms' contain the atoms to be used
568       for distance/angle origins and atoms to be used as targets.
569       Item 'SGData' has the space group information (see :ref:`Space Group object<SGData_table>`)
570
571    :returns: AtomLabels,DistArray,AngArray where:
572
573      **AtomLabels** is a dict of atom labels, keys are the atom number
574
575      **DistArray** is a dict keyed by the origin atom number where the value is a list
576      of distance entries. The value for each distance is a list containing:
577
578        0) the target atom number (int);
579        1) the unit cell offsets added to x,y & z (tuple of int values)
580        2) the symmetry operator number (which may be modified to indicate centering and center of symmetry)
581        3) an interatomic distance in A (float)
582        4) an uncertainty on the distance in A or 0.0 (float)
583
584      **AngArray** is a dict keyed by the origin (central) atom number where
585      the value is a list of
586      angle entries. The value for each angle entry consists of three values:
587
588        0) a distance item reference for one neighbor (int)
589        1) a distance item reference for a second neighbor (int)
590        2) a angle, uncertainty pair; the s.u. may be zero (tuple of two floats)
591
592      The AngArray distance reference items refer directly to the index of the items in the
593      DistArray item for the list of distances for the central atom.
594    '''
595    import numpy.ma as ma
596
597    SGData = DisAglData['SGData']
598    Cell = DisAglData['Cell']
599
600    Amat,Bmat = G2lat.cell2AB(Cell[:6])
601    covData = {}
602    if 'covData' in DisAglData:
603        covData = DisAglData['covData']
604        covMatrix = covData['covMatrix']
605        varyList = covData['varyList']
606        pfx = str(DisAglData['pId'])+'::'
607
608    Factor = DisAglCtls['Factors']
609    Radii = dict(zip(DisAglCtls['AtomTypes'],zip(DisAglCtls['BondRadii'],DisAglCtls['AngleRadii'])))
610    indices = (-2,-1,0,1,2)
611    Units = np.array([[h,k,l] for h in indices for k in indices for l in indices])
612    origAtoms = DisAglData['OrigAtoms']
613    targAtoms = DisAglData['TargAtoms']
614    AtomLabels = {}
615    for Oatom in origAtoms:
616        AtomLabels[Oatom[0]] = Oatom[1]
617    for Oatom in targAtoms:
618        AtomLabels[Oatom[0]] = Oatom[1]
619    DistArray = {}
620    AngArray = {}
621    for Oatom in origAtoms:
622        DistArray[Oatom[0]] = []
623        AngArray[Oatom[0]] = []
624        OxyzNames = ''
625        IndBlist = []
626        Dist = []
627        Vect = []
628        VectA = []
629        angles = []
630        for Tatom in targAtoms:
631            Xvcov = []
632            TxyzNames = ''
633            if 'covData' in DisAglData:
634                OxyzNames = [pfx+'dAx:%d'%(Oatom[0]),pfx+'dAy:%d'%(Oatom[0]),pfx+'dAz:%d'%(Oatom[0])]
635                TxyzNames = [pfx+'dAx:%d'%(Tatom[0]),pfx+'dAy:%d'%(Tatom[0]),pfx+'dAz:%d'%(Tatom[0])]
636                Xvcov = G2mth.getVCov(OxyzNames+TxyzNames,varyList,covMatrix)
637            result = G2spc.GenAtom(Tatom[3:6],SGData,False,Move=False)
638            BsumR = (Radii[Oatom[2]][0]+Radii[Tatom[2]][0])*Factor[0]
639            AsumR = (Radii[Oatom[2]][1]+Radii[Tatom[2]][1])*Factor[1]
640            for [Txyz,Top,Tunit,Spn] in result:
641                Dx = (Txyz-np.array(Oatom[3:6]))+Units
642                dx = np.inner(Amat,Dx)
643                dist = ma.masked_less(np.sqrt(np.sum(dx**2,axis=0)),0.5)
644                IndB = ma.nonzero(ma.masked_greater(dist-BsumR,0.))
645                if np.any(IndB):
646                    for indb in IndB:
647                        for i in range(len(indb)):
648                            if str(dx.T[indb][i]) not in IndBlist:
649                                IndBlist.append(str(dx.T[indb][i]))
650                                unit = Units[indb][i]
651                                tunit = (unit[0]+Tunit[0],unit[1]+Tunit[1],unit[2]+Tunit[2])
652                                pdpx = G2mth.getDistDerv(Oatom[3:6],Tatom[3:6],Amat,unit,Top,SGData)
653                                sig = 0.0
654                                if len(Xvcov):
655                                    sig = np.sqrt(np.inner(pdpx,np.inner(pdpx,Xvcov)))
656                                Dist.append([Oatom[0],Tatom[0],tunit,Top,ma.getdata(dist[indb])[i],sig])
657                                if (Dist[-1][-2]-AsumR) <= 0.:
658                                    Vect.append(dx.T[indb][i]/Dist[-1][-2])
659                                    VectA.append([OxyzNames,np.array(Oatom[3:6]),TxyzNames,np.array(Tatom[3:6]),unit,Top])
660                                else:
661                                    Vect.append([0.,0.,0.])
662                                    VectA.append([])
663        for D in Dist:
664            DistArray[Oatom[0]].append(D[1:])
665        Vect = np.array(Vect)
666        angles = np.zeros((len(Vect),len(Vect)))
667        angsig = np.zeros((len(Vect),len(Vect)))
668        for i,veca in enumerate(Vect):
669            if np.any(veca):
670                for j,vecb in enumerate(Vect):
671                    if np.any(vecb):
672                        angles[i][j],angsig[i][j] = G2mth.getAngSig(VectA[i],VectA[j],Amat,SGData,covData)
673                        if i <= j: continue
674                        AngArray[Oatom[0]].append((i,j,
675                            G2mth.getAngSig(VectA[i],VectA[j],Amat,SGData,covData)))
676    return AtomLabels,DistArray,AngArray
677
678def PrintDistAngle(DisAglCtls,DisAglData,out=sys.stdout):
679    '''Print distances and angles
680
681    :param dict DisAglCtls: contains distance/angle radii usually defined using
682       :func:`GSASIIctrlGUI.DisAglDialog`
683    :param dict DisAglData: contains phase data:
684       Items 'OrigAtoms' and 'TargAtoms' contain the atoms to be used
685       for distance/angle origins and atoms to be used as targets.
686       Item 'SGData' has the space group information (see :ref:`Space Group object<SGData_table>`)
687    :param file out: file object for output. Defaults to sys.stdout.
688    '''
689    def MyPrint(s):
690        out.write(s+'\n')
691        # print(s,file=out) # use in Python 3
692
693    def ShowBanner(name):
694        MyPrint(80*'*')
695        MyPrint('   Interatomic Distances and Angles for phase '+name)
696        MyPrint((80*'*')+'\n')
697
698    ShowBanner(DisAglCtls['Name'])
699    SGData = DisAglData['SGData']
700    SGtext,SGtable = G2spc.SGPrint(SGData)
701    for line in SGtext: MyPrint(line)
702    if len(SGtable) > 1:
703        for i,item in enumerate(SGtable[::2]):
704            if 2*i+1 == len(SGtable):
705                line = ' %s'%(item.ljust(30))
706            else:
707                line = ' %s %s'%(item.ljust(30),SGtable[2*i+1].ljust(30))
708            MyPrint(line)
709    else:
710        MyPrint(' ( 1)    %s'%(SGtable[0])) #triclinic case
711    Cell = DisAglData['Cell']
712
713    Amat,Bmat = G2lat.cell2AB(Cell[:6])
714    covData = {}
715    if 'covData' in DisAglData:
716        covData = DisAglData['covData']
717        pfx = str(DisAglData['pId'])+'::'
718        A = G2lat.cell2A(Cell[:6])
719        cellSig = G2stIO.getCellEsd(pfx,SGData,A,covData)
720        names = [' a = ',' b = ',' c = ',' alpha = ',' beta = ',' gamma = ',' Volume = ']
721        valEsd = [G2mth.ValEsd(Cell[i],cellSig[i],True) for i in range(7)]
722        line = '\n Unit cell:'
723        for name,vals in zip(names,valEsd):
724            line += name+vals
725        MyPrint(line)
726    else:
727        MyPrint('\n Unit cell: a = '+('%.5f'%Cell[0])+' b = '+('%.5f'%Cell[1])+' c = '+('%.5f'%Cell[2])+
728            ' alpha = '+('%.3f'%Cell[3])+' beta = '+('%.3f'%Cell[4])+' gamma = '+
729            ('%.3f'%Cell[5])+' Volume = '+('%.3f'%Cell[6]))
730
731    AtomLabels,DistArray,AngArray = RetDistAngle(DisAglCtls,DisAglData)
732    origAtoms = DisAglData['OrigAtoms']
733    for Oatom in origAtoms:
734        i = Oatom[0]
735        Dist = DistArray[i]
736        nDist = len(Dist)
737        angles = np.zeros((nDist,nDist))
738        angsig = np.zeros((nDist,nDist))
739        for k,j,tup in AngArray[i]:
740            angles[k][j],angsig[k][j] = angles[j][k],angsig[j][k] = tup
741        line = ''
742        for i,x in enumerate(Oatom[3:6]):
743            line += ('%12.5f'%x).rstrip('0')
744        MyPrint('\n Distances & angles for '+Oatom[1]+' at '+line.rstrip())
745        MyPrint(80*'*')
746        line = ''
747        for dist in Dist[:-1]:
748            line += '%12s'%(AtomLabels[dist[0]].center(12))
749        MyPrint('  To       cell +(sym. op.)      dist.  '+line.rstrip())
750        for i,dist in enumerate(Dist):
751            line = ''
752            for j,angle in enumerate(angles[i][0:i]):
753                sig = angsig[i][j]
754                if angle:
755                    if sig:
756                        line += '%12s'%(G2mth.ValEsd(angle,sig,True).center(12))
757                    else:
758                        val = '%.3f'%(angle)
759                        line += '%12s'%(val.center(12))
760                else:
761                    line += 12*' '
762            if dist[4]:            #sig exists!
763                val = G2mth.ValEsd(dist[3],dist[4])
764            else:
765                val = '%8.4f'%(dist[3])
766            tunit = '[%2d%2d%2d]'% dist[1]
767            MyPrint((%8s%10s+(%4d) %12s'%(AtomLabels[dist[0]].ljust(8),tunit.ljust(10),dist[2],val.center(12)))+line.rstrip())
768
769def DisAglTor(DATData):
770    'Needs a doc string'
771    SGData = DATData['SGData']
772    Cell = DATData['Cell']
773
774    Amat,Bmat = G2lat.cell2AB(Cell[:6])
775    covData = {}
776    pfx = ''
777    if 'covData' in DATData:
778        covData = DATData['covData']
779        pfx = str(DATData['pId'])+'::'
780    Datoms = []
781    Oatoms = []
782    for i,atom in enumerate(DATData['Datoms']):
783        symop = atom[-1].split('+')
784        if len(symop) == 1:
785            symop.append('0,0,0')
786        symop[0] = int(symop[0])
787        symop[1] = eval(symop[1])
788        atom.append(symop)
789        Datoms.append(atom)
790        oatom = DATData['Oatoms'][i]
791        names = ['','','']
792        if pfx:
793            names = [pfx+'dAx:'+str(oatom[0]),pfx+'dAy:'+str(oatom[0]),pfx+'dAz:'+str(oatom[0])]
794        oatom += [names,]
795        Oatoms.append(oatom)
796    atmSeq = [atom[1]+'('+atom[-2]+')' for atom in Datoms]
797    if DATData['Natoms'] == 4:  #torsion
798        Tors,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData)
799        print (' Torsion angle for %s atom sequence: %s = %s'%(DATData['Name'],str(atmSeq).replace("'","")[1:-1],G2mth.ValEsd(Tors,sig)))
800        print (' NB: Atom sequence determined by selection order')
801        return      # done with torsion
802    elif DATData['Natoms'] == 3:  #angle
803        Ang,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData)
804        print (' Angle in %s for atom sequence: %s = %s'%(DATData['Name'],str(atmSeq).replace("'","")[1:-1],G2mth.ValEsd(Ang,sig)))
805        print (' NB: Atom sequence determined by selection order')
806    else:   #2 atoms - distance
807        Dist,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData)
808        print (' Distance in %s for atom sequence: %s = %s'%(DATData['Name'],str(atmSeq).replace("'","")[1:-1],G2mth.ValEsd(Dist,sig)))
809
810def BestPlane(PlaneData):
811    'Needs a doc string'
812
813    def ShowBanner(name):
814        print (80*'*')
815        print ('   Best plane result for phase '+name)
816        print (80*'*','\n')
817
818    ShowBanner(PlaneData['Name'])
819
820    Cell = PlaneData['Cell']
821    Amat,Bmat = G2lat.cell2AB(Cell[:6])
822    Atoms = PlaneData['Atoms']
823    sumXYZ = np.zeros(3)
824    XYZ = []
825    Natoms = len(Atoms)
826    for atom in Atoms:
827        xyz = np.array(atom[3:6])
828        XYZ.append(xyz)
829        sumXYZ += xyz
830    sumXYZ /= Natoms
831    XYZ = np.array(XYZ)-sumXYZ
832    XYZ = np.inner(Amat,XYZ).T
833    Zmat = np.zeros((3,3))
834    for i,xyz in enumerate(XYZ):
835        Zmat += np.outer(xyz.T,xyz)
836    print (' Selected atoms centered at %10.5f %10.5f %10.5f'%(sumXYZ[0],sumXYZ[1],sumXYZ[2]))
837    Evec,Emat = nl.eig(Zmat)
838    Evec = np.sqrt(Evec)/(Natoms-3)
839    Order = np.argsort(Evec)
840    XYZ = np.inner(XYZ,Emat.T).T
841    XYZ = np.array([XYZ[Order[2]],XYZ[Order[1]],XYZ[Order[0]]]).T
842    print (' Atoms in Cartesian best plane coordinates:')
843    print (' Name         X         Y         Z')
844    for i,xyz in enumerate(XYZ):
845        print (' %6s%10.3f%10.3f%10.3f'%(Atoms[i][1].ljust(6),xyz[0],xyz[1],xyz[2]))
846    print ('\n Best plane RMS X =%8.3f, Y =%8.3f, Z =%8.3f'%(Evec[Order[2]],Evec[Order[1]],Evec[Order[0]]))
847
848def main():
849    'Called to run a refinement when this module is executed '
850    starttime = time.time()
851    arg = sys.argv
852    if len(arg) > 1:
853        GPXfile = arg[1]
854        if not ospath.exists(GPXfile):
855            print ('ERROR - '+GPXfile+" doesn't exist!")
856            exit()
857    else:
858        print ('ERROR - missing filename')
859        exit()
860    # TODO: figure out if this is a sequential refinement and call SeqRefine(GPXfile,None)
861    Refine(GPXfile,None)
862    print("Done. Execution time {:.2f} sec.".format(time.time()-starttime))
863
864if __name__ == '__main__':
865    GSASIIpath.InvokeDebugOpts()
866    main()
Note: See TracBrowser for help on using the repository browser.