Merge pull request #120 from marscher/rdl_discard_imag_for_reversible_input

[analyis/RDL] discard imaginary part for reversible input/norm

Author: marscher

msmtools/analysis/api.py

@@ -718,7 +718,7 @@ def mfpt(T, target, origin=None, tau=1, mu=None):
 
 
 def hitting_probability(T, target):
-    """
+    r"""
     Computes the hitting probabilities for all states to the target states.
 
     The hitting probability of state i to the target set A is defined as the minimal,
@@ -732,7 +732,7 @@ def hitting_probability(T, target):
     ----------
     T : (M, M) ndarray or scipy.sparse matrix
         Transition matrix
-    B : array_like
+    target: array_like
         List of integer state labels for the target set
 
     Returns

msmtools/analysis/dense/decomposition.py

@@ -74,7 +74,7 @@ def eigenvalues(T, k=None, reversible=False, mu=None):
     if reversible:
         try:
             evals = eigenvalues_rev(T, k=k, mu=mu)
-        except:
+        except ValueError:
             evals = eigvals(T).real  # use fallback code but cast to real
     else:
         evals = eigvals(T)  # nonreversible
@@ -314,6 +314,9 @@ def rdl_decomposition(T, k=None, reversible=False, norm='standard', mu=None):
     else:
         R, D, L = rdl_decomposition_nrev(T, norm=norm)
 
+    if reversible or norm == 'reversible':
+        D = D.real
+
     if k is None:
         return R, D, L
     else:

msmtools/analysis/sparse/decomposition.py

@@ -296,9 +296,13 @@ def rdl_decomposition(T, k=None, norm='auto', ncv=None, reversible=False, mu=Non
         else:
             norm = 'standard'
     if reversible:
-        return rdl_decomposition_rev(T, k, norm=norm, ncv=ncv, mu=mu)
+        R, D, L = rdl_decomposition_rev(T, k, norm=norm, ncv=ncv, mu=mu)
     else:
-        return rdl_decomposition_nrev(T, k, norm=norm, ncv=ncv)
+        R, D, L = rdl_decomposition_nrev(T, k, norm=norm, ncv=ncv)
+
+    if reversible or norm == 'reversible':
+        D = D.real
+    return R, D, L
 
 
 def rdl_decomposition_nrev(T, k, norm='standard', ncv=None):

msmtools/estimation/dense/__init__.py

@@ -21,5 +21,5 @@
 from . import transition_matrix
 from . import covariance
 
-#from .mle import mle_trev, mle_trev_given_pi
+from .mle import mle_trev, mle_trev_given_pi
 from .tmat_sampling import tmatrix_sampler

msmtools/estimation/dense/tmat_sampling/tmatrix_sampler.py

@@ -62,8 +62,7 @@ def __init__(self, C, reversible=False, mu=None, P0=None, nsteps=1, prior='spars
                 nsteps = 1  # just force to 1, because this is independent sampling
                 self.sampler = SamplerNonRev(C-1.0)
             else:
-                msg = "Non reversible sampling with fixed stationary vector not implemented"
-                raise ValueError(msg)
+                raise ValueError('Non reversible sampling with fixed stationary vector not implemented')
 
         # remember number of steps to decorrelate between samples
         self.nsteps = nsteps

msmtools/estimation/sparse/__init__.py

@@ -28,3 +28,5 @@
 from . import likelihood
 from . import transition_matrix
 from . import prior
+
+from .mle import mle_trev, mle_trev_given_pi

tests/analysis/test_decomposition.py

@@ -32,7 +32,7 @@
 from msmtools.util.exceptions import SpectralWarning, ImaginaryEigenValueWarning
 from msmtools.util.birth_death_chain import BirthDeathChain
 
-from msmtools.analysis import stationary_distribution, eigenvalues, eigenvectors
+from msmtools.analysis import stationary_distribution, eigenvalues, eigenvectors, is_reversible
 from msmtools.analysis import rdl_decomposition, timescales
 
 ################################################################################
@@ -162,6 +162,7 @@ def test_eigenvectors_reversible(self):
 
     def test_rdl_decomposition(self):
         P = self.bdc.transition_matrix()
+        assert is_reversible(P)
         mu = self.bdc.stationary_distribution()
 
         """Non-reversible"""
@@ -194,6 +195,7 @@ def test_rdl_decomposition(self):
 
         """k=None"""
         Rn, Dn, Ln = rdl_decomposition(P, norm='reversible')
+        assert Dn.dtype in (np.float32, np.float64)
         Xn = np.dot(Ln, Rn)
         """Right-eigenvectors"""
         assert_allclose(np.dot(P, Rn), np.dot(Rn, Dn))