source: trunk/GSASIIstrMain.py @ 4133

Last change on this file since 4133 was 4133, checked in by toby, 2 years ago

seq ref: copynext for phase parms or get from SeqRef? table

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