def test_covariances_and_eigenvalues(self):
reader = FeatureReader(self.trajnames, self.temppdb)
for tau in [1, 10, 100, 1000, 2000]:
trans = tica(lag=tau, dim=self.dim, kinetic_map=False)
trans.data_producer = reader
log.info('number of trajectories reported by tica %d' % trans.number_of_trajectories())
trans.parametrize()
data = trans.get_output()
log.info('max. eigenvalue: %f' % np.max(trans.eigenvalues))
self.assertTrue(np.all(trans.eigenvalues <= 1.0))
# check ICs
check = tica(data=data, lag=tau, dim=self.dim)
np.testing.assert_allclose(np.eye(self.dim), check.cov, atol=1e-8)
np.testing.assert_allclose(check.mean, 0.0, atol=1e-8)
ic_cov_tau = np.zeros((self.dim, self.dim))
ic_cov_tau[np.diag_indices(self.dim)] = trans.eigenvalues
np.testing.assert_allclose(ic_cov_tau, check.cov_tau, atol=1e-8)
python类diag_indices()的实例源码
test_featurereader_and_tica_projection.py 文件源码
项目：coordinates
作者: markovmodel
项目源码
文件源码
阅读 48
收藏 0
点赞 0
评论 0
def update_hypers(self, params):
for i in range(self.nolayers):
layeri = self.layers[i]
Mi = layeri.M
Dini = layeri.Din
Douti = layeri.Dout
layeri.ls.set_value(params['ls'][i])
layeri.sf.set_value(params['sf'][i])
layeri.sn.set_value(params['sn'][i])
triu_ind = np.triu_indices(Mi)
diag_ind = np.diag_indices(Mi)
for d in range(Douti):
layeri.zu[d].set_value(params['zu'][i][d])
theta_m_d = params['eta2'][i][d]
theta_R_d = params['eta1_R'][i][d]
R = np.zeros((Mi, Mi))
R[triu_ind] = theta_R_d.reshape(theta_R_d.shape[0], )
R[diag_ind] = np.exp(R[diag_ind])
layeri.theta_1_R[d] = R
layeri.theta_1[d] = np.dot(R.T, R)
layeri.theta_2[d] = theta_m_d
def _passthrough_delay(m, c):
"""Analytically derived state-space when p = q = m.
We use this because it is numerically stable for high m.
"""
j = np.arange(1, m+1, dtype=np.float64)
u = (m + j) * (m - j + 1) / (c * j)
A = np.zeros((m, m))
B = np.zeros((m, 1))
C = np.zeros((1, m))
D = np.zeros((1,))
A[0, :] = B[0, 0] = -u[0]
A[1:, :-1][np.diag_indices(m-1)] = u[1:]
D[0] = (-1)**m
C[0, np.arange(m) % 2 == 0] = 2*D[0]
return LinearSystem((A, B, C, D), analog=True)
def _proper_delay(q, c):
"""Analytically derived state-space when p = q - 1.
We use this because it is numerically stable for high q
and doesn't have a passthrough.
"""
j = np.arange(q, dtype=np.float64)
u = (q + j) * (q - j) / (c * (j + 1))
A = np.zeros((q, q))
B = np.zeros((q, 1))
C = np.zeros((1, q))
D = np.zeros((1,))
A[0, :] = -u[0]
B[0, 0] = u[0]
A[1:, :-1][np.diag_indices(q-1)] = u[1:]
C[0, :] = (j + 1) / float(q) * (-1) ** (q - 1 - j)
return LinearSystem((A, B, C, D), analog=True)
def convert_solution(z, Cs):
if issparse(Cs):
Cs = Cs.toarray()
M = Cs.shape[0]
x = z[0:M]
y = z[M:]
w=np.exp(y)
pi=w/w.sum()
X=pi[:,np.newaxis]*x[np.newaxis,:]
Y=X+np.transpose(X)
denom=Y
enum=Cs*np.transpose(pi)
P=enum/denom
ind=np.diag_indices(Cs.shape[0])
P[ind]=0.0
rowsums=P.sum(axis=1)
P[ind]=1.0-rowsums
return pi, P
###############################################################################
# Objective, Gradient, and Hessian
###############################################################################
def get_connectivity_matrix_nodiag(self):
"""
Returns a similar matrix as in Ontology.get_connectivity_matrix(),
but the diagonal of the matrix is 0.
Note: !!!!!!!!!!!!!!!!!!!!!!!!
d[a, a] == 0 instead of 1
"""
if not hasattr(self, 'd_nodiag'):
d = self.get_connectivity_matrix()
self.d_nodiag = d.copy()
self.d_nodiag[np.diag_indices(self.d_nodiag.shape[0])] = 0
assert not np.isfortran(self.d_nodiag)
return self.d_nodiag
def corrsort(features, use_tsp=False):
"""
Given a features array, one feature per axis=0 entry, return axis=0 indices
such that adjacent indices correspond to features that are correlated.
cf. Traveling Salesman Problem. Not an optimal solution.
use_tsp:
Use tsp solver. See tsp_solver.greedy module that is used for this.
Slows run-time considerably: O(N^4) computation, O(N^2) memory.
Without use_tsp, both computation and memory are O(N^2).
"""
correlations = np.ma.corrcoef(arrays.plane(features))
if use_tsp: return tsp.solve_tsp(-correlations)
size = features.shape[0]
correlations.mask[np.diag_indices(size)] = True
# initialize results with the pair with the highest correlations.
largest = np.argmax(correlations)
results = [int(largest / size), largest % size]
correlations.mask[:, results[0]] = True
while len(results) < size:
correlations.mask[:, results[-1]] = True
results.append(np.argmax(correlations[results[-1]]))
return results
def _get_batch_diagonal_cpu(array):
batch_size, m, n = array.shape
assert m == n
rows, cols = np.diag_indices(n)
return array[:, rows, cols]
def _set_batch_diagonal_cpu(array, diag_val):
batch_size, m, n = array.shape
assert m == n
rows, cols = np.diag_indices(n)
array[:, rows, cols] = diag_val
def _diagonal_idx_array(batch_size, n):
idx_offsets = np.arange(
start=0, stop=batch_size * n * n, step=n * n, dtype=np.int32).reshape(
(batch_size, 1))
idx = np.ravel_multi_index(
np.diag_indices(n), (n, n)).reshape((1, n)).astype(np.int32)
return cuda.to_gpu(idx + idx_offsets)
def check_forward(self, diag_data, non_diag_data):
diag = chainer.Variable(diag_data)
non_diag = chainer.Variable(non_diag_data)
y = lower_triangular_matrix(diag, non_diag)
correct_y = numpy.zeros(
(self.batch_size, self.n, self.n), dtype=numpy.float32)
tril_rows, tril_cols = numpy.tril_indices(self.n, -1)
correct_y[:, tril_rows, tril_cols] = cuda.to_cpu(non_diag_data)
diag_rows, diag_cols = numpy.diag_indices(self.n)
correct_y[:, diag_rows, diag_cols] = cuda.to_cpu(diag_data)
gradient_check.assert_allclose(correct_y, cuda.to_cpu(y.data))
def kth_diag_indices(n, k):
"""
Return indices of bins k steps away from the diagonal.
(from http://stackoverflow.com/questions/10925671/numpy-k-th-diagonal-indices)
"""
rows, cols = np.diag_indices(n)
if k < 0:
return rows[:k], cols[-k:]
elif k > 0:
return rows[k:], cols[:-k]
else:
return rows, cols
def compute_posterior_grad_u(self, dmu, dSu):
# return grads wrt u params and Kuuinv
triu_ind = np.triu_indices(self.M)
diag_ind = np.diag_indices(self.M)
if self.nat_param:
dSu_via_m = np.einsum('da,db->dab', dmu, self.theta_2)
dSu += dSu_via_m
dSuinv = - np.einsum('dab,dbc,dce->dae', self.Su, dSu, self.Su)
dKuuinv = np.sum(dSuinv, axis=0)
dtheta1 = dSuinv
deta2 = np.einsum('dab,db->da', self.Su, dmu)
else:
deta2 = dmu
dtheta1 = dSu
dKuuinv = 0
dtheta1T = np.transpose(dtheta1, [0, 2, 1])
dtheta1_R = np.einsum('dab,dbc->dac', self.theta_1_R, dtheta1 + dtheta1T)
deta1_R = np.zeros([self.Dout, self.M * (self.M + 1) / 2])
for d in range(self.Dout):
dtheta1_R_d = dtheta1_R[d, :, :]
theta1_R_d = self.theta_1_R[d, :, :]
dtheta1_R_d[diag_ind] = dtheta1_R_d[diag_ind] * theta1_R_d[diag_ind]
dtheta1_R_d = dtheta1_R_d[triu_ind]
deta1_R[d, :] = dtheta1_R_d.reshape((dtheta1_R_d.shape[0], ))
return deta1_R, deta2, dKuuinv
def get_hypers(self, key_suffix=''):
"""Summary
Args:
key_suffix (str, optional): Description
Returns:
TYPE: Description
"""
params = {}
M = self.M
Din = self.Din
Dout = self.Dout
params['ls' + key_suffix] = self.ls
params['sf' + key_suffix] = self.sf
triu_ind = np.triu_indices(M)
diag_ind = np.diag_indices(M)
params_eta2 = self.theta_2
params_eta1_R = np.zeros((Dout, M * (M + 1) / 2))
params_zu_i = self.zu
for d in range(Dout):
Rd = np.copy(self.theta_1_R[d, :, :])
Rd[diag_ind] = np.log(Rd[diag_ind])
params_eta1_R[d, :] = np.copy(Rd[triu_ind])
params['zu' + key_suffix] = self.zu
params['eta1_R' + key_suffix] = params_eta1_R
params['eta2' + key_suffix] = params_eta2
return params
def update_hypers(self, params, key_suffix=''):
"""Summary
Args:
params (TYPE): Description
key_suffix (str, optional): Description
"""
M = self.M
self.ls = params['ls' + key_suffix]
self.sf = params['sf' + key_suffix]
triu_ind = np.triu_indices(M)
diag_ind = np.diag_indices(M)
zu = params['zu' + key_suffix]
self.zu = zu
for d in range(self.Dout):
theta_m_d = params['eta2' + key_suffix][d, :]
theta_R_d = params['eta1_R' + key_suffix][d, :]
R = np.zeros((M, M))
R[triu_ind] = np.copy(theta_R_d.reshape(theta_R_d.shape[0], ))
R[diag_ind] = np.exp(R[diag_ind])
self.theta_1_R[d, :, :] = R
self.theta_1[d, :, :] = np.dot(R.T, R)
self.theta_2[d, :] = theta_m_d
# update Kuu given new hypers
self.compute_kuu()
# compute mu and Su for each layer
self.update_posterior()
def compute_cav_grad_u(self, dmu, dSu, alpha):
# return grads wrt u params and Kuuinv
triu_ind = np.triu_indices(self.M)
diag_ind = np.diag_indices(self.M)
beta = (self.N - alpha) * 1.0 / self.N
if self.nat_param:
dSu_via_m = np.einsum('da,db->dab', dmu, beta * self.theta_2)
dSu += dSu_via_m
dSuinv = - np.einsum('dab,dbc,dce->dae', self.Suhat, dSu, self.Suhat)
dKuuinv = np.sum(dSuinv, axis=0)
dtheta1 = beta * dSuinv
deta2 = beta * np.einsum('dab,db->da', self.Suhat, dmu)
else:
Suhat = self.Suhat
Suinv = self.Suinv
mu = self.mu
data_f_2 = np.einsum('dab,db->da', Suinv, mu)
dSuhat_via_mhat = np.einsum('da,db->dab', dmu, beta * data_f_2)
dSuhat = dSu + dSuhat_via_mhat
dmuhat = dmu
dSuhatinv = - np.einsum('dab,dbc,dce->dae', Suhat, dSuhat, Suhat)
dSuinv_1 = beta * dSuhatinv
Suhatdmu = np.einsum('dab,db->da', Suhat, dmuhat)
dSuinv = dSuinv_1 + beta * np.einsum('da,db->dab', Suhatdmu, mu)
dtheta1 = - np.einsum('dab,dbc,dce->dae', Suinv, dSuinv, Suinv)
deta2 = beta * np.einsum('dab,db->da', Suinv, Suhatdmu)
dKuuinv = (1 - beta) / beta * np.sum(dSuinv_1, axis=0)
dtheta1T = np.transpose(dtheta1, [0, 2, 1])
dtheta1_R = np.einsum('dab,dbc->dac', self.theta_1_R, dtheta1 + dtheta1T)
deta1_R = np.zeros([self.Dout, self.M * (self.M + 1) / 2])
for d in range(self.Dout):
dtheta1_R_d = dtheta1_R[d, :, :]
theta1_R_d = self.theta_1_R[d, :, :]
dtheta1_R_d[diag_ind] = dtheta1_R_d[diag_ind] * theta1_R_d[diag_ind]
dtheta1_R_d = dtheta1_R_d[triu_ind]
deta1_R[d, :] = dtheta1_R_d.reshape((dtheta1_R_d.shape[0], ))
return deta1_R, deta2, dKuuinv
def test_integration_adaptive_graph_lasso(self, params_in):
'''
Just tests inputs/outputs (not validity of result).
'''
n_features = 20
n_samples = 25
cov, prec, adj = ClusterGraph(
n_blocks=1,
chain_blocks=False,
seed=1,
).create(n_features, 0.8)
prng = np.random.RandomState(2)
X = prng.multivariate_normal(np.zeros(n_features), cov, size=n_samples)
model = AdaptiveGraphLasso(**params_in)
model.fit(X)
assert model.estimator_ is not None
assert model.lam_ is not None
assert np.sum(model.lam_[np.diag_indices(n_features)]) == 0
if params_in['method'] == 'binary':
uvals = set(model.lam_.flat)
assert len(uvals) == 2
assert 0 in uvals
assert 1 in uvals
elif params_in['method'] == 'inverse' or\
params_in['method'] == 'inverse_squared':
uvals = set(model.lam_.flat[model.lam_.flat != 0])
assert len(uvals) > 0
def _nonzero_intersection(m, m_hat):
'''Count the number of nonzeros in and between m and m_hat.
Returns
----------
m_nnz : number of nonzeros in m (w/o diagonal)
m_hat_nnz : number of nonzeros in m_hat (w/o diagonal)
intersection_nnz : number of nonzeros in intersection of m/m_hat
(w/o diagonal)
'''
n_features, _ = m.shape
m_no_diag = m.copy()
m_no_diag[np.diag_indices(n_features)] = 0
m_hat_no_diag = m_hat.copy()
m_hat_no_diag[np.diag_indices(n_features)] = 0
m_hat_nnz = len(np.nonzero(m_hat_no_diag.flat)[0])
m_nnz = len(np.nonzero(m_no_diag.flat)[0])
intersection_nnz = len(np.intersect1d(
np.nonzero(m_no_diag.flat)[0],
np.nonzero(m_hat_no_diag.flat)[0]
))
return m_nnz, m_hat_nnz, intersection_nnz
def _binary_weights(self, estimator):
n_features, _ = estimator.precision_.shape
lam = np.zeros((n_features, n_features))
lam[estimator.precision_ == 0] = 1
lam[np.diag_indices(n_features)] = 0
return lam
def _inverse_squared_weights(self, estimator):
n_features, _ = estimator.precision_.shape
lam = np.zeros((n_features, n_features))
mask = estimator.precision_ != 0
lam[mask] = 1. / (np.abs(estimator.precision_[mask]) ** 2)
mask_0 = estimator.precision_ == 0
lam[mask_0] = np.max(lam[mask].flat) # non-zero in appropriate range
lam[np.diag_indices(n_features)] = 0
return lam
def _inverse_weights(self, estimator):
n_features, _ = estimator.precision_.shape
lam = np.zeros((n_features, n_features))
mask = estimator.precision_ != 0
lam[mask] = 1. / np.abs(estimator.precision_[mask])
mask_0 = estimator.precision_ == 0
lam[mask_0] = np.max(lam[mask].flat) # non-zero in appropriate range
lam[np.diag_indices(n_features)] = 0
return lam
def _count_support_diff(m, m_hat):
n_features, _ = m.shape
m_no_diag = m.copy()
m_no_diag[np.diag_indices(n_features)] = 0
m_hat_no_diag = m_hat.copy()
m_hat_no_diag[np.diag_indices(n_features)] = 0
m_nnz = len(np.nonzero(m_no_diag.flat)[0])
m_hat_nnz = len(np.nonzero(m_hat_no_diag.flat)[0])
nnz_intersect = len(np.intersect1d(np.nonzero(m_no_diag.flat)[0],
np.nonzero(m_hat_no_diag.flat)[0]))
return m_nnz + m_hat_nnz - (2 * nnz_intersect)
def dist_diff(self, geom):
diff = self.coords[:, None, :]-geom.coords[None, :, :]
dist = np.sqrt(np.sum(diff**2, 2))
dist[np.diag_indices(len(self))] = np.inf
return dist, diff
def forward(self, x):
# create diagonal matrices
m = np.zeros((x.size * self.dim)).reshape(-1, self.dim, self.dim)
x = x.reshape(-1, self.dim)
m[(np.s_[:],) + np.diag_indices(x.shape[1])] = x
return m
def __init__(self, n=10, l=None):
self.n = n
self.l = l
self.W = randfilt(self.W if hasattr(self, 'W') else None, (self.l, self.n))
self.C = np.zeros((n, n), dtype='float32')
self.C[np.diag_indices(n)] = 1
def __init__(self, n=10, l=None):
self.n = n
self.l = l
self.W = randfilt(self.W if hasattr(self, 'W') else None, (self.l, self.n))
self.C = np.zeros((n, n), dtype='float32') - 1
self.C[np.diag_indices(n)] = 1
self.SII = None
self.SIO = None
def test_naf_layer_full():
batch_size = 2
for nb_actions in (1, 3):
# Construct single model with NAF as the only layer, hence it is fully deterministic
# since no weights are used, which would be randomly initialized.
L_flat_input = Input(shape=((nb_actions * nb_actions + nb_actions) // 2,))
mu_input = Input(shape=(nb_actions,))
action_input = Input(shape=(nb_actions,))
x = NAFLayer(nb_actions, mode='full')([L_flat_input, mu_input, action_input])
model = Model(inputs=[L_flat_input, mu_input, action_input], outputs=x)
model.compile(loss='mse', optimizer='sgd')
# Create random test data.
L_flat = np.random.random((batch_size, (nb_actions * nb_actions + nb_actions) // 2)).astype('float32')
mu = np.random.random((batch_size, nb_actions)).astype('float32')
action = np.random.random((batch_size, nb_actions)).astype('float32')
# Perform reference computations in numpy since these are much easier to verify.
L = np.zeros((batch_size, nb_actions, nb_actions)).astype('float32')
LT = np.copy(L)
for l, l_T, l_flat in zip(L, LT, L_flat):
l[np.tril_indices(nb_actions)] = l_flat
l[np.diag_indices(nb_actions)] = np.exp(l[np.diag_indices(nb_actions)])
l_T[:, :] = l.T
P = np.array([np.dot(l, l_T) for l, l_T in zip(L, LT)]).astype('float32')
A_ref = np.array([np.dot(np.dot(a - m, p), a - m) for a, m, p in zip(action, mu, P)]).astype('float32')
A_ref *= -.5
# Finally, compute the output of the net, which should be identical to the previously
# computed reference.
A_net = model.predict([L_flat, mu, action]).flatten()
assert_allclose(A_net, A_ref, rtol=1e-5)
def test_naf_layer_diag():
batch_size = 2
for nb_actions in (1, 3):
# Construct single model with NAF as the only layer, hence it is fully deterministic
# since no weights are used, which would be randomly initialized.
L_flat_input = Input(shape=(nb_actions,))
mu_input = Input(shape=(nb_actions,))
action_input = Input(shape=(nb_actions,))
x = NAFLayer(nb_actions, mode='diag')([L_flat_input, mu_input, action_input])
model = Model(inputs=[L_flat_input, mu_input, action_input], outputs=x)
model.compile(loss='mse', optimizer='sgd')
# Create random test data.
L_flat = np.random.random((batch_size, nb_actions)).astype('float32')
mu = np.random.random((batch_size, nb_actions)).astype('float32')
action = np.random.random((batch_size, nb_actions)).astype('float32')
# Perform reference computations in numpy since these are much easier to verify.
P = np.zeros((batch_size, nb_actions, nb_actions)).astype('float32')
for p, l_flat in zip(P, L_flat):
p[np.diag_indices(nb_actions)] = l_flat
print(P, L_flat)
A_ref = np.array([np.dot(np.dot(a - m, p), a - m) for a, m, p in zip(action, mu, P)]).astype('float32')
A_ref *= -.5
# Finally, compute the output of the net, which should be identical to the previously
# computed reference.
A_net = model.predict([L_flat, mu, action]).flatten()
assert_allclose(A_net, A_ref, rtol=1e-5)
def setupMatrices(dim, mult=0.05):
di = np.diag_indices(dim)
bt = np.random.randn(dim,)*mult
Wt = np.zeros((dim,dim))
Wt[di] = np.random.randn(dim,)*mult
We = np.zeros((dim,dim))
We[di] = np.random.randn(dim,)*mult
return Wt,bt,We
def toPartialCorr(Prec):
D=Prec[np.diag_indices(Prec.shape[0])]
P=Prec.copy()
D=np.outer(D,D)
return -P/np.sqrt(D)
