source: trunk/GSASIIstrMain.py @ 3764

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

fixes to nonFourier modulation & derivatives - zigzag works

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