def infExact_scipy_post(self, K, covars, y, sig2e, fixedEffects):
n = y.shape[0]
#mean vector
m = covars.dot(fixedEffects)
if (K.shape[1] < K.shape[0]): K_true = K.dot(K.T)
else: K_true = K
if sig2e<1e-6:
L = la.cholesky(K_true + sig2e*np.eye(n), overwrite_a=True, check_finite=False) #Cholesky factor of covariance with noise
sl = 1
pL = -self.solveChol(L, np.eye(n)) #L = -inv(K+inv(sW^2))
else:
L = la.cholesky(K_true/sig2e + np.eye(n), overwrite_a=True, check_finite=False) #Cholesky factor of B
sl = sig2e
pL = L #L = chol(eye(n)+sW*sW'.*K)
alpha = self.solveChol(L, y-m, overwrite_b=False) / sl
post = dict([])
post['alpha'] = alpha #return the posterior parameters
post['sW'] = np.ones(n) / np.sqrt(sig2e) #sqrt of noise precision vector
post['L'] = pL
return post
python类inv()的实例源码
def transfer_color(content, style):
import scipy.linalg as sl
# Mean and covariance of content
content_mean = np.mean(content, axis = (0, 1))
content_diff = content - content_mean
content_diff = np.reshape(content_diff, (-1, content_diff.shape[2]))
content_covariance = np.matmul(content_diff.T, content_diff) / (content_diff.shape[0])
# Mean and covariance of style
style_mean = np.mean(style, axis = (0, 1))
style_diff = style - style_mean
style_diff = np.reshape(style_diff, (-1, style_diff.shape[2]))
style_covariance = np.matmul(style_diff.T, style_diff) / (style_diff.shape[0])
# Calculate A and b
A = np.matmul(sl.sqrtm(content_covariance), sl.inv(sl.sqrtm(style_covariance)))
b = content_mean - np.matmul(A, style_mean)
# Construct new style
new_style = np.reshape(style, (-1, style.shape[2])).T
new_style = np.matmul(A, new_style).T
new_style = np.reshape(new_style, style.shape)
new_style = new_style + b
return new_style
def updateGamma(self):
task_matrices = np.zeros((self.n_tasks, self.num_feats, self.num_feats))
for m in range(self.n_tasks):
rho_vector = rhoFunction(np.array(self.xi[m]))
rho_vector = rho_vector.reshape((1,-1)) # Make rho vector 2D
task_X = self.task_dict[m]['X']
# Note that the transposing doesn't exactly match the paper because our data format is slightly different
rho_matrix = abs(rho_vector) * task_X.T
task_matrices[m,:,:] = np.dot(rho_matrix, task_X)
for k in range(self.K):
inner_sum = np.zeros((self.num_feats,self.num_feats))
for m in range(self.n_tasks):
inner_sum = inner_sum + self.phi[m,k] * task_matrices[m,:,:]
self.gamma[k] = la.inv(la.inv(self.sigma) + 2*inner_sum)
if self.debug:
print "gamma computation {0}".format(k), la.det(la.inv(self.sigma) + 2*inner_sum)
def normal_density_at_zero(m, c):
""" Compute the normal density with given mean and covariance at zero.
Parameters
----------
m : vector
Mean.
c : ndarray
Covariance matrix. Assumption: c is square matrix and its size is
compatible with that of m.
Returns
-------
g : float
Computed density value.
"""
dim = len(m)
g = 1 / ((2 * pi)**(dim / 2) * sqrt(absolute(det(c)))) *\
exp(-1/2 * dot(dot(m, inv(c)), m))
return g
def fit(self, X, y):
"""
Fits a t-Student Process regressor
Parameters
----------
X: np.ndarray, shape=(nsamples, nfeatures)
Training instances to fit the GP.
y: np.ndarray, shape=(nsamples,)
Corresponding continuous target values to `X`.
"""
self.X = X
self.y = y
self.n1 = X.shape[0]
if self.optimize:
self.optHyp(param_key=self.covfunc.parameters, param_bounds=self.covfunc.bounds)
self.K11 = self.covfunc.K(self.X, self.X)
self.beta1 = np.dot(np.dot(self.y.T, inv(self.K11)), self.y)
self.logp = logpdf(self.y, self.nu, mu=np.zeros(self.n1), Sigma=self.K11)
def predict(self, a_hist, t):
"""
This function implements the prediction formula discussed is section 6 (1.59)
It takes a realization for a^N, and the period in which the prediciton is formed
Output: E[abar | a_t, a_{t-1}, ..., a_1, a_0]
"""
N = np.asarray(a_hist).shape[0] - 1
a_hist = np.asarray(a_hist).reshape(N + 1, 1)
V = self.construct_V(N + 1)
aux_matrix = np.zeros((N + 1, N + 1))
aux_matrix[:(t + 1), :(t + 1)] = np.eye(t + 1)
L = la.cholesky(V).T
Ea_hist = la.inv(L) @ aux_matrix @ L @ a_hist
return Ea_hist
def plot4():
# Density 1
Z = gen_gaussian_plot_vals(x_hat, ?)
cs1 = ax.contour(X, Y, Z, 6, colors="black")
ax.clabel(cs1, inline=1, fontsize=10)
# Density 2
M = ? * G.T * linalg.inv(G * ? * G.T + R)
x_hat_F = x_hat + M * (y - G * x_hat)
?_F = ? - M * G * ?
Z_F = gen_gaussian_plot_vals(x_hat_F, ?_F)
cs2 = ax.contour(X, Y, Z_F, 6, colors="black")
ax.clabel(cs2, inline=1, fontsize=10)
# Density 3
new_x_hat = A * x_hat_F
new_? = A * ?_F * A.T + Q
new_Z = gen_gaussian_plot_vals(new_x_hat, new_?)
cs3 = ax.contour(X, Y, new_Z, 6, colors="black")
ax.clabel(cs3, inline=1, fontsize=10)
ax.contourf(X, Y, new_Z, 6, alpha=0.6, cmap=cm.jet)
ax.text(float(y[0]), float(y[1]), r"$y$", fontsize=20, color="black")
# == Choose a plot to generate == #
def factor(self):
""" Factorize the camera matrix into K,R,t as P = K[R|t]. """
# factor first 3*3 part
K,R = linalg.rq(self.P[:,:3])
# make diagonal of K positive
T = diag(sign(diag(K)))
if linalg.det(T) < 0:
T[1,1] *= -1
self.K = dot(K,T)
self.R = dot(T,R) # T is its own inverse
self.t = dot(linalg.inv(self.K),self.P[:,3])
return self.K, self.R, self.t
def __init__(self, suvi, *args):
if list_or_tuple(suvi):
if len(suvi) == 3:
s, u, v = suvi
inv = Identity(s.shape[0])
else:
s, u, v, inv = suvi
else:
s = suvi
u = args[0]
v = args[1]
if len(args)>2:
inv = args[2]
else:
inv = Identity(s.shape[0])
self.s = s
self.u = u
self.v = v
self.inv = inv
def __init__(self, suvi, *args):
if list_or_tuple(suvi):
if len(suvi) == 3:
s, u, v = suvi
inv = Identity(s.shape[0])
else:
s, u, v, inv = suvi
else:
s = suvi
u = args[0]
v = args[1]
if len(args)>2:
inv = args[2]
else:
inv = Identity(s.shape[0])
self.s = s
self.u = u
self.v = v
self.inv = inv
def __init__(self, suvi, *args):
if list_or_tuple(suvi):
if len(suvi) == 3:
s, u, v = suvi
inv = Identity(s.shape[0])
else:
s, u, v, inv = suvi
else:
s = suvi
u = args[0]
v = args[1]
if len(args)>2:
inv = args[2]
else:
inv = Identity(s.shape[0])
self.s = s
self.u = u
self.v = v
self.inv = inv
def __init__(self, suvi, *args):
if list_or_tuple(suvi):
if len(suvi) == 3:
s, u, v = suvi
inv = Identity(s.shape[0])
else:
s, u, v, inv = suvi
else:
s = suvi
u = args[0]
v = args[1]
if len(args)>2:
inv = args[2]
else:
inv = Identity(s.shape[0])
self.s = s
self.u = u
self.v = v
self.inv = inv
def test_optimized_A_and_?_w_make_eye(cls):
""" I = (V^T)(?_w)(V), where A = inv(V^T)
NOTES:
(1) There are two ways to compute ?_w:
?_w = (A)(A^T)
?_w = n / (n-1) * S_w
(2) *** COMPUTE PHI WITH S_w: ?_w = n/(n-1) * S_w. ***
(3) Do NOT use ?_w = (A)(A^T) because that is trivially true:
(V^T)(?_w)(V), where V = inv(A^T), which gives
(inv(A))(A)(A^T)(inv(A^T)) = (I)(I) = I.
"""
tolerance = 1e-13 # Should be smaller than n / (n - 1).
S_w = cls.model.S_w
n = cls.model.n_avg
V = inv(cls.model.A.T)
cls.assertTrue(tolerance < (n / (n - 1)))
?_w = n / (n - 1) * S_w
result = np.matmul(np.matmul(V.T, ?_w), V)
cls.assert_same(result, np.eye(cls.dims), tolerance=tolerance)
def test_calc_A(self):
tolerance = 1e-100
dims = self.dims
?_w = np.diag(np.ones(dims))
W = np.random.randint(0, 9, self.dims ** 2).reshape(dims, dims) + \
np.eye(dims)
n_avg = 9
A_truth = ?_w * n_avg / (n_avg - 1)
A_truth = np.sqrt(A_truth)
A_truth = np.matmul(np.linalg.inv(W).T, A_truth)
A_model = self.model.calc_A(n_avg, ?_w, W)
self.assert_same(A_model, A_truth, tolerance=tolerance)
self.assert_invertible(self.model.W)
def compute(self,signal):
check_MD(signal.X)
N=signal.N
M,L=signal.X.shape
sigma2=signal.sigma2
H=np.matrix(signal.H)
X=np.matrix(signal.X)
Rx=X*X.H/L
PH=np.eye(N)-H*lg.pinv(H)
D=np.matrix(np.diag(1j*np.arange(N)))
M=H.H*D.H*PH*D*H
CRB=(sigma2/2)*lg.inv(np.real(np.multiply(M,Rx)))
self.w=np.diag(CRB)
return self
def compute_criterion(self,y):
self.N=len(y)
#construct matrices
A=np.matrix(self.A)
b=np.matrix(self.b).T
H=np.matrix(self.H)
x=lg.inv(H.T*H)*H.T*y #estimation of x
if self.estimate_sigma2 is True:
r,p=self.A.shape
coef=(self.N-p)/r
den=lg.norm(y-H*x)**2
else:
den=self.sigma2
coef=1
term1=A*x-b
num=term1.T*lg.inv(A*lg.inv(H.T*H)*A.T)*term1
self.criterion=coef*num/den ## See page 274 / 345
def test_func_realize(radii):
sys = Alpha(0.1)
T = np.asarray([[1., 2.], [0, -1.]])
for Tinv in (None, inv(T)):
realize_result = _realize(sys, radii, T, Tinv)
assert realize_result.sys is sys
assert np.allclose(inv(realize_result.T), realize_result.Tinv)
rsys = realize_result.realization
assert ss_equal(rsys, sys.transform(realize_result.T))
# Check that the state vector are related by T
length = 1000
dt = 0.001
x_old = np.asarray([sub.impulse(length, dt) for sub in sys])
x_new = np.asarray([sub.impulse(length, dt) for sub in rsys])
r = np.atleast_2d(np.asarray(radii).T).T
assert np.allclose(np.dot(T, x_new * r), x_old)
def test_identity(radii):
sys = Alpha(0.1)
identity = Identity()
assert repr(identity) == "Identity()"
I = np.eye(len(sys))
realize_result = identity(sys, radii)
assert realize_result.sys is sys
assert np.allclose(realize_result.T, I * radii)
assert np.allclose(realize_result.Tinv, inv(I * radii))
rsys = realize_result.realization
assert ss_equal(rsys, sys.transform(realize_result.T))
# Check that it's still the same system, even though different matrices
assert sys_equal(sys, rsys)
if radii == 1:
assert ss_equal(sys, rsys)
else:
assert not np.allclose(sys.B, rsys.B)
assert not np.allclose(sys.C, rsys.C)
# Check that the state vectors have scaled power
assert np.allclose(state_norm(sys) / radii, state_norm(rsys))
def _realize(sys, radii, T, Tinv=None):
"""Helper function for producing a RealizerResult."""
sys = LinearSystem(sys)
r = np.asarray(radii, dtype=np.float64)
if r.ndim == 0:
r = np.ones(len(sys)) * r
elif r.ndim > 1:
raise ValueError("radii (%s) must be a 1-dim array or scalar" % (
radii,))
elif len(r) != len(sys):
raise ValueError("radii (%s) length must match state dimension %d" % (
radii, len(sys)))
T = T * r[None, :]
if Tinv is None: # this needs to be computed eventually anyways
Tinv = inv(T)
else:
Tinv = Tinv / r[:, None]
return RealizerResult(sys, T, Tinv, sys.transform(T, Tinv))
def decompose(P):
M = P[:3, :3]
T = P[:3, 3]
K, R = scipy.linalg.rq(M)
for i in range(2):
if K[i,i] < 0:
K[:, i] *= -1
R[i, :] *= -1
if K[2,2] > 0:
K[:, 2] *= -1
R[2, :] *= -1
if det(R) < 0:
R *= -1
T = linalg.inv(dot(K, -R)).dot(T.reshape((3, 1)))
K /= -K[2,2]
return K, R, T
def _fwd_bem_multi_solution(solids, gamma, nps):
"""Do multi surface solution
* Invert I - solids/(2*M_PI)
* Take deflation into account
* The matrix is destroyed after inversion
* This is the general multilayer case
"""
pi2 = 1.0 / (2 * np.pi)
n_tot = np.sum(nps)
assert solids.shape == (n_tot, n_tot)
nsurf = len(nps)
defl = 1.0 / n_tot
# Modify the matrix
offsets = np.cumsum(np.concatenate(([0], nps)))
for si_1 in range(nsurf):
for si_2 in range(nsurf):
mult = pi2 if gamma is None else pi2 * gamma[si_1, si_2]
slice_j = slice(offsets[si_1], offsets[si_1 + 1])
slice_k = slice(offsets[si_2], offsets[si_2 + 1])
solids[slice_j, slice_k] = defl - solids[slice_j, slice_k] * mult
solids += np.eye(n_tot)
return linalg.inv(solids, overwrite_a=True)
def _check_infos_trans(infos):
"""XXX this goes to tests later, currently not used"""
chan_max_idx = np.argmax([c['nchan'] for c in infos])
chan_template = infos[chan_max_idx]['ch_names']
channels = [c['ch_names'] for c in infos]
common_channels = set(chan_template).intersection(*channels)
common_chs = [[c['chs'][c['ch_names'].index(ch)] for ch in common_channels]
for c in infos]
dev_ctf_trans = [i['dev_ctf_t']['trans'] for i in infos]
cns = [[c['ch_name'] for c in cc] for cc in common_chs]
for cn1, cn2 in itt.combinations(cns, 2):
assert cn1 == cn2
# BTI stores data in head coords, as a consequence the coordinates
# change across run, we apply the ctf->ctf_head transform here
# to check that all transforms are correct.
cts = [np.array([linalg.inv(_loc_to_coil_trans(c['loc'])).dot(t)
for c in cc])
for t, cc in zip(dev_ctf_trans, common_chs)]
for ct1, ct2 in itt.combinations(cts, 2):
np.testing.assert_array_almost_equal(ct1, ct2, 12)
def test_graph_lasso_cv(random_state=1):
# Sample data from a sparse multivariate normal
dim = 5
n_samples = 6
random_state = check_random_state(random_state)
prec = make_sparse_spd_matrix(dim, alpha=.96,
random_state=random_state)
cov = linalg.inv(prec)
X = random_state.multivariate_normal(np.zeros(dim), cov, size=n_samples)
# Capture stdout, to smoke test the verbose mode
orig_stdout = sys.stdout
try:
sys.stdout = StringIO()
# We need verbose very high so that Parallel prints on stdout
GraphLassoCV(verbose=100, alphas=5, tol=1e-1).fit(X)
finally:
sys.stdout = orig_stdout
# Smoke test with specified alphas
GraphLassoCV(alphas=[0.8, 0.5], tol=1e-1, n_jobs=1).fit(X)
def get_precision(self):
"""Compute data precision matrix with the generative model.
Equals the inverse of the covariance but computed with
the matrix inversion lemma for efficiency.
Returns
-------
precision : array, shape=(n_features, n_features)
Estimated precision of data.
"""
n_features = self.components_.shape[1]
# handle corner cases first
if self.n_components_ == 0:
return np.eye(n_features) / self.noise_variance_
if self.n_components_ == n_features:
return linalg.inv(self.get_covariance())
# Get precision using matrix inversion lemma
components_ = self.components_
exp_var = self.explained_variance_
exp_var_diff = np.maximum(exp_var - self.noise_variance_, 0.)
precision = np.dot(components_, components_.T) / self.noise_variance_
precision.flat[::len(precision) + 1] += 1. / exp_var_diff
precision = np.dot(components_.T,
np.dot(linalg.inv(precision), components_))
precision /= -(self.noise_variance_ ** 2)
precision.flat[::len(precision) + 1] += 1. / self.noise_variance_
return precision
def updateTheta(self):
for k in range(self.K):
inner_sum = np.zeros((1,self.num_feats))
for m in range(self.n_tasks):
inner_sum = inner_sum + self.phi[m,k] * np.atleast_2d(self.task_vectors[m,:])
self.theta[k,:] = (np.dot(self.gamma[k],(np.dot(la.inv(self.sigma),self.mu.T) + inner_sum.T) )).T
LSFIR.py 文件源码
项目：Least-Squared-Error-Based-FIR-Filters
作者: fourier-being
项目源码
文件源码
阅读 28
收藏 0
点赞 0
评论 0
def lpfls2notch(N,wp,ws,wn1,wn2,W):
M = (N-1)/2
nq = np.arange(0,2*M+1)
nb = np.arange(0,M+1)
q = (wp/np.pi)*np.sinc((wp/np.pi)*nq) - W*(ws/np.pi)*np.sinc((ws/np.pi)*nq)
b = (wp/np.pi)*np.sinc((wp/np.pi)*nb)
q[0] = wp/np.pi + W*(1-ws/np.pi) # since sin(pi*n)/pi*n = 1, not 0
b = np.asmatrix(b)
b = b.transpose()
Q1 = ln.toeplitz(q[0:M+1])
Q2 = ln.hankel(q[0:M+1],q[M:])
Q = Q1+Q2
G1 = np.cos(wn1*nb)
G2 = np.cos(wn2*nb)
G = np.matrix([G1,G2])
d = np.array([0,0])
d = np.asmatrix(d)
d = d.transpose()
c = np.asmatrix(ln.solve(Q,b))
mu = ln.solve(G*ln.inv(Q)*G.transpose(),G*c - d)
a = c - ln.solve(Q,G.transpose()*mu)
h = np.zeros(N)
for i in nb:
h[i] = 0.5*a[M-i]
h[N-1-i] = h[i]
h[M] = 2*h[M]
hmax = max(np.absolute(h))
for i in nq:
h[i] = (8191/hmax)*h[i]
return h
LSFIR.py 文件源码
项目：Least-Squared-Error-Based-FIR-Filters
作者: fourier-being
项目源码
文件源码
阅读 26
收藏 0
点赞 0
评论 0
def lpfls1notch(N,wp,ws,wn1,W):
M = (N-1)/2
nq = np.arange(0,2*M+1)
nb = np.arange(0,M+1)
q = (wp/np.pi)*np.sinc((wp/np.pi)*nq) - W*(ws/np.pi)*np.sinc((ws/np.pi)*nq)
b = (wp/np.pi)*np.sinc((wp/np.pi)*nb)
q[0] = wp/np.pi + W*(1-ws/np.pi) # since sin(pi*n)/pi*n = 1, not 0
b = np.asmatrix(b)
b = b.transpose()
Q1 = ln.toeplitz(q[0:M+1])
Q2 = ln.hankel(q[0:M+1],q[M:])
Q = Q1+Q2
G1 = np.cos(wn1*nb)
G = np.matrix([G1])
d = np.array([0])
d = np.asmatrix(d)
c = np.asmatrix(ln.solve(Q,b))
mu = ln.solve(G*ln.inv(Q)*G.transpose(),G*c - d)
a = c - ln.solve(Q,G.transpose()*mu)
h = np.zeros(N)
for i in nb:
h[i] = 0.5*a[M-i]
h[N-1-i] = h[i]
h[M] = 2*h[M]
hmax = max(np.absolute(h))
for i in nq:
h[i] = (8191/hmax)*h[i]
return h
def __init__(self, K,R,t):
self.K = K
self.R = R
self.t = t
self.invR = self.R.T
self.P = np.matmul(self.K,
np.concatenate((self.R, self.t.reshape((3,1))), axis=1))
self.invP3 = la.inv(self.P[:3,:3])
self.center = np.matmul(-self.invP3, self.P[:3,3])
def triangulate(cam1, cam2, p1, p2):
p1 = np.asarray([p1[0],p1[1],1])
p2 = np.asarray([p2[0],p2[1],1])
c1, v1 = cam_center_vector(cam1, p1)
c2, v2 = cam_center_vector(cam2, p2)
t = c2 - c1
v3 = np.cross(v1, v2)
X = np.stack((v1, v3, -v2), axis=1)
alpha = np.matmul(la.inv(X), t)
output = c1 + v1 * alpha[0] + alpha[1]*0.5*v3
return output
def calc_deconv_matrix(matrix_vector_dab_he):
"""
Custom calculated matrix of lab's stains DAB + Hematoxylin
The raw matrix was moved to the global scope before main() function as a constant
"""
matrix_vector_dab_he[2, :] = np.cross(matrix_vector_dab_he[0, :], matrix_vector_dab_he[1, :])
matrix_dh = linalg.inv(matrix_vector_dab_he)
return matrix_dh
