bug fix in analytical transition matrix calculation

aristoteleo · Sep 13, 2019 · 27ca92b · 27ca92b
1 parent a1ae881
commit 27ca92b
Show file tree

Hide file tree

Showing 2 changed files with 10 additions and 7 deletions.
diff --git 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
@@ -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