source: trunk/GSASIIstrMain.py @ 3732

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

more fixes for constraint refactor

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