markovmodel / PyEMMA

đźš‚ Python API for Emma's Markov Model Algorithms đźš‚
http://pyemma.org
GNU Lesser General Public License v3.0
311 stars 119 forks source link

[msm.its] fails if sparse estimation fails #403

Closed gph82 closed 9 years ago

gph82 commented 9 years ago

I have a count-matrix for which sparse-MLE-estimation fails but dense-MLE-estimation succeeds: eml2

Apart from looking into why one converges and the other does not, I am interested in the behavior of [msm.its]. If parse the dtrajs and lagtime behind said count-matrix to [msm/its], it fails (see error below). Why? Because [analysis/dense/decomposition] is trying to operate on a matrix full of zeroes (which is what the sparse estimator returns after a non-converged run).

How can one handle best this (and similar) exceptions within [msm/its]?

pyemma.msm.its(dtrajs, lags=[10])
---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
<ipython-input-137-23ddf9cb551b> in <module>()
----> 1 pyemma.msm.its(dtrajs, lags=[10])

/home/mi/gph82/programs/PyEmma/pyemma/msm/api.pyc in its(dtrajs, lags, nits, reversible, connected)
    104 
    105     """
--> 106     itsobj = ImpliedTimescales(dtrajs, lags=lags, nits=nits, reversible=reversible, connected=connected)
    107     return itsobj
    108 

/home/mi/gph82/programs/PyEmma/pyemma/msm/ui/timescales.pyc in __init__(self, dtrajs, lags, nits, connected, reversible, failfast)
    107 
    108         # estimate
--> 109         self._estimate()
    110 
    111     def _generate_lags(self, maxlag, multiplier):

/home/mi/gph82/programs/PyEmma/pyemma/msm/ui/timescales.pyc in _estimate(self)
    155             C = cmatrix(self._dtrajs, tau)
    156             # estimate timescales
--> 157             ts = self._estimate_ts_tau(C, tau)
    158             if (ts is None):
    159                 maxnlags = i

/home/mi/gph82/programs/PyEmma/pyemma/msm/ui/timescales.pyc in _estimate_ts_tau(self, C, tau)
    136             T = T.toarray()
    137             # timescales
--> 138             ts = timescales(T, tau, k=min(self._nits, len(T)) + 1, reversible=self._reversible)[1:]
    139             return ts
    140         else:

/home/mi/gph82/programs/PyEmma/pyemma/msm/analysis/api.pyc in timescales(T, tau, k, ncv, reversible, mu)
    489         return sparse.decomposition.timescales(T, tau=tau, k=k, ncv=ncv)
    490     elif _isdense(T):
--> 491         return dense.decomposition.timescales(T, tau=tau, k=k, reversible=reversible, mu=mu)
    492     else:
    493         raise _type_not_supported

/home/mi/gph82/programs/PyEmma/pyemma/msm/analysis/dense/decomposition.pyc in timescales(T, tau, k, reversible, mu)
    354 
    355     """
--> 356     values = eigenvalues(T, reversible=reversible, mu=mu)
    357 
    358     """Sort by absolute value"""

/home/mi/gph82/programs/PyEmma/pyemma/msm/analysis/dense/decomposition.pyc in eigenvalues(T, k, reversible, mu)
    166         # compute stationary distribution if not given
    167         if mu is None:
--> 168             mu = stationary_distribution_from_backward_iteration(T)
    169         # symmetrize T
    170         smu = np.sqrt(mu)

/home/mi/gph82/programs/PyEmma/pyemma/msm/analysis/dense/decomposition.pyc in stationary_distribution_from_backward_iteration(P, eps)
    101     mu = 1.0 - eps
    102     x0 = np.ones(P.shape[0])
--> 103     y = backward_iteration(A, mu, x0)
    104     pi = y / y.sum()
    105     return pi

/home/mi/gph82/programs/PyEmma/pyemma/msm/analysis/dense/decomposition.pyc in backward_iteration(A, mu, x0, tol, maxiter)
     78             return y
     79     msg = "Failed to converge after %d iterations, residuum is %e" % (maxiter, r)
---> 80     raise RuntimeError(msg)
     81 
     82 

RuntimeError: Failed to converge after 100 iterations, residuum is 1.000000e+00
fabian-paul commented 9 years ago

I wouldn't expect that the estimation returns a matrix filled with zeros. We have to investigate this.

gph82 commented 9 years ago

probably it gets preallocated and never filled?

fabian-paul commented 9 years ago

Which would mean that during the iterations (of the estimator) the matrix is never changed? That sound unlikely.

fabian-paul commented 9 years ago

Ah, that was actually my mistake. If the sparse estimator fails with a warning, it always returns a zero matrix. I will fix this.

gph82 commented 9 years ago

Thanks! i think we also need a fallback-to-dense behaviour if sparse-MLE does not converge.

The reasons for not converging in the sparse mode are another issue

fabian-paul commented 9 years ago

I think the criteria for convergence in sparse and dense mode are still different. I have to adapt the sparse implementation. I don't think that a fall-back will be necessary then.

gph82 commented 9 years ago

ok, thanks!

On 07/10/2015 08:15 AM, fabian-paul wrote:

I think the criteria for convergence in sparse and dense mode are still different. I have to adapt the sparse implementation. I don't think that a fall-back will be necessary then.

— Reply to this email directly or view it on GitHub https://github.com/markovmodel/PyEMMA/issues/403#issuecomment-120244267.

Dr. Guillermo Pérez-Hernández Freie Universität Berlin Institute for Mathematics Arnimallee 6 D-14195 Berlin tel 0049 30 838 75776

franknoe commented 9 years ago

This didn't get closed automatically.