source: trunk/GSASIIstrMain.py @ 4111

Last change on this file since 4111 was 4111, checked in by toby, 4 years ago

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