source: trunk/GSASIIstrMain.py @ 3855

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

major speed up with very big sequential fits

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