Terminal states D-orthogonality error

le-ander commented 3 years ago

We have briefly discussed this error before and from what I got, this should be about to be fixed anyways. I just wanted to leave this here for tracking and maybe to get a time estimate from you guys when you think this would be fixed? Thanks a lot!

cr.tl.terminal_states(adata, cluster_key='louvain_separate_fine', n_states=None, weight_connectivities=0.2)

[[ 1.00000000e+00 -2.04305480e-03  1.71689486e-13]
 [-2.04305480e-03  1.00000000e+00  1.76311602e-16]
 [ 1.71687968e-13  1.76555548e-16  1.00000000e+00]]
WARNING: Using `3` components would split a block of complex conjugates. Increasing `n_components` to `4`
[[ 1.00000000e+00 -2.04305480e-03  1.01769287e-13  3.91856686e-14]
 [-2.04305480e-03  1.00000000e+00 -3.79125171e-16  3.71041088e-16]
 [ 1.01772757e-13 -3.79064185e-16  1.00000000e+00 -4.63171168e-16]
 [ 3.91856686e-14  3.70986878e-16 -4.61436445e-16  1.00000000e+00]]

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
/opt/python/lib/python3.8/site-packages/cellrank/tl/estimators/_decomposition.py in compute_schur(self, n_components, initial_distribution, method, which, alpha)
    447         try:
--> 448             self._gpcca._do_schur_helper(n_components)
    449         except ValueError:

/opt/python/lib/python3.8/site-packages/cellrank/_vendor/msmtools/analysis/dense/gpcca.py in _do_schur_helper(self, m)
   1080         else:
-> 1081             self.X, self.R, self.eigenvalues = _do_schur(self.P.copy(), self.eta, m, self.z, self.method)
   1082 

/opt/python/lib/python3.8/site-packages/cellrank/_vendor/msmtools/analysis/dense/gpcca.py in _do_schur(P, eta, m, z, method, tol_krylov)
    257         print(X.T.dot(X*eta[:, None]))
--> 258         raise ValueError("Schur vectors appear to not be D-orthogonal.")
    259 

ValueError: Schur vectors appear to not be D-orthogonal.

During handling of the above exception, another exception occurred:

ValueError                                Traceback (most recent call last)
<ipython-input-18-e6c35a0187c4> in <module>
----> 1 cr.tl.terminal_states(adata, cluster_key='louvain_separate_fine', n_states=None, weight_connectivities=0.2)
      2 cr.pl.terminal_states(adata, discrete=True)

/opt/python/lib/python3.8/site-packages/cellrank/tl/_init_term_states.py in terminal_states(adata, estimator, mode, n_states, cluster_key, key, show_plots, copy, return_estimator, fit_kwargs, **kwargs)
    216 ) -> Optional[AnnData]:  # noqa
    217 
--> 218     return _initial_terminal(
    219         adata,
    220         estimator=estimator,

/opt/python/lib/python3.8/site-packages/cellrank/tl/_init_term_states.py in _initial_terminal(adata, estimator, backward, mode, backward_mode, n_states, cluster_key, key, show_plots, copy, return_estimator, fit_kwargs, **kwargs)
    127         )
    128 
--> 129     mc.fit(
    130         n_lineages=n_states,
    131         cluster_key=cluster_key,

/opt/python/lib/python3.8/site-packages/cellrank/tl/estimators/_gppca.py in fit(self, n_lineages, cluster_key, keys, method, compute_absorption_probabilities, **kwargs)
   1081         """
   1082 
-> 1083         super().fit(
   1084             n_lineages=n_lineages,
   1085             cluster_key=cluster_key,

/opt/python/lib/python3.8/site-packages/cellrank/tl/estimators/_base_estimator.py in fit(self, keys, compute_absorption_probabilities, **kwargs)
    983                 - :paramref:`{dp}`
    984         """
--> 985         self._fit_terminal_states(**kwargs)
    986         if compute_absorption_probabilities:
    987             self.compute_absorption_probabilities(keys=keys)

/opt/python/lib/python3.8/site-packages/cellrank/tl/estimators/_gppca.py in _fit_terminal_states(self, n_lineages, cluster_key, method, **kwargs)
    996 
    997         if n_lineages > 1:
--> 998             self.compute_schur(n_lineages + 1, method=method)
    999 
   1000         try:

/opt/python/lib/python3.8/site-packages/cellrank/tl/estimators/_decomposition.py in compute_schur(self, n_components, initial_distribution, method, which, alpha)
    452                 f"Increasing `n_components` to `{n_components + 1}`"
    453             )
--> 454             self._gpcca._do_schur_helper(n_components + 1)
    455 
    456         # make it available for pl

/opt/python/lib/python3.8/site-packages/cellrank/_vendor/msmtools/analysis/dense/gpcca.py in _do_schur_helper(self, m)
   1079                         print('INFO: Using pre-computed schur decomposition')
   1080         else:
-> 1081             self.X, self.R, self.eigenvalues = _do_schur(self.P.copy(), self.eta, m, self.z, self.method)
   1082 
   1083     def minChi(self, m_min, m_max):

/opt/python/lib/python3.8/site-packages/cellrank/_vendor/msmtools/analysis/dense/gpcca.py in _do_schur(P, eta, m, z, method, tol_krylov)
    256     if not np.allclose(X.T.dot(X*eta[:, None]), np.eye(X.shape[1]), atol=1e-6, rtol=1e-5):
    257         print(X.T.dot(X*eta[:, None]))
--> 258         raise ValueError("Schur vectors appear to not be D-orthogonal.")
    259 
    260     # Raise, if X doesn't fullfill the invariant subspace condition!

ValueError: Schur vectors appear to not be D-orthogonal.

Versions:

cellrank==1.1.0 scanpy==1.6.0 anndata==0.7.5 numpy==1.19.4 numba==0.51.2 scipy==1.5.4 pandas==1.1.4 scikit-learn==0.23.2 statsmodels==0.12.1 python-igraph==0.8.3 scvelo==0.2.2 pygam==0.8.0 matplotlib==3.3.3 seaborn==0.11.0

Marius1311 commented 3 years ago

Thaks for posting this @le-ander! We hope to merge to PR into msmtools later this week or early next week.

le-ander commented 3 years ago

That's great news, thanks!

michalk8 commented 3 years ago

If you want to try this sooner, here's a how-to:

git clone https://github.com/theislab/cellrank && cd cellrank
modify the URL in vedorize.toml to git+git://github.com/michalk8/msmtools@c88f4a5c7c464a38cef4264b2fdecafb8f7f6ffd
run python-vendorize
pip install -e.

le-ander commented 3 years ago

Thanks a lot michal!

Marius1311 commented 3 years ago

Sorry for the delay, we're still actively working on this and hope to merge the PR in the next couple of days.

le-ander commented 3 years ago

I actually found my outlier cells and kicked them out. so not an issue for me anymore at the moment :)

Marius1311 commented 3 years ago

fabs, I will leave this open anyways as a reminder until the PR is closed.

linasieverling commented 3 years ago

Hi,

I am getting the same error ("ValueError: Schur vectors appear to not be D-orthogonal."). I tried to modify the URL in vendorize.toml as described above, but also get an error:

$ python-vendorize Collecting git+git://github.com/michalk8/msmtools@c88f4a5c7c464a38cef4264b2fdecafb8f7f6ffd Cloning git://github.com/michalk8/msmtools (to revision c88f4a5c7c464a38cef4264b2fdecafb8f7f6ffd) to /tmp/pip-req-build-iavunliv ERROR: Command errored out with exit status 128: git clone -q git://github.com/michalk8/msmtools /tmp/pip-req-build-iavunliv Check the logs for full command output. Traceback (most recent call last): File "/home/sieverli/.local/bin/python-vendorize", line 20, in main() File "/home/sieverli/.local/bin/python-vendorize", line 10, in main vendorize_requirements(path="vendorize.toml") File "/home/sieverli/.local/lib/python3.6/site-packages/vendorize/init.py", line 20, in vendorize_requirements target_directory=target_directory, File "/home/sieverli/.local/lib/python3.6/site-packages/vendorize/init.py", line 28, in vendorize_requirement cwd=cwd) File "/usr/lib64/python3.6/subprocess.py", line 311, in check_call raise CalledProcessError(retcode, cmd) subprocess.CalledProcessError: Command '['pip', 'install', '--no-dependencies', '--target', 'cellrank/_vendor', 'git+git://github.com/michalk8/msmtools@c88f4a5c7c464a38cef4264b2fdecafb8f7f6ffd']' returned non-zero exit status 1.

Do you have any suggestions how to solve this? Thanks!

michalk8 commented 3 years ago

Hi @linasieverling , tried it just now and didn't experience any issues. I think the possible culprits are:

not writeable /tmp directory - in that case, you can try running TMPDIR=<add here a path to a writeable temp directory> python-vendorize
some issue with the git protocol - you can try the https as pip install git+https://github.com/michalk8/msmtools@c88f4a5c7c464a38cef4264b2fdecafb8f7f6ffd --no-dependencies --target cellrank/_vendor/ --upgrade (--upgrade means that we overwrite stuff in cellrank/_vendor so you don't have to do *)
not writeable cellrank/_vendor (highly unlikely)

* Also for the update to take place, you need to remove the cellrank/_vendor, if present.

linasieverling commented 3 years ago

Thank you, the https solution worked!

Marius1311 commented 3 years ago

amazing, thanks @michalk8 for being fast as always!

Marius1311 commented 3 years ago

Still working on a principled solution to this, we hope to be done by the end of the week. basically, we're re-structuring the underlying GPCCA implementation to make it more easily maintainable.

Marius1311 commented 3 years ago

This will be fixed with #465.

michalk8 commented 3 years ago

closed via #472

Marius1311 commented 3 years ago

@le-ander, this is finally fixed! We completely refactored one of the "engines" of CellRank, the GPCCA method and included it in a new packaged called pyGPCCA. In this context, we also resolved the D-orthogonality bug.

le-ander commented 3 years ago

amazing stuff! I guess there will be a new cellrank release with these things soon? :)

Marius1311 commented 3 years ago

yes! We're about to release.

theislab / cellrank

Terminal states D-orthogonality error #453

Versions: