source: trunk/GSASIImath.py @ 657

Last change on this file since 657 was 657, checked in by vondreele, 10 years ago

fixups of findOffset & mapsearch

File size: 30.2 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    if SGData['SpGrp'] == 'P 1':
614        return [0,0,0]   
615    mapShape = rho.shape
616    steps = np.array(mapShape)
617    hklShape = Fhkl.shape
618    mapHalf = np.array(mapShape)/2
619    Fmax = np.max(np.absolute(Fhkl))
620    hklHalf = np.array(hklShape)/2
621    sortHKL = np.argsort(Fhkl.flatten())
622    Fdict = {}
623    for hkl in sortHKL:
624        HKL = np.unravel_index(hkl,hklShape)
625        F = Fhkl[HKL[0]][HKL[1]][HKL[2]]
626        if F == 0.:
627            break
628        Fdict['%.6f'%(np.absolute(F))] = hkl
629    Flist = np.flipud(np.sort(Fdict.keys()))
630    F = str(1.e6)
631    i = 0
632    DH = []
633    Dphi = []
634    while i < 20:
635        F = Flist[i]
636        hkl = np.unravel_index(Fdict[F],hklShape)
637        iabsnt,mulp,Uniq,Phi = G2spc.GenHKLf(list(hkl-hklHalf),SGData)
638        Uniq = np.array(Uniq,dtype='i')
639        Phi = np.array(Phi)
640        Uniq = np.concatenate((Uniq,-Uniq))+hklHalf         # put in Friedel pairs & make as index to Farray
641        Phi = np.concatenate((Phi,-Phi))                      # and their phase shifts
642        Fh0 = Fhkl[hkl[0],hkl[1],hkl[2]]
643        ang0 = np.angle(Fh0,deg=True)/360.
644        for j,H in enumerate(Uniq[1:]):
645            ang = (np.angle(Fhkl[H[0],H[1],H[2]],deg=True)/360.-Phi[j+1])
646            dH = H-hkl
647            dang = ang-ang0
648            if np.any(np.abs(dH)-8 > 0):    #keep low order DHs
649                continue
650            DH.append(dH)
651            Dphi.append((dang+0.5) % 1.0)
652        i += 1
653    DH = np.array(DH)
654    Dphi = np.array(Dphi)
655    X,Y,Z = np.mgrid[0:1:1./steps[0],0:1:1./steps[1],0:1:1./steps[2]]
656    XYZ = np.array(zip(X.flatten(),Y.flatten(),Z.flatten()))
657    Mmap = np.reshape(np.sum(((np.dot(XYZ,DH.T)+.5)%1.-Dphi)**2,axis=1),newshape=steps)
658    chisq = np.min(Mmap)
659    DX = -np.array(np.unravel_index(np.argmin(Mmap),Mmap.shape))
660    print ' map offset chi**2: %.3f, map offset: %d %d %d'%(chisq,DX[0],DX[1],DX[2])
661    return DX
662   
663def ChargeFlip(data,reflData,pgbar):
664    generalData = data['General']
665    mapData = generalData['Map']
666    flipData = generalData['Flip']
667    FFtable = {}
668    if 'None' not in flipData['Norm element']:
669        normElem = flipData['Norm element'].upper()
670        FFs = G2el.GetFormFactorCoeff(normElem.split('+')[0].split('-')[0])
671        for ff in FFs:
672            if ff['Symbol'] == normElem:
673                FFtable.update(ff)
674    dmin = flipData['Resolution']
675    SGData = generalData['SGData']
676    cell = generalData['Cell'][1:8]       
677    A = G2lat.cell2A(cell[:6])
678    Vol = cell[6]
679    Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
680    Ehkl = np.zeros(shape=2*Hmax,dtype='c16')       #2X64bits per complex no.
681    time0 = time.time()
682    for ref in reflData:
683        dsp = ref[4]
684        if dsp >= dmin:
685            ff = 0.1*Vol    #est. no. atoms for ~10A**3/atom
686            if FFtable:
687                SQ = 0.25/dsp**2
688                ff *= G2el.ScatFac(FFtable,SQ)[0]
689            if ref[8] > 0.:
690                E = np.sqrt(ref[8])/ff
691            else:
692                E = 0.
693            ph = ref[10]
694            ph = rn.uniform(0.,360.)
695            for i,hkl in enumerate(ref[11]):
696                hkl = np.asarray(hkl,dtype='i')
697                dp = 360.*ref[12][i]
698                a = cosd(ph+dp)
699                b = sind(ph+dp)
700                phasep = complex(a,b)
701                phasem = complex(a,-b)
702                h,k,l = hkl+Hmax
703                Ehkl[h,k,l] = E*phasep
704                h,k,l = -hkl+Hmax       #Friedel pair refl.
705                Ehkl[h,k,l] = E*phasem
706#    Ehkl[Hmax] = 0.00001           #this to preserve F[0,0,0]
707    CEhkl = copy.copy(Ehkl)
708    MEhkl = ma.array(Ehkl,mask=(Ehkl==0.0))
709    Emask = ma.getmask(MEhkl)
710    sumE = np.sum(ma.array(np.absolute(CEhkl),mask=Emask))
711    Ncyc = 0
712    old = np.seterr(all='raise')
713    while True:       
714        try:
715            CErho = np.real(fft.fftn(fft.fftshift(CEhkl)))*(1.+0j)
716            CEsig = np.std(CErho)
717            CFrho = np.where(np.real(CErho) >= flipData['k-factor']*CEsig,CErho,-CErho)
718            CFhkl = fft.ifftshift(fft.ifftn(CFrho))
719            phase = CFhkl/np.absolute(CFhkl)
720            CEhkl = np.absolute(Ehkl)*phase
721            Ncyc += 1
722            sumCF = np.sum(ma.array(np.absolute(CFhkl),mask=Emask))
723            DEhkl = np.absolute(np.absolute(Ehkl)/sumE-np.absolute(CFhkl)/sumCF)
724            Rcf = min(100.,np.sum(ma.array(DEhkl,mask=Emask)*100.))
725            if Rcf < 5.:
726                break
727            GoOn = pgbar.Update(Rcf,newmsg='%s%8.3f%s\n%s %d'%('Residual Rcf =',Rcf,'%','No.cycles = ',Ncyc))[0]
728            if not GoOn or Ncyc > 10000:
729                break
730        except FloatingPointError:
731            Rcf = 100.
732            break
733    np.seterr(**old)
734    print ' Charge flip time: %.4f'%(time.time()-time0),'no. elements: %d'%(Ehkl.size)
735    print ' No.cycles = ',Ncyc,'Residual Rcf =%8.3f%s'%(Rcf,'%')+' Map size:',CErho.shape
736    CErho = np.real(fft.fftn(fft.fftshift(CEhkl)))
737    roll = findOffset(SGData,CErho,CEhkl)
738       
739    mapData['Rcf'] = Rcf
740    mapData['rho'] = np.roll(np.roll(np.roll(CErho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
741    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
742    mapData['rollMap'] = [0,0,0]
743    return mapData
744   
745def SearchMap(data,keepDup=False,Pgbar=None):
746   
747    norm = 1./(np.sqrt(3.)*np.sqrt(2.*np.pi)**3)
748   
749    def noDuplicate(xyz,peaks,SGData,incre):                  #be careful where this is used - it's slow
750        equivs = G2spc.GenAtom(xyz,SGData)
751        xyzs = [equiv[0] for equiv in equivs]
752        for x in xyzs:
753            if True in [np.allclose(x,peak,atol=0.02) for peak in peaks]:
754                return False
755        return True
756           
757    def findRoll(rhoMask,mapHalf):
758        indx = np.array(ma.nonzero(rhoMask)).T
759        rhoList = np.array([rho[i,j,k] for i,j,k in indx])
760        rhoMax = np.max(rhoList)
761        return mapHalf-indx[np.argmax(rhoList)]
762       
763    def rhoCalc(parms,rX,rY,rZ,res,SGLaue):
764        Mag,x0,y0,z0,sig = parms
765#        if SGLaue in ['3','3m1','31m','6/m','6/mmm']:
766#            return norm*Mag*np.exp(-((x0-rX)**2+(y0-rY)**2+(x0-rX)*(y0-rY)+(z0-rZ)**2)/(2.*sig**2))/(sig*res**3)
767#        else:
768#            return norm*Mag*np.exp(-((x0-rX)**2+(y0-rY)**2+(z0-rZ)**2)/(2.*sig**2))/(sig*res**3)
769        return norm*Mag*np.exp(-((x0-rX)**2+(y0-rY)**2+(z0-rZ)**2)/(2.*sig**2))/(sig*res**3)
770       
771    def peakFunc(parms,rX,rY,rZ,rho,res,SGLaue):
772        Mag,x0,y0,z0,sig = parms
773        M = rho-rhoCalc(parms,rX,rY,rZ,res,SGLaue)
774        return M
775       
776    def peakHess(parms,rX,rY,rZ,rho,res,SGLaue):
777        Mag,x0,y0,z0,sig = parms
778        dMdv = np.zeros(([5,]+list(rX.shape)))
779        delt = .01
780        for i in range(5):
781            parms[i] -= delt
782            rhoCm = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
783            parms[i] += 2.*delt
784            rhoCp = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
785            parms[i] -= delt
786            dMdv[i] = (rhoCp-rhoCm)/(2.*delt)
787        rhoC = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
788        Vec = np.sum(np.sum(np.sum(dMdv*(rho-rhoC),axis=3),axis=2),axis=1)
789        dMdv = np.reshape(dMdv,(5,rX.size))
790        Hess = np.inner(dMdv,dMdv)
791       
792        return Vec,Hess
793       
794    generalData = data['General']
795    phaseName = generalData['Name']
796    SGData = generalData['SGData']
797    cell = generalData['Cell'][1:8]       
798    A = G2lat.cell2A(cell[:6])
799    drawingData = data['Drawing']
800    peaks = []
801    mags = []
802    try:
803        mapData = generalData['Map']
804        contLevel = mapData['cutOff']*mapData['rhoMax']/100.
805        rho = copy.copy(mapData['rho'])     #don't mess up original
806        mapHalf = np.array(rho.shape)/2
807        res = mapData['Resolution']
808        incre = 1./np.array(rho.shape)
809        step = max(1.0,1./res)+1
810        steps = np.array(3*[step,])
811    except KeyError:
812        print '**** ERROR - Fourier map not defined'
813        return peaks,mags
814    while True:
815        rhoMask = ma.array(rho,mask=(rho<contLevel))
816        if not ma.count(rhoMask):
817            break
818        rMI = findRoll(rhoMask,mapHalf)
819        rho = np.roll(np.roll(np.roll(rho,rMI[0],axis=0),rMI[1],axis=1),rMI[2],axis=2)
820        rMM = mapHalf-steps
821        rMP = mapHalf+steps+1
822        rhoPeak = rho[rMM[0]:rMP[0],rMM[1]:rMP[1],rMM[2]:rMP[2]]
823        peakInt = np.sum(rhoPeak)*res**3
824        rX,rY,rZ = np.mgrid[rMM[0]:rMP[0],rMM[1]:rMP[1],rMM[2]:rMP[2]]
825        x0 = [peakInt,mapHalf[0],mapHalf[1],mapHalf[2],2.0]          #magnitude, position & width(sig)
826        result = HessianLSQ(peakFunc,x0,Hess=peakHess,
827            args=(rX,rY,rZ,rhoPeak,res,SGData['SGLaue']),ftol=.01,maxcyc=10)
828        x1 = result[0]
829        if np.any(x1 < 0):
830            break
831        peak = (np.array(x1[1:4])-rMI)*incre
832        if not len(peaks):
833            peaks.append(peak)
834            mags.append(x1[0])
835        else:
836            if keepDup:
837                peaks.append(peak)
838                mags.append(x1[0])
839            elif noDuplicate(peak,peaks,SGData,incre) and x1[0] > 0.:
840                peaks.append(peak)
841                mags.append(x1[0])
842            GoOn = Pgbar.Update(len(peaks),newmsg='%s%d'%('No. Peaks found =',len(peaks)))[0]
843            if not GoOn or len(peaks) > 300:
844                break
845        rho[rMM[0]:rMP[0],rMM[1]:rMP[1],rMM[2]:rMP[2]] = peakFunc(x1,rX,rY,rZ,rhoPeak,res,SGData['SGLaue'])
846        rho = np.roll(np.roll(np.roll(rho,-rMI[2],axis=2),-rMI[1],axis=1),-rMI[0],axis=0)
847    return np.array(peaks),np.array([mags,]).T
Note: See TracBrowser for help on using the repository browser.