/usr/share/pyshared/mdp/utils/_symeig.py is in python-mdp 3.3-1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 | import mdp
from mdp import numx, numx_linalg
class SymeigException(mdp.MDPException):
pass
# the following functions and classes were part of the scipy_emulation.py file
_type_keys = ['f', 'd', 'F', 'D']
_type_conv = {('f','d'): 'd', ('f','F'): 'F', ('f','D'): 'D',
('d','F'): 'D', ('d','D'): 'D',
('F','d'): 'D', ('F','D'): 'D'}
def _greatest_common_dtype(alist):
"""
Apply conversion rules to find the common conversion type
dtype 'd' is default for 'i' or unknown types
(known types: 'f','d','F','D').
"""
dtype = 'f'
for array in alist:
if array is None:
continue
tc = array.dtype.char
if tc not in _type_keys:
tc = 'd'
transition = (dtype, tc)
if transition in _type_conv:
dtype = _type_conv[transition]
return dtype
def _assert_eigenvalues_real_and_positive(w, dtype):
tol = numx.finfo(dtype.type).eps * 100
if abs(w.imag).max() > tol:
err = "Some eigenvalues have significant imaginary part: %s " % str(w)
raise mdp.SymeigException(err)
#if w.real.min() < 0:
# err = "Got negative eigenvalues: %s" % str(w)
# raise SymeigException(err)
def wrap_eigh(A, B = None, eigenvectors = True, turbo = "on", range = None,
type = 1, overwrite = False):
"""Wrapper for scipy.linalg.eigh for scipy version > 0.7"""
args = {}
args['a'] = A
args['b'] = B
args['eigvals_only'] = not eigenvectors
args['overwrite_a'] = overwrite
args['overwrite_b'] = overwrite
if turbo == "on":
args['turbo'] = True
else:
args['turbo'] = False
args['type'] = type
if range is not None:
n = A.shape[0]
lo, hi = range
if lo < 1:
lo = 1
if lo > n:
lo = n
if hi > n:
hi = n
if lo > hi:
lo, hi = hi, lo
# in scipy.linalg.eigh the range starts from 0
lo -= 1
hi -= 1
range = (lo, hi)
args['eigvals'] = range
try:
return numx_linalg.eigh(**args)
except numx_linalg.LinAlgError, exception:
raise SymeigException(str(exception))
def _symeig_fake(A, B = None, eigenvectors = True, turbo = "on", range = None,
type = 1, overwrite = False):
"""Solve standard and generalized eigenvalue problem for symmetric
(hermitian) definite positive matrices.
This function is a wrapper of LinearAlgebra.eigenvectors or
numarray.linear_algebra.eigenvectors with an interface compatible with symeig.
Syntax:
w,Z = symeig(A)
w = symeig(A,eigenvectors=0)
w,Z = symeig(A,range=(lo,hi))
w,Z = symeig(A,B,range=(lo,hi))
Inputs:
A -- An N x N matrix.
B -- An N x N matrix.
eigenvectors -- if set return eigenvalues and eigenvectors, otherwise
only eigenvalues
turbo -- not implemented
range -- the tuple (lo,hi) represent the indexes of the smallest and
largest (in ascending order) eigenvalues to be returned.
1 <= lo < hi <= N
if range = None, returns all eigenvalues and eigenvectors.
type -- not implemented, always solve A*x = (lambda)*B*x
overwrite -- not implemented
Outputs:
w -- (selected) eigenvalues in ascending order.
Z -- if range = None, Z contains the matrix of eigenvectors,
normalized as follows:
Z^H * A * Z = lambda and Z^H * B * Z = I
where ^H means conjugate transpose.
if range, an N x M matrix containing the orthonormal
eigenvectors of the matrix A corresponding to the selected
eigenvalues, with the i-th column of Z holding the eigenvector
associated with w[i]. The eigenvectors are normalized as above.
"""
dtype = numx.dtype(_greatest_common_dtype([A, B]))
try:
if B is None:
w, Z = numx_linalg.eigh(A)
else:
# make B the identity matrix
wB, ZB = numx_linalg.eigh(B)
_assert_eigenvalues_real_and_positive(wB, dtype)
ZB = ZB.real / numx.sqrt(wB.real)
# transform A in the new basis: A = ZB^T * A * ZB
A = mdp.utils.mult(mdp.utils.mult(ZB.T, A), ZB)
# diagonalize A
w, ZA = numx_linalg.eigh(A)
Z = mdp.utils.mult(ZB, ZA)
except numx_linalg.LinAlgError, exception:
raise SymeigException(str(exception))
_assert_eigenvalues_real_and_positive(w, dtype)
w = w.real
Z = Z.real
idx = w.argsort()
w = w.take(idx)
Z = Z.take(idx, axis=1)
# sanitize range:
n = A.shape[0]
if range is not None:
lo, hi = range
if lo < 1:
lo = 1
if lo > n:
lo = n
if hi > n:
hi = n
if lo > hi:
lo, hi = hi, lo
Z = Z[:, lo-1:hi]
w = w[lo-1:hi]
# the final call to refcast is necessary because of a bug in the casting
# behavior of Numeric and numarray: eigenvector does not wrap the LAPACK
# single precision routines
if eigenvectors:
return mdp.utils.refcast(w, dtype), mdp.utils.refcast(Z, dtype)
else:
return mdp.utils.refcast(w, dtype)
|