Page MenuHomec4science
Files F92030793

mathtools.py
No OneTemporary
Actions

File Metadata

Created
Sat, Nov 16, 18:45

mathtools.py
View Options

"""Contains simple algorithms.
Copyright (C) 2013, Joshua More and Michele Ceriotti
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program. If not, see <http.//www.gnu.org/licenses/>.
Functions:
matrix_exp: Computes the exponential of a square matrix via a Taylor series.
stab_cholesky: A numerically stable version of the Cholesky decomposition.
h2abc: Takes the representation of the system box in terms of an upper
triangular matrix of column vectors, and returns the representation in
terms of the lattice vector lengths and the angles between them
in radians.
h2abc_deg: Takes the representation of the system box in terms of an upper
triangular matrix of column vectors, and returns the representation in
terms of the lattice vector lengths and the angles between them in
degrees.
abc2h: Takes the representation of the system box in terms of the lattice
vector lengths and the angles between them, and returns the
representation in terms of an upper triangular lattice vector matrix.
invert_ut3x3: Inverts a 3*3 upper triangular matrix.
det_ut3x3(h): Finds the determinant of a 3*3 upper triangular matrix.
eigensystem_ut3x3: Finds the eigenvector matrix and eigenvalues of a 3*3
upper triangular matrix
exp_ut3x3: Computes the exponential of a 3*3 upper triangular matrix.
root_herm: Computes the square root of a positive-definite hermitian
matrix.
logsumlog: Routine to accumulate the logarithm of a sum
"""
__all__ = ['matrix_exp', 'stab_cholesky', 'h2abc', 'h2abc_deg', 'abc2h',
'invert_ut3x3', 'det_ut3x3', 'eigensystem_ut3x3', 'exp_ut3x3',
'root_herm', 'logsumlog' ]
import numpy as np
import math
from ipi.utils.messages import verbosity, warning
def logsumlog(lasa, lbsb):
"""Computes log(|A+B|) and sign(A+B) given log(|A|), log(|B|),
sign(A), sign(B).
Args:
lasa: (log(|A|), sign(A)) as a tuple
lbsb: (log(|B|), sign(B)) as a tuple
Returns:
(log(|A+B|), sign(A+B)) as a tuple
"""
(la,sa) = lasa
(lb,sb) = lbsb
if (la > lb):
sr = sa
lr = la + np.log(1.0 + sb*np.exp(lb-la))
else:
sr = sb
lr = lb + np.log(1.0 + sa*np.exp(la-lb))
return (lr,sr)
def matrix_exp(M, ntaylor=15, nsquare=15):
"""Computes the exponential of a square matrix via a Taylor series.
Calculates the matrix exponential by first calculating exp(M/(2**nsquare)),
then squaring the result the appropriate number of times.
Args:
M: Matrix to be exponentiated.
ntaylor: Optional integer giving the number of terms in the Taylor series.
Defaults to 15.
nsquare: Optional integer giving how many times the original matrix will
be halved. Defaults to 15.
Returns:
The matrix exponential of M.
"""
n = M.shape[1]
tc = np.zeros(ntaylor+1)
tc[0] = 1.0
for i in range(ntaylor):
tc[i+1] = tc[i]/(i+1)
SM = np.copy(M)/2.0**nsquare
EM = np.identity(n,float)*tc[ntaylor]
for i in range(ntaylor-1,-1,-1):
EM = np.dot(SM,EM)
EM += np.identity(n)*tc[i]
for i in range(nsquare):
EM = np.dot(EM,EM)
return EM
def stab_cholesky(M):
""" A numerically stable version of the Cholesky decomposition.
Used in the GLE implementation. Since many of the matrices used in this
algorithm have very large and very small numbers in at once, to handle a
wide range of frequencies, a naive algorithm can end up having to calculate
the square root of a negative number, which breaks the algorithm. This is
due to numerical precision errors turning a very tiny positive eigenvalue
into a tiny negative value.
Instead of this, an LDU decomposition is used, and any small negative numbers
in the diagonal D matrix are assumed to be due to numerical precision errors,
and so are replaced with zero.
Args:
M: The matrix to be decomposed.
"""
n = M.shape[1]
D = np.zeros(n,float)
L = np.zeros(M.shape,float)
for i in range(n):
L[i,i] = 1.
for j in range(i):
L[i,j] = M[i,j]
for k in range(j):
L[i,j] -= L[i,k]*L[j,k]*D[k]
if (not D[j] == 0.0):
L[i,j] = L[i,j]/D[j]
D[i] = M[i,i]
for k in range(i):
D[i] -= L[i,k]*L[i,k]*D[k]
S = np.zeros(M.shape,float)
for i in range(n):
if (D[i]>0):
D[i] = math.sqrt(D[i])
else:
warning("Zeroing negative element in stab-cholesky decomposition: " + str(D[i]), verbosity.low)
D[i] = 0
for j in range(i+1):
S[i,j] += L[i,j]*D[j]
return S
def h2abc(h):
"""Returns a description of the cell in terms of the length of the
lattice vectors and the angles between them in radians.
Args:
h: Cell matrix in upper triangular column vector form.
Returns:
A list containing the lattice vector lengths and the angles between them.
"""
a = float(h[0,0])
b = math.sqrt(h[0,1]**2 + h[1,1]**2)
c = math.sqrt(h[0,2]**2 + h[1,2]**2 + h[2,2]**2)
gamma = math.acos(h[0,1]/b)
beta = math.acos(h[0,2]/c)
alpha = math.acos(np.dot(h[:,1], h[:,2])/(b*c))
return a, b, c, alpha, beta, gamma
def h2abc_deg(h):
"""Returns a description of the cell in terms of the length of the
lattice vectors and the angles between them in degrees.
Args:
h: Cell matrix in upper triangular column vector form.
Returns:
A list containing the lattice vector lengths and the angles between them
in degrees.
"""
(a, b, c, alpha, beta, gamma) = h2abc(h)
return a, b, c, alpha*180/math.pi, beta*180/math.pi, gamma*180/math.pi
def abc2h(a, b, c, alpha, beta, gamma):
"""Returns a lattice vector matrix given a description in terms of the
lattice vector lengths and the angles in between.
Args:
a: First cell vector length.
b: Second cell vector length.
c: Third cell vector length.
alpha: Angle between sides b and c in radians.
beta: Angle between sides a and c in radians.
gamma: Angle between sides b and a in radians.
Returns:
An array giving the lattice vector matrix in upper triangular form.
"""
h = np.zeros((3,3) ,float)
h[0,0] = a
h[0,1] = b*math.cos(gamma)
h[0,2] = c*math.cos(beta)
h[1,1] = b*math.sin(gamma)
h[1,2] = (b*c*math.cos(alpha) - h[0,1]*h[0,2])/h[1,1]
h[2,2] = math.sqrt(c**2 - h[0,2]**2 - h[1,2]**2)
return h
def invert_ut3x3(h):
"""Inverts a 3*3 upper triangular matrix.
Args:
h: An upper triangular 3*3 matrix.
Returns:
The inverse matrix of h.
"""
ih = np.zeros((3,3), float)
for i in range(3):
ih[i,i] = 1.0/h[i,i]
ih[0,1] = -ih[0,0]*h[0,1]*ih[1,1]
ih[1,2] = -ih[1,1]*h[1,2]*ih[2,2]
ih[0,2] = -ih[1,2]*h[0,1]*ih[0,0] - ih[0,0]*h[0,2]*ih[2,2]
return ih
def eigensystem_ut3x3(p):
"""Finds the eigenvector matrix of a 3*3 upper-triangular matrix.
Args:
p: An upper triangular 3*3 matrix.
Returns:
An array giving the 3 eigenvalues of p, and the eigenvector matrix of p.
"""
eigp = np.zeros((3,3), float)
eigvals = np.zeros(3, float)
for i in range(3):
eigp[i,i] = 1
eigp[0,1] = -p[0,1]/(p[0,0] - p[1,1])
eigp[1,2] = -p[1,2]/(p[1,1] - p[2,2])
eigp[0,2] = -(p[0,1]*p[1,2] - p[0,2]*p[1,1] + p[0,2]*p[2,2])/((p[0,0] - p[2,2])*(p[2,2] - p[1,1]))
for i in range(3):
eigvals[i] = p[i,i]
return eigp, eigvals
def det_ut3x3(h):
"""Calculates the determinant of a 3*3 upper triangular matrix.
Note that the volume of the system box when the lattice vector matrix is
expressed as a 3*3 upper triangular matrix is given by the determinant of
this matrix.
Args:
h: An upper triangular 3*3 matrix.
Returns:
The determinant of h.
"""
return h[0,0]*h[1,1]*h[2,2]
MINSERIES=1e-8
def exp_ut3x3(h):
"""Computes the matrix exponential for a 3x3 upper-triangular matrix.
Note that for 3*3 upper triangular matrices this is the best method, as
it is stable. This is terms which become unstable as the
denominator tends to zero are calculated via a Taylor series in this limit.
Args:
h: An upper triangular 3*3 matrix.
Returns:
The matrix exponential of h.
"""
eh = np.zeros((3,3), float)
e00 = math.exp(h[0,0])
e11 = math.exp(h[1,1])
e22 = math.exp(h[2,2])
eh[0,0] = e00
eh[1,1] = e11
eh[2,2] = e22
if (abs((h[0,0] - h[1,1])/h[0,0])>MINSERIES):
r01 = (e00 - e11)/(h[0,0] - h[1,1])
else:
r01 = e00*(1 + (h[0,0] - h[1,1])*(0.5 + (h[0,0] - h[1,1])/6.0))
if (abs((h[1,1] - h[2,2])/h[1,1])>MINSERIES):
r12 = (e11 - e22)/(h[1,1] - h[2,2])
else:
r12 = e11*(1 + (h[1,1] - h[2,2])*(0.5 + (h[1,1] - h[2,2])/6.0))
if (abs((h[2,2] - h[0,0])/h[2,2])>MINSERIES):
r02 = (e22 - e00)/(h[2,2] - h[0,0])
else:
r02 = e22*(1 + (h[2,2] - h[0,0])*(0.5 + (h[2,2] - h[0,0])/6.0))
eh[0,1] = h[0,1]*r01
eh[1,2] = h[1,2]*r12
eh[0,2] = h[0,2]*r02
if (abs((h[2,2] - h[0,0])/h[2,2])>MINSERIES):
eh[0,2] += h[0,1]*h[0,2]*(r01 - r12)/(h[0,0] - h[2,2])
elif (abs((h[1,1] - h[0,0])/h[1,1])>MINSERIES):
eh[0,2] += h[0,1]*h[0,2]*(r12 - r02)/(h[1,1] - h[0,0])
elif (abs((h[1,1]-h[2,2])/h[1,1])>MINSERIES):
eh[0,2] += h[0,1]*h[0,2]*(r02 - r01)/(h[2,2] - h[1,1])
else:
eh[0,2] += h[0,1]*h[0,2]*e00/24.0*(12.0 + 4*(h[1,1] + h[2,2] - 2*h[0,0]) + (h[1,1] - h[0,0])*(h[2,2] - h[0,0]))
return eh
def root_herm(A):
"""Gives the square root of a hermitian matrix with real eigenvalues.
Args:
A: A Hermitian matrix.
Returns:
A matrix such that itself matrix multiplied by its transpose gives the
original matrix.
"""
if not (abs(A.T - A) < 1e-10).all():
raise ValueError("Non-Hermitian matrix passed to root_herm function")
eigvals, eigvecs = np.linalg.eigh(A)
ndgrs = len(eigvals)
diag = np.zeros((ndgrs,ndgrs))
for i in range(ndgrs):
if eigvals[i] >= 0:
diag[i,i] = math.sqrt(eigvals[i])
else:
warning("Zeroing negative element in matrix square root: " + str(eigvals[i]), verbosity.low)
diag[i,i] = 0
return np.dot(eigvecs, np.dot(diag, eigvecs.T))

Event Timeline

c4science · Help