def run_test_matmul_aa_correlator_kernel(self, ntime, nstand, nchan, misalign=0):
x_shape = (ntime, nchan, nstand*2)
perm = [1,0,2]
x8 = ((np.random.random(size=x_shape+(2,))*2-1)*127).astype(np.int8)
x = x8.astype(np.float32).view(np.complex64).reshape(x_shape)
x = x.transpose(perm)
x = x[..., misalign:]
b_gold = np.matmul(H(x), x)
triu = np.triu_indices(x.shape[-1], 1)
b_gold[..., triu[0], triu[1]] = 0
x = x8.view(bf.DataType.ci8).reshape(x_shape)
x = bf.asarray(x, space='cuda')
x = x.transpose(perm)
x = x[..., misalign:]
b = bf.zeros_like(b_gold, space='cuda')
self.linalg.matmul(1, None, x, 0, b)
b = b.copy('system')
np.testing.assert_allclose(b, b_gold, RTOL*10, ATOL)
python类triu_indices()的实例源码
def get_close_markers(markers,centroids=None, min_distance=20):
if centroids is None:
centroids = [m['centroid']for m in markers]
centroids = np.array(centroids)
ti = np.triu_indices(centroids.shape[0], 1)
def full_idx(i):
#get the pair from condensed matrix index
#defindend inline because ti changes every time
return np.array([ti[0][i], ti[1][i]])
#calculate pairwise distance, return dense distace matrix (upper triangle)
distances = pdist(centroids,'euclidean')
close_pairs = np.where(distances<min_distance)
return full_idx(close_pairs)
def doTotalOrderExperiment(N, binaryWeights = False):
I, J = np.meshgrid(np.arange(N), np.arange(N))
I = I[np.triu_indices(N, 1)]
J = J[np.triu_indices(N, 1)]
#[I, J] = [I[0:N-1], J[0:N-1]]
NEdges = len(I)
R = np.zeros((NEdges, 2))
R[:, 0] = J
R[:, 1] = I
#W = np.random.rand(NEdges)
W = np.ones(NEdges)
#Note: When using binary weights, Y is not necessarily a cocycle
Y = I - J
if binaryWeights:
Y = np.ones(NEdges)
(s, I, H) = doHodge(R, W, Y, verbose = True)
printConsistencyRatios(Y, I, H, W)
#Random flip experiment with linear order
def from_tri_2_sym(tri, dim):
"""convert a upper triangular matrix in 1D format
to 2D symmetric matrix
Parameters
----------
tri: 1D array
Contains elements of upper triangular matrix
dim : int
The dimension of target matrix.
Returns
-------
symm : 2D array
Symmetric matrix in shape=[dim, dim]
"""
symm = np.zeros((dim, dim))
symm[np.triu_indices(dim)] = tri
return symm
def run_test_matmul_aa_ci8_shape(self, shape, transpose=False):
# **TODO: This currently never triggers the transpose path in the backend
shape_complex = shape[:-1] + (shape[-1] * 2,)
# Note: The xGPU-like correlation kernel does not support input values of -128 (only [-127:127])
a8 = ((np.random.random(size=shape_complex) * 2 - 1) * 127).astype(np.int8)
a_gold = a8.astype(np.float32).view(np.complex64)
if transpose:
a_gold = H(a_gold)
# Note: np.matmul seems to be slow and inaccurate when there are batch dims
c_gold = np.matmul(a_gold, H(a_gold))
triu = np.triu_indices(shape[-2] if not transpose else shape[-1], 1)
c_gold[..., triu[0], triu[1]] = 0
a = a8.view(bf.DataType.ci8)
a = bf.asarray(a, space='cuda')
if transpose:
a = H(a)
c = bf.zeros_like(c_gold, space='cuda')
self.linalg.matmul(1, a, None, 0, c)
c = c.copy('system')
np.testing.assert_allclose(c, c_gold, RTOL, ATOL)
def run_test_matmul_aa_dtype_shape(self, shape, dtype, axes=None, conj=False):
a = ((np.random.random(size=shape)) * 127).astype(dtype)
if axes is None:
axes = range(len(shape))
aa = a.transpose(axes)
if conj:
aa = aa.conj()
c_gold = np.matmul(aa, H(aa))
triu = np.triu_indices(shape[axes[-2]], 1)
c_gold[..., triu[0], triu[1]] = 0
a = bf.asarray(a, space='cuda')
aa = a.transpose(axes)
if conj:
aa = aa.conj()
c = bf.zeros_like(c_gold, space='cuda')
self.linalg.matmul(1, aa, None, 0, c)
c = c.copy('system')
np.testing.assert_allclose(c, c_gold, RTOL, ATOL)
def run_benchmark_matmul_aa_correlator_kernel(self, ntime, nstand, nchan):
x_shape = (ntime, nchan, nstand*2)
perm = [1,0,2]
x8 = ((np.random.random(size=x_shape+(2,))*2-1)*127).astype(np.int8)
x = x8.astype(np.float32).view(np.complex64).reshape(x_shape)
x = x.transpose(perm)
b_gold = np.matmul(H(x[:,[0],:]), x[:,[0],:])
triu = np.triu_indices(x_shape[-1], 1)
b_gold[..., triu[0], triu[1]] = 0
x = x8.view(bf.DataType.ci8).reshape(x_shape)
x = bf.asarray(x, space='cuda')
x = x.transpose(perm)
b = bf.zeros_like(b_gold, space='cuda')
bf.device.stream_synchronize();
t0 = time.time()
nrep = 200
for _ in xrange(nrep):
self.linalg.matmul(1, None, x, 0, b)
bf.device.stream_synchronize();
dt = time.time() - t0
nflop = nrep * nchan * ntime * nstand*(nstand+1)/2 * 2*2 * 8
print nstand, '\t', nflop / dt / 1e9, 'GFLOP/s'
print '\t\t', nrep*ntime*nchan / dt / 1e6, 'MHz'
def _compute_dispersion_matrix(X, labels):
n = len(np.unique(labels))
dist = np.zeros((n, n))
ITR = list(itertools.combinations_with_replacement(range(n), 2))
for i, j in tqdm(ITR):
if i == j:
d = pdist(X[labels == i], metric='cosine')
else:
d = cdist(X[labels == i], X[labels == j], metric='cosine')
# Only take upper diagonal (+diagonal elements)
d = d[np.triu_indices(n=d.shape[0], m=d.shape[1], k=0)]
dist[i, j] = dist[j, i] = d.mean()
return dist
def get_gcovmat(h2, rg):
"""
Args: h2: vector with SNP heritabilities
rg: vector with genetic correlations
Returns: numpy trait by trait array with h2 on diagonal and genetic covariance on offdiagnoals
"""
mat = numpy.zeros((len(h2), len(h2)))
mat[numpy.triu_indices(len(h2), 1)] = rg
mat = mat + mat.T
mat = mat * numpy.sqrt(numpy.outer(h2, h2))
numpy.fill_diagonal(mat, h2)
return numpy.array(mat)
# When input files are score files, not beta files, mtot may be unknown.
# Here mtot=1e6 is assumed. The absolute value of the expected variances for each trait is irrelevant for the multi-trait weighting, so it doesn't matter too much what this value is, expecially if M > N.
def applyNNBig(Xin,model,msize=500,start=150):
#Returns an adjacency matrix
n_features=Xin.shape[1]
X=Xin.copy()
#Center
X -= X.mean(axis=0)
std = X.std(axis=0)
std[std == 0] = 1
X /= std
larger=np.zeros((msize,msize))
larger[start:start+n_features,start:start+n_features]=X.T.dot(X)/X.shape[0]
emp_cov_matrix=np.expand_dims(larger,0)
pred=model.predict(np.expand_dims(emp_cov_matrix,0))
pred=pred.reshape(msize,msize)[start:start+n_features,start:start+n_features]
C=np.zeros((X.shape[1],X.shape[1]))
C[np.triu_indices(n_features,k=1)]=pred[np.triu_indices(n_features,k=1)]
C=C+C.T
return C
def linkage(df, n_groups):
# create the distance matrix based on the forbenius norm: |A-B|_F where A is
# a 24 x N matrix with N the number of timeseries inside the dataframe df
# TODO: We can save have time as we only need the upper triangle once as the
# distance matrix is symmetric
if True:
Y = np.empty((n_groups, n_groups,))
Y[:] = np.NAN
for i in range(len(Y)):
for j in range(len(Y[i,:])):
A = df.loc[i+1].values
B = df.loc[j+1].values
#print('Computing distance of:{},{}'.format(i,j))
Y[i,j] = norm(A-B, ord='fro')
# condensed distance matrix as vector for linkage (upper triangle as a vector)
y = Y[np.triu_indices(n_groups, 1)]
# create linkage matrix with wards algorithm an euclidean norm
Z = hac.linkage(y, method='ward', metric='euclidean')
# R = hac.inconsistent(Z, d=10)
return Z
def estimate_ls(self, X):
Ntrain = X.shape[0]
if Ntrain < 10000:
X1 = np.copy(X)
else:
randind = np.random.permutation(Ntrain)
X1 = X[randind[1:10000], :]
d2 = compute_distance_matrix(X1)
D = X1.shape[1]
N = X1.shape[0]
triu_ind = np.triu_indices(N)
ls = np.zeros((D, ))
for i in range(D):
d2i = d2[:, :, i]
d2imed = np.median(d2i[triu_ind])
ls[i] = np.log(d2imed)
return ls
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 estimate_ls(self, X):
Ntrain = X.shape[0]
if Ntrain < 10000:
X1 = np.copy(X)
else:
randind = np.random.permutation(Ntrain)
X1 = X[randind[0:(5*self.no_pseudos[0])], :]
# d2 = compute_distance_matrix(X1)
D = X1.shape[1]
N = X1.shape[0]
triu_ind = np.triu_indices(N)
ls = np.zeros((D, ))
for i in range(D):
X1i = np.reshape(X1[:, i], (N, 1))
d2i = cdist(X1i, X1i, 'euclidean')
# d2i = d2[:, :, i]
d2imed = np.median(d2i[triu_ind])
# d2imed = 0.01
# print d2imed,
ls[i] = np.log(d2imed + 1e-16)
return ls
def array2tree(d, names, outbase="", method="ward"):
"""Return tree representation for array"""
import ete3
Z = fastcluster.linkage(d[np.triu_indices(d.shape[0], 1)], method=method)
tree = sch.to_tree(Z, False)
t = ete3.Tree(getNewick(tree, "", tree.dist, names))
# save tree & newick
if outbase:
pdf, nw = outbase+".nw.pdf", outbase+".nw"
with open(nw, "w") as out:
out.write(t.write())
ts = ete3.TreeStyle()
ts.show_leaf_name = False
ts.layout_fn = mylayout
t.render(pdf, tree_style=ts)
return t
def xyz2U(x, y, z):
"""
compute the potential
"""
dsq = 0.
indices = np.triu_indices(x.size, k=1)
for coord in [x, y, z]:
coord_i = np.tile(coord, (coord.size, 1))
coord_j = coord_i.T
dsq += (coord_i[indices]-coord_j[indices])**2
d = np.sqrt(dsq)
U = np.sum(1./d)
return U
def ang_potential(x0):
"""
If distance is computed along sphere rather than through 3-space.
"""
theta = x0[0:x0.size/2]
phi = np.pi/2-x0[x0.size/2:]
indices = np.triu_indices(theta.size, k=1)
theta_i = np.tile(theta, (theta.size, 1))
theta_j = theta_i.T
phi_i = np.tile(phi, (phi.size, 1))
phi_j = phi_i.T
d = _angularSeparation(theta_i[indices], phi_i[indices], theta_j[indices], phi_j[indices])
U = np.sum(1./d)
return U
def elec_potential_xyz(x0):
x0 = x0.reshape(3, x0.size/3)
x = x0[0, :]
y = x0[1, :]
z = x0[2, :]
dsq = 0.
r = np.sqrt(x**2 + y**2 + z**2)
x = x/r
y = y/r
z = z/r
indices = np.triu_indices(x.size, k=1)
for coord in [x, y, z]:
coord_i = np.tile(coord, (coord.size, 1))
coord_j = coord_i.T
d = (coord_i[indices]-coord_j[indices])**2
dsq += d
U = np.sum(1./np.sqrt(dsq))
return U
def prepare_pairs(X, y, site, indices):
""" Prepare the graph pairs before feeding them to the network """
N, M, F = X.shape
n_pairs = int(len(indices) * (len(indices) - 1) / 2)
triu_pairs = np.triu_indices(len(indices), 1)
X_pairs = np.ones((n_pairs, M, F, 2))
X_pairs[:, :, :, 0] = X[indices][triu_pairs[0]]
X_pairs[:, :, :, 1] = X[indices][triu_pairs[1]]
site_pairs = np.ones(int(n_pairs))
site_pairs[site[indices][triu_pairs[0]] != site[indices][triu_pairs[1]]] = 0
y_pairs = np.ones(int(n_pairs))
y_pairs[y[indices][triu_pairs[0]] != y[indices][triu_pairs[1]]] = 0 # -1
print(n_pairs)
return X_pairs, y_pairs, site_pairs
def plot(self, fig=None, ax=None):
if fig is None:
fig = plt.figure()
if ax is None:
ax = fig.gca()
fmax = self.fs / 2.0
bicoherence = np.copy(self.bicoherence)
n_freq = bicoherence.shape[0]
np.flipud(bicoherence)[np.triu_indices(n_freq, 1)] = 0
bicoherence[np.triu_indices(n_freq, 1)] = 0
bicoherence = bicoherence[:, :n_freq // 2 + 1]
vmin, vmax = compute_vmin_vmax(bicoherence, tick=1e-15, percentile=1)
ax.imshow(bicoherence, cmap=plt.cm.viridis, aspect='auto', vmin=vmin,
vmax=vmax, origin='lower', extent=(0, fmax // 2, 0, fmax),
interpolation='none')
ax.set_title('Bicoherence (%s)' % self.method)
ax.set_xlabel('Frequency (Hz)')
ax.set_ylabel('Frequency (Hz)')
# add_colorbar(fig, cax, vmin, vmax, unit='', ax=ax)
return ax
def cluster_words(words, service_name, size):
stopwords = ["GET", "POST", "total", "http-requests", service_name, "-", "_"]
cleaned_words = []
for word in words:
for stopword in stopwords:
word = word.replace(stopword, "")
cleaned_words.append(word)
def distance(coord):
i, j = coord
return 1 - jaro_distance(cleaned_words[i], cleaned_words[j])
indices = np.triu_indices(len(words), 1)
distances = np.apply_along_axis(distance, 0, indices)
return cluster_of_size(linkage(distances), size)
def doRandomFlipExperiment(N, PercentFlips):
I, J = np.meshgrid(np.arange(N), np.arange(N))
I = I[np.triu_indices(N, 1)]
J = J[np.triu_indices(N, 1)]
NEdges = len(I)
R = np.zeros((NEdges, 2))
R[:, 0] = J
R[:, 1] = I
# toKeep = int(NEdges/200)
# R = R[np.random.permutation(NEdges)[0:toKeep], :]
# NEdges = toKeep
W = np.random.rand(NEdges)
#W = np.ones(NEdges)
Y = np.ones(NEdges)
NFlips = int(PercentFlips*len(Y))
Y[np.random.permutation(NEdges)[0:NFlips]] = -1
(s, I, H) = doHodge(R, W, Y)
[INorm, HNorm] = [getWNorm(I, W), getWNorm(H, W)]
return (INorm, HNorm)
#Do a bunch of random flip experiments, varying the percentage
#of flips, and take the mean consistency norms for each percentage
#over a number of trials
def doBatchExperiments(N, omissions, flips, NTrials = 10):
I, J = np.meshgrid(np.arange(N), np.arange(N))
I = I[np.triu_indices(N, 1)]
J = J[np.triu_indices(N, 1)]
Rankings = np.zeros((len(omissions), len(flips), NTrials, N))
INorms = np.zeros((len(omissions), len(flips), NTrials))
HNorms = np.zeros((len(omissions), len(flips), NTrials))
for t in range(NTrials):
for i in range(0, len(omissions)):
om = omissions[i]
NEdges = int(len(I)*(1-om))
for j in range(len(flips)):
print("Doing trial %i of %i, omission %i of %i, flip %i of %i"%(t, NTrials, i, len(omissions), j, len(flips)))
R = np.zeros((NEdges, 2))
idx = np.random.permutation(len(I))[0:NEdges]
R[:, 0] = I[idx]
R[:, 1] = J[idx]
f = flips[j]
Y = np.ones(NEdges)
W = 0.5 + 0.5*np.random.rand(NEdges)
#W = np.ones(NEdges)
NFlips = int(f*len(Y))
Y[np.random.permutation(NEdges)[0:NFlips]] = -1
(s, IM, H) = doHodge(R, W, Y, verbose = True)
Rankings[i, j, t, :] = s
[INorms[i, j, t], HNorms[i, j, t]] = [getWNorm(IM, W), getWNorm(H, W)]
KTScores = scoreBatchRankings(Rankings)
sio.savemat("BatchResults.mat", {"Rankings":Rankings, "INorms":INorms, "HNorms":HNorms, "omissions":omissions, "flips":flips, "KTScores":KTScores})
def triu_indices(n, k=0):
m = numpy.ones((n, n), int)
a = triu(m, k)
return numpy.where(a != 0)
def __call__(self, row):
'''
Compute partition function stats over each document.
'''
text = row['text']
stat_names = [
'Z_mu', 'Z_std', 'Z_skew', 'Z_kurtosis',
'I_mu', 'I_std', 'I_skew', 'I_kurtosis',
]
stats = {}
for key in stat_names:
stats[key] = 0.0
# Only keep words that are defined in the embedding
valid_tokens = [w for w in text.split() if w in self.Z]
# Take only the unique words in the document
all_tokens = np.array(list(set(valid_tokens)))
if len(all_tokens) > 3:
# Possibly clip the values here as very large Z don't contribute
doc_z = np.array([self.Z[w] for w in all_tokens])
compute_stats(doc_z, stats, "Z")
# Take top x% most descriptive words
z_sort_idx = np.argsort(doc_z)[::-1]
z_cut = max(int(self.intra_document_cutoff * len(doc_z)), 3)
important_index = z_sort_idx[:z_cut]
sub_tokens = all_tokens[important_index]
doc_v = np.array([self.model[w] for w in sub_tokens])
upper_idx = np.triu_indices(doc_v.shape[0], k=1)
dist = np.dot(doc_v, doc_v.T)[upper_idx]
compute_stats(dist, stats, "I")
stats['_ref'] = row['_ref']
return stats
def init_hypers(self, y_train):
"""Summary
Returns:
TYPE: Description
Args:
y_train (TYPE): Description
"""
N = self.N
M = self.M
Din = self.Din
Dout = self.Dout
x_train = self.x_train
if N < 10000:
centroids, label = kmeans2(x_train, M, minit='points')
else:
randind = np.random.permutation(N)
centroids = x_train[randind[0:M], :]
zu = centroids
if N < 10000:
X1 = np.copy(x_train)
else:
randind = np.random.permutation(N)
X1 = X[randind[:5000], :]
x_dist = cdist(X1, X1, 'euclidean')
triu_ind = np.triu_indices(N)
ls = np.zeros((Din, ))
d2imed = np.median(x_dist[triu_ind])
for i in range(Din):
ls[i] = 2 * np.log(d2imed + 1e-16)
sf = np.log(np.array([0.5]))
params = dict()
params['sf'] = sf
params['ls'] = ls
params['zu'] = zu
params['sn'] = np.log(0.01)
return params
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
