source: trunk/GSASIIstrMain.py @ 4370

Last change on this file since 4370 was 4370, checked in by vondreele, 3 years ago

Add new parameter to Data tab - Layer DispSizerThis? is intended to cover multilayered materials of different phases with layers perpendicular to beam - it might have other uses.

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