def __call__(self, params):
print '???', params
sd1 = params[0]
sd2 = params[1]
cor = params[2]
if sd1 < 0. or sd1 > 10. or sd2 < 0. or sd2 > 10. or cor < -1. or cor > 1.:
return np.inf
bandwidth = maths.stats.choleskysqrt2d(sd1, sd2, cor)
bandwidthdet = la.det(bandwidth)
bandwidthinv = la.inv(bandwidth)
diff = sample[self.__iidx] - sample[self.__jidx]
temp = diff.dot(bandwidthinv.T)
temp *= temp
e = np.exp(np.sum(temp, axis=1))
s = np.sum(e**(-.25) - 4 * e**(-.5))
cost = self.__n / bandwidthdet + (2. / bandwidthdet) * s
print '!!!', cost
return cost / 10000.
python类det()的实例源码
def do(self, a, b):
d = linalg.det(a)
(s, ld) = linalg.slogdet(a)
if asarray(a).dtype.type in (single, double):
ad = asarray(a).astype(double)
else:
ad = asarray(a).astype(cdouble)
ev = linalg.eigvals(ad)
assert_almost_equal(d, multiply.reduce(ev, axis=-1))
assert_almost_equal(s * np.exp(ld), multiply.reduce(ev, axis=-1))
s = np.atleast_1d(s)
ld = np.atleast_1d(ld)
m = (s != 0)
assert_almost_equal(np.abs(s[m]), 1)
assert_equal(ld[~m], -inf)
def test_byteorder_check():
# Byte order check should pass for native order
if sys.byteorder == 'little':
native = '<'
else:
native = '>'
for dtt in (np.float32, np.float64):
arr = np.eye(4, dtype=dtt)
n_arr = arr.newbyteorder(native)
sw_arr = arr.newbyteorder('S').byteswap()
assert_equal(arr.dtype.byteorder, '=')
for routine in (linalg.inv, linalg.det, linalg.pinv):
# Normal call
res = routine(arr)
# Native but not '='
assert_array_equal(res, routine(n_arr))
# Swapped
assert_array_equal(res, routine(sw_arr))
def do(self, a, b):
d = linalg.det(a)
(s, ld) = linalg.slogdet(a)
if asarray(a).dtype.type in (single, double):
ad = asarray(a).astype(double)
else:
ad = asarray(a).astype(cdouble)
ev = linalg.eigvals(ad)
assert_almost_equal(d, multiply.reduce(ev, axis=-1))
assert_almost_equal(s * np.exp(ld), multiply.reduce(ev, axis=-1))
s = np.atleast_1d(s)
ld = np.atleast_1d(ld)
m = (s != 0)
assert_almost_equal(np.abs(s[m]), 1)
assert_equal(ld[~m], -inf)
test_linalg.py 文件源码
项目:PyDataLondon29-EmbarrassinglyParallelDAWithAWSLambda
作者: SignalMedia
项目源码
文件源码
阅读 21
收藏 0
点赞 0
评论 0
def do(self, a, b):
d = linalg.det(a)
(s, ld) = linalg.slogdet(a)
if asarray(a).dtype.type in (single, double):
ad = asarray(a).astype(double)
else:
ad = asarray(a).astype(cdouble)
ev = linalg.eigvals(ad)
assert_almost_equal(d, multiply.reduce(ev, axis=-1))
assert_almost_equal(s * np.exp(ld), multiply.reduce(ev, axis=-1))
s = np.atleast_1d(s)
ld = np.atleast_1d(ld)
m = (s != 0)
assert_almost_equal(np.abs(s[m]), 1)
assert_equal(ld[~m], -inf)
test_linalg.py 文件源码
项目:PyDataLondon29-EmbarrassinglyParallelDAWithAWSLambda
作者: SignalMedia
项目源码
文件源码
阅读 21
收藏 0
点赞 0
评论 0
def test_byteorder_check():
# Byte order check should pass for native order
if sys.byteorder == 'little':
native = '<'
else:
native = '>'
for dtt in (np.float32, np.float64):
arr = np.eye(4, dtype=dtt)
n_arr = arr.newbyteorder(native)
sw_arr = arr.newbyteorder('S').byteswap()
assert_equal(arr.dtype.byteorder, '=')
for routine in (linalg.inv, linalg.det, linalg.pinv):
# Normal call
res = routine(arr)
# Native but not '='
assert_array_equal(res, routine(n_arr))
# Swapped
assert_array_equal(res, routine(sw_arr))
def do(self, a, b):
d = linalg.det(a)
(s, ld) = linalg.slogdet(a)
if asarray(a).dtype.type in (single, double):
ad = asarray(a).astype(double)
else:
ad = asarray(a).astype(cdouble)
ev = linalg.eigvals(ad)
assert_almost_equal(d, multiply.reduce(ev, axis=-1))
assert_almost_equal(s * np.exp(ld), multiply.reduce(ev, axis=-1))
s = np.atleast_1d(s)
ld = np.atleast_1d(ld)
m = (s != 0)
assert_almost_equal(np.abs(s[m]), 1)
assert_equal(ld[~m], -inf)
def test_byteorder_check():
# Byte order check should pass for native order
if sys.byteorder == 'little':
native = '<'
else:
native = '>'
for dtt in (np.float32, np.float64):
arr = np.eye(4, dtype=dtt)
n_arr = arr.newbyteorder(native)
sw_arr = arr.newbyteorder('S').byteswap()
assert_equal(arr.dtype.byteorder, '=')
for routine in (linalg.inv, linalg.det, linalg.pinv):
# Normal call
res = routine(arr)
# Native but not '='
assert_array_equal(res, routine(n_arr))
# Swapped
assert_array_equal(res, routine(sw_arr))
def do(self, a, b):
d = linalg.det(a)
(s, ld) = linalg.slogdet(a)
if asarray(a).dtype.type in (single, double):
ad = asarray(a).astype(double)
else:
ad = asarray(a).astype(cdouble)
ev = linalg.eigvals(ad)
assert_almost_equal(d, multiply.reduce(ev, axis=-1))
assert_almost_equal(s * np.exp(ld), multiply.reduce(ev, axis=-1))
s = np.atleast_1d(s)
ld = np.atleast_1d(ld)
m = (s != 0)
assert_almost_equal(np.abs(s[m]), 1)
assert_equal(ld[~m], -inf)
def test_byteorder_check():
# Byte order check should pass for native order
if sys.byteorder == 'little':
native = '<'
else:
native = '>'
for dtt in (np.float32, np.float64):
arr = np.eye(4, dtype=dtt)
n_arr = arr.newbyteorder(native)
sw_arr = arr.newbyteorder('S').byteswap()
assert_equal(arr.dtype.byteorder, '=')
for routine in (linalg.inv, linalg.det, linalg.pinv):
# Normal call
res = routine(arr)
# Native but not '='
assert_array_equal(res, routine(n_arr))
# Swapped
assert_array_equal(res, routine(sw_arr))
def do(self, a, b):
d = linalg.det(a)
(s, ld) = linalg.slogdet(a)
if asarray(a).dtype.type in (single, double):
ad = asarray(a).astype(double)
else:
ad = asarray(a).astype(cdouble)
ev = linalg.eigvals(ad)
assert_almost_equal(d, multiply.reduce(ev, axis=-1))
assert_almost_equal(s * np.exp(ld), multiply.reduce(ev, axis=-1))
s = np.atleast_1d(s)
ld = np.atleast_1d(ld)
m = (s != 0)
assert_almost_equal(np.abs(s[m]), 1)
assert_equal(ld[~m], -inf)
def test_byteorder_check():
# Byte order check should pass for native order
if sys.byteorder == 'little':
native = '<'
else:
native = '>'
for dtt in (np.float32, np.float64):
arr = np.eye(4, dtype=dtt)
n_arr = arr.newbyteorder(native)
sw_arr = arr.newbyteorder('S').byteswap()
assert_equal(arr.dtype.byteorder, '=')
for routine in (linalg.inv, linalg.det, linalg.pinv):
# Normal call
res = routine(arr)
# Native but not '='
assert_array_equal(res, routine(n_arr))
# Swapped
assert_array_equal(res, routine(sw_arr))
def F_value_multivariate(ER, EF, dfnum, dfden):
"""
Returns an F-statistic given the following:
ER = error associated with the null hypothesis (the Restricted model)
EF = error associated with the alternate hypothesis (the Full model)
dfR = degrees of freedom the Restricted model
dfF = degrees of freedom associated with the Restricted model
where ER and EF are matrices from a multivariate F calculation.
"""
if type(ER) in [IntType, FloatType]:
ER = N.array([[ER]])
if type(EF) in [IntType, FloatType]:
EF = N.array([[EF]])
n_um = (LA.det(ER) - LA.det(EF)) / float(dfnum)
d_en = LA.det(EF) / float(dfden)
return n_um / d_en
#####################################
####### ASUPPORT FUNCTIONS ########
#####################################
def _int_var_rbf(self, X, hyp, jitter=1e-8):
"""
Posterior integral variance of the Gaussian Process quadrature.
X - vector (1, 2*xdim**2+xdim)
hyp - kernel hyperparameters [s2, el_1, ... el_d]
"""
# reshape X to SP matrix
X = np.reshape(X, (self.n, self.d))
# set kernel hyper-parameters
s2, el = hyp[0], hyp[1:]
self.kern.param_array[0] = s2 # variance
self.kern.param_array[1:] = el # lengthscale
K = self.kern.K(X)
L = np.diag(el ** 2)
# posterior variance of the integral
ks = s2 * np.sqrt(det(L + np.eye(self.d))) * multivariate_normal(mean=np.zeros(self.d), cov=L).pdf(X)
postvar = -ks.dot(solve(K + jitter * np.eye(self.n), ks.T))
return postvar
def _int_var_rbf_hyp(self, hyp, X, jitter=1e-8):
"""
Posterior integral variance as a function of hyper-parameters
:param hyp: RBF kernel hyper-parameters [s2, el_1, ..., el_d]
:param X: sigma-points
:param jitter: numerical jitter (for stabilizing computations)
:return: posterior integral variance
"""
# reshape X to SP matrix
X = np.reshape(X, (self.n, self.d))
# set kernel hyper-parameters
s2, el = 1, hyp # sig_var hyper always set to 1
self.kern.param_array[0] = s2 # variance
self.kern.param_array[1:] = el # lengthscale
K = self.kern.K(X)
L = np.diag(el ** 2)
# posterior variance of the integral
ks = s2 * np.sqrt(det(L + np.eye(self.d))) * multivariate_normal(mean=np.zeros(self.d), cov=L).pdf(X)
postvar = s2 * np.sqrt(det(2 * inv(L) + np.eye(self.d))) ** -1 - ks.dot(
solve(K + jitter * np.eye(self.n), ks.T))
return postvar
def genCovariace(size):
MaxIter = 1e+6
S = np.zeros((size,size))
itn = 0
while(alg.det(S) <= 1e-3 and itn < MaxIter):
itn = itn + 1
#print int(numpy.log2(size))*size
G6 = GenRndGnm(PUNGraph, size, int((size*(size-1))*0.05))
S = np.zeros((size,size))
for EI in G6.Edges():
S[EI.GetSrcNId(), EI.GetDstNId()] = 0.6
S = S + S.T + S.max()*np.matrix(np.eye(size))
if itn == MaxIter:
print 'fail to find an invertible sparse inverse covariance matrix'
S = np.asarray(S)
return S
def genCovariace(size):
MaxIter = 1e+6
S = np.zeros((size,size))
itn = 0
while(alg.det(S) <= 1e-3 and itn < MaxIter):
itn = itn + 1
#print int(numpy.log2(size))*size
G6 = GenRndGnm(PUNGraph, size, int((size*(size-1))*0.05))
#G6 = snap.GenRndGnm(snap.PUNGraph, 5, 5)
S = np.zeros((size,size))
for EI in G6.Edges():
S[EI.GetSrcNId(), EI.GetDstNId()] = 0.6
S = S + S.T + S.max()*np.matrix(np.eye(size))
if itn == MaxIter:
print 'fail to find an invertible sparse inverse covariance matrix'
S = np.asarray(S)
return S
def do(self, a, b):
d = linalg.det(a)
(s, ld) = linalg.slogdet(a)
if asarray(a).dtype.type in (single, double):
ad = asarray(a).astype(double)
else:
ad = asarray(a).astype(cdouble)
ev = linalg.eigvals(ad)
assert_almost_equal(d, multiply.reduce(ev, axis=-1))
assert_almost_equal(s * np.exp(ld), multiply.reduce(ev, axis=-1))
s = np.atleast_1d(s)
ld = np.atleast_1d(ld)
m = (s != 0)
assert_almost_equal(np.abs(s[m]), 1)
assert_equal(ld[~m], -inf)
def test_byteorder_check():
# Byte order check should pass for native order
if sys.byteorder == 'little':
native = '<'
else:
native = '>'
for dtt in (np.float32, np.float64):
arr = np.eye(4, dtype=dtt)
n_arr = arr.newbyteorder(native)
sw_arr = arr.newbyteorder('S').byteswap()
assert_equal(arr.dtype.byteorder, '=')
for routine in (linalg.inv, linalg.det, linalg.pinv):
# Normal call
res = routine(arr)
# Native but not '='
assert_array_equal(res, routine(n_arr))
# Swapped
assert_array_equal(res, routine(sw_arr))
def F(self, other):
"""
Compute the fundamental matrix with respect to other camera
http://www.robots.ox.ac.uk/~vgg/hzbook/code/vgg_multiview/vgg_F_from_P.m
The computed fundamental matrix, given by the formula 17.3 (p. 412) in
Hartley & Zisserman book (2nd ed.).
Use as:
F_10 = poses[0].F(poses[1])
l_1 = F_10 * x_0
"""
X1 = self.P[[1,2],:]
X2 = self.P[[2,0],:]
X3 = self.P[[0,1],:]
Y1 = other.P[[1,2],:]
Y2 = other.P[[2,0],:]
Y3 = other.P[[0,1],:]
F = np.float64([
[det(np.vstack([X1, Y1])), det(np.vstack([X2, Y1])), det(np.vstack([X3, Y1]))],
[det(np.vstack([X1, Y2])), det(np.vstack([X2, Y2])), det(np.vstack([X3, Y2]))],
[det(np.vstack([X1, Y3])), det(np.vstack([X2, Y3])), det(np.vstack([X3, Y3]))]
])
return F # / F[2,2]
def __init__(self, density, bandwidth):
self.__density = density
self.__bandwidth = bandwidth
self.__bandwidthinv = la.inv(bandwidth)
self.__bandwidthdet = la.det(bandwidth)
def test_zero(self):
assert_equal(linalg.det([[0.0]]), 0.0)
assert_equal(type(linalg.det([[0.0]])), double)
assert_equal(linalg.det([[0.0j]]), 0.0)
assert_equal(type(linalg.det([[0.0j]])), cdouble)
assert_equal(linalg.slogdet([[0.0]]), (0.0, -inf))
assert_equal(type(linalg.slogdet([[0.0]])[0]), double)
assert_equal(type(linalg.slogdet([[0.0]])[1]), double)
assert_equal(linalg.slogdet([[0.0j]]), (0.0j, -inf))
assert_equal(type(linalg.slogdet([[0.0j]])[0]), cdouble)
assert_equal(type(linalg.slogdet([[0.0j]])[1]), double)
def test_types(self):
def check(dtype):
x = np.array([[1, 0.5], [0.5, 1]], dtype=dtype)
assert_equal(np.linalg.det(x).dtype, dtype)
ph, s = np.linalg.slogdet(x)
assert_equal(s.dtype, get_real_dtype(dtype))
assert_equal(ph.dtype, dtype)
for dtype in [single, double, csingle, cdouble]:
yield check, dtype
def is_left(x0, x1, x2):
"""Returns True if x0 is left of the line between x1 and x2,
False otherwise. Ref: https://stackoverflow.com/questions/1560492"""
assert x1.shape == x2.shape == (2,)
matrix = array([x1-x0, x2-x0])
if len(x0.shape) == 2:
matrix = matrix.transpose((1, 2, 0))
return det(matrix) > 0
def log_prob(self, x):
from numpy.linalg import det
mu = self.mu
S = self.sigma
D = len(mu)
q = self.__q(x)
return -0.5 * (D * log(2 * pi) + log(abs(det(S)))) - 0.5 * q ** 2
def rmsd(X, Y):
"""
Calculate the root mean squared deviation (RMSD) using Kabsch' formula.
@param X: (n, d) input vector
@type X: numpy array
@param Y: (n, d) input vector
@type Y: numpy array
@return: rmsd value between the input vectors
@rtype: float
"""
from numpy import sum, dot, sqrt, clip, average
from numpy.linalg import svd, det
X = X - X.mean(0)
Y = Y - Y.mean(0)
R_x = sum(X ** 2)
R_y = sum(Y ** 2)
V, L, U = svd(dot(Y.T, X))
if det(dot(V, U)) < 0.:
L[-1] *= -1
return sqrt(clip(R_x + R_y - 2 * sum(L), 0., 1e300) / len(X))
def gauss2D(x, mean, cov):
# x, mean, cov??numpy.array??
z = -np.dot(np.dot((x-mean).T,inv(cov)),(x-mean))/2.0
temp = pow(sqrt(2.0*pi),len(x))*sqrt(det(cov))
return (1.0/temp)*exp(z)
def emit_prob_updated(self, X, post_state): # ??????
for k in range(self.n_state):
for j in range(self.x_size):
self.emit_means[k][j] = np.sum(post_state[:,k] *X[:,j]) / np.sum(post_state[:,k])
X_cov = np.dot((X-self.emit_means[k]).T, (post_state[:,k]*(X-self.emit_means[k]).T).T)
self.emit_covars[k] = X_cov / np.sum(post_state[:,k])
if det(self.emit_covars[k]) == 0: # ????????
self.emit_covars[k] = self.emit_covars[k] + 0.01*np.eye(len(X[0]))
def test_zero(self):
assert_equal(linalg.det([[0.0]]), 0.0)
assert_equal(type(linalg.det([[0.0]])), double)
assert_equal(linalg.det([[0.0j]]), 0.0)
assert_equal(type(linalg.det([[0.0j]])), cdouble)
assert_equal(linalg.slogdet([[0.0]]), (0.0, -inf))
assert_equal(type(linalg.slogdet([[0.0]])[0]), double)
assert_equal(type(linalg.slogdet([[0.0]])[1]), double)
assert_equal(linalg.slogdet([[0.0j]]), (0.0j, -inf))
assert_equal(type(linalg.slogdet([[0.0j]])[0]), cdouble)
assert_equal(type(linalg.slogdet([[0.0j]])[1]), double)
def test_types(self):
def check(dtype):
x = np.array([[1, 0.5], [0.5, 1]], dtype=dtype)
assert_equal(np.linalg.det(x).dtype, dtype)
ph, s = np.linalg.slogdet(x)
assert_equal(s.dtype, get_real_dtype(dtype))
assert_equal(ph.dtype, dtype)
for dtype in [single, double, csingle, cdouble]:
yield check, dtype