Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
def test_invsqrtm():
"""Test matrix inverse square root"""
C = 2*np.eye(3)
Ctrue = (1.0/np.sqrt(2))*np.eye(3)
assert_array_almost_equal(invsqrtm(C), Ctrue)
def transport(Covs, Cref, metric='riemann'):
"""Parallel transport of two set of covariance matrix.
"""
C = mean_covariance(Covs, metric=metric)
iC = invsqrtm(C)
E = sqrtm(numpy.dot(numpy.dot(iC, Cref), iC))
out = numpy.array([numpy.dot(numpy.dot(E, c), E.T) for c in Covs])
return out
Nt, Ne, Ne = covmats.shape
crit = numpy.inf
k = 0
# init with AJD
B, _ = ajd_pham(covmats)
while (crit > tol) and (k < maxiter):
k += 1
J = numpy.zeros((Ne, Ne))
for index, Ci in enumerate(covmats):
tmp = logm(numpy.dot(numpy.dot(B.T, Ci), B))
J += sample_weight[index] * tmp
update = numpy.diag(numpy.diag(expm(J)))
B = numpy.dot(B, invsqrtm(update))
crit = distance_riemann(numpy.eye(Ne), update)
A = numpy.linalg.inv(B)
J = numpy.zeros((Ne, Ne))
for index, Ci in enumerate(covmats):
tmp = logm(numpy.dot(numpy.dot(B.T, Ci), B))
J += sample_weight[index] * tmp
C = numpy.dot(numpy.dot(A.T, expm(J)), A)
return C
# init
sample_weight = _get_sample_weight(sample_weight, covmats)
Nt, Ne, Ne = covmats.shape
if init is None:
C = numpy.mean(covmats, axis=0)
else:
C = init
k = 0
nu = 1.0
tau = numpy.finfo(numpy.float64).max
crit = numpy.finfo(numpy.float64).max
# stop when J<10^-9 or max iteration = 50
while (crit > tol) and (k < maxiter) and (nu > tol):
k = k + 1
C12 = sqrtm(C)
Cm12 = invsqrtm(C)
J = numpy.zeros((Ne, Ne))
for index in range(Nt):
tmp = numpy.dot(numpy.dot(Cm12, covmats[index, :, :]), Cm12)
J += sample_weight[index] * logm(tmp)
crit = numpy.linalg.norm(J, ord='fro')
h = nu * crit
C = numpy.dot(numpy.dot(C12, expm(nu * J)), C12)
if h < tau:
nu = 0.95 * nu
tau = h
else:
nu = 0.5 * nu
return C
def tangent_space(covmats, Cref):
"""Project a set of covariance matrices in the tangent space. according to
the reference point Cref
:param covmats: np.ndarray
Covariance matrices set, Ntrials X Nchannels X Nchannels
:param Cref: np.ndarray
The reference covariance matrix
:returns: np.ndarray
the Tangent space , a matrix of Ntrials X (Nchannels*(Nchannels+1)/2)
"""
Nt, Ne, Ne = covmats.shape
Cm12 = invsqrtm(Cref)
idx = numpy.triu_indices_from(Cref)
Nf = int(Ne * (Ne + 1) / 2)
T = numpy.empty((Nt, Nf))
coeffs = (numpy.sqrt(2) * numpy.triu(numpy.ones((Ne, Ne)), 1) +
numpy.eye(Ne))[idx]
for index in range(Nt):
tmp = numpy.dot(numpy.dot(Cm12, covmats[index, :, :]), Cm12)
tmp = logm(tmp)
T[index, :] = numpy.multiply(coeffs, tmp[idx])
return T
def geodesic_riemann(A, B, alpha=0.5):
"""Return the matrix at the position alpha on the riemannian geodesic between A and B :
.. math::
\mathbf{C} = \mathbf{A}^{1/2} \left( \mathbf{A}^{-1/2} \mathbf{B} \mathbf{A}^{-1/2} \\right)^\\alpha \mathbf{A}^{1/2}
C is equal to A if alpha = 0 and B if alpha = 1
:param A: the first coavriance matrix
:param B: the second coavriance matrix
:param alpha: the position on the geodesic
:returns: the covariance matrix
"""
sA = sqrtm(A)
isA = invsqrtm(A)
C = numpy.dot(numpy.dot(isA, B), isA)
D = powm(C, alpha)
E = numpy.dot(numpy.dot(sA, D), sA)
return E