def get_inversion(self):
diag_blocks = self.blocks
inverted_blocks = np.zeros(len(diag_blocks), dtype=object)
for i in xrange(len(diag_blocks)):
inverted_blocks[i] = np.linalg.inv(diag_blocks[i])
inverted_diag = block_diag(*inverted_blocks)
return inverted_diag
python类block_diag()的实例源码
def calculate_expanded_base_transformation_matrix(src_base, dst_base, src_order, dst_order, use_eye=False):
"""
constructs a transformation matrix from basis given by 'src_base' to basis given by 'dst_base' that also
transforms all temporal derivatives of the given weights.
:param src_base: the source basis, given by an array of BaseFractions
:param dst_base: the destination basis, given by an array of BaseFractions
:param src_order: temporal derivative order available in src
:param dst_order: temporal derivative order needed in dst
:param use_eye: use identity as base transformation matrix
:return: transformation matrix as 2d np.ndarray
"""
if src_order < dst_order:
raise ValueError("higher derivative order needed than provided!")
# build core transformation
if use_eye:
core_transformation = np.eye(src_base.size)
else:
core_transformation = calculate_base_transformation_matrix(src_base, dst_base)
# build block matrix
part_transformation = block_diag(*[core_transformation for i in range(dst_order + 1)])
complete_transformation = np.hstack([part_transformation] + [np.zeros((part_transformation.shape[0], src_base.size))
for i in range(src_order - dst_order)])
return complete_transformation
def Identified_Model(y, t, library, estimator) :
'''
Simulates the model from Sparse identification.
Inputs
------
library: library object used in the sparse identification
(e.g. poly_lib = PolynomialFeatures(degree=3) )
estimator: estimator object obtained from the sparse identification
Output
------
dy : numpy array object containing the derivatives evaluated using the
model identified from sparse regression.
'''
dy = np.zeros_like(y)
lib = library.fit_transform(y.reshape(1,-1))
Theta = block_diag(lib, lib, lib)
dy = Theta.dot(estimator.coef_)
return dy
def mtx_fri2visi_ri(M, p_mic_x, p_mic_y, D1, D2):
"""
build the matrix that maps the Fourier series to the visibility in terms of
REAL-VALUED entries only. (matrix size double)
:param M: the Fourier series expansion is limited from -M to M
:param p_mic_x: a vector that contains microphones x coordinates
:param p_mic_y: a vector that contains microphones y coordinates
:param D1: expansion matrix for the real-part
:param D2: expansion matrix for the imaginary-part
:return:
"""
return np.dot(cpx_mtx2real(mtx_freq2visi(M, p_mic_x, p_mic_y)),
linalg.block_diag(D1, D2))
def mtx_updated_G_multiband(phi_recon, M, mtx_amp2visi_ri,
mtx_fri2visi_ri, num_bands):
"""
Update the linear transformation matrix that links the FRI sequence to the
visibilities by using the reconstructed Dirac locations.
:param phi_recon: the reconstructed Dirac locations (azimuths)
:param M: the Fourier series expansion is between -M to M
:param p_mic_x: a vector that contains microphones' x-coordinates
:param p_mic_y: a vector that contains microphones' y-coordinates
:param mtx_freq2visi: the linear mapping from Fourier series to visibilities
:return:
"""
L = 2 * M + 1
ms_half = np.reshape(np.arange(-M, 1, step=1), (-1, 1), order='F')
phi_recon = np.reshape(phi_recon, (1, -1), order='F')
mtx_amp2freq = np.exp(-1j * ms_half * phi_recon) # size: (M + 1) x K
mtx_amp2freq_ri = np.vstack((mtx_amp2freq.real, mtx_amp2freq.imag[:-1, :])) # size: (2M + 1) x K
mtx_fri2amp_ri = linalg.lstsq(mtx_amp2freq_ri, np.eye(L))[0]
# projection mtx_freq2visi to the null space of mtx_fri2amp
mtx_null_proj = np.eye(L) - np.dot(mtx_fri2amp_ri.T,
linalg.lstsq(mtx_fri2amp_ri.T, np.eye(L))[0])
G_updated = np.dot(mtx_amp2visi_ri,
linalg.block_diag(*([mtx_fri2amp_ri] * num_bands))
) + \
np.dot(mtx_fri2visi_ri,
linalg.block_diag(*([mtx_null_proj] * num_bands))
)
return G_updated
def getClique(N=100, K=4):
from scipy.linalg import block_diag
b = []
for k in range(K):
n = N // K
b.append(np.ones((n,n), int))
C = block_diag(*b)
return C
### @Issue42: fronteNetwork should be imported fron frontend
### =====> : resolve this with @class_method (from_hardrive etc...)
def __init__(self, systems, dt=None, elementwise=False, method='zoh'):
if not is_iterable(systems) or isinstance(systems, LinearSystem):
systems = [systems]
self.systems = systems
self.dt = dt
self.elementwise = elementwise
self.A = []
self.B = []
self.C = []
self.D = []
for sys in systems:
sys = LinearSystem(sys)
if dt is not None:
sys = cont2discrete(sys, dt, method=method)
elif sys.analog:
raise ValueError(
"system (%s) must be digital if not given dt" % sys)
A, B, C, D = sys.ss
self.A.append(A)
self.B.append(B)
self.C.append(C)
self.D.append(D)
# TODO: If all of the synapses are single order, than A is diagonal
# and so np.dot(self.A, self._x) is trivial. But perhaps
# block_diag is already optimized for this.
# Note: ideally we could put this into CCF to reduce the A mapping
# to a single dot product and a shift operation. But in general
# since this is MIMO it is not controllable from a single input.
# Instead we might want to consider balanced reduction to
# improve efficiency.
self.A = block_diag(*self.A)
self.B = block_diag(*self.B) if elementwise else np.vstack(self.B)
self.C = block_diag(*self.C)
self.D = block_diag(*self.D) if elementwise else np.vstack(self.D)
# TODO: shape validation
self._x = np.zeros(len(self.A))[:, None]
def _build_sparse_mtx():
# Create 3 topics and each topic has 3 distinct words.
# (Each word only belongs to a single topic.)
n_topics = 3
block = n_topics * np.ones((3, 3))
blocks = [block] * n_topics
X = block_diag(*blocks)
X = csr_matrix(X)
return (n_topics, X)
def stiffness_global(x, y, z, E, G, Ax, Iz=0, Iy=0, Ay=0, Az=0, theta=0, J=0):
t = block_diag(*(transformMatrix(x, y, z, theta),) * 2)
r = t.transpose()
for m in stiffness_local(norm((x, y, z)), E, G, Ax, Iz, Iy, Ay, Az, J):
yield dot(dot(t, m), r)
def test_kernel_backend(self):
# set up signal parameter
selection = Selection(selections.by_backend)
log10_sigma = parameter.Uniform(-10, -5)
log10_lam = parameter.Uniform(np.log10(86400), np.log10(1500*86400))
basis = create_quant_matrix(dt=7*86400)
prior = se_kernel(log10_sigma=log10_sigma, log10_lam=log10_lam)
se = gs.BasisGP(prior, basis, selection=selection, name='se')
sem = se(self.psr)
# parameters
log10_sigmas = [-7, -6, -6.4, -8.5]
log10_lams = [8.3, 7.4, 6.8, 5.6]
params = {'B1855+09_se_430_ASP_log10_lam': log10_lams[0],
'B1855+09_se_430_ASP_log10_sigma': log10_sigmas[0],
'B1855+09_se_430_PUPPI_log10_lam': log10_lams[1],
'B1855+09_se_430_PUPPI_log10_sigma': log10_sigmas[1],
'B1855+09_se_L-wide_ASP_log10_lam': log10_lams[2],
'B1855+09_se_L-wide_ASP_log10_sigma': log10_sigmas[2],
'B1855+09_se_L-wide_PUPPI_log10_lam': log10_lams[3],
'B1855+09_se_L-wide_PUPPI_log10_sigma': log10_sigmas[3]}
# get the basis
bflags = self.psr.backend_flags
Fmats, fs, phis = [], [], []
for ct, flag in enumerate(np.unique(bflags)):
mask = bflags == flag
U, avetoas = create_quant_matrix(self.psr.toas[mask], dt=7*86400)
Fmats.append(U)
fs.append(avetoas)
phis.append(se_kernel(avetoas, log10_sigma=log10_sigmas[ct],
log10_lam=log10_lams[ct]))
nf = sum(F.shape[1] for F in Fmats)
U = np.zeros((len(self.psr.toas), nf))
K = sl.block_diag(*phis)
Kinv = np.linalg.inv(K)
nftot = 0
for ct, flag in enumerate(np.unique(bflags)):
mask = bflags == flag
nn = Fmats[ct].shape[1]
U[mask, nftot:nn+nftot] = Fmats[ct]
nftot += nn
msg = 'Kernel basis incorrect for backend signal.'
assert np.allclose(U, sem.get_basis(params)), msg
# spectrum test
msg = 'Kernel incorrect for backend signal.'
assert np.allclose(sem.get_phi(params), K), msg
# inverse spectrum test
msg = 'Kernel inverse incorrect for backend signal.'
assert np.allclose(sem.get_phiinv(params), Kinv), msg
def get_phi(self, params, cliques=False):
phis = [signalcollection.get_phi(params) for
signalcollection in self._signalcollections]
# if we found common signals, we'll return a big phivec matrix,
# otherwise a list of phivec vectors (some of which possibly None)
if self._commonsignals:
if np.any([phi.ndim == 2 for phi in phis if phi is not None]):
# if we have any dense matrices,
Phi = sl.block_diag(*[np.diag(phi) if phi.ndim == 1 else phi
for phi in phis
if phi is not None])
else:
Phi = np.diag(np.concatenate([phi for phi in phis
if phi is not None]))
# get a dictionary of slices locating each pulsar in Phi matrix
slices = self._get_slices(phis)
# self._cliques is a vector of the same size as the Phi matrix
# for each Phi index i, self._cliques[i] is -1 if row/column
# belong to no clique, or it gives the clique number otherwise
if cliques:
self._resetcliques(Phi.shape[0])
self._setpulsarcliques(slices, phis)
# iterate over all common signal classes
for csclass, csdict in self._commonsignals.items():
# first figure out which indices are used in this common signal
# and update the clique index
if cliques:
self._setcliques(slices, csdict)
# now iterate over all pairs of common signal instances
pairs = itertools.combinations(csdict.items(),2)
for (cs1, csc1), (cs2, csc2) in pairs:
crossdiag = csclass.get_phicross(cs1, cs2, params)
block1, idx1 = slices[csc1], csc1._idx[cs1]
block2, idx2 = slices[csc2], csc2._idx[cs2]
Phi[block1,block2][idx1,idx2] += crossdiag
Phi[block2,block1][idx2,idx1] += crossdiag
return Phi
else:
return phis
def model_training(K, y, C, n):
"""
train the model
:param K: kernel matrix
:param y: labels which has the same length as all functions f
:param C: num of categories
:param n: num of training data
:return:
"""
tolerance = 0.0001
step_size = 0.00001
s = 0.0005
# initialization
f = np.zeros((C*n,)) # initialize f=0(unbiased) which is an constant=0 function and means no GP prior in this case
block_K = [K[i * n:(i + 1) * n, i * n:(i + 1) * n] for i in range(K.shape[0] / n)]
# Newton iteration
for j in range(100):
pi_vector, pi_matrix = compute_pi(f, C, n)
D = np.zeros((C * n, C * n))
np.fill_diagonal(D, pi_vector)
savetxt_compact('D.txt',D)
block_D = [D[i * n:(i + 1) * n, i * n:(i + 1) * n] for i in range(D.shape[0] / n)]
savetxt_compact('block_D0.txt', block_D[0])
savetxt_compact('block_D2.txt', block_D[2])
E_c_sum = np.zeros((n, n))
for c in range(C):
L = np.linalg.cholesky(np.eye(n) + np.dot(np.sqrt(block_D[c]), np.dot(block_K[c], np.sqrt(block_D[c]))))
L_inv = np.linalg.inv(L)
E_c_part = np.dot(np.sqrt(block_D[c]), np.dot(L_inv.T, np.dot(L_inv, np.sqrt(block_D[c]))))
# create a block diagonal matrix E
if c == 0:
E = E_c_part
else:
E = block_diag(E, E_c_part)
E_c_sum += E_c_part
L_whole = np.linalg.cholesky(np.eye(C*n) + np.dot(np.sqrt(D), np.dot(K, np.sqrt(D))))
L_whole_inv = np.linalg.inv(L_whole)
E = np.dot(np.sqrt(D), np.dot(L_whole_inv.T, np.dot(L_whole_inv, np.sqrt(D))))
#E = np.dot(np.sqrt(D), np.dot(np.linalg.inv(np.eye(C*n) + np.dot(np.sqrt(D), np.dot(K, np.sqrt(D)))), np.sqrt(D)))
R = np.dot(np.linalg.inv(D), pi_matrix)
#M = np.linalg.cholesky(E_c_sum)
M = np.linalg.cholesky(np.dot(R.T, np.dot(E, R)))
W = D - np.dot(pi_matrix, pi_matrix.T)
L_K = np.linalg.cholesky(s * np.eye(C*n) + K)
L_K_inv = np.linalg.inv(L_K)
b = np.dot((1-step_size) * np.dot(L_K_inv.T, L_K_inv) + W, f) + y - pi_vector
c = np.dot(E, np.dot(K, b))
M_inv = np.linalg.inv(M)
a = b - c + np.dot(E, np.dot(R, np.dot(M_inv.T, np.dot(M_inv, np.dot(R.T, c)))))
f_new = np.dot(K, a)
error = np.sqrt(np.sum((f_new - f) ** 2))
f = f_new
print `j + 1` + "th iteration, error:" + `error`
if error <= tolerance:
print "The function has already converged after " + `j + 1` + " iterations!"
print "The error is " + `error`
print "training end!"
break
def householder(self):
r"""
Implementation of Householder reflections method to performing QR
decomposition.
Returns
-------
qr : tuple
Returns a tuple containing the orthogonal matrix Q and the upper-triangular
matrix R resulting from QR decomposition.
Notes
-----
The Householder reflection approach to QR decomposition is the more common approach
due to its numerical stability compared to Gram-Schmidt and its relative speed to
Givens rotations. The orthogonal matrix :math:`Q` is defined as successive Householder
matrices :math:`H_1 \cdots H_n` while :math:`R` is upper triangular, defined as
:math:`R = Q^T A`.
Householder matrices :math:`H` are defined as:
.. math::
H = I - 2vv^T
References
----------
Golub, G., & Van Loan, C. (2013). Matrix computations (3rd ed.). Baltimore (MD): Johns Hopkins U.P.
Householder transformation. (2017, March 19). In Wikipedia, The Free Encyclopedia.
From https://en.wikipedia.org/w/index.php?title=Householder_transformation&oldid=771169379
Trefethen, L., & Bau, D. (1997). Numerical linear algebra (1st ed.). Philadelphia: SIAM.
"""
h = []
r = self.x.copy()
if self.m > self.n:
c = self.n
else:
c = self.m
for j in np.arange(c):
hj = _householder_mat(r[j:self.m, j])
if j > 0:
hj = block_diag(np.eye(j), hj)
r = np.dot(hj, r)
h.append(hj)
self.q = reduce(np.dot, reversed(h))[0:self.n].T
r = np.array(r)[0:self.n]
qr = (self.q, r)
return qr
