Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fixed randomness in tl.diffmap - compute_eigen (v0 argument) #1858

Merged
merged 7 commits into from Jun 22, 2021
Merged
Show file tree
Hide file tree
Changes from 1 commit
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
11 changes: 10 additions & 1 deletion scanpy/neighbors/__init__.py
Expand Up @@ -937,6 +937,7 @@ def compute_eigen(
n_comps: int = 15,
sym: Optional[bool] = None,
sort: Literal['decrease', 'increase'] = 'decrease',
random_state: AnyRandom = 0,
):
"""\
Compute eigen decomposition of transition matrix.
Expand All @@ -950,6 +951,8 @@ def compute_eigen(
Instead of computing the eigendecomposition of the assymetric
transition matrix, computed the eigendecomposition of the symmetric
Ktilde matrix.
random_state
A numpy random seed

Returns
-------
Expand Down Expand Up @@ -978,8 +981,14 @@ def compute_eigen(
which = 'LM' if sort == 'decrease' else 'SM'
# it pays off to increase the stability with a bit more precision
matrix = matrix.astype(np.float64)

# Setting the random initial vector
random_state = check_random_state(random_state)
np.random.set_state(random_state.get_state())
ivirshup marked this conversation as resolved.
Show resolved Hide resolved
v0 = np.random.randn((matrix.shape[0]))

evals, evecs = scipy.sparse.linalg.eigsh(
matrix, k=n_comps, which=which, ncv=ncv
matrix, k=n_comps, which=which, ncv=ncv, v0=v0
)
evals, evecs = evals.astype(np.float32), evecs.astype(np.float32)
if sort == 'decrease':
Expand Down
8 changes: 7 additions & 1 deletion scanpy/tools/_diffmap.py
Expand Up @@ -2,12 +2,14 @@
from typing import Optional

from ._dpt import _diffmap
from .._utils import AnyRandom


def diffmap(
adata: AnnData,
n_comps: int = 15,
neighbors_key: Optional[str] = None,
random_state: AnyRandom = 0,
copy: bool = False,
):
"""\
Expand Down Expand Up @@ -39,6 +41,8 @@ def diffmap(
.obsp[.uns[neighbors_key]['connectivities_key']],
.obsp[.uns[neighbors_key]['distances_key']] for connectivities and distances
respectively.
random_state
A numpy random seed
Iwo-K marked this conversation as resolved.
Show resolved Hide resolved
copy
Return a copy instead of writing to adata.

Expand Down Expand Up @@ -70,5 +74,7 @@ def diffmap(
if n_comps <= 2:
raise ValueError('Provide any value greater than 2 for `n_comps`. ')
adata = adata.copy() if copy else adata
_diffmap(adata, n_comps=n_comps, neighbors_key=neighbors_key)
_diffmap(
adata, n_comps=n_comps, neighbors_key=neighbors_key, random_state=random_state
)
return adata if copy else None
4 changes: 2 additions & 2 deletions scanpy/tools/_dpt.py
Expand Up @@ -10,11 +10,11 @@
from ..neighbors import Neighbors, OnFlySymMatrix


def _diffmap(adata, n_comps=15, neighbors_key=None):
def _diffmap(adata, n_comps=15, neighbors_key=None, random_state=0):
start = logg.info(f'computing Diffusion Maps using n_comps={n_comps}(=n_dcs)')
dpt = DPT(adata, neighbors_key=neighbors_key)
dpt.compute_transitions()
dpt.compute_eigen(n_comps=n_comps)
dpt.compute_eigen(n_comps=n_comps, random_state=random_state)
adata.obsm['X_diffmap'] = dpt.eigen_basis
adata.uns['diffmap_evals'] = dpt.eigen_values
logg.info(
Expand Down