source: trunk/GSASIIstrMain.py @ 4188

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

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