Implied timescales plot

CMakeLists.txt

@@ -9,7 +9,7 @@ if(DEFINED PROJECT_NAME)
 endif()
 project(deeptime LANGUAGES C CXX VERSION 0.0.0)
 
-find_package(OpenMP)
+find_package(OpenMP 4)
 if(OpenMP_FOUND)
     add_definitions(-DUSE_OPENMP)
 endif()

deeptime/__init__.py

@@ -14,6 +14,11 @@
 from . import decomposition
 from . import markov
 
+from .util.platform import module_available
+if module_available("matplotlib"):
+    from . import plots
+del module_available
+
 
 def capi_includes(inc_clustering: bool = False, inc_markov: bool = False, inc_markov_hmm: bool = False,
                   inc_data: bool = False):

deeptime/data/_drunkards_walk_simulator.py

@@ -3,7 +3,7 @@
 
 import numpy as np
 
-from deeptime.util.decorators import plotting_function
+from deeptime.util.decorators import plotting_function_with_networkx
 
 Coordinate = Tuple[int, int]
 
@@ -245,7 +245,7 @@ def walk(self, start: Coordinate, n_steps: int, stop: bool = True, return_states
             return states
 
     @staticmethod
-    @plotting_function
+    @plotting_function_with_networkx
     def plot_path(ax, path, intermediates: bool = True, color_lerp: bool = True, **plot_kw):  # pragma: no cover
         r""" Plots a path onto a drunkard's walk map. The path is interpolated with splines and, if desired, has a
         linearly interpolated color along its path.
@@ -286,7 +286,7 @@ def plot_path(ax, path, intermediates: bool = True, color_lerp: bool = True, **p
         else:
             ax.plot(xint, yint, **plot_kw)
 
-    @plotting_function
+    @plotting_function_with_networkx
     def plot_network(self, ax, F, cmap=None, connection_threshold: float = 0.):  # pragma: no cover
         r""" Plots a flux network onto a matplotlib axis assuming that the states are ordered according to a iota
         range (i.e., {0, 1, ..., n_states-1}) on the grid.
@@ -335,7 +335,7 @@ def plot_network(self, ax, F, cmap=None, connection_threshold: float = 0.):  # p
                                connectionstyle='arc3, rad=0.1')
         return edge_colors
 
-    @plotting_function
+    @plotting_function_with_networkx
     def plot_2d_map(self, ax, barriers: bool = True, barrier_mode: str = 'filled'):  # pragma: no cover
         r""" Plots the drunkard's map onto an already existing matplotlib axis. One can select whether to draw barriers
         and if barriers should be drawn, whether they have a 'filled' face or be 'hollow' with just the border

deeptime/decomposition/_koopman.py

@@ -642,6 +642,12 @@ def timescales(self, k=None, lagtime: Optional[int] = None) -> np.ndarray:
             lagtime = self._cov.lagtime
         return - lagtime / np.log(np.abs(self.singular_values[:k]))
 
+    @property
+    def lagtime(self):
+        r""" The lagtime corresponding to this model. See also
+        :meth:`CovarianceModel.lagtime <deeptime.covariance.CovarianceModel.lagtime>`. """
+        return self._cov.lagtime
+
     @property
     def feature_component_correlation(self):
         r"""Instantaneous correlation matrix between mean-free input features and projection components.

deeptime/decomposition/_vamp.py

@@ -513,7 +513,7 @@ def chapman_kolmogorov_validator(self, mlags, test_model: CovarianceKoopmanModel
         The Champan-Kolmogorov test is to compare the predictions to the estimates.
         """
         test_model = self.fetch_model() if test_model is None else test_model
-        assert test_model is not None, "We need a test model via argument or an estimator which was already" \
+        assert test_model is not None, "We need a test model via argument or an estimator which was already " \
                                        "fit to data."
         lagtime = self.lagtime
         if n_observables is not None:

deeptime/markov/_base.py

@@ -168,6 +168,26 @@ def evaluate_samples(self, quantity, delimiter='/', *args, **kwargs):
         from deeptime.util.stats import evaluate_samples as _eval
         return _eval(self.samples, quantity=quantity, delimiter=delimiter, *args, **kwargs)
 
+    def timescales(self, k=None):
+        r""" Relaxation timescales corresponding to the eigenvalues.
+
+        Parameters
+        ----------
+        k : int, optional, default=None
+            The number of timescales (excluding the stationary process).
+
+        Returns
+        -------
+        timescales : tuple(iterable, iterable)
+            Timescales of the prior and timescales of the samples.
+        """
+        return self.prior.timescales(k=k), self.evaluate_samples('timescales', k=k)
+
+    @property
+    def lagtime(self):
+        r"""Lagtime of the models."""
+        return self.prior.lagtime
+
 
 class MembershipsChapmanKolmogorovValidator(LaggedModelValidator):
     r""" Validates a model estimated at lag time tau by testing its predictions for longer lag times.

deeptime/markov/hmm/_hidden_markov_model.py

@@ -303,6 +303,12 @@ def transition_counts(self) -> Optional[np.ndarray]:
         """
         return self.transition_model.count_model.count_matrix if self.transition_model.count_model is not None else None
 
+    def timescales(self, k=None):
+        r""" Yields the timescales of the hidden transition model.
+        See :meth:`MarkovStateModel.timescales <deeptime.markov.msm.MarkovStateModel.timescales>`.
+        """
+        return self.transition_model.timescales(k=k)
+
     @property
     def hidden_state_trajectories(self) -> Optional[List[np.ndarray]]:
         r"""

deeptime/markov/msm/_bayesian_msm.py

@@ -173,7 +173,7 @@ def fetch_model(self) -> Optional[BayesianPosterior]:
         """
         return self._model
 
-    def fit(self, data, callback: Callable = None):
+    def fit(self, data, callback: Callable = None, **kw):
         """
         Performs the estimation on either a count matrix or a previously estimated TransitionCountModel.
 
@@ -185,23 +185,23 @@ def fit(self, data, callback: Callable = None):
         callback: callable, optional, default=None
             Function to be called to indicate progress of sampling.
 
+        Other Parameters
+        ----------------
+        ignore_counting_mode : bool, optional, default=False
+            Method does not raise if counting mode isn't of the "effective" family. Use with caution.
+
         Returns
         -------
         self : BayesianMSM
             Reference to self.
         """
-        from deeptime.markov import TransitionCountModel
-        if isinstance(data, TransitionCountModel) and data.counting_mode is not None \
-                and "effective" not in data.counting_mode:
-            raise ValueError("The transition count model was not estimated using an effective counting method, "
-                             "therefore counts are likely to be strongly correlated yielding wrong confidences.")
-
         if isinstance(data, Estimator):
             if data.has_model:
                 data = data.fetch_model()
             else:
                 raise ValueError("Can only use estimators as input if they have been fit previously.")
 
+        from deeptime.markov import TransitionCountModel
         if isinstance(data, TransitionCountModel) or is_square_matrix(data):
             msm = MaximumLikelihoodMSM(
                 reversible=self.reversible, stationary_distribution_constraint=self.stationary_distribution_constraint,
@@ -214,7 +214,7 @@ def fit(self, data, callback: Callable = None):
                              "TransitionCountEstimator) or a MarkovStateModel instance or an estimator producing "
                              "Markov state models.")
 
-        return self.fit_from_msm(msm, callback=callback)
+        return self.fit_from_msm(msm, callback=callback, **kw)
 
     def sample(self, prior: MarkovStateModel, n_samples: int, n_steps: Optional[int] = None, callback=None):
         r""" Performs sampling based on a prior.
@@ -275,7 +275,7 @@ def sample(self, prior: MarkovStateModel, n_samples: int, n_steps: Optional[int]
         ]
         return samples
 
-    def fit_from_msm(self, msm: MarkovStateModel, callback=None):
+    def fit_from_msm(self, msm: MarkovStateModel, callback=None, **kw):
         r""" Fits a bayesian posterior from a given Markov state model. The MSM must contain a count model to be able
         to produce confidences. Note that the count model should be produced using effective counting, otherwise
         counts are correlated and computed confidences are wrong.
@@ -287,6 +287,11 @@ def fit_from_msm(self, msm: MarkovStateModel, callback=None):
         callback : callable, optional, default=None
             Function to be called to indicate progress of sampling.
 
+        Other Parameters
+        ----------------
+        ignore_counting_mode : bool, optional, default=False
+            Method does not raise if counting mode isn't of the "effective" family. Use with caution.
+
         Returns
         -------
         self : BayesianMSM
@@ -295,6 +300,11 @@ def fit_from_msm(self, msm: MarkovStateModel, callback=None):
         if not msm.has_count_model:
             raise ValueError("Can only sample confidences with a count model. The counting mode should be 'effective'"
                              " to avoid correlations between counts and therefore wrong confidences.")
+        if not kw.get("ignore_counting_mode", False) \
+                and msm.count_model.counting_mode is not None and "effective" not in msm.count_model.counting_mode:
+            raise ValueError("The transition count model was not estimated using an effective counting method, "
+                             "therefore counts are likely to be strongly correlated yielding wrong confidences. "
+                             "To ignore this, set `ignore_counting_mode` to True in the call to `fit`.")
         # use the same count matrix as the MLE. This is why we have effective as a default
         samples = self.sample(msm, self.n_samples, self.n_steps, callback)
         self._model = BayesianPosterior(prior=msm, samples=samples)
@@ -336,4 +346,4 @@ def _ck_estimate_model_for_lag(estimator, model, data, lagtime):
 class BayesianMSMChapmanKolmogorovValidator(MembershipsChapmanKolmogorovValidator):
 
     def fit(self, data, n_jobs=None, progress=None, **kw):
-        return super().fit(data, n_jobs, progress, estimate_model_for_lag=_ck_estimate_model_for_lag, **kw)
+        return super().fit(data, n_jobs, progress, estimate_model_for_lag=_ck_estimate_model_for_lag, **kw)

deeptime/markov/msm/_markov_state_model.py

@@ -48,6 +48,10 @@ class MarkovStateModel(Model):
     transition_matrix_tolerance : float, default=1e-8
         The tolerance under which a matrix is still considered a transition matrix (only non-negative elements and
         row sums of 1).
+    lagtime : int, optional, default=None
+        The lagtime of this MSM. If there is a count model, the MSM assumes the lagtime of the count model, otherwise
+        falls back to the lagtime set via this constructor argument or a lagtime of `1` if no lagtime is provided
+        at all.
 
     See Also
     --------
@@ -64,11 +68,12 @@ class MarkovStateModel(Model):
     """
 
     def __init__(self, transition_matrix, stationary_distribution=None, reversible=None,
-                 n_eigenvalues=None, ncv=None, count_model=None, transition_matrix_tolerance=1e-8):
+                 n_eigenvalues=None, ncv=None, count_model=None, transition_matrix_tolerance=1e-8, lagtime=None):
         super().__init__()
         self._is_reversible = reversible
         self._ncv = ncv
         self._transition_matrix_tolerance = transition_matrix_tolerance
+        self._lagtime = lagtime
 
         self.update_transition_matrix(transition_matrix)
 
@@ -204,7 +209,7 @@ def to_koopman_model(self, empirical: bool = True, epsilon: float = 1e-10):
     @property
     def lagtime(self) -> int:
         r""" The lagtime this model was estimated at. In case no count model was provided, this property defaults
-        to a lagtime of `1`.
+        to the lagtime set in the constructor or a lagtime of :math:`\tau = 1` if it was left `None`.
 
         Returns
         -------
@@ -213,7 +218,7 @@ def lagtime(self) -> int:
         """
         if self.count_model is not None:
             return self.count_model.lagtime
-        return 1
+        return 1 if self._lagtime is None else self._lagtime
 
     @property
     def transition_matrix(self):
@@ -474,20 +479,19 @@ def eigenvectors_right(self, k=None):
 
     def timescales(self, k=None):
         r"""
-        The relaxation timescales corresponding to the eigenvalues
+        Relaxation timescales corresponding to the eigenvalues.
 
         Parameters
         ----------
-        k : int
-            number of timescales to be returned. By default all available
-            eigenvalues, minus 1.
+        k : int, optional, default=None
+            Number of timescales to be returned. By default, this uses all available eigenvalues except for the
+            first (stationary) eigenvalue.
 
         Returns
         -------
         ts : ndarray(m)
-            relaxation timescales in units of the input trajectory time step,
-            defined by :math:`-\tau / ln | \lambda_i |, i = 2,...,k+1`.
-
+            Relaxation timescales in units of the input trajectory time step,
+            defined by :math:`-\tau / \mathrm{ln} | \lambda_i |, i = 2,...,k+1`.
         """
         if k is None:
             self._ensure_eigenvalues()

deeptime/markov/tools/estimation/dense/tmat_sampling/tmatrix_sampler.py

@@ -55,16 +55,13 @@ def sample(self, nsamples=1, return_statdist=False, callback=None):
         else:
             n = self.count_matrix.shape[0]
             P_samples = np.zeros((nsamples, n, n))
-            if return_statdist:
-                pi_samples = np.zeros((nsamples, n))
-                for i in range(nsamples):
-                    P_samples[i, :, :], pi_samples[i, :] = self.sampler.sample(N=self.n_steps, return_statdist=True)
-                    if callback is not None:
-                        callback()
-                return P_samples, pi_samples
-            else:
-                for i in range(nsamples):
-                    P_samples[i, :, :] = self.sampler.sample(N=self.n_steps, return_statdist=False)
-                    if callback is not None:
-                        callback()
-                return P_samples
+            pi_samples = np.zeros((nsamples, n)) if return_statdist else None
+            for i in range(nsamples):
+                out = self.sampler.sample(N=self.n_steps, return_statdist=return_statdist)
+                if return_statdist:
+                    P_samples[i, :, :], pi_samples[i, :] = out
+                else:
+                    P_samples[i, :, :] = out
+                if callback is not None:
+                    callback()
+            return P_samples if not return_statdist else (P_samples, pi_samples)

deeptime/numeric/_utils.py

@@ -46,13 +46,13 @@ def is_diagonal_matrix(matrix: _np.ndarray) -> bool:
     return _np.all(matrix == _np.diag(_np.diagonal(matrix)))
 
 
-def is_square_matrix(arr: _np.ndarray) -> bool:
+def is_square_matrix(arr) -> bool:
     r""" Determines whether an array is a square matrix. This means that ndim must be 2 and shape[0] must be equal
     to shape[1].
 
     Parameters
     ----------
-    arr : ndarray or sparse array
+    arr : object
         The array to check.
 
     Returns

deeptime/plots/__init__.py

@@ -0,0 +1,12 @@
+r"""
+.. currentmodule: deeptime.plots
+
+.. autosummary::
+    :toctree: generated/
+    :template: class_nomodule.rst
+
+    ImpliedTimescalesData
+    plot_implied_timescales
+"""
+
+from .implied_timescales import plot_implied_timescales, ImpliedTimescalesData

deeptime/plots/chapman_kolmogorov.py

@@ -0,0 +1,8 @@
+
+
+class ChapmanKolmogorovData:
+    pass
+
+
+def plot_chapman_kolmogorov_test(ax, data: ChapmanKolmogorovData):
+    pass