source: trunk/GSASIImath.py @ 655

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

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