source: trunk/GSASIImath.py @ 658

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

use steps/4 or /6 for fourier & charge flipping
more mods to findOffset
remove a stray debug print

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