source: trunk/GSASIImath.py @ 569

Last change on this file since 569 was 569, checked in by vondreele, 11 years ago

change image default to 'Paired'
fix mismatch for 'rotation'

File size: 27.9 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        if 0.0 in Adiag:                #hard singularity in matrix
106            psing = list(np.where(Adiag == 0.)[0])
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                psing = list(np.where(np.diag(nl.gr(Amatlam)[1]) < 1.e-14)[0])
118                return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':psing}]
119            Xvec /= Adiag
120            M2 = func(x0+Xvec,*args)
121            nfev += 1
122            chisq1 = np.sum(M2**2)
123            if chisq1 > chisq0:
124                lam *= 10.
125            else:
126                x0 += Xvec
127                lam /= 10.
128                break
129        if (chisq0-chisq1)/chisq0 < ftol:
130            break
131        icycle += 1
132    M = func(x0,*args)
133    nfev += 1
134    Yvec,Amat = Hess(x0,*args)
135    try:
136        Bmat = nl.inv(Amat)
137        return [x0,Bmat,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':[]}]
138    except nl.LinAlgError:
139        psing = []
140        if maxcyc:
141            psing = list(np.where(np.diag(nl.qr(Amat)[1]) < 1.e-14)[0])
142        return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':psing}]
143   
144def getVCov(varyNames,varyList,covMatrix):
145    vcov = np.zeros((len(varyNames),len(varyNames)))
146    for i1,name1 in enumerate(varyNames):
147        for i2,name2 in enumerate(varyNames):
148            try:
149                vcov[i1][i2] = covMatrix[varyList.index(name1)][varyList.index(name2)]
150            except ValueError:
151                vcov[i1][i2] = 0.0
152    return vcov
153   
154def getDistDerv(Oxyz,Txyz,Amat,Tunit,Top,SGData):
155   
156    def calcDist(Ox,Tx,U,inv,C,M,T,Amat):
157        TxT = inv*(np.inner(M,Tx)+T)+C+U
158        return np.sqrt(np.sum(np.inner(Amat,(TxT-Ox))**2))
159       
160    inv = Top/abs(Top)
161    cent = abs(Top)/100
162    op = abs(Top)%100-1
163    M,T = SGData['SGOps'][op]
164    C = SGData['SGCen'][cent]
165    dx = .00001
166    deriv = np.zeros(6)
167    for i in [0,1,2]:
168        Oxyz[i] += dx
169        d0 = calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)
170        Oxyz[i] -= 2*dx
171        deriv[i] = (calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)-d0)/(2.*dx)
172        Oxyz[i] += dx
173        Txyz[i] += dx
174        d0 = calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)
175        Txyz[i] -= 2*dx
176        deriv[i+3] = (calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)-d0)/(2.*dx)
177        Txyz[i] += dx
178    return deriv
179   
180def getAngSig(VA,VB,Amat,SGData,covData={}):
181   
182    def calcVec(Ox,Tx,U,inv,C,M,T,Amat):
183        TxT = inv*(np.inner(M,Tx)+T)+C
184        TxT = G2spc.MoveToUnitCell(TxT)+U
185        return np.inner(Amat,(TxT-Ox))
186       
187    def calcAngle(Ox,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat):
188        VecA = calcVec(Ox,TxA,unitA,invA,CA,MA,TA,Amat)
189        VecA /= np.sqrt(np.sum(VecA**2))
190        VecB = calcVec(Ox,TxB,unitB,invB,CB,MB,TB,Amat)
191        VecB /= np.sqrt(np.sum(VecB**2))
192        edge = VecB-VecA
193        edge = np.sum(edge**2)
194        angle = (2.-edge)/2.
195        angle = max(angle,-1.)
196        return acosd(angle)
197       
198    OxAN,OxA,TxAN,TxA,unitA,TopA = VA
199    OxBN,OxB,TxBN,TxB,unitB,TopB = VB
200    invA = invB = 1
201    invA = TopA/abs(TopA)
202    invB = TopB/abs(TopB)
203    centA = abs(TopA)/100
204    centB = abs(TopB)/100
205    opA = abs(TopA)%100-1
206    opB = abs(TopB)%100-1
207    MA,TA = SGData['SGOps'][opA]
208    MB,TB = SGData['SGOps'][opB]
209    CA = SGData['SGCen'][centA]
210    CB = SGData['SGCen'][centB]
211    if 'covMatrix' in covData:
212        covMatrix = covData['covMatrix']
213        varyList = covData['varyList']
214        AngVcov = getVCov(OxAN+TxAN+TxBN,varyList,covMatrix)
215        dx = .00001
216        dadx = np.zeros(9)
217        Ang = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
218        for i in [0,1,2]:
219            OxA[i] += dx
220            a0 = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
221            OxA[i] -= 2*dx
222            dadx[i] = (calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)-a0)/dx
223            OxA[i] += dx
224           
225            TxA[i] += dx
226            a0 = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
227            TxA[i] -= 2*dx
228            dadx[i+3] = (calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)-a0)/dx
229            TxA[i] += dx
230           
231            TxB[i] += dx
232            a0 = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
233            TxB[i] -= 2*dx
234            dadx[i+6] = (calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)-a0)/dx
235            TxB[i] += dx
236           
237        sigAng = np.sqrt(np.inner(dadx,np.inner(AngVcov,dadx)))
238        if sigAng < 0.01:
239            sigAng = 0.0
240        return Ang,sigAng
241    else:
242        return calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat),0.0
243
244def GetDistSig(Oatoms,Atoms,Amat,SGData,covData={}):
245
246    def calcDist(Atoms,SyOps,Amat):
247        XYZ = []
248        for i,atom in enumerate(Atoms):
249            Inv,M,T,C,U = SyOps[i]
250            XYZ.append(np.array(atom[1:4]))
251            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
252            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
253        V1 = XYZ[1]-XYZ[0]
254        return np.sqrt(np.sum(V1**2))
255       
256    Inv = []
257    SyOps = []
258    names = []
259    for i,atom in enumerate(Oatoms):
260        names += atom[-1]
261        Op,unit = Atoms[i][-1]
262        inv = Op/abs(Op)
263        m,t = SGData['SGOps'][abs(Op)%100-1]
264        c = SGData['SGCen'][abs(Op)/100]
265        SyOps.append([inv,m,t,c,unit])
266    Dist = calcDist(Oatoms,SyOps,Amat)
267   
268    sig = -0.001
269    if 'covMatrix' in covData:
270        parmNames = []
271        dx = .00001
272        dadx = np.zeros(6)
273        for i in range(6):
274            ia = i/3
275            ix = i%3
276            Oatoms[ia][ix+1] += dx
277            a0 = calcDist(Oatoms,SyOps,Amat)
278            Oatoms[ia][ix+1] -= 2*dx
279            dadx[i] = (calcDist(Oatoms,SyOps,Amat)-a0)/(2.*dx)
280        covMatrix = covData['covMatrix']
281        varyList = covData['varyList']
282        DistVcov = getVCov(names,varyList,covMatrix)
283        sig = np.sqrt(np.inner(dadx,np.inner(DistVcov,dadx)))
284        if sig < 0.001:
285            sig = -0.001
286   
287    return Dist,sig
288
289def GetAngleSig(Oatoms,Atoms,Amat,SGData,covData={}):
290       
291    def calcAngle(Atoms,SyOps,Amat):
292        XYZ = []
293        for i,atom in enumerate(Atoms):
294            Inv,M,T,C,U = SyOps[i]
295            XYZ.append(np.array(atom[1:4]))
296            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
297            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
298        V1 = XYZ[1]-XYZ[0]
299        V1 /= np.sqrt(np.sum(V1**2))
300        V2 = XYZ[1]-XYZ[2]
301        V2 /= np.sqrt(np.sum(V2**2))
302        V3 = V2-V1
303        cang = min(1.,max((2.-np.sum(V3**2))/2.,-1.))
304        return acosd(cang)
305
306    Inv = []
307    SyOps = []
308    names = []
309    for i,atom in enumerate(Oatoms):
310        names += atom[-1]
311        Op,unit = Atoms[i][-1]
312        inv = Op/abs(Op)
313        m,t = SGData['SGOps'][abs(Op)%100-1]
314        c = SGData['SGCen'][abs(Op)/100]
315        SyOps.append([inv,m,t,c,unit])
316    Angle = calcAngle(Oatoms,SyOps,Amat)
317   
318    sig = -0.01
319    if 'covMatrix' in covData:
320        parmNames = []
321        dx = .00001
322        dadx = np.zeros(9)
323        for i in range(9):
324            ia = i/3
325            ix = i%3
326            Oatoms[ia][ix+1] += dx
327            a0 = calcAngle(Oatoms,SyOps,Amat)
328            Oatoms[ia][ix+1] -= 2*dx
329            dadx[i] = (calcAngle(Oatoms,SyOps,Amat)-a0)/(2.*dx)
330        covMatrix = covData['covMatrix']
331        varyList = covData['varyList']
332        AngVcov = getVCov(names,varyList,covMatrix)
333        sig = np.sqrt(np.inner(dadx,np.inner(AngVcov,dadx)))
334        if sig < 0.01:
335            sig = -0.01
336   
337    return Angle,sig
338
339def GetTorsionSig(Oatoms,Atoms,Amat,SGData,covData={}):
340   
341    def calcTorsion(Atoms,SyOps,Amat):
342       
343        XYZ = []
344        for i,atom in enumerate(Atoms):
345            Inv,M,T,C,U = SyOps[i]
346            XYZ.append(np.array(atom[1:4]))
347            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
348            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
349        V1 = XYZ[1]-XYZ[0]
350        V2 = XYZ[2]-XYZ[1]
351        V3 = XYZ[3]-XYZ[2]
352        V1 /= np.sqrt(np.sum(V1**2))
353        V2 /= np.sqrt(np.sum(V2**2))
354        V3 /= np.sqrt(np.sum(V3**2))
355        M = np.array([V1,V2,V3])
356        D = nl.det(M)
357        Ang = 1.0
358        P12 = np.dot(V1,V2)
359        P13 = np.dot(V1,V3)
360        P23 = np.dot(V2,V3)
361        Tors = acosd((P12*P23-P13)/(np.sqrt(1.-P12**2)*np.sqrt(1.-P23**2)))*D/abs(D)
362        return Tors
363           
364    Inv = []
365    SyOps = []
366    names = []
367    for i,atom in enumerate(Oatoms):
368        names += atom[-1]
369        Op,unit = Atoms[i][-1]
370        inv = Op/abs(Op)
371        m,t = SGData['SGOps'][abs(Op)%100-1]
372        c = SGData['SGCen'][abs(Op)/100]
373        SyOps.append([inv,m,t,c,unit])
374    Tors = calcTorsion(Oatoms,SyOps,Amat)
375   
376    sig = -0.01
377    if 'covMatrix' in covData:
378        parmNames = []
379        dx = .00001
380        dadx = np.zeros(12)
381        for i in range(12):
382            ia = i/3
383            ix = i%3
384            Oatoms[ia][ix+1] += dx
385            a0 = calcTorsion(Oatoms,SyOps,Amat)
386            Oatoms[ia][ix+1] -= 2*dx
387            dadx[i] = (calcTorsion(Oatoms,SyOps,Amat)-a0)/(2.*dx)
388        covMatrix = covData['covMatrix']
389        varyList = covData['varyList']
390        TorVcov = getVCov(names,varyList,covMatrix)
391        sig = np.sqrt(np.inner(dadx,np.inner(TorVcov,dadx)))
392        if sig < 0.01:
393            sig = -0.01
394   
395    return Tors,sig
396       
397def GetDATSig(Oatoms,Atoms,Amat,SGData,covData={}):
398   
399    def calcDist(Atoms,SyOps,Amat):
400        XYZ = []
401        for i,atom in enumerate(Atoms):
402            Inv,M,T,C,U = SyOps[i]
403            XYZ.append(np.array(atom[1:4]))
404            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
405            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
406        V1 = XYZ[1]-XYZ[0]
407        return np.sqrt(np.sum(V1**2))
408       
409    def calcAngle(Atoms,SyOps,Amat):
410        XYZ = []
411        for i,atom in enumerate(Atoms):
412            Inv,M,T,C,U = SyOps[i]
413            XYZ.append(np.array(atom[1:4]))
414            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
415            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
416        V1 = XYZ[1]-XYZ[0]
417        V1 /= np.sqrt(np.sum(V1**2))
418        V2 = XYZ[1]-XYZ[2]
419        V2 /= np.sqrt(np.sum(V2**2))
420        V3 = V2-V1
421        cang = min(1.,max((2.-np.sum(V3**2))/2.,-1.))
422        return acosd(cang)
423
424    def calcTorsion(Atoms,SyOps,Amat):
425       
426        XYZ = []
427        for i,atom in enumerate(Atoms):
428            Inv,M,T,C,U = SyOps[i]
429            XYZ.append(np.array(atom[1:4]))
430            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
431            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
432        V1 = XYZ[1]-XYZ[0]
433        V2 = XYZ[2]-XYZ[1]
434        V3 = XYZ[3]-XYZ[2]
435        V1 /= np.sqrt(np.sum(V1**2))
436        V2 /= np.sqrt(np.sum(V2**2))
437        V3 /= np.sqrt(np.sum(V3**2))
438        M = np.array([V1,V2,V3])
439        D = nl.det(M)
440        Ang = 1.0
441        P12 = np.dot(V1,V2)
442        P13 = np.dot(V1,V3)
443        P23 = np.dot(V2,V3)
444        Tors = acosd((P12*P23-P13)/(np.sqrt(1.-P12**2)*np.sqrt(1.-P23**2)))*D/abs(D)
445        return Tors
446           
447    Inv = []
448    SyOps = []
449    names = []
450    for i,atom in enumerate(Oatoms):
451        names += atom[-1]
452        Op,unit = Atoms[i][-1]
453        inv = Op/abs(Op)
454        m,t = SGData['SGOps'][abs(Op)%100-1]
455        c = SGData['SGCen'][abs(Op)/100]
456        SyOps.append([inv,m,t,c,unit])
457    M = len(Oatoms)
458    if M == 2:
459        Val = calcDist(Oatoms,SyOps,Amat)
460    elif M == 3:
461        Val = calcAngle(Oatoms,SyOps,Amat)
462    else:
463        Val = calcTorsion(Oatoms,SyOps,Amat)
464   
465    sigVals = [-0.001,-0.01,-0.01]
466    sig = sigVals[M-3]
467    if 'covMatrix' in covData:
468        parmNames = []
469        dx = .00001
470        N = M*3
471        dadx = np.zeros(N)
472        for i in range(N):
473            ia = i/3
474            ix = i%3
475            Oatoms[ia][ix+1] += dx
476            if M == 2:
477                a0 = calcDist(Oatoms,SyOps,Amat)
478            elif M == 3:
479                a0 = calcAngle(Oatoms,SyOps,Amat)
480            else:
481                a0 = calcTorsion(Oatoms,SyOps,Amat)
482            Oatoms[ia][ix+1] -= 2*dx
483            if M == 2:
484                dadx[i] = (calcDist(Oatoms,SyOps,Amat)-a0)/(2.*dx)               
485            elif M == 3:
486                dadx[i] = (calcAngle(Oatoms,SyOps,Amat)-a0)/(2.*dx)               
487            else:
488                dadx[i] = (calcTorsion(Oatoms,SyOps,Amat)-a0)/(2.*dx)
489        covMatrix = covData['covMatrix']
490        varyList = covData['varyList']
491        Vcov = getVCov(names,varyList,covMatrix)
492        sig = np.sqrt(np.inner(dadx,np.inner(Vcov,dadx)))
493        if sig < sigVals[M-3]:
494            sig = sigVals[M-3]
495   
496    return Val,sig
497       
498   
499def ValEsd(value,esd=0,nTZ=False):                  #NOT complete - don't use
500    # returns value(esd) string; nTZ=True for no trailing zeros
501    # use esd < 0 for level of precision shown e.g. esd=-0.01 gives 2 places beyond decimal
502    #get the 2 significant digits in the esd
503    edig = lambda esd: int(round(10**(math.log10(esd) % 1+1)))
504    #get the number of digits to represent them
505    epl = lambda esd: 2+int(1.545-math.log10(10*edig(esd)))
506   
507    mdec = lambda esd: -int(round(math.log10(abs(esd))))+1
508    ndec = lambda esd: int(1.545-math.log10(abs(esd)))
509    if esd > 0:
510        fmt = '"%.'+str(ndec(esd))+'f(%d)"'
511        return str(fmt%(value,int(round(esd*10**(mdec(esd)))))).strip('"')
512    elif esd < 0:
513         return str(round(value,mdec(esd)-1))
514    else:
515        text = str("%f"%(value))
516        if nTZ:
517            return text.rstrip('0')
518        else:
519            return text
520
521def FourierMap(data,reflData):
522   
523    import scipy.fftpack as fft
524    generalData = data['General']
525    if not generalData['Map']['MapType']:
526        print '**** ERROR - Fourier map not defined'
527        return
528    mapData = generalData['Map']
529    dmin = mapData['Resolution']
530    SGData = generalData['SGData']
531    cell = generalData['Cell'][1:8]       
532    A = G2lat.cell2A(cell[:6])
533    Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
534    Fhkl = np.zeros(shape=2*Hmax,dtype='c16')
535#    Fhkl[0,0,0] = generalData['F000X']
536    time0 = time.time()
537    for ref in reflData:
538        if ref[4] >= dmin:
539            Fosq,Fcsq,ph = ref[8:11]
540            for i,hkl in enumerate(ref[11]):
541                hkl = np.asarray(hkl,dtype='i')
542                dp = 360.*ref[12][i]
543                a = cosd(ph+dp)
544                b = sind(ph+dp)
545                phasep = complex(a,b)
546                phasem = complex(a,-b)
547                if 'Fobs' in mapData['MapType']:
548                    F = np.sqrt(Fosq)
549                    h,k,l = hkl+Hmax
550                    Fhkl[h,k,l] = F*phasep
551                    h,k,l = -hkl+Hmax
552                    Fhkl[h,k,l] = F*phasem
553                elif 'Fcalc' in mapData['MapType']:
554                    F = np.sqrt(Fcsq)
555                    h,k,l = hkl+Hmax
556                    Fhkl[h,k,l] = F*phasep
557                    h,k,l = -hkl+Hmax
558                    Fhkl[h,k,l] = F*phasem
559                elif 'delt-F' in mapData['MapType']:
560                    dF = np.sqrt(Fosq)-np.sqrt(Fcsq)
561                    h,k,l = hkl+Hmax
562                    Fhkl[h,k,l] = dF*phasep
563                    h,k,l = -hkl+Hmax
564                    Fhkl[h,k,l] = dF*phasem
565                elif '2*Fo-Fc' in mapData['MapType']:
566                    F = 2.*np.sqrt(Fosq)-np.sqrt(Fcsq)
567                    h,k,l = hkl+Hmax
568                    Fhkl[h,k,l] = F*phasep
569                    h,k,l = -hkl+Hmax
570                    Fhkl[h,k,l] = F*phasem
571                elif 'Patterson' in mapData['MapType']:
572                    h,k,l = hkl+Hmax
573                    Fhkl[h,k,l] = complex(Fosq,0.)
574                    h,k,l = -hkl+Hmax
575                    Fhkl[h,k,l] = complex(Fosq,0.)
576    Fhkl = fft.fftshift(Fhkl)
577    rho = fft.fftn(Fhkl)/cell[6]
578    print 'Fourier map time: %.4f'%(time.time()-time0),'no. elements: %d'%(Fhkl.size)
579    mapData['rho'] = np.real(rho)
580    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
581    return mapData
582   
583def findOffset(SGData,rho,Fhkl):
584    mapShape = rho.shape
585    hklShape = Fhkl.shape
586    mapHalf = np.array(mapShape)/2
587    Fmax = np.max(np.absolute(Fhkl))
588    hklHalf = np.array(hklShape)/2
589    sortHKL = np.argsort(Fhkl.flatten())
590    Fdict = {}
591    for hkl in sortHKL:
592        HKL = np.unravel_index(hkl,hklShape)
593        F = Fhkl[HKL[0]][HKL[1]][HKL[2]]
594        if F == 0.:
595            break
596        Fdict['%.6f'%(np.absolute(F))] = hkl
597    Flist = np.flipud(np.sort(Fdict.keys()))
598    F = str(1.e6)
599    i = 0
600    while float(F) > 0.5*Fmax:
601        F = Flist[i]
602        hkl = np.unravel_index(Fdict[F],hklShape)
603        iabsnt,mulp,Uniq,phi = G2spc.GenHKLf(list(hkl-hklHalf),SGData)
604        Uniq = np.array(Uniq,dtype='i')
605        Uniq = np.concatenate((Uniq,-Uniq))+hklHalf
606        phi = np.concatenate((phi,-phi))
607        print hkl-hklHalf
608        for j,H in enumerate(Uniq):
609            Fh = Fhkl[H[0],H[1],H[2]]
610            h,k,l = H-hklHalf
611            print '(%3d,%3d,%3d) %5.2f %9.5f'%(h,k,l,phi[j],np.angle(Fh,deg=True))       
612        i += 1
613       
614       
615   
616    return [0,0,0]
617
618def ChargeFlip(data,reflData,pgbar):
619#    import scipy.fftpack as fft
620    import numpy.fft as fft
621    generalData = data['General']
622    mapData = generalData['Map']
623    flipData = generalData['Flip']
624    FFtable = {}
625    if 'None' not in flipData['Norm element']:
626        normElem = flipData['Norm element'].upper()
627        FFs = G2el.GetFormFactorCoeff(normElem.split('+')[0].split('-')[0])
628        for ff in FFs:
629            if ff['Symbol'] == normElem:
630                FFtable.update(ff)
631    dmin = flipData['Resolution']
632    SGData = generalData['SGData']
633    cell = generalData['Cell'][1:8]       
634    A = G2lat.cell2A(cell[:6])
635    Vol = cell[6]
636    Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
637    Ehkl = np.zeros(shape=2*Hmax,dtype='c16')       #2X64bits per complex no.
638    time0 = time.time()
639    for ref in reflData:
640        dsp = ref[4]
641        if dsp >= dmin:
642            ff = 0.1*Vol    #est. no. atoms for ~10A**3/atom
643            if FFtable:
644                SQ = 0.25/dsp**2
645                ff *= G2el.ScatFac(FFtable,SQ)[0]
646            E = np.sqrt(ref[8])/ff
647            ph = ref[10]
648            if SGData['SGInv']:
649                ph = rn.randint(0,1)*180.
650            else:
651                ph = rn.uniform(0.,360.)
652            for i,hkl in enumerate(ref[11]):
653                hkl = np.asarray(hkl,dtype='i')
654                dp = 360.*ref[12][i]
655                a = cosd(ph+dp)
656                b = sind(ph+dp)
657                phasep = complex(a,b)
658                phasem = complex(a,-b)
659                h,k,l = hkl+Hmax
660                Ehkl[h,k,l] = E*phasep
661                h,k,l = -hkl+Hmax       #Friedel pair refl.
662                Ehkl[h,k,l] = E*phasem
663#    Ehkl[Hmax] = 0.00001           #this to preserve F[0,0,0]
664    CEhkl = copy.copy(Ehkl)
665    MEhkl = ma.array(Ehkl,mask=(Ehkl==0.0))
666    Emask = ma.getmask(MEhkl)
667    sumE = np.sum(ma.array(np.absolute(CEhkl),mask=Emask))
668    Ncyc = 0
669    while True:       
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    print 'Charge flip time: %.4f'%(time.time()-time0),'no. elements: %d'%(Ehkl.size)
686    print 'No.cycles = ',Ncyc,'Residual Rcf =%8.3f%s'%(Rcf,'%')
687    CErho = np.real(fft.fftn(fft.fftshift(CEhkl)))
688    roll = findOffset(SGData,CErho,CEhkl)
689    mapData['rho'] = np.roll(np.roll(np.roll(CErho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
690    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
691    return mapData
692   
693def SearchMap(data,keepDup=False):
694   
695    norm = 1./(np.sqrt(3.)*np.sqrt(2.*np.pi)**3)
696   
697    def noDuplicate(xyz,peaks,SGData):                  #be careful where this is used - it's slow
698        equivs = G2spc.GenAtom(xyz,SGData)
699        xyzs = [equiv[0] for equiv in equivs]
700        for x in xyzs:
701            if True in [np.allclose(x,peak,atol=0.002) for peak in peaks]:
702                return False
703        return True
704           
705    def findRoll(rhoMask,mapHalf):
706        indx = np.array(ma.nonzero(rhoMask)).T
707        rhoList = np.array([rho[i,j,k] for i,j,k in indx])
708        rhoMax = np.max(rhoList)
709        return mapHalf-indx[np.argmax(rhoList)]
710       
711    def rhoCalc(parms,rX,rY,rZ,res,SGLaue):
712        Mag,x0,y0,z0,sig = parms
713        if SGLaue in ['3','3m1','31m','6/m','6/mmm']:
714            return norm*Mag*np.exp(-((x0-rX)**2+(y0-rY)**2+(x0-rX)*(y0-rY)+(z0-rZ)**2)/(2.*sig**2))/(sig*res**3)
715        else:
716            return norm*Mag*np.exp(-((x0-rX)**2+(y0-rY)**2+(z0-rZ)**2)/(2.*sig**2))/(sig*res**3)
717       
718    def peakFunc(parms,rX,rY,rZ,rho,res,SGLaue):
719        Mag,x0,y0,z0,sig = parms
720        M = rho-rhoCalc(parms,rX,rY,rZ,res,SGLaue)
721        return M
722       
723    def peakHess(parms,rX,rY,rZ,rho,res,SGLaue):
724        Mag,x0,y0,z0,sig = parms
725        dMdv = np.zeros(([5,]+list(rX.shape)))
726        delt = .01
727        for i in range(5):
728            parms[i] -= delt
729            rhoCm = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
730            parms[i] += 2.*delt
731            rhoCp = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
732            parms[i] -= delt
733            dMdv[i] = (rhoCp-rhoCm)/(2.*delt)
734        rhoC = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
735        Vec = np.sum(np.sum(np.sum(dMdv*(rho-rhoC),axis=3),axis=2),axis=1)
736        dMdv = np.reshape(dMdv,(5,rX.size))
737        Hess = np.inner(dMdv,dMdv)
738       
739        return Vec,Hess
740       
741    generalData = data['General']
742    phaseName = generalData['Name']
743    SGData = generalData['SGData']
744    cell = generalData['Cell'][1:8]       
745    A = G2lat.cell2A(cell[:6])
746    drawingData = data['Drawing']
747    peaks = []
748    mags = []
749    try:
750        mapData = generalData['Map']
751        contLevel = mapData['cutOff']*mapData['rhoMax']/100.
752        rho = copy.copy(mapData['rho'])     #don't mess up original
753        mapHalf = np.array(rho.shape)/2
754        res = mapData['Resolution']
755        incre = 1./np.array(rho.shape)
756        step = max(1.0,1./res)+1
757        steps = np.array(3*[step,])
758    except KeyError:
759        print '**** ERROR - Fourier map not defined'
760        return peaks,mags
761    while True:
762        rhoMask = ma.array(rho,mask=(rho<contLevel))
763        if not ma.count(rhoMask):
764            break
765        rMI = findRoll(rhoMask,mapHalf)
766        rho = np.roll(np.roll(np.roll(rho,rMI[0],axis=0),rMI[1],axis=1),rMI[2],axis=2)
767        rMM = mapHalf-steps
768        rMP = mapHalf+steps+1
769        rhoPeak = rho[rMM[0]:rMP[0],rMM[1]:rMP[1],rMM[2]:rMP[2]]
770        peakInt = np.sum(rhoPeak)*res**3
771        rX,rY,rZ = np.mgrid[rMM[0]:rMP[0],rMM[1]:rMP[1],rMM[2]:rMP[2]]
772        x0 = [peakInt,mapHalf[0],mapHalf[1],mapHalf[2],2.0]          #magnitude, position & width(sig)
773        result = HessianLSQ(peakFunc,x0,Hess=peakHess,
774            args=(rX,rY,rZ,rhoPeak,res,SGData['SGLaue']),ftol=.0001,maxcyc=100)
775        x1 = result[0]
776        if np.any(x1 < 0):
777            break
778        peak = (np.array(x1[1:4])-rMI)*incre
779        if not len(peaks):
780            peaks.append(peak)
781            mags.append(x1[0])
782        else:
783            if keepDup:
784                peaks.append(peak)
785                mags.append(x1[0])
786            elif noDuplicate(peak,peaks,SGData) and x1[0] > 0.:
787                peaks.append(peak)
788                mags.append(x1[0])
789        rho[rMM[0]:rMP[0],rMM[1]:rMP[1],rMM[2]:rMP[2]] = peakFunc(result[0],rX,rY,rZ,rhoPeak,res,SGData['SGLaue'])
790        rho = np.roll(np.roll(np.roll(rho,-rMI[2],axis=2),-rMI[1],axis=1),-rMI[0],axis=0)
791    return np.array(peaks),np.array([mags,]).T
Note: See TracBrowser for help on using the repository browser.