From 27ca92b4beae92634a288d34a6ea14ab8e5e80d5 Mon Sep 17 00:00:00 2001 From: Xiaojie Qiu Date: Fri, 13 Sep 2019 01:19:43 -0700 Subject: [PATCH] bug fix in analytical transition matrix calculation --- dynamo/tools/cell_velocities.py | 9 ++++++--- dynamo/tools/dimension_reduction.py | 8 ++++---- 2 files changed, 10 insertions(+), 7 deletions(-) diff --git a/dynamo/tools/cell_velocities.py b/dynamo/tools/cell_velocities.py index f4c8d3a97..b00c67297 100755 --- a/dynamo/tools/cell_velocities.py +++ b/dynamo/tools/cell_velocities.py @@ -50,7 +50,7 @@ def cell_velocities(adata, vkey='pca', basis='umap', method='analytical', neg_ce Y = X_pca[indices[i, 1:]] q, u = markov_combination(y, v, Y) Q[i] = q.T - U[i] = u.T + U[i] = (X_embedding[indices[i, 1:]] - X_embedding[i]).T.dot(np.array(q)).T # project in two dimension delta_X[i, :] = X_embedding[i, :] + U[i] @@ -117,13 +117,15 @@ def markov_combination(x, v, X): from cvxopt import matrix, solvers n = X.shape[0] - R = matrix(X - x).T + R = matrix((X - x).astype('double')).T H = R.T * R - f = matrix(v).T * R + f = matrix((v).astype('double')).T * R G = np.vstack((-np.eye(n), np.ones(n))) h = np.zeros(n+1) h[-1] = 1 + solvers.options['show_progress'] = False + p = solvers.qp(H, -f.T, G=matrix(G), h=matrix(h))['x'] u = R * p return p, u @@ -204,3 +206,4 @@ def expected_return_time(M, backward=False): T = 1 / steady_state return T + diff --git a/dynamo/tools/dimension_reduction.py b/dynamo/tools/dimension_reduction.py index 94c13c8ae..01ed8d79f 100755 --- a/dynamo/tools/dimension_reduction.py +++ b/dynamo/tools/dimension_reduction.py @@ -39,14 +39,14 @@ def extract_indices_dist_from_graph(graph, n_neighbors): ind_mat[cur_cell, 0] = cur_cell dist_mat[cur_cell, 0] = 0 - # there could be more than n_neighbors because of an approximate search - if len(cur_neighbors[1]) > n_neighbors - 1: + # there could be more or less than n_neighbors because of an approximate search + if len(cur_neighbors[1]) != n_neighbors - 1: sorted_indices = np.argsort(graph[cur_cell][:, cur_neighbors[1]].A)[0][:(n_neighbors - 1)] ind_mat[cur_cell, 1:] = cur_neighbors[1][sorted_indices] dist_mat[cur_cell, 1:] = graph[cur_cell][0, cur_neighbors[1][sorted_indices]].A else: - ind_mat[cur_cell, 1:] = cur_neighbors[1][1:] # could not broadcast input array from shape (13) into shape (14) - dist_mat[cur_cell, 1:] = graph[cur_cell][:, cur_neighbors[1]].A[1:] + ind_mat[cur_cell, 1:] = cur_neighbors[1] # could not broadcast input array from shape (13) into shape (14) + dist_mat[cur_cell, 1:] = graph[cur_cell][:, cur_neighbors[1]].A return ind_mat, dist_mat