source: trunk/GSASIImath.py @ 486

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

edit help & fix references to it
begin implementation of Fourier calcs.
File size: 17.2 KB
Line 
1#GSASIImath - major mathematics routines
2########### SVN repository information ###################
3# $Date: 2012-01-13 11:48:53 -0600 (Fri, 13 Jan 2012) $
4# $Author: vondreele & toby $
5# $Revision: 451 $
6# $URL: https://subversion.xor.aps.anl.gov/pyGSAS/trunk/GSASIImath.py $
7# $Id: GSASIImath.py 451 2012-01-13 17:48:53Z vondreele $
8########### SVN repository information ###################
9import sys
10import os
11import os.path as ospath
12import numpy as np
13import numpy.linalg as nl
14import cPickle
15import time
16import math
17import GSASIIpath
18import GSASIIspc as G2spc
19import scipy.optimize as so
20import scipy.linalg as sl
21
22sind = lambda x: np.sin(x*np.pi/180.)
23cosd = lambda x: np.cos(x*np.pi/180.)
24tand = lambda x: np.tan(x*np.pi/180.)
25asind = lambda x: 180.*np.arcsin(x)/np.pi
26acosd = lambda x: 180.*np.arccos(x)/np.pi
27atan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
28
29def HessianLSQ(func,x0,Hess,args=(),ftol=1.49012e-8,xtol=1.49012e-8, maxcyc=0):
30   
31    """
32    Minimize the sum of squares of a set of equations.
33
34    ::
35   
36                    Nobs
37        x = arg min(sum(func(y)**2,axis=0))
38                    y=0
39
40    Parameters
41    ----------
42    func : callable
43        should take at least one (possibly length N vector) argument and
44        returns M floating point numbers.
45    x0 : ndarray
46        The starting estimate for the minimization of length N
47    Hess : callable
48        A required function or method to compute the weighted vector and Hessian for func.
49        It must be a symmetric NxN array
50    args : tuple
51        Any extra arguments to func are placed in this tuple.
52    ftol : float
53        Relative error desired in the sum of squares.
54    xtol : float
55        Relative error desired in the approximate solution.
56    maxcyc : int
57        The maximum number of cycles of refinement to execute, if -1 refine
58        until other limits are met (ftol, xtol)
59
60    Returns
61    -------
62    x : ndarray
63        The solution (or the result of the last iteration for an unsuccessful
64        call).
65    cov_x : ndarray
66        Uses the fjac and ipvt optional outputs to construct an
67        estimate of the jacobian around the solution.  ``None`` if a
68        singular matrix encountered (indicates very flat curvature in
69        some direction).  This matrix must be multiplied by the
70        residual standard deviation to get the covariance of the
71        parameter estimates -- see curve_fit.
72    infodict : dict
73        a dictionary of optional outputs with the key s::
74
75            - 'fvec' : the function evaluated at the output
76
77
78    Notes
79    -----
80
81    """
82               
83    x0 = np.array(x0, ndmin=1)      #might be redundant?
84    n = len(x0)
85    if type(args) != type(()):
86        args = (args,)
87       
88    icycle = 0
89    One = np.ones((n,n))
90    lam = 0.001
91    lamMax = lam
92    nfev = 0
93    while icycle <= maxcyc:
94        lamMax = max(lamMax,lam)
95        M = func(x0,*args)
96        nfev += 1
97        chisq0 = np.sum(M**2)
98        Yvec,Amat = Hess(x0,*args)
99        Adiag = np.sqrt(np.diag(Amat))
100        if 0.0 in Adiag:                #hard singularity in matrix
101            psing = list(np.where(Adiag == 0.)[0])
102            return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':psing}]
103        Anorm = np.outer(Adiag,Adiag)
104        Yvec /= Adiag
105        Amat /= Anorm       
106        while True:
107            Lam = np.eye(Amat.shape[0])*lam
108            Amatlam = Amat*(One+Lam)
109            try:
110                Xvec = nl.solve(Amatlam,Yvec)
111            except LinAlgError:
112                psing = list(np.where(np.diag(nl.gr(Amatlam)[1]) < 1.e-14)[0])
113                return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':psing}]
114            Xvec /= Adiag
115            M2 = func(x0+Xvec,*args)
116            nfev += 1
117            chisq1 = np.sum(M2**2)
118            if chisq1 > chisq0:
119                lam *= 10.
120            else:
121                x0 += Xvec
122                lam /= 10.
123                break
124        if (chisq0-chisq1)/chisq0 < ftol:
125            break
126        icycle += 1
127    M = func(x0,*args)
128    nfev += 1
129    Yvec,Amat = Hess(x0,*args)
130    try:
131        Bmat = nl.inv(Amat)
132        return [x0,Bmat,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':[]}]
133    except LinAlgError:
134        psing = list(np.where(np.diag(nl.gr(Amat)[1]) < 1.e-14)[0])
135        return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':psing}]
136
137def calcFouriermap():
138    print 'Calculate Fourier map'
139
140   
141def getVCov(varyNames,varyList,covMatrix):
142    vcov = np.zeros((len(varyNames),len(varyNames)))
143    for i1,name1 in enumerate(varyNames):
144        for i2,name2 in enumerate(varyNames):
145            try:
146                vcov[i1][i2] = covMatrix[varyList.index(name1)][varyList.index(name2)]
147            except ValueError:
148                vcov[i1][i2] = 0.0
149    return vcov
150   
151def getDistDerv(Oxyz,Txyz,Amat,Tunit,Top,SGData):
152   
153    def calcDist(Ox,Tx,U,inv,C,M,T,Amat):
154        TxT = inv*(np.inner(M,Tx)+T)+C+U
155        return np.sqrt(np.sum(np.inner(Amat,(TxT-Ox))**2))
156       
157    inv = Top/abs(Top)
158    cent = abs(Top)/100
159    op = abs(Top)%100-1
160    M,T = SGData['SGOps'][op]
161    C = SGData['SGCen'][cent]
162    dx = .00001
163    deriv = np.zeros(6)
164    for i in [0,1,2]:
165        Oxyz[i] += dx
166        d0 = calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)
167        Oxyz[i] -= 2*dx
168        deriv[i] = (calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)-d0)/(2.*dx)
169        Oxyz[i] += dx
170        Txyz[i] += dx
171        d0 = calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)
172        Txyz[i] -= 2*dx
173        deriv[i+3] = (calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)-d0)/(2.*dx)
174        Txyz[i] += dx
175    return deriv
176   
177def getAngSig(VA,VB,Amat,SGData,covData={}):
178   
179    def calcVec(Ox,Tx,U,inv,C,M,T,Amat):
180        TxT = inv*(np.inner(M,Tx)+T)+C
181        TxT = G2spc.MoveToUnitCell(TxT)+U
182        return np.inner(Amat,(TxT-Ox))
183       
184    def calcAngle(Ox,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat):
185        VecA = calcVec(Ox,TxA,unitA,invA,CA,MA,TA,Amat)
186        VecA /= np.sqrt(np.sum(VecA**2))
187        VecB = calcVec(Ox,TxB,unitB,invB,CB,MB,TB,Amat)
188        VecB /= np.sqrt(np.sum(VecB**2))
189        edge = VecB-VecA
190        edge = np.sum(edge**2)
191        angle = (2.-edge)/2.
192        angle = max(angle,-1.)
193        return acosd(angle)
194       
195    OxAN,OxA,TxAN,TxA,unitA,TopA = VA
196    OxBN,OxB,TxBN,TxB,unitB,TopB = VB
197    invA = invB = 1
198    invA = TopA/abs(TopA)
199    invB = TopB/abs(TopB)
200    centA = abs(TopA)/100
201    centB = abs(TopB)/100
202    opA = abs(TopA)%100-1
203    opB = abs(TopB)%100-1
204    MA,TA = SGData['SGOps'][opA]
205    MB,TB = SGData['SGOps'][opB]
206    CA = SGData['SGCen'][centA]
207    CB = SGData['SGCen'][centB]
208    if 'covMatrix' in covData:
209        covMatrix = covData['covMatrix']
210        varyList = covData['varyList']
211        AngVcov = getVCov(OxAN+TxAN+TxBN,varyList,covMatrix)
212        dx = .00001
213        dadx = np.zeros(9)
214        Ang = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
215        for i in [0,1,2]:
216            OxA[i] += dx
217            a0 = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
218            OxA[i] -= 2*dx
219            dadx[i] = (calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)-a0)/dx
220            OxA[i] += dx
221           
222            TxA[i] += dx
223            a0 = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
224            TxA[i] -= 2*dx
225            dadx[i+3] = (calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)-a0)/dx
226            TxA[i] += dx
227           
228            TxB[i] += dx
229            a0 = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
230            TxB[i] -= 2*dx
231            dadx[i+6] = (calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)-a0)/dx
232            TxB[i] += dx
233           
234        sigAng = np.sqrt(np.inner(dadx,np.inner(AngVcov,dadx)))
235        if sigAng < 0.01:
236            sigAng = 0.0
237        return Ang,sigAng
238    else:
239        return calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat),0.0
240
241def GetDistSig(Oatoms,Atoms,Amat,SGData,covData={}):
242
243    def calcDist(Atoms,SyOps,Amat):
244        XYZ = []
245        for i,atom in enumerate(Atoms):
246            Inv,M,T,C,U = SyOps[i]
247            XYZ.append(np.array(atom[1:4]))
248            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
249            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
250        V1 = XYZ[1]-XYZ[0]
251        return np.sqrt(np.sum(V1**2))
252       
253    Inv = []
254    SyOps = []
255    names = []
256    for i,atom in enumerate(Oatoms):
257        names += atom[-1]
258        Op,unit = Atoms[i][-1]
259        inv = Op/abs(Op)
260        m,t = SGData['SGOps'][abs(Op)%100-1]
261        c = SGData['SGCen'][abs(Op)/100]
262        SyOps.append([inv,m,t,c,unit])
263    Dist = calcDist(Oatoms,SyOps,Amat)
264   
265    sig = -0.001
266    if 'covMatrix' in covData:
267        parmNames = []
268        dx = .00001
269        dadx = np.zeros(6)
270        for i in range(6):
271            ia = i/3
272            ix = i%3
273            Oatoms[ia][ix+1] += dx
274            a0 = calcDist(Oatoms,SyOps,Amat)
275            Oatoms[ia][ix+1] -= 2*dx
276            dadx[i] = (calcDist(Oatoms,SyOps,Amat)-a0)/(2.*dx)
277        covMatrix = covData['covMatrix']
278        varyList = covData['varyList']
279        DistVcov = getVCov(names,varyList,covMatrix)
280        sig = np.sqrt(np.inner(dadx,np.inner(DistVcov,dadx)))
281        if sig < 0.001:
282            sig = -0.001
283   
284    return Dist,sig
285
286def GetAngleSig(Oatoms,Atoms,Amat,SGData,covData={}):
287       
288    def calcAngle(Atoms,SyOps,Amat):
289        XYZ = []
290        for i,atom in enumerate(Atoms):
291            Inv,M,T,C,U = SyOps[i]
292            XYZ.append(np.array(atom[1:4]))
293            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
294            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
295        V1 = XYZ[1]-XYZ[0]
296        V1 /= np.sqrt(np.sum(V1**2))
297        V2 = XYZ[1]-XYZ[2]
298        V2 /= np.sqrt(np.sum(V2**2))
299        V3 = V2-V1
300        cang = min(1.,max((2.-np.sum(V3**2))/2.,-1.))
301        return acosd(cang)
302
303    Inv = []
304    SyOps = []
305    names = []
306    for i,atom in enumerate(Oatoms):
307        names += atom[-1]
308        Op,unit = Atoms[i][-1]
309        inv = Op/abs(Op)
310        m,t = SGData['SGOps'][abs(Op)%100-1]
311        c = SGData['SGCen'][abs(Op)/100]
312        SyOps.append([inv,m,t,c,unit])
313    Angle = calcAngle(Oatoms,SyOps,Amat)
314   
315    sig = -0.01
316    if 'covMatrix' in covData:
317        parmNames = []
318        dx = .00001
319        dadx = np.zeros(9)
320        for i in range(9):
321            ia = i/3
322            ix = i%3
323            Oatoms[ia][ix+1] += dx
324            a0 = calcAngle(Oatoms,SyOps,Amat)
325            Oatoms[ia][ix+1] -= 2*dx
326            dadx[i] = (calcAngle(Oatoms,SyOps,Amat)-a0)/(2.*dx)
327        covMatrix = covData['covMatrix']
328        varyList = covData['varyList']
329        AngVcov = getVCov(names,varyList,covMatrix)
330        sig = np.sqrt(np.inner(dadx,np.inner(AngVcov,dadx)))
331        if sig < 0.01:
332            sig = -0.01
333   
334    return Angle,sig
335
336def GetTorsionSig(Oatoms,Atoms,Amat,SGData,covData={}):
337   
338    def calcTorsion(Atoms,SyOps,Amat):
339       
340        XYZ = []
341        for i,atom in enumerate(Atoms):
342            Inv,M,T,C,U = SyOps[i]
343            XYZ.append(np.array(atom[1:4]))
344            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
345            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
346        V1 = XYZ[1]-XYZ[0]
347        V2 = XYZ[2]-XYZ[1]
348        V3 = XYZ[3]-XYZ[2]
349        V1 /= np.sqrt(np.sum(V1**2))
350        V2 /= np.sqrt(np.sum(V2**2))
351        V3 /= np.sqrt(np.sum(V3**2))
352        M = np.array([V1,V2,V3])
353        D = nl.det(M)
354        Ang = 1.0
355        P12 = np.dot(V1,V2)
356        P13 = np.dot(V1,V3)
357        P23 = np.dot(V2,V3)
358        Tors = acosd((P12*P23-P13)/(np.sqrt(1.-P12**2)*np.sqrt(1.-P23**2)))*D/abs(D)
359        return Tors
360           
361    Inv = []
362    SyOps = []
363    names = []
364    for i,atom in enumerate(Oatoms):
365        names += atom[-1]
366        Op,unit = Atoms[i][-1]
367        inv = Op/abs(Op)
368        m,t = SGData['SGOps'][abs(Op)%100-1]
369        c = SGData['SGCen'][abs(Op)/100]
370        SyOps.append([inv,m,t,c,unit])
371    Tors = calcTorsion(Oatoms,SyOps,Amat)
372   
373    sig = -0.01
374    if 'covMatrix' in covData:
375        parmNames = []
376        dx = .00001
377        dadx = np.zeros(12)
378        for i in range(12):
379            ia = i/3
380            ix = i%3
381            Oatoms[ia][ix+1] += dx
382            a0 = calcTorsion(Oatoms,SyOps,Amat)
383            Oatoms[ia][ix+1] -= 2*dx
384            dadx[i] = (calcTorsion(Oatoms,SyOps,Amat)-a0)/(2.*dx)
385        covMatrix = covData['covMatrix']
386        varyList = covData['varyList']
387        TorVcov = getVCov(names,varyList,covMatrix)
388        sig = np.sqrt(np.inner(dadx,np.inner(TorVcov,dadx)))
389        if sig < 0.01:
390            sig = -0.01
391   
392    return Tors,sig
393       
394def GetDATSig(Oatoms,Atoms,Amat,SGData,covData={}):
395   
396    def calcDist(Atoms,SyOps,Amat):
397        XYZ = []
398        for i,atom in enumerate(Atoms):
399            Inv,M,T,C,U = SyOps[i]
400            XYZ.append(np.array(atom[1:4]))
401            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
402            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
403        V1 = XYZ[1]-XYZ[0]
404        return np.sqrt(np.sum(V1**2))
405       
406    def calcAngle(Atoms,SyOps,Amat):
407        XYZ = []
408        for i,atom in enumerate(Atoms):
409            Inv,M,T,C,U = SyOps[i]
410            XYZ.append(np.array(atom[1:4]))
411            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
412            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
413        V1 = XYZ[1]-XYZ[0]
414        V1 /= np.sqrt(np.sum(V1**2))
415        V2 = XYZ[1]-XYZ[2]
416        V2 /= np.sqrt(np.sum(V2**2))
417        V3 = V2-V1
418        cang = min(1.,max((2.-np.sum(V3**2))/2.,-1.))
419        return acosd(cang)
420
421    def calcTorsion(Atoms,SyOps,Amat):
422       
423        XYZ = []
424        for i,atom in enumerate(Atoms):
425            Inv,M,T,C,U = SyOps[i]
426            XYZ.append(np.array(atom[1:4]))
427            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
428            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
429        V1 = XYZ[1]-XYZ[0]
430        V2 = XYZ[2]-XYZ[1]
431        V3 = XYZ[3]-XYZ[2]
432        V1 /= np.sqrt(np.sum(V1**2))
433        V2 /= np.sqrt(np.sum(V2**2))
434        V3 /= np.sqrt(np.sum(V3**2))
435        M = np.array([V1,V2,V3])
436        D = nl.det(M)
437        Ang = 1.0
438        P12 = np.dot(V1,V2)
439        P13 = np.dot(V1,V3)
440        P23 = np.dot(V2,V3)
441        Tors = acosd((P12*P23-P13)/(np.sqrt(1.-P12**2)*np.sqrt(1.-P23**2)))*D/abs(D)
442        return Tors
443           
444    Inv = []
445    SyOps = []
446    names = []
447    for i,atom in enumerate(Oatoms):
448        names += atom[-1]
449        Op,unit = Atoms[i][-1]
450        inv = Op/abs(Op)
451        m,t = SGData['SGOps'][abs(Op)%100-1]
452        c = SGData['SGCen'][abs(Op)/100]
453        SyOps.append([inv,m,t,c,unit])
454    M = len(Oatoms)
455    if M == 2:
456        Val = calcDist(Oatoms,SyOps,Amat)
457    elif M == 3:
458        Val = calcAngle(Oatoms,SyOps,Amat)
459    else:
460        Val = calcTorsion(Oatoms,SyOps,Amat)
461   
462    sigVals = [-0.001,-0.01,-0.01]
463    sig = sigVals[M-3]
464    if 'covMatrix' in covData:
465        parmNames = []
466        dx = .00001
467        N = M*3
468        dadx = np.zeros(N)
469        for i in range(N):
470            ia = i/3
471            ix = i%3
472            Oatoms[ia][ix+1] += dx
473            if M == 2:
474                a0 = calcDist(Oatoms,SyOps,Amat)
475            elif M == 3:
476                a0 = calcAngle(Oatoms,SyOps,Amat)
477            else:
478                a0 = calcTorsion(Oatoms,SyOps,Amat)
479            Oatoms[ia][ix+1] -= 2*dx
480            if M == 2:
481                dadx[i] = (calcDist(Oatoms,SyOps,Amat)-a0)/(2.*dx)               
482            elif M == 3:
483                dadx[i] = (calcAngle(Oatoms,SyOps,Amat)-a0)/(2.*dx)               
484            else:
485                dadx[i] = (calcTorsion(Oatoms,SyOps,Amat)-a0)/(2.*dx)
486        covMatrix = covData['covMatrix']
487        varyList = covData['varyList']
488        Vcov = getVCov(names,varyList,covMatrix)
489        sig = np.sqrt(np.inner(dadx,np.inner(Vcov,dadx)))
490        if sig < sigVals[M-3]:
491            sig = sigVals[M-3]
492   
493    return Val,sig
494       
495   
496def ValEsd(value,esd=0,nTZ=False):                  #NOT complete - don't use
497    # returns value(esd) string; nTZ=True for no trailing zeros
498    # use esd < 0 for level of precision shown e.g. esd=-0.01 gives 2 places beyond decimal
499    #get the 2 significant digits in the esd
500    edig = lambda esd: int(round(10**(math.log10(esd) % 1+1)))
501    #get the number of digits to represent them
502    epl = lambda esd: 2+int(1.545-math.log10(10*edig(esd)))
503   
504    mdec = lambda esd: -int(round(math.log10(abs(esd))))+1
505    ndec = lambda esd: int(1.545-math.log10(abs(esd)))
506    if esd > 0:
507        fmt = '"%.'+str(ndec(esd))+'f(%d)"'
508        return str(fmt%(value,int(round(esd*10**(mdec(esd)))))).strip('"')
509    elif esd < 0:
510         return str(round(value,mdec(esd)-1))
511    else:
512        text = str("%f"%(value))
513        if nTZ:
514            return text.rstrip('0')
515        else:
516            return text
517
518   
Note: See TracBrowser for help on using the repository browser.