source: trunk/GSASIImath.py @ 578

Last change on this file since 578 was 578, checked in by vondreele, 12 years ago

fix "hard" Hessian singular matrix problems - now deletes parameters
fix Pawley overlap reflection constraint problem
File size: 28.0 KB
Line 
1#GSASIImath - major mathematics routines
2########### SVN repository information ###################
3# $Date: 2012-01-13 11:48:53 -0600 (Fri, 13 Jan 2012) $
4# $Author: vondreele & toby $
5# $Revision: 451 $
6# $URL: https://subversion.xor.aps.anl.gov/pyGSAS/trunk/GSASIImath.py $
7# $Id: GSASIImath.py 451 2012-01-13 17:48:53Z vondreele $
8########### SVN repository information ###################
9import sys
10import os
11import os.path as ospath
12import random as rn
13import numpy as np
14import numpy.linalg as nl
15import numpy.ma as ma
16import cPickle
17import time
18import math
19import copy
20import GSASIIpath
21import GSASIIElem as G2el
22import GSASIIlattice as G2lat
23import GSASIIspc as G2spc
24import scipy.optimize as so
25import scipy.linalg as sl
26
27sind = lambda x: np.sin(x*np.pi/180.)
28cosd = lambda x: np.cos(x*np.pi/180.)
29tand = lambda x: np.tan(x*np.pi/180.)
30asind = lambda x: 180.*np.arcsin(x)/np.pi
31acosd = lambda x: 180.*np.arccos(x)/np.pi
32atan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
33
34def HessianLSQ(func,x0,Hess,args=(),ftol=1.49012e-8,xtol=1.49012e-8, maxcyc=0):
35   
36    """
37    Minimize the sum of squares of a set of equations.
38
39    ::
40   
41                    Nobs
42        x = arg min(sum(func(y)**2,axis=0))
43                    y=0
44
45    Parameters
46    ----------
47    func : callable
48        should take at least one (possibly length N vector) argument and
49        returns M floating point numbers.
50    x0 : ndarray
51        The starting estimate for the minimization of length N
52    Hess : callable
53        A required function or method to compute the weighted vector and Hessian for func.
54        It must be a symmetric NxN array
55    args : tuple
56        Any extra arguments to func are placed in this tuple.
57    ftol : float
58        Relative error desired in the sum of squares.
59    xtol : float
60        Relative error desired in the approximate solution.
61    maxcyc : int
62        The maximum number of cycles of refinement to execute, if -1 refine
63        until other limits are met (ftol, xtol)
64
65    Returns
66    -------
67    x : ndarray
68        The solution (or the result of the last iteration for an unsuccessful
69        call).
70    cov_x : ndarray
71        Uses the fjac and ipvt optional outputs to construct an
72        estimate of the jacobian around the solution.  ``None`` if a
73        singular matrix encountered (indicates very flat curvature in
74        some direction).  This matrix must be multiplied by the
75        residual standard deviation to get the covariance of the
76        parameter estimates -- see curve_fit.
77    infodict : dict
78        a dictionary of optional outputs with the key s::
79
80            - 'fvec' : the function evaluated at the output
81
82
83    Notes
84    -----
85
86    """
87               
88    x0 = np.array(x0, ndmin=1)      #might be redundant?
89    n = len(x0)
90    if type(args) != type(()):
91        args = (args,)
92       
93    icycle = 0
94    One = np.ones((n,n))
95    lam = 0.001
96    lamMax = lam
97    nfev = 0
98    while icycle < maxcyc:
99        lamMax = max(lamMax,lam)
100        M = func(x0,*args)
101        nfev += 1
102        chisq0 = np.sum(M**2)
103        Yvec,Amat = Hess(x0,*args)
104        Adiag = np.sqrt(np.diag(Amat))
105        psing = np.where(np.abs(Adiag) < 1.e-14,True,False)
106        if np.any(psing):                #hard singularity in matrix
107            return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':psing}]
108        Anorm = np.outer(Adiag,Adiag)
109        Yvec /= Adiag
110        Amat /= Anorm       
111        while True:
112            Lam = np.eye(Amat.shape[0])*lam
113            Amatlam = Amat*(One+Lam)
114            try:
115                Xvec = nl.solve(Amatlam,Yvec)
116            except LinAlgError:
117                print 'ouch #1'
118                psing = list(np.where(np.diag(nl.gr(Amatlam)[1]) < 1.e-14)[0])
119                return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':psing}]
120            Xvec /= Adiag
121            M2 = func(x0+Xvec,*args)
122            nfev += 1
123            chisq1 = np.sum(M2**2)
124            if chisq1 > chisq0:
125                lam *= 10.
126            else:
127                x0 += Xvec
128                lam /= 10.
129                break
130        if (chisq0-chisq1)/chisq0 < ftol:
131            break
132        icycle += 1
133    M = func(x0,*args)
134    nfev += 1
135    Yvec,Amat = Hess(x0,*args)
136    try:
137        Bmat = nl.inv(Amat)
138        return [x0,Bmat,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':[]}]
139    except nl.LinAlgError:
140        print 'ouch #2'
141        psing = []
142        if maxcyc:
143            psing = list(np.where(np.diag(nl.qr(Amat)[1]) < 1.e-14)[0])
144        return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':psing}]
145   
146def getVCov(varyNames,varyList,covMatrix):
147    vcov = np.zeros((len(varyNames),len(varyNames)))
148    for i1,name1 in enumerate(varyNames):
149        for i2,name2 in enumerate(varyNames):
150            try:
151                vcov[i1][i2] = covMatrix[varyList.index(name1)][varyList.index(name2)]
152            except ValueError:
153                vcov[i1][i2] = 0.0
154    return vcov
155   
156def getDistDerv(Oxyz,Txyz,Amat,Tunit,Top,SGData):
157   
158    def calcDist(Ox,Tx,U,inv,C,M,T,Amat):
159        TxT = inv*(np.inner(M,Tx)+T)+C+U
160        return np.sqrt(np.sum(np.inner(Amat,(TxT-Ox))**2))
161       
162    inv = Top/abs(Top)
163    cent = abs(Top)/100
164    op = abs(Top)%100-1
165    M,T = SGData['SGOps'][op]
166    C = SGData['SGCen'][cent]
167    dx = .00001
168    deriv = np.zeros(6)
169    for i in [0,1,2]:
170        Oxyz[i] += dx
171        d0 = calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)
172        Oxyz[i] -= 2*dx
173        deriv[i] = (calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)-d0)/(2.*dx)
174        Oxyz[i] += dx
175        Txyz[i] += dx
176        d0 = calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)
177        Txyz[i] -= 2*dx
178        deriv[i+3] = (calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)-d0)/(2.*dx)
179        Txyz[i] += dx
180    return deriv
181   
182def getAngSig(VA,VB,Amat,SGData,covData={}):
183   
184    def calcVec(Ox,Tx,U,inv,C,M,T,Amat):
185        TxT = inv*(np.inner(M,Tx)+T)+C
186        TxT = G2spc.MoveToUnitCell(TxT)+U
187        return np.inner(Amat,(TxT-Ox))
188       
189    def calcAngle(Ox,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat):
190        VecA = calcVec(Ox,TxA,unitA,invA,CA,MA,TA,Amat)
191        VecA /= np.sqrt(np.sum(VecA**2))
192        VecB = calcVec(Ox,TxB,unitB,invB,CB,MB,TB,Amat)
193        VecB /= np.sqrt(np.sum(VecB**2))
194        edge = VecB-VecA
195        edge = np.sum(edge**2)
196        angle = (2.-edge)/2.
197        angle = max(angle,-1.)
198        return acosd(angle)
199       
200    OxAN,OxA,TxAN,TxA,unitA,TopA = VA
201    OxBN,OxB,TxBN,TxB,unitB,TopB = VB
202    invA = invB = 1
203    invA = TopA/abs(TopA)
204    invB = TopB/abs(TopB)
205    centA = abs(TopA)/100
206    centB = abs(TopB)/100
207    opA = abs(TopA)%100-1
208    opB = abs(TopB)%100-1
209    MA,TA = SGData['SGOps'][opA]
210    MB,TB = SGData['SGOps'][opB]
211    CA = SGData['SGCen'][centA]
212    CB = SGData['SGCen'][centB]
213    if 'covMatrix' in covData:
214        covMatrix = covData['covMatrix']
215        varyList = covData['varyList']
216        AngVcov = getVCov(OxAN+TxAN+TxBN,varyList,covMatrix)
217        dx = .00001
218        dadx = np.zeros(9)
219        Ang = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
220        for i in [0,1,2]:
221            OxA[i] += dx
222            a0 = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
223            OxA[i] -= 2*dx
224            dadx[i] = (calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)-a0)/dx
225            OxA[i] += dx
226           
227            TxA[i] += dx
228            a0 = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
229            TxA[i] -= 2*dx
230            dadx[i+3] = (calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)-a0)/dx
231            TxA[i] += dx
232           
233            TxB[i] += dx
234            a0 = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
235            TxB[i] -= 2*dx
236            dadx[i+6] = (calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)-a0)/dx
237            TxB[i] += dx
238           
239        sigAng = np.sqrt(np.inner(dadx,np.inner(AngVcov,dadx)))
240        if sigAng < 0.01:
241            sigAng = 0.0
242        return Ang,sigAng
243    else:
244        return calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat),0.0
245
246def GetDistSig(Oatoms,Atoms,Amat,SGData,covData={}):
247
248    def calcDist(Atoms,SyOps,Amat):
249        XYZ = []
250        for i,atom in enumerate(Atoms):
251            Inv,M,T,C,U = SyOps[i]
252            XYZ.append(np.array(atom[1:4]))
253            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
254            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
255        V1 = XYZ[1]-XYZ[0]
256        return np.sqrt(np.sum(V1**2))
257       
258    Inv = []
259    SyOps = []
260    names = []
261    for i,atom in enumerate(Oatoms):
262        names += atom[-1]
263        Op,unit = Atoms[i][-1]
264        inv = Op/abs(Op)
265        m,t = SGData['SGOps'][abs(Op)%100-1]
266        c = SGData['SGCen'][abs(Op)/100]
267        SyOps.append([inv,m,t,c,unit])
268    Dist = calcDist(Oatoms,SyOps,Amat)
269   
270    sig = -0.001
271    if 'covMatrix' in covData:
272        parmNames = []
273        dx = .00001
274        dadx = np.zeros(6)
275        for i in range(6):
276            ia = i/3
277            ix = i%3
278            Oatoms[ia][ix+1] += dx
279            a0 = calcDist(Oatoms,SyOps,Amat)
280            Oatoms[ia][ix+1] -= 2*dx
281            dadx[i] = (calcDist(Oatoms,SyOps,Amat)-a0)/(2.*dx)
282        covMatrix = covData['covMatrix']
283        varyList = covData['varyList']
284        DistVcov = getVCov(names,varyList,covMatrix)
285        sig = np.sqrt(np.inner(dadx,np.inner(DistVcov,dadx)))
286        if sig < 0.001:
287            sig = -0.001
288   
289    return Dist,sig
290
291def GetAngleSig(Oatoms,Atoms,Amat,SGData,covData={}):
292       
293    def calcAngle(Atoms,SyOps,Amat):
294        XYZ = []
295        for i,atom in enumerate(Atoms):
296            Inv,M,T,C,U = SyOps[i]
297            XYZ.append(np.array(atom[1:4]))
298            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
299            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
300        V1 = XYZ[1]-XYZ[0]
301        V1 /= np.sqrt(np.sum(V1**2))
302        V2 = XYZ[1]-XYZ[2]
303        V2 /= np.sqrt(np.sum(V2**2))
304        V3 = V2-V1
305        cang = min(1.,max((2.-np.sum(V3**2))/2.,-1.))
306        return acosd(cang)
307
308    Inv = []
309    SyOps = []
310    names = []
311    for i,atom in enumerate(Oatoms):
312        names += atom[-1]
313        Op,unit = Atoms[i][-1]
314        inv = Op/abs(Op)
315        m,t = SGData['SGOps'][abs(Op)%100-1]
316        c = SGData['SGCen'][abs(Op)/100]
317        SyOps.append([inv,m,t,c,unit])
318    Angle = calcAngle(Oatoms,SyOps,Amat)
319   
320    sig = -0.01
321    if 'covMatrix' in covData:
322        parmNames = []
323        dx = .00001
324        dadx = np.zeros(9)
325        for i in range(9):
326            ia = i/3
327            ix = i%3
328            Oatoms[ia][ix+1] += dx
329            a0 = calcAngle(Oatoms,SyOps,Amat)
330            Oatoms[ia][ix+1] -= 2*dx
331            dadx[i] = (calcAngle(Oatoms,SyOps,Amat)-a0)/(2.*dx)
332        covMatrix = covData['covMatrix']
333        varyList = covData['varyList']
334        AngVcov = getVCov(names,varyList,covMatrix)
335        sig = np.sqrt(np.inner(dadx,np.inner(AngVcov,dadx)))
336        if sig < 0.01:
337            sig = -0.01
338   
339    return Angle,sig
340
341def GetTorsionSig(Oatoms,Atoms,Amat,SGData,covData={}):
342   
343    def calcTorsion(Atoms,SyOps,Amat):
344       
345        XYZ = []
346        for i,atom in enumerate(Atoms):
347            Inv,M,T,C,U = SyOps[i]
348            XYZ.append(np.array(atom[1:4]))
349            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
350            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
351        V1 = XYZ[1]-XYZ[0]
352        V2 = XYZ[2]-XYZ[1]
353        V3 = XYZ[3]-XYZ[2]
354        V1 /= np.sqrt(np.sum(V1**2))
355        V2 /= np.sqrt(np.sum(V2**2))
356        V3 /= np.sqrt(np.sum(V3**2))
357        M = np.array([V1,V2,V3])
358        D = nl.det(M)
359        Ang = 1.0
360        P12 = np.dot(V1,V2)
361        P13 = np.dot(V1,V3)
362        P23 = np.dot(V2,V3)
363        Tors = acosd((P12*P23-P13)/(np.sqrt(1.-P12**2)*np.sqrt(1.-P23**2)))*D/abs(D)
364        return Tors
365           
366    Inv = []
367    SyOps = []
368    names = []
369    for i,atom in enumerate(Oatoms):
370        names += atom[-1]
371        Op,unit = Atoms[i][-1]
372        inv = Op/abs(Op)
373        m,t = SGData['SGOps'][abs(Op)%100-1]
374        c = SGData['SGCen'][abs(Op)/100]
375        SyOps.append([inv,m,t,c,unit])
376    Tors = calcTorsion(Oatoms,SyOps,Amat)
377   
378    sig = -0.01
379    if 'covMatrix' in covData:
380        parmNames = []
381        dx = .00001
382        dadx = np.zeros(12)
383        for i in range(12):
384            ia = i/3
385            ix = i%3
386            Oatoms[ia][ix+1] += dx
387            a0 = calcTorsion(Oatoms,SyOps,Amat)
388            Oatoms[ia][ix+1] -= 2*dx
389            dadx[i] = (calcTorsion(Oatoms,SyOps,Amat)-a0)/(2.*dx)
390        covMatrix = covData['covMatrix']
391        varyList = covData['varyList']
392        TorVcov = getVCov(names,varyList,covMatrix)
393        sig = np.sqrt(np.inner(dadx,np.inner(TorVcov,dadx)))
394        if sig < 0.01:
395            sig = -0.01
396   
397    return Tors,sig
398       
399def GetDATSig(Oatoms,Atoms,Amat,SGData,covData={}):
400   
401    def calcDist(Atoms,SyOps,Amat):
402        XYZ = []
403        for i,atom in enumerate(Atoms):
404            Inv,M,T,C,U = SyOps[i]
405            XYZ.append(np.array(atom[1:4]))
406            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
407            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
408        V1 = XYZ[1]-XYZ[0]
409        return np.sqrt(np.sum(V1**2))
410       
411    def calcAngle(Atoms,SyOps,Amat):
412        XYZ = []
413        for i,atom in enumerate(Atoms):
414            Inv,M,T,C,U = SyOps[i]
415            XYZ.append(np.array(atom[1:4]))
416            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
417            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
418        V1 = XYZ[1]-XYZ[0]
419        V1 /= np.sqrt(np.sum(V1**2))
420        V2 = XYZ[1]-XYZ[2]
421        V2 /= np.sqrt(np.sum(V2**2))
422        V3 = V2-V1
423        cang = min(1.,max((2.-np.sum(V3**2))/2.,-1.))
424        return acosd(cang)
425
426    def calcTorsion(Atoms,SyOps,Amat):
427       
428        XYZ = []
429        for i,atom in enumerate(Atoms):
430            Inv,M,T,C,U = SyOps[i]
431            XYZ.append(np.array(atom[1:4]))
432            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
433            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
434        V1 = XYZ[1]-XYZ[0]
435        V2 = XYZ[2]-XYZ[1]
436        V3 = XYZ[3]-XYZ[2]
437        V1 /= np.sqrt(np.sum(V1**2))
438        V2 /= np.sqrt(np.sum(V2**2))
439        V3 /= np.sqrt(np.sum(V3**2))
440        M = np.array([V1,V2,V3])
441        D = nl.det(M)
442        Ang = 1.0
443        P12 = np.dot(V1,V2)
444        P13 = np.dot(V1,V3)
445        P23 = np.dot(V2,V3)
446        Tors = acosd((P12*P23-P13)/(np.sqrt(1.-P12**2)*np.sqrt(1.-P23**2)))*D/abs(D)
447        return Tors
448           
449    Inv = []
450    SyOps = []
451    names = []
452    for i,atom in enumerate(Oatoms):
453        names += atom[-1]
454        Op,unit = Atoms[i][-1]
455        inv = Op/abs(Op)
456        m,t = SGData['SGOps'][abs(Op)%100-1]
457        c = SGData['SGCen'][abs(Op)/100]
458        SyOps.append([inv,m,t,c,unit])
459    M = len(Oatoms)
460    if M == 2:
461        Val = calcDist(Oatoms,SyOps,Amat)
462    elif M == 3:
463        Val = calcAngle(Oatoms,SyOps,Amat)
464    else:
465        Val = calcTorsion(Oatoms,SyOps,Amat)
466   
467    sigVals = [-0.001,-0.01,-0.01]
468    sig = sigVals[M-3]
469    if 'covMatrix' in covData:
470        parmNames = []
471        dx = .00001
472        N = M*3
473        dadx = np.zeros(N)
474        for i in range(N):
475            ia = i/3
476            ix = i%3
477            Oatoms[ia][ix+1] += dx
478            if M == 2:
479                a0 = calcDist(Oatoms,SyOps,Amat)
480            elif M == 3:
481                a0 = calcAngle(Oatoms,SyOps,Amat)
482            else:
483                a0 = calcTorsion(Oatoms,SyOps,Amat)
484            Oatoms[ia][ix+1] -= 2*dx
485            if M == 2:
486                dadx[i] = (calcDist(Oatoms,SyOps,Amat)-a0)/(2.*dx)               
487            elif M == 3:
488                dadx[i] = (calcAngle(Oatoms,SyOps,Amat)-a0)/(2.*dx)               
489            else:
490                dadx[i] = (calcTorsion(Oatoms,SyOps,Amat)-a0)/(2.*dx)
491        covMatrix = covData['covMatrix']
492        varyList = covData['varyList']
493        Vcov = getVCov(names,varyList,covMatrix)
494        sig = np.sqrt(np.inner(dadx,np.inner(Vcov,dadx)))
495        if sig < sigVals[M-3]:
496            sig = sigVals[M-3]
497   
498    return Val,sig
499       
500   
501def ValEsd(value,esd=0,nTZ=False):                  #NOT complete - don't use
502    # returns value(esd) string; nTZ=True for no trailing zeros
503    # use esd < 0 for level of precision shown e.g. esd=-0.01 gives 2 places beyond decimal
504    #get the 2 significant digits in the esd
505    edig = lambda esd: int(round(10**(math.log10(esd) % 1+1)))
506    #get the number of digits to represent them
507    epl = lambda esd: 2+int(1.545-math.log10(10*edig(esd)))
508   
509    mdec = lambda esd: -int(round(math.log10(abs(esd))))+1
510    ndec = lambda esd: int(1.545-math.log10(abs(esd)))
511    if esd > 0:
512        fmt = '"%.'+str(ndec(esd))+'f(%d)"'
513        return str(fmt%(value,int(round(esd*10**(mdec(esd)))))).strip('"')
514    elif esd < 0:
515         return str(round(value,mdec(esd)-1))
516    else:
517        text = str("%f"%(value))
518        if nTZ:
519            return text.rstrip('0')
520        else:
521            return text
522
523def FourierMap(data,reflData):
524   
525#    import scipy.fftpack as fft
526    import numpy.fft as fft
527    generalData = data['General']
528    if not generalData['Map']['MapType']:
529        print '**** ERROR - Fourier map not defined'
530        return
531    mapData = generalData['Map']
532    dmin = mapData['Resolution']
533    SGData = generalData['SGData']
534    cell = generalData['Cell'][1:8]       
535    A = G2lat.cell2A(cell[:6])
536    Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
537    Fhkl = np.zeros(shape=2*Hmax,dtype='c16')
538#    Fhkl[0,0,0] = generalData['F000X']
539    time0 = time.time()
540    for ref in reflData:
541        if ref[4] >= dmin:
542            Fosq,Fcsq,ph = ref[8:11]
543            for i,hkl in enumerate(ref[11]):
544                hkl = np.asarray(hkl,dtype='i')
545                dp = 360.*ref[12][i]
546                a = cosd(ph+dp)
547                b = sind(ph+dp)
548                phasep = complex(a,b)
549                phasem = complex(a,-b)
550                if 'Fobs' in mapData['MapType']:
551                    F = np.sqrt(Fosq)
552                    h,k,l = hkl+Hmax
553                    Fhkl[h,k,l] = F*phasep
554                    h,k,l = -hkl+Hmax
555                    Fhkl[h,k,l] = F*phasem
556                elif 'Fcalc' in mapData['MapType']:
557                    F = np.sqrt(Fcsq)
558                    h,k,l = hkl+Hmax
559                    Fhkl[h,k,l] = F*phasep
560                    h,k,l = -hkl+Hmax
561                    Fhkl[h,k,l] = F*phasem
562                elif 'delt-F' in mapData['MapType']:
563                    dF = np.sqrt(Fosq)-np.sqrt(Fcsq)
564                    h,k,l = hkl+Hmax
565                    Fhkl[h,k,l] = dF*phasep
566                    h,k,l = -hkl+Hmax
567                    Fhkl[h,k,l] = dF*phasem
568                elif '2*Fo-Fc' in mapData['MapType']:
569                    F = 2.*np.sqrt(Fosq)-np.sqrt(Fcsq)
570                    h,k,l = hkl+Hmax
571                    Fhkl[h,k,l] = F*phasep
572                    h,k,l = -hkl+Hmax
573                    Fhkl[h,k,l] = F*phasem
574                elif 'Patterson' in mapData['MapType']:
575                    h,k,l = hkl+Hmax
576                    Fhkl[h,k,l] = complex(Fosq,0.)
577                    h,k,l = -hkl+Hmax
578                    Fhkl[h,k,l] = complex(Fosq,0.)
579    Fhkl = fft.fftshift(Fhkl)
580    rho = fft.fftn(Fhkl)/cell[6]
581    print 'Fourier map time: %.4f'%(time.time()-time0),'no. elements: %d'%(Fhkl.size)
582    mapData['rho'] = np.real(rho)
583    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
584    return mapData
585   
586def findOffset(SGData,rho,Fhkl):
587    mapShape = rho.shape
588    hklShape = Fhkl.shape
589    mapHalf = np.array(mapShape)/2
590    Fmax = np.max(np.absolute(Fhkl))
591    hklHalf = np.array(hklShape)/2
592    sortHKL = np.argsort(Fhkl.flatten())
593    Fdict = {}
594    for hkl in sortHKL:
595        HKL = np.unravel_index(hkl,hklShape)
596        F = Fhkl[HKL[0]][HKL[1]][HKL[2]]
597        if F == 0.:
598            break
599        Fdict['%.6f'%(np.absolute(F))] = hkl
600    Flist = np.flipud(np.sort(Fdict.keys()))
601    F = str(1.e6)
602    i = 0
603    while float(F) > 0.5*Fmax:
604        F = Flist[i]
605        hkl = np.unravel_index(Fdict[F],hklShape)
606        iabsnt,mulp,Uniq,phi = G2spc.GenHKLf(list(hkl-hklHalf),SGData)
607        Uniq = np.array(Uniq,dtype='i')+hklHalf
608        print hkl-hklHalf
609        for j,H in enumerate(Uniq):
610            Fh = Fhkl[H[0],H[1],H[2]]
611            h,k,l = H-hklHalf
612            print '(%3d,%3d,%3d) %5.2f %9.5f'%(h,k,l,phi[j],np.angle(Fh,deg=True))       
613        i += 1
614       
615       
616   
617    return [0,0,0]
618
619def ChargeFlip(data,reflData,pgbar):
620#    import scipy.fftpack as fft
621    import numpy.fft as fft
622    generalData = data['General']
623    mapData = generalData['Map']
624    flipData = generalData['Flip']
625    FFtable = {}
626    if 'None' not in flipData['Norm element']:
627        normElem = flipData['Norm element'].upper()
628        FFs = G2el.GetFormFactorCoeff(normElem.split('+')[0].split('-')[0])
629        for ff in FFs:
630            if ff['Symbol'] == normElem:
631                FFtable.update(ff)
632    dmin = flipData['Resolution']
633    SGData = generalData['SGData']
634    cell = generalData['Cell'][1:8]       
635    A = G2lat.cell2A(cell[:6])
636    Vol = cell[6]
637    Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
638    Ehkl = np.zeros(shape=2*Hmax,dtype='c16')       #2X64bits per complex no.
639    time0 = time.time()
640    for ref in reflData:
641        dsp = ref[4]
642        if dsp >= dmin:
643            ff = 0.1*Vol    #est. no. atoms for ~10A**3/atom
644            if FFtable:
645                SQ = 0.25/dsp**2
646                ff *= G2el.ScatFac(FFtable,SQ)[0]
647            E = np.sqrt(ref[8])/ff
648            ph = ref[10]
649#            ph = rn.uniform(0.,360.)
650            for i,hkl in enumerate(ref[11]):
651                hkl = np.asarray(hkl,dtype='i')
652                dp = 360.*ref[12][i]
653                a = cosd(ph+dp)
654                b = sind(ph+dp)
655                phasep = complex(a,b)
656                phasem = complex(a,-b)
657                h,k,l = hkl+Hmax
658                Ehkl[h,k,l] = E*phasep
659                h,k,l = -hkl+Hmax       #Friedel pair refl.
660                Ehkl[h,k,l] = E*phasem
661#    Ehkl[Hmax] = 0.00001           #this to preserve F[0,0,0]
662    CEhkl = copy.copy(Ehkl)
663    MEhkl = ma.array(Ehkl,mask=(Ehkl==0.0))
664    Emask = ma.getmask(MEhkl)
665    sumE = np.sum(ma.array(np.absolute(CEhkl),mask=Emask))
666    Ncyc = 0
667    old = np.seterr(all='raise')
668    while True:       
669        try:
670            CErho = np.real(fft.fftn(fft.fftshift(CEhkl)))*(1.+0j)
671            CEsig = np.std(CErho)
672            CFrho = np.where(np.real(CErho) >= flipData['k-factor']*CEsig,CErho,-CErho)
673            CFhkl = fft.ifftshift(fft.ifftn(CFrho))
674            phase = CFhkl/np.absolute(CFhkl)
675            CEhkl = np.absolute(Ehkl)*phase
676            Ncyc += 1
677            sumCF = np.sum(ma.array(np.absolute(CFhkl),mask=Emask))
678            DEhkl = np.absolute(np.absolute(Ehkl)/sumE-np.absolute(CFhkl)/sumCF)
679            Rcf = min(100.,np.sum(ma.array(DEhkl,mask=Emask)*100.))
680            if Rcf < 5.:
681                break
682            GoOn = pgbar.Update(Rcf,newmsg='%s%8.3f%s\n%s %d'%('Residual Rcf =',Rcf,'%','No.cycles = ',Ncyc))[0]
683            if not GoOn:
684                break
685        except FloatingPointError:
686            Rcf = 100.
687            break
688    np.seterr(**old)
689    print 'Charge flip time: %.4f'%(time.time()-time0),'no. elements: %d'%(Ehkl.size)
690    print 'No.cycles = ',Ncyc,'Residual Rcf =%8.3f%s'%(Rcf,'%')
691    CErho = np.real(fft.fftn(fft.fftshift(CEhkl)))
692    roll = findOffset(SGData,CErho,CEhkl)
693    mapData['Rcf'] = Rcf
694    mapData['rho'] = np.roll(np.roll(np.roll(CErho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
695    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
696    mapData['rollMap'] = [0,0,0]
697    return mapData
698   
699def SearchMap(data,keepDup=False):
700   
701    norm = 1./(np.sqrt(3.)*np.sqrt(2.*np.pi)**3)
702   
703    def noDuplicate(xyz,peaks,SGData):                  #be careful where this is used - it's slow
704        equivs = G2spc.GenAtom(xyz,SGData)
705        xyzs = [equiv[0] for equiv in equivs]
706        for x in xyzs:
707            if True in [np.allclose(x,peak,atol=0.002) for peak in peaks]:
708                return False
709        return True
710           
711    def findRoll(rhoMask,mapHalf):
712        indx = np.array(ma.nonzero(rhoMask)).T
713        rhoList = np.array([rho[i,j,k] for i,j,k in indx])
714        rhoMax = np.max(rhoList)
715        return mapHalf-indx[np.argmax(rhoList)]
716       
717    def rhoCalc(parms,rX,rY,rZ,res,SGLaue):
718        Mag,x0,y0,z0,sig = parms
719        return norm*Mag*np.exp(-((x0-rX)**2+(y0-rY)**2+(z0-rZ)**2)/(2.*sig**2))/(sig*res**3)
720       
721    def peakFunc(parms,rX,rY,rZ,rho,res,SGLaue):
722        Mag,x0,y0,z0,sig = parms
723        M = rho-rhoCalc(parms,rX,rY,rZ,res,SGLaue)
724        return M
725       
726    def peakHess(parms,rX,rY,rZ,rho,res,SGLaue):
727        Mag,x0,y0,z0,sig = parms
728        dMdv = np.zeros(([5,]+list(rX.shape)))
729        delt = .01
730        for i in range(5):
731            parms[i] -= delt
732            rhoCm = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
733            parms[i] += 2.*delt
734            rhoCp = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
735            parms[i] -= delt
736            dMdv[i] = (rhoCp-rhoCm)/(2.*delt)
737        rhoC = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
738        Vec = np.sum(np.sum(np.sum(dMdv*(rho-rhoC),axis=3),axis=2),axis=1)
739        dMdv = np.reshape(dMdv,(5,rX.size))
740        Hess = np.inner(dMdv,dMdv)
741       
742        return Vec,Hess
743       
744    generalData = data['General']
745    phaseName = generalData['Name']
746    SGData = generalData['SGData']
747    cell = generalData['Cell'][1:8]       
748    A = G2lat.cell2A(cell[:6])
749    drawingData = data['Drawing']
750    peaks = []
751    mags = []
752    try:
753        mapData = generalData['Map']
754        contLevel = mapData['cutOff']*mapData['rhoMax']/100.
755        rho = copy.copy(mapData['rho'])     #don't mess up original
756        mapHalf = np.array(rho.shape)/2
757        res = mapData['Resolution']
758        incre = 1./np.array(rho.shape)
759        step = max(1.0,1./res)+1
760        steps = np.array(3*[step,])
761    except KeyError:
762        print '**** ERROR - Fourier map not defined'
763        return peaks,mags
764    while True:
765        rhoMask = ma.array(rho,mask=(rho<contLevel))
766        if not ma.count(rhoMask):
767            break
768        rMI = findRoll(rhoMask,mapHalf)
769        rho = np.roll(np.roll(np.roll(rho,rMI[0],axis=0),rMI[1],axis=1),rMI[2],axis=2)
770        rMM = mapHalf-steps
771        rMP = mapHalf+steps+1
772        rhoPeak = rho[rMM[0]:rMP[0],rMM[1]:rMP[1],rMM[2]:rMP[2]]
773        peakInt = np.sum(rhoPeak)*res**3
774        rX,rY,rZ = np.mgrid[rMM[0]:rMP[0],rMM[1]:rMP[1],rMM[2]:rMP[2]]
775        x0 = [peakInt,mapHalf[0],mapHalf[1],mapHalf[2],2.0]          #magnitude, position & width(sig)
776        result = HessianLSQ(peakFunc,x0,Hess=peakHess,
777            args=(rX,rY,rZ,rhoPeak,res,SGData['SGLaue']),ftol=.0001,maxcyc=100)
778        x1 = result[0]
779        if np.any(x1 < 0):
780            break
781        peak = (np.array(x1[1:4])-rMI)*incre
782        if not len(peaks):
783            peaks.append(peak)
784            mags.append(x1[0])
785        else:
786            if keepDup:
787                peaks.append(peak)
788                mags.append(x1[0])
789            elif noDuplicate(peak,peaks,SGData) and x1[0] > 0.:
790                peaks.append(peak)
791                mags.append(x1[0])
792            if len(peaks) > 100:
793                break
794        rho[rMM[0]:rMP[0],rMM[1]:rMP[1],rMM[2]:rMP[2]] = peakFunc(result[0],rX,rY,rZ,rhoPeak,res,SGData['SGLaue'])
795        rho = np.roll(np.roll(np.roll(rho,-rMI[2],axis=2),-rMI[1],axis=1),-rMI[0],axis=0)
796    return np.array(peaks),np.array([mags,]).T
Note: See TracBrowser for help on using the repository browser.