The following are code examples for showing how to use . They are extracted from open source Python projects. You can vote up the examples you like or vote down the exmaples you don’t like. You can also save this page to your account.
Example 1
def get_metrics_per_contact_type(predicted_cmap, target_cmap):
L = predicted_cmap.shape[0]
short_range = np.zeros((L, L), dtype = np.bool)
medium_range = np.zeros((L, L), dtype = np.bool)
long_range = np.zeros((L, L), dtype = np.bool)
more_than_6 = np.tril_indices(L, -6)
more_than_12 = np.tril_indices(L, -12)
more_than_24 = np.tril_indices(L, -24)
short_range[more_than_6] = True
short_range[more_than_12] = False
medium_range[more_than_12] = True
medium_range[more_than_24] = False
long_range[more_than_24] = True
short_range = np.logical_or(short_range, short_range.T)
medium_range = np.logical_or(medium_range, medium_range.T)
long_range = np.logical_or(long_range, long_range.T)
short_range_metrics = get_metrics(predicted_cmap[short_range], target_cmap[short_range])
medium_range_metrics = get_metrics(predicted_cmap[medium_range], target_cmap[medium_range])
long_range_metrics = get_metrics(predicted_cmap[long_range], target_cmap[long_range])
return short_range_metrics, medium_range_metrics, long_range_metrics Example 2
def _chol_idx(self, n_C, rank):
l_idx = np.tril_indices(n_C)
if rank is not None:
# The rank of covariance matrix is specified
idx_rank = np.where(l_idx[1] < rank)
l_idx = (l_idx[0][idx_rank], l_idx[1][idx_rank])
logger.info('Using the rank specified by the user: '
'{}'.format(rank))
else:
rank = n_C
# if not specified, we assume you want to
# estimate a full rank matrix
logger.warning('Please be aware that you did not specify the'
' rank of covariance matrix to estimate.'
'I will assume that the covariance matrix '
'shared among voxels is of full rank.'
'Rank = {}'.format(rank))
logger.warning('Please be aware that estimating a matrix of '
'high rank can be very slow.'
'If you have a good reason to specify a rank '
'lower than the number of experiment conditions,'
' do so.')
return l_idx, rank Example 3
def forward(self, x):
"""
Transforms from the free state to the variable.
Args:
x: Free state vector. Must have length of `self.num_matrices` *
triangular_number.
Returns:
Reconstructed variable.
"""
L = self._validate_vector_length(len(x))
matsize = int((L * 8 + 1) ** 0.5 * 0.5 - 0.5)
xr = np.reshape(x, (self.num_matrices, -1))
var = np.zeros((matsize, matsize, self.num_matrices), settings.float_type)
for i in range(self.num_matrices):
indices = np.tril_indices(matsize, 0)
var[indices + (np.zeros(len(indices[0])).astype(int) + i,)] = xr[i, :]
return var.squeeze() if self.squeeze else var Example 4
def vec_to_tri(vectors, N):
"""
Takes a D x M tensor `vectors' and maps it to a D x matrix_size X matrix_sizetensor
where the where the lower triangle of each matrix_size x matrix_size matrix is
constructed by unpacking each M-vector.
Native TensorFlow version of Custom Op by Mark van der Wilk.
def int_shape(x):
return list(map(int, x.get_shape()))
D, M = int_shape(vectors)
N = int( np.floor( 0.5 * np.sqrt( M * 8. + 1. ) - 0.5 ) )
# Check M is a valid triangle number
assert((matrix * (N + 1)) == (2 * M))
"""
indices = list(zip(*np.tril_indices(N)))
indices = tf.constant([list(i) for i in indices], dtype=tf.int64)
def vec_to_tri_vector(vector):
return tf.scatter_nd(indices=indices, shape=[N, N], updates=vector)
return tf.map_fn(vec_to_tri_vector, vectors) Example 5
def makeImgPatchPrototype(D, compID):
''' Create image patch prototype for specific component
Returns
--------
Xprototype : sqrt(D) x sqrt(D) matrix
'''
# Create a "prototype" image patch of PxP pixels
P = np.sqrt(D)
Xprototype = np.zeros((P, P))
if compID % 4 == 0:
Xprototype[:P / 2] = 1.0
Xprototype = np.rot90(Xprototype, compID / 4)
if compID % 4 == 1:
Xprototype[np.tril_indices(P)] = 1
Xprototype = np.rot90(Xprototype, (compID - 1) / 4)
if compID % 4 == 2:
Xprototype[np.tril_indices(P, 2)] = 1
Xprototype = np.rot90(Xprototype, (compID - 2) / 4)
if compID % 4 == 3:
Xprototype[np.tril_indices(P, -2)] = 1
Xprototype = np.rot90(Xprototype, (compID - 3) / 4)
return Xprototype Example 6
def extract_patches(self, all_seqdata, with_targets=True, inplace=False):
print("\tConverting dataset to Numpy array...")
L0_X, L0_Y = list(), list()
for seqdata in all_seqdata:
if not inplace:
seqdata = copy.deepcopy(seqdata)
L = len(seqdata)
tril_indices = np.tril_indices(L, -1) # Lower triangle without diagonal
cpreds = seqdata.cpreds
for cpred in cpreds:
assert(len(cpred.shape) == 2)
if cpred is not None:
cpred_shape = cpred.shape
# Create contact maps from file data
print(len(cpreds))
for i in range(len(cpreds)):
if cpreds[i] is None:
cpreds[i] = np.full(cpred_shape, Params.DISTANCE_NAN, dtype=Params.FLOAT_DTYPE)
else:
cpreds[i][np.isnan(cpreds[i])] = Params.DISTANCE_NAN
cpreds[i] = np.asarray(cpreds[i][tril_indices], dtype=Params.FLOAT_DTYPE)
cpreds = np.asarray(cpreds, dtype=Params.FLOAT_DTYPE).T
L0_X.append(cpreds)
if with_targets:
contacts = distances_to_contacts(seqdata.distances)
contacts = contacts[tril_indices]
L0_Y.append(contacts)
L0_X = np.concatenate(L0_X, axis=0)
if with_targets:
L0_Y = np.concatenate(L0_Y, axis=0)
return (L0_X, L0_Y) if with_targets else L0_X Example 7
def separation_under_k_aa(cmap):
L = cmap.shape[0]
mask = np.ones((L, L), dtype=np.bool)
more_than_k = np.tril_indices(L, -Params.MIN_AA_SEPARATION)
mask[more_than_k] = False
mask = np.logical_and(mask, mask.T)
return mask Example 8
def _non_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.tril_indices(n, -1), (n, n)).reshape((1, -1)).astype(np.int32)
return cuda.to_gpu(idx + idx_offsets) Example 9
def _get_batch_non_diagonal_cpu(array):
batch_size, m, n = array.shape
assert m == n
rows, cols = np.tril_indices(n, -1)
return array[:, rows, cols] Example 10
def _set_batch_non_diagonal_cpu(array, non_diag_val):
batch_size, m, n = array.shape
assert m == n
rows, cols = np.tril_indices(n, -1)
array[:, rows, cols] = non_diag_val Example 11
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)) Example 12
def getImageDescriptors_HOG_cdist(self, all_emb, ref_emb, ref_mask):
# unnormalized cosine distance for HOG
dist = numpy.dot(all_emb, ref_emb.T)
# normalize by length of query descriptor projected on reference
norm = numpy.sqrt(numpy.dot(numpy.square(all_emb), ref_mask.T))
dist /= norm
dist[numpy.isinf(dist)] = 0.
dist[numpy.isnan(dist)] = 0.
# dist[numpy.triu_indices(dist.shape[0], 1)] = numpy.maximum(dist[numpy.triu_indices(dist.shape[0], 1)],
# dist.T[numpy.triu_indices(dist.shape[0], 1)])
# dist[numpy.tril_indices(dist.shape[0], -1)] = 0.
# dist += dist.T
return dist Example 13
def mk_test(x, alpha = 0.05):
n = len(x)
# calculate S
listMa = np.matrix(x) # convert input List to 1D matrix
subMa = np.sign(listMa.T - listMa) # calculate all possible differences in matrix
# with itself and save only sign of difference (-1,0,1)
s = np.sum( subMa[np.tril_indices(n,-1)] ) # sum lower left triangle of matrix
# calculate the unique data
# return_counts=True returns a second array that is equivalent to tp in old version
unique_x = np.unique(x, return_counts=True)
g = len(unique_x[0])
# calculate the var(s)
if n == g: # there is no tie
var_s = (n*(n-1)*(2*n+5))/18
else: # there are some ties in data
tp = unique_x[1]
var_s = (n*(n-1)*(2*n+5) + np.sum(tp*(tp-1)*(2*tp+5)))/18
if s>0:
z = (s - 1)/np.sqrt(var_s)
elif s == 0:
z = 0
elif s<0:
z = (s + 1)/np.sqrt(var_s)
# calculate the p_value
p = 2*(1-scipy.stats.norm.cdf(abs(z))) # two tail test
h = abs(z) > scipy.stats.norm.ppf(1-alpha/2)
return h, p
#########################################################################################
#hdf to tiff
######################################################################################### Example 14
def reference_inverse(self, matrices):
# This is the inverse operation of the vec_to_tri op being tested.
D, N, _ = matrices.shape
M = (N * (N + 1)) // 2
tril_indices = np.tril_indices(N)
output = np.zeros((D, M))
for vector_index in range(D):
matrix = matrices[vector_index, :]
output[vector_index, :] = matrix[tril_indices]
return output Example 15
def backward(self, y):
"""
Transforms from the variable to the free state.
Args:
y: Variable representation.
Returns:
Free state.
"""
N = int(np.sqrt(y.size / self.num_matrices))
reshaped = np.reshape(y, (N, N, self.num_matrices))
size = len(reshaped)
triangular = reshaped[np.tril_indices(size, 0)].T
return triangular Example 16
def unvech_kh(v, cols_stacked=True, norm_conserved=False):
# Restore matrix dimension and add triangular
v = v.flatten()
d = int(0.5 * (np.sqrt(8 * len(v) + 1) - 1))
X = np.empty((d, d))
triu_idx_r = []
triu_idx_c = []
# Find appropriate indexes
if cols_stacked:
for c in range(0, d):
for r in range(0, c+1):
triu_idx_r.append(r)
triu_idx_c.append(c)
else:
for r in range(0, d):
for c in range(r, d):
triu_idx_r.append(r)
triu_idx_c.append(c)
# Restore original matrix
triu_idx = (triu_idx_r, triu_idx_c)
X[triu_idx] = v
X[np.tril_indices(d)] = X.T[np.tril_indices(d)]
# Undo rescaling on off diagonal elements
if norm_conserved:
X[np.triu_indices(d, 1)] /= np.sqrt(2)
X[np.tril_indices(d, -1)] /= np.sqrt(2)
return X Example 17
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) Example 18
def __init__(self, tv_dim, ndim, nmix):
self.tv_dim = tv_dim
self.ndim = ndim
self.nmix = nmix
self.itril = np.tril_indices(tv_dim)
self.Sigma = np.empty((self.ndim * self.nmix, 1))
self.T_iS = None # np.empty((self.tv_dim, self.ndim * self.nmix))
self.T_iS_Tt = None # np.empty((self.nmix, self.tv_dim * (self.tv_dim+1)/2))
self.Tm = np.empty((self.tv_dim, self.ndim * self.nmix))
self.Im = np.eye(self.tv_dim) Example 19
def corrcoef_matrix(matrix):
# Code originating from http://stackoverflow.com/a/24547964 by http://stackoverflow.com/users/2455058/jingchao
r = np.corrcoef(matrix)
rf = r[np.triu_indices(r.shape[0], 1)]
df = matrix.shape[1] - 2
ts = rf * rf * (df / (1 - rf * rf))
pf = betai(0.5 * df, 0.5, df / (df + ts))
p = np.zeros(shape=r.shape)
p[np.triu_indices(p.shape[0], 1)] = pf
p[np.tril_indices(p.shape[0], -1)] = pf
p[np.diag_indices(p.shape[0])] = np.ones(p.shape[0])
return r, p Example 20
def vector_from_array(array):
"""Get triangle of output in vector from a correlation-type array
Parameters
----------
array : np.array
Returns
-------
vector (np.array)
Notes
-----
Old Matlab code indexes by column (aka.Fortran-style), so to get the indices
of the top triangle, we have to do some reshaping.
Otherwise, if the vector made up by rows is OK, then simply :
triangle = np.triu_indices(array.size, k=1), out = array[triangle]
"""
triangle_lower = np.tril_indices(array.shape[0], k=-1)
flatten_idx = np.arange(array.size).reshape(array.shape)[triangle_lower]
triangle = np.unravel_index(flatten_idx, array.shape, order='F')
# triangle = np.triu_indices(array.size, k=1)
# out = array[triangle]
return array[triangle] Example 21
def feature_extractor_2(state, action, normalize=True):
# features are
# - bias
# - the elements of the state vector
# - all second order terms (mixed and pure)
# - sines of every first and second order term
# - cosines of every first and second order term
# all of that times two (i.e. for each of the actions)
s_prime = np.r_[1, state] # dim(state) + 1
# dim(state) (dim(state) + 1) / 2
quadratic = np.outer(s_prime, s_prime)[np.tril_indices(s_prime.shape[0])].reshape(-1)
# (dim(state) (dim(state) + 1) / 2) - 1
sines = np.sin(quadratic[1:])
cosines = np.cos(quadratic[1:])
# dim(state)(dim(state) + 1) - 2 + dim(state)(dim(state) + 1) / 2
state_feats = np.r_[quadratic, sines, cosines]
# dim(state_feats) * 2
features = np.outer(state_feats, np.array([0, 1]) == action).T.reshape(-1)
if normalize:
# normalize everything but the bias.
norm = features / (np.linalg.norm(features))
norm[0] = int(action == 0) * 1.0
norm[state_feats.shape[0]] = int(action == 1) * 1.0
return norm
else:
return features Example 22
def _get_lowTriIDs(K):
if K in lowTriIDsDict:
return lowTriIDsDict[K]
else:
ltIDs = np.tril_indices(K)
lowTriIDsDict[K] = ltIDs
return ltIDs Example 23
def _get_lowTriIDs(K):
if K in lowTriIDsDict:
return lowTriIDsDict[K]
else:
ltIDs = np.tril_indices(K)
lowTriIDsDict[K] = ltIDs
return ltIDs Example 24
def _scorematrix2rankedlist_greedy(M, nPairs, doKeepZeros=False):
''' Return the nPairs highest-ranked pairs in score matrix M
Args
-------
M : score matrix, K x K
should have only entries kA,kB where kA <= kB
Returns
--------
aList : list of integer ids for rows of M
bList : list of integer ids for cols of M
Example
---------
_scorematrix2rankedlist( [0 2 3], [0 0 1], [0 0 0], 3)
>> [ (0,2), (0,1), (1,2)]
'''
M = M.copy()
M[np.tril_indices(M.shape[0])] = - np.inf
Mflat = M.flatten()
sortIDs = np.argsort(-1 * Mflat)
# Remove any entries that are -Inf
sortIDs = sortIDs[Mflat[sortIDs] != -np.inf]
if not doKeepZeros:
# Remove any entries that are zero
sortIDs = sortIDs[Mflat[sortIDs] != 0]
bestrs, bestcs = np.unravel_index(sortIDs, M.shape)
return bestrs[:nPairs].tolist(), bestcs[:nPairs].tolist() Example 25
def _tril_indices(n, k=0):
"""Replacement for tril_indices that is provided for numpy >= 1.4"""
mask = np.greater_equal(np.subtract.outer(np.arange(n), np.arange(n)), -k)
indices = np.where(mask)
return indices Example 26
def _get_strength_matrix(self):
assert self.pmi.shape == self.ts_correlation.shape
assert is_square(self.pmi)
a = np.multiply(self.pmi, self.ts_correlation)
lower_indexes = np.tril_indices(self.pmi.shape[0])
a[lower_indexes] = np.nan
return a Example 27
def tril_elements(M):
'''
Somewhat like matlab's "diag" function, but for lower triangular matrices
tril_elements(randn(D*(D+1)//2))
'''
if len(M.shape)==2:
# M is a matrix
if not M.shape[0]==M.shape[1]:
raise ValueError("Extracting upper triangle elements supported only on square arrays")
# Extract upper trianglular elements
i = np.tril_indices(M.shape[0])
return M[i]
if len(M.shape)==1:
# M is a vector
# N(N+1)/2 = K
# N(N+1) = 2K
# NN+N = 2K
# NN+N-2K=0
# A x^2 + Bx + C
# -1 +- sqrt(1-4*1*(-2K))
# -----------------------
# 2
#
# (sqrt(1+8*K)-1)/2
K = M.shape[0]
N = (np.sqrt(1+8*K)-1)/2
if N!=round(N):
raise ValueError('Cannot pack %d elements into a square triangular matrix'%K)
N = int(N)
result = np.zeros((N,N))
result[np.tril_indices(N)] = M
return result
raise ValueError("Must be 2D matrix or 1D vector") Example 28
def givens_rotation(A):
"""Perform QR decomposition of matrix A using Givens rotation."""
(num_rows, num_cols) = np.shape(A)
# Initialize orthogonal matrix Q and upper triangular matrix R.
Q = np.identity(num_rows)
R = np.copy(A)
# Iterate over lower triangular matrix.
(rows, cols) = np.tril_indices(num_rows, -1, num_cols)
for (row, col) in zip(rows, cols):
# Compute Givens rotation matrix and
# zero-out lower triangular matrix entries.
if R[row, col] != 0:
(c, s) = _givens_rotation_matrix_entries(R[col, col], R[row, col])
G = np.identity(num_rows)
G[[col, row], [col, row]] = c
G[row, col] = s
G[col, row] = -s
R = np.dot(G, R)
Q = np.dot(Q, G.T)
return (Q, R) Example 29
def _fill_lower_triangular(self, x):
"""Numpy implementation of `fill_lower_triangular`."""
x = np.asarray(x)
d = x.shape[-1]
# d = n(n+1)/2 implies n is:
n = int(0.5 * (math.sqrt(1. + 8. * d) - 1.))
ids = np.tril_indices(n)
y = np.zeros(list(x.shape[:-1]) + [n, n], dtype=x.dtype)
y[..., ids[0], ids[1]] = x
return y Example 30
def _binary_sim(matrix):
"""Compute a jaccard similarity matrix.
Vecorization based on: https://stackoverflow.com/a/40579567
Parameters
----------
matrix : np.array
Returns
-------
np.array
matrix of shape (n_rows, n_rows) giving the similarity
between rows of the input matrix.
"""
# Generate the indices of the lower triangle of our result matrix.
# The diagonal is offset by -1 because the identity in a similarity
# matrix is always 1.
r, c = np.tril_indices(matrix.shape[0], -1)
# Particularly large groups can blow out memory usage. Chunk the calculation
# into steps that require no more than ~100MB of memory at a time.
# We have 4 2d intermediate arrays in memory at a given time, plus the
# input and output:
# p1 = max_rows * matrix.shape[1]
# p2 = max_rows * matrix.shape[1]
# intersection = max_rows * matrix.shape[1] * 4
# union = max_rows * matrix.shape[1] * 8
# This adds up to:
# memory usage = max_rows * matrix.shape[1] * 14
mem_limit = 100 * pow(2, 20)
max_rows = mem_limit / (14 * matrix.shape[1])
out = np.eye(matrix.shape[0])
for c_batch, r_batch in _batch(c, r, max_rows):
# Use those indices to build two matrices that contains all
# the rows we need to do a similarity comparison on
p1 = matrix[c_batch]
p2 = matrix[r_batch]
# Run the main jaccard calculation
intersection = np.logical_and(p1, p2).sum(1)
union = np.logical_or(p1, p2).sum(1).astype(np.float64)
# Build the result matrix with 1's on the diagonal
# Insert the result of our similarity calculation at their original indices
out[c_batch, r_batch] = intersection / union
# Above only populated half of the matrix, the mirrored diagonal contains
# only zeros. Fix that up by transposing. Adding the transposed matrix double
# counts the diagonal so subtract that back out. We could skip this step and
# leave half the matrix empty, but it takes a fraction of a ms to be correct
# even on mid-sized inputs (~200 queries).
return out + out.T - np.diag(np.diag(out)) Example 31
def gaus_posterior(dim, std0):
"""Initialise a posterior Gaussian distribution with a diagonal covariance.
Even though this is initialised with a diagonal covariance, a full
covariance will be learned, using a lower triangular Cholesky
parameterisation.
Parameters
----------
dim : tuple or list
the dimension of this distribution.
std0 : float
the initial (unoptimized) diagonal standard deviation of this
distribution.
Returns
-------
Q : tf.contrib.distributions.MultivariateNormalTriL
the initialised posterior Gaussian object.
Note
----
This will make tf.Variables on the randomly initialised mean and covariance
of the posterior. The initialisation of the mean is from a Normal with zero
mean, and ``std0`` standard deviation, and the initialisation of the (lower
triangular of the) covariance is from a gamma distribution with an alpha of
``std0`` and a beta of 1.
"""
o, i = dim
# Optimize only values in lower triangular
u, v = np.tril_indices(i)
indices = (u * i + v)[:, np.newaxis]
l0 = np.tile(np.eye(i), [o, 1, 1])[:, u, v].T
l0 = l0 * tf.random_gamma(alpha=std0, shape=l0.shape, seed=next(seedgen))
lflat = tf.Variable(l0, name="W_cov_q")
Lt = tf.transpose(tf.scatter_nd(indices, lflat, shape=(i * i, o)))
L = tf.reshape(Lt, (o, i, i))
mu_0 = tf.random_normal((o, i), stddev=std0, seed=next(seedgen))
mu = tf.Variable(mu_0, name="W_mu_q")
Q = MultivariateNormalTriL(mu, L)
return Q
#
# KL divergence calculations
# Example 32
def impute_missing_bins(hic_matrix, regions=None, per_chromosome=True, stat=np.ma.mean):
"""
Impute missing contacts in a Hi-C matrix.
For inter-chromosomal data uses the mean of all inter-chromosomal contacts,
for intra-chromosomal data uses the mean of intra-chromosomal counts at the corresponding diagonal.
:param hic_matrix: A square numpy array
:param regions: A list of :class:`~GenomicRegion`s - if omitted, will create a dummy list
:param per_chromosome: Do imputation on a per-chromosome basis (recommended)
:param stat: The aggregation statistic to be used for imputation, defaults to the mean.
"""
if regions is None:
for i in range(hic_matrix.shape[0]):
regions.append(GenomicRegion(chromosome='', start=i, end=i))
chr_bins = dict()
for i, region in enumerate(regions):
if region.chromosome not in chr_bins:
chr_bins[region.chromosome] = [i, i]
else:
chr_bins[region.chromosome][1] = i
n = len(regions)
if not hasattr(hic_matrix, "mask"):
hic_matrix = masked_matrix(hic_matrix)
imputed = hic_matrix.copy()
if per_chromosome:
for c_start, c_end in chr_bins.itervalues():
# Correcting intrachromoc_startmal contacts by mean contact count at each diagonal
for i in range(c_end - c_start):
ind = kth_diag_indices(c_end - c_start, -i)
diag = imputed[c_start:c_end, c_start:c_end][ind]
diag[diag.mask] = stat(diag)
imputed[c_start:c_end, c_start:c_end][ind] = diag
# Correcting interchromoc_startmal contacts by mean of all contact counts between
# each set of chromoc_startmes
for other_start, other_end in chr_bins.itervalues():
# Only correct upper triangle
if other_start <= c_start:
continue
inter = imputed[c_start:c_end, other_start:other_end]
inter[inter.mask] = stat(inter)
imputed[c_start:c_end, other_start:other_end] = inter
else:
for i in range(n):
diag = imputed[kth_diag_indices(n, -i)]
diag[diag.mask] = stat(diag)
imputed[kth_diag_indices(n, -i)] = diag
# Copying upper triangle to lower triangle
imputed[np.tril_indices(n)] = imputed.T[np.tril_indices(n)]
return imputed Example 33
def _plot(self, region=None, cax=None):
if region is None:
raise ValueError("Cannot plot triangle plot for whole genome.")
hm, sr = sub_matrix_regions(self.hic_matrix, self.regions, region)
hm[np.tril_indices(hm.shape[0])] = np.nan
# Remove part of matrix further away than max_dist
if self.max_dist is not None:
for i in range(hm.shape[0]):
i_region = sr[i]
for j in range(hm.shape[1]):
j_region = sr[j]
if j_region.start-i_region.end > self.max_dist:
hm[i, j] = np.nan
hm_masked = np.ma.MaskedArray(hm, mask=np.isnan(hm))
# prepare an array of the corner coordinates of the Hic-matrix
# Distances have to be scaled by sqrt(2), because the diagonals of the bins
# are sqrt(2)*len(bin_size)
sqrt2 = math.sqrt(2)
bin_coords = np.r_[[(x.start - 1) for x in sr], sr[-1].end]/sqrt2
X, Y = np.meshgrid(bin_coords, bin_coords)
# rotatate coordinate matrix 45 degrees
sin45 = math.sin(math.radians(45))
X_, Y_ = X*sin45 + Y*sin45, X*sin45 - Y*sin45
# shift x coords to correct start coordinate and center the diagonal directly on the
# x-axis
X_ -= X_[1, 0] - (sr[0].start - 1)
Y_ -= .5*np.min(Y_) + .5*np.max(Y_)
# create plot
self.hicmesh = self.ax.pcolormesh(X_, Y_, hm_masked, cmap=self.colormap, norm=self.norm)
# set limits and aspect ratio
#self.ax.set_aspect(aspect="equal")
ylim_max = 0.5*(region.end-region.start)
if self.max_dist is not None and self.max_dist/2 < ylim_max:
ylim_max = self.max_dist/2
self.ax.set_ylim(0, ylim_max)
# remove y ticks
self.ax.set_yticks([])
# hide background patch
self.ax.patch.set_visible(False)
if self.show_colorbar:
self.add_colorbar(cax) Example 34
def mixtureMV(self, mu, sig, ps):
"""LearningRules.mixtureMV
Sample from a multivariate gaussian mixture with \mu = means, \Sigma = cov matrix
"""
# print "%s.mixtureMV: Implement me" % (self.__class__.__name__)
# print "mixtureMV", "mu", mu.shape, "sig", sig.shape, "ps", ps.shape
mixcomps = self.mixcomps
d = self.dim
assert mixcomps == ps.shape[0]
assert not np.isnan(mu).any()
assert not np.isnan(sig).any()
assert not np.isnan(ps).any()
# check flat inputs
if mu.shape[1] == 1:
mu = mu.reshape((self.mixcomps, self.dim))
# mu_ = mu.reshape((mixcomps, d))
if sig.shape[1] == 1:
# S_diag = sig[:(mixcomps * self.dim)]
# S_triu = sig[(mixcomps * self.dim):]
S_diag = sig[:(mixcomps * d)].reshape((mixcomps, d))
S_triu = sig[(mixcomps * d):].reshape((mixcomps, -1))
S = []
for i in range(mixcomps):
S_ = np.diag(S_diag[i])
S_[np.triu_indices(d, 1)] = S_triu[i]
S_[np.tril_indices(d, -1)] = S_triu[i] # check if thats correct
S.append(S_)
sig = np.array(S)
# print "mixtureMV", "mu", mu.shape, "sig", sig.shape, "ps", ps.shape
# d = mu.shape[0]
compidx = self.mixture_sample_prior(ps)
mu_ = mu[compidx]
S_ = sig[compidx]
# print "compidx", compidx, "mu_", mu_, "S_", S_
y = np.random.multivariate_normal(mu_, S_)
return y
# mixture, mdn_loss, and softmax are taken from karpathy's MixtureDensityNets.py Example 35
def test_mixtureMV(args):
"""test_mixtureMV
Test mixtureMV() by sampling from dummy statistics mu, S, p
"""
mixcomps = 3
d = 2
mu = np.random.uniform(-1, 1, size = (mixcomps, d))
S = np.array([np.eye(d) * 0.01 for c in range(mixcomps)])
for s in S:
s += 0.05 * np.random.uniform(-1, 1, size = s.shape)
s[np.tril_indices(d, -1)] = s[np.triu_indices(d, 1)]
pi = np.ones((mixcomps, )) * 1.0/mixcomps
lr = LearningRules(ndim_out = d, dim = d)
lr.learnFORCEmdn_setup(mixcomps = mixcomps)
fig = plt.figure()
gs = gridspec.GridSpec(2, mixcomps)
for g in gs:
fig.add_subplot(g)
ax = fig.axes[0]
ax.plot(mu[:,0], mu[:,1], "ko", alpha = 0.5)
# ax.set_xlim((-1.1, 1.1))
# ax.set_ylim((-1.1, 1.1))
ax = fig.axes[0] # (2 * mixcomps) + c]
for i in range(1000):
y = lr.mixtureMV(mu, S, pi)
print "y", y
ax.plot([y[0]], [y[1]], "ro", alpha = 0.2, markersize = 3.0)
ax.set_aspect(1)
for c in range(mixcomps):
ax = fig.axes[mixcomps + c]
ax.imshow(S[c], cmap = plt.get_cmap('PiYG'))
fig.show()
plt.show()
print "mu", mu.shape, "S", S.shape Example 36
def updateScoreMat_wholeELBO(ScoreMat, curModel, SS, doAllPairs=0):
''' Calculate upper-tri matrix of exact ELBO gap for each candidate pair
Returns
---------
Mraw : 2D array, size K x K. Uppert tri entries carry content.
Mraw[j,k] gives the scalar ELBO gap for the potential merge of j,k
'''
K = SS.K
if doAllPairs:
AGap = curModel.allocModel.calcHardMergeGap_AllPairs(SS)
OGap = curModel.obsModel.calcHardMergeGap_AllPairs(SS)
ScoreMat = AGap + OGap
ScoreMat[np.tril_indices(SS.K)] = -np.inf
for k, uID in enumerate(SS.uIDs):
CountTracker[uID] = SS.getCountVec()[k]
nUpdated = SS.K * (SS.K - 1) / 2
else:
ScoreMat[np.tril_indices(SS.K)] = -np.inf
# Rescore only specific pairs that are positive
redoMask = ScoreMat > -1 * ELBO_GAP_ACCEPT_TOL
for k, uID in enumerate(SS.uIDs):
if CountTracker[uID] == 0:
# Always precompute for brand-new comps
redoMask[k, :] = 1
redoMask[:, k] = 1
else:
absDiff = np.abs(SS.getCountVec()[k] - CountTracker[uID])
percDiff = absDiff / (CountTracker[uID] + 1e-10)
if percDiff > 0.25:
redoMask[k, :] = 1
redoMask[:, k] = 1
CountTracker[uID] = SS.getCountVec()[k]
redoMask[np.tril_indices(SS.K)] = 0
aList, bList = np.unravel_index(np.flatnonzero(redoMask), (SS.K, SS.K))
if len(aList) > 0:
mPairIDs = zip(aList, bList)
AGap = curModel.allocModel.calcHardMergeGap_SpecificPairs(
SS, mPairIDs)
OGap = curModel.obsModel.calcHardMergeGap_SpecificPairs(
SS, mPairIDs)
ScoreMat[aList, bList] = AGap + OGap
nUpdated = len(aList)
MergeLogger.log('MERGE ScoreMat Updates: %d entries.' % (nUpdated),
level='debug')
return ScoreMat Example 37
def posthoc_tukey_hsd(x, g, alpha = 0.05):
'''Pairwise comparisons with TukeyHSD confidence intervals. This is a
convenience function to make statsmodels pairwise_tukeyhsd() method more
applicable for further use.
Parameters
----------
x : array_like or pandas Series object, 1d
An array, any object exposing the array interface, containing
the response variable. NaN values will cause an error. Please
handle manually.
g : array_like or pandas Series object, 1d
An array, any object exposing the array interface, containing
the groups' names. Can be string or integers.
alpha : float, optional
Significance level for the test. Default is 0.05.
Returns
-------
Numpy ndarray where 0 is False (not significant), 1 is True (significant),
and -1 is for diagonal elements.
Examples
--------
>>> x = [[1,2,3,4,5], [35,31,75,40,21], [10,6,9,6,1]]
>>> g = [['a'] * 5, ['b'] * 5, ['c'] * 5]
>>> sp.posthoc_tukey_hsd(np.concatenate(x), np.concatenate(g))
array([[-1, 1, 0],
[ 1, -1, 1],
[ 0, 1, -1]])
'''
result = pairwise_tukeyhsd(x, g, alpha=0.05)
groups = np.array(result.groupsunique, dtype=np.str)
groups_len = len(groups)
vs = np.arange(groups_len, dtype=np.int)[:,None].T.repeat(groups_len, axis=0)
for a in result.summary()[1:]:
a0 = str(a[0])
a1 = str(a[1])
a0i = np.where(groups == a0)[0][0]
a1i = np.where(groups == a1)[0][0]
vs[a0i, a1i] = 1 if str(a[5]) == 'True' else 0
vs = np.triu(vs)
np.fill_diagonal(vs, -1)
tri_lower = np.tril_indices(vs.shape[0], -1)
vs[tri_lower] = vs.T[tri_lower]
return vs 