source: trunk/GSASIIstrMain.py @ 4399

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

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