Misc (#153)

Author: clonker

deeptime/data/_datasets.py

@@ -796,6 +796,8 @@ def prinz_potential(h=1e-5, n_steps=500, temperature_factor=1., mass=1., damping
 
     where :math:`m` is the mass, :math:`d` the damping factor, and :math:`\eta_t \sim \mathcal{N}(0, 1)`.
 
+    The locations of the minima can be accessed via the `minima` attribute.
+
     .. plot:: datasets/plot_prinz.py
 
     Parameters
@@ -844,11 +846,13 @@ def prinz_potential(h=1e-5, n_steps=500, temperature_factor=1., mass=1., damping
     .. footbibliography::
     """
     from ._data_bindings import Prinz
-    return TimeIndependentSystem(Prinz(), h, n_steps, props={
+    system = TimeIndependentSystem(Prinz(), h, n_steps, props={
         'kT': temperature_factor,
         'mass': mass,
         'damping': damping
     })
+    system.minima = [-0.73943018, -0.22373758, 0.26914935, 0.67329636]
+    return system
 
 
 def triple_well_1d(h=1e-3, n_steps=500):

deeptime/markov/hmm/_bindings/include/OutputModelUtils.h

@@ -302,11 +302,12 @@ std::tuple<np_array<dtype>, np_array<dtype>> fit(std::size_t nHiddenStates, cons
         for (decltype(nObsTrajs) k = 0; k < nObsTrajs; ++k, ++weightsIt, ++obsIt) {
             const auto &w = py::cast<np_array<dtype>>(*weightsIt);
             const auto &obs = py::cast<np_array<dtype>>(*obsIt);
+            const auto* obsPtr = obs.data();
             for (decltype(nHiddenStates) i = 0; i < nHiddenStates; ++i) {
                 dtype dot = 0;
                 dtype wStateSum = 0;
                 for (ssize_t t = 0; t < obs.shape(0); ++t) {
-                    dot += w.at(t, i) * obs.at(t);
+                    dot += w.at(t, i) * obsPtr[t];
                     wStateSum += w.at(t, i);
                 }
                 // update nominator
@@ -329,12 +330,13 @@ std::tuple<np_array<dtype>, np_array<dtype>> fit(std::size_t nHiddenStates, cons
         for (decltype(nObsTrajs) k = 0; k < nObsTrajs; ++k, ++weightsIt, ++obsIt) {
             const auto &w = py::cast<np_array<dtype>>(*weightsIt);
             const auto &obs = py::cast<np_array<dtype>>(*obsIt);
+            const auto *obsPtr = obs.data();
 
             for (decltype(nHiddenStates) i = 0; i < nHiddenStates; ++i) {
                 dtype wStateSum = 0;
                 dtype sigmaUpdate = 0;
                 for (ssize_t t = 0; t < obs.shape(0); ++t) {
-                    auto sqrty = static_cast<dtype>(obs.at(t)) - static_cast<dtype>(means.at(i));
+                    auto sqrty = static_cast<dtype>(obsPtr[t]) - static_cast<dtype>(means.at(i));
                     sigmaUpdate += w.at(t, i) * sqrty*sqrty;
                     wStateSum += w.at(t, i);
                 }

deeptime/markov/hmm/init/gaussian/_init_gaussian_impl.py

@@ -1,14 +1,14 @@
 import numpy as np
 
 
-def from_data(dtrajs, n_hidden_states, reversible):
+def from_data(trajs, n_hidden_states, reversible):
     r""" Makes an initial guess :class:`HMM <HiddenMarkovModel>` with Gaussian output model.
 
     To this end, a Gaussian mixture model is estimated using `scikit-learn <https://scikit-learn.org/>`_.
 
     Parameters
     ----------
-    dtrajs : array_like or list of array_like
+    trajs : array_like or list of array_like
         Trajectories which are used for making the initial guess.
     n_hidden_states : int
         Number of hidden states.
@@ -32,15 +32,17 @@ def from_data(dtrajs, n_hidden_states, reversible):
     import deeptime.markov.tools.analysis as msmana
     from deeptime.util.types import ensure_timeseries_data
 
-    dtrajs = ensure_timeseries_data(dtrajs)
-    collected_observations = np.concatenate(dtrajs)
+    trajs = ensure_timeseries_data(trajs)
+    collected_observations = np.concatenate(trajs)
+    if collected_observations.ndim == 1:
+        collected_observations = collected_observations[..., None]
     gmm = GaussianMixture(n_components=n_hidden_states)
-    gmm.fit(collected_observations[:, None])
+    gmm.fit(collected_observations)
     output_model = GaussianOutputModel(n_hidden_states, means=gmm.means_[:, 0], sigmas=np.sqrt(gmm.covariances_[:, 0]))
 
     # Compute fractional state memberships.
     Nij = np.zeros((n_hidden_states, n_hidden_states))
-    for o_t in dtrajs:
+    for o_t in trajs:
         # length of trajectory
         T = o_t.shape[0]
         # output probability

deeptime/markov/tools/analysis/_api.py

@@ -16,11 +16,12 @@
 from ._decomposition import timescales_from_eigenvalues
 
 from . import dense
-from . import sparse
 from . import _assessment
 from . import _mean_first_passage_time
 from . import _decomposition
 from . import _fingerprints
+from . import _committor
+from . import _expectations
 
 __docformat__ = "restructuredtext en"
 __authors__ = __author__ = "Benjamin Trendelkamp-Schroer, Martin Scherer, Jan-Hendrik Prinz, Frank Noe"
@@ -739,26 +740,15 @@ def committor(T, A, B, forward=True, mu=None):
     T = ensure_number_array(T, ndim=2)
     A = ensure_integer_array(A, ndim=1)
     B = ensure_integer_array(B, ndim=1)
-    if _issparse(T):
-        if forward:
-            return sparse.committor.forward_committor(T, A, B)
-        else:
-            """ if P is time reversible backward commitor is equal 1 - q+"""
-            if is_reversible(T, mu=mu):
-                return 1.0 - sparse.committor.forward_committor(T, A, B)
-
-            else:
-                return sparse.committor.backward_committor(T, A, B)
-
+    if forward:
+        return _committor.forward_committor(T, A, B)
     else:
-        if forward:
-            return dense.committor.forward_committor(T, A, B)
+        """ if P is time reversible backward commitor is equal 1 - q+"""
+        if is_reversible(T, mu=mu):
+            return 1.0 - _committor.forward_committor(T, A, B)
+
         else:
-            """ if P is time reversible backward commitor is equal 1 - q+"""
-            if is_reversible(T, mu=mu):
-                return 1.0 - dense.committor.forward_committor(T, A, B)
-            else:
-                return dense.committor.backward_committor(T, A, B)
+            return _committor.backward_committor(T, A, B)
 
 
 ################################################################################
@@ -811,10 +801,7 @@ def expected_counts(T, p0, N):
     T = ensure_number_array(T, ndim=2)
     p0 = ensure_floating_array(p0, ndim=1)
     # go
-    if _issparse(T):
-        return sparse.expectations.expected_counts(p0, T, N)
-    else:
-        return dense.expectations.expected_counts(p0, T, N)
+    return _expectations.expected_counts(p0, T, N)
 
 
 def expected_counts_stationary(T, N, mu=None):
@@ -867,10 +854,7 @@ def expected_counts_stationary(T, N, mu=None):
     if mu is not None:
         mu = ensure_floating_array(mu, ndim=1)
     # go
-    if _issparse(T):
-        return sparse.expectations.expected_counts_stationary(T, N, mu=mu)
-    else:
-        return dense.expectations.expected_counts_stationary(T, N, mu=mu)
+    return _expectations.expected_counts_stationary(T, N, mu=mu)
 
 
 ################################################################################

deeptime/markov/tools/analysis/dense/_committor.py deeptime/markov/tools/analysis/_committor.py

@@ -1,14 +1,35 @@
-r"""This module provides functions for the computation of forward and
-backward comittors using dense linear algebra.
+import numpy as np
+from scipy.linalg import solve
+from scipy.sparse.linalg import spsolve
 
-.. moduleauthor:: B.Trendelkamp-Schroer <benjamin DOT trendelkamp-schroer AT fu-berlin DOT de>
+import scipy.sparse as sparse
 
-"""
+from ._stationary_vector import stationary_distribution
 
-import numpy as np
-from scipy.linalg import solve
 
-from .._stationary_vector import stationary_distribution
+def _set_up_linear_system(K, A, B):
+    """Assemble left-hand side W for linear system"""
+    """Equation (I)"""
+    W = 1.0 * K
+    """Equation (II)"""
+    if sparse.issparse(W):
+        W = W.todok()
+        W[list(A), :] = 0.0
+        W.tocsr()
+        W = W + sparse.coo_matrix((np.ones(len(A)), (list(A), list(A))), shape=W.shape).tocsr()
+    else:
+        W[list(A), :] = 0.0
+        W[list(A), list(A)] = 1.0
+    """Equation (III)"""
+    if sparse.issparse(W):
+        W = W.todok()
+        W[list(B), :] = 0.0
+        W.tocsr()
+        W = W + sparse.coo_matrix((np.ones(len(B)), (list(B), list(B))), shape=W.shape).tocsr()
+    else:
+        W[list(B), :] = 0.0
+        W[list(B), list(B)] = 1.0
+    return W
 
 
 def forward_committor(T, A, B):
@@ -49,28 +70,18 @@ def forward_committor(T, A, B):
     A = set(A)
     B = set(B)
     AB = A.intersection(B)
-    notAB = X.difference(A).difference(B)
     if len(AB) > 0:
         raise ValueError("Sets A and B have to be disjoint")
     L = T - np.eye(T.shape[0])  # Generator matrix
 
-    """Assemble left hand-side W for linear system"""
-    """Equation (I)"""
-    W = 1.0 * L
-    """Equation (II)"""
-    W[list(A), :] = 0.0
-    W[list(A), list(A)] = 1.0
-    """Equation (III)"""
-    W[list(B), :] = 0.0
-    W[list(B), list(B)] = 1.0
-
+    W = _set_up_linear_system(L, A, B)
     """Assemble right hand side r for linear system"""
     """Equation (I+II)"""
     r = np.zeros(T.shape[0])
     """Equation (III)"""
     r[list(B)] = 1.0
 
-    u = solve(W, r)
+    u = solve(W, r) if not sparse.issparse(W) else spsolve(W, r)
     return u
 
 
@@ -115,29 +126,25 @@ def backward_committor(T, A, B, mu=None):
     A = set(A)
     B = set(B)
     AB = A.intersection(B)
-    notAB = X.difference(A).difference(B)
     if len(AB) > 0:
         raise ValueError("Sets A and B have to be disjoint")
     if mu is None:
         mu = stationary_distribution(T)
-    K = np.transpose(mu[:, np.newaxis] * (T - np.eye(T.shape[0])))
+    if sparse.issparse(T):
+        L = T - sparse.eye(T.shape[0], T.shape[0])
+        D = sparse.diags([mu, ], [0, ])
+        K = (D.dot(L)).T
+    else:
+        K = np.transpose(mu[:, np.newaxis] * (T - np.eye(T.shape[0])))
 
     """Assemble left-hand side W for linear system"""
-    """Equation (I)"""
-    W = 1.0 * K
-    """Equation (II)"""
-    W[list(A), :] = 0.0
-    W[list(A), list(A)] = 1.0
-    """Equation (III)"""
-    W[list(B), :] = 0.0
-    W[list(B), list(B)] = 1.0
-
+    W = _set_up_linear_system(K, A, B)
     """Assemble right-hand side r for linear system"""
     """Equation (I)+(III)"""
     r = np.zeros(T.shape[0])
     """Equation (II)"""
     r[list(A)] = 1.0
 
-    u = solve(W, r)
+    u = solve(W, r) if not sparse.issparse(W) else spsolve(W, r)
 
     return u

deeptime/markov/tools/analysis/dense/_expectations.py deeptime/markov/tools/analysis/_expectations.py

@@ -6,9 +6,9 @@
 """
 
 import numpy as np
-
-from .._decomposition import rdl_decomposition
-from .._stationary_vector import stationary_distribution
+import scipy.sparse as sparse
+from ._decomposition import rdl_decomposition
+from ._stationary_vector import stationary_distribution
 
 
 def expected_counts(p0, T, n):
@@ -40,7 +40,7 @@ def expected_counts(p0, T, n):
 
     """
     M = T.shape[0]
-    if n <= M:
+    if sparse.issparse(T) or n <= M:
         return ec_matrix_vector(p0, T, n)
     else:
         return ec_geometric_series(p0, T, n)
@@ -72,12 +72,19 @@ def expected_counts_stationary(T, n, mu=None):
 
     """
     if n <= 0:
-        EC = np.zeros(T.shape)
+        if sparse.issparse(T):
+            EC = sparse.coo_matrix(T.shape, dtype=float)
+        else:
+            EC = np.zeros(T.shape)
         return EC
     else:
         if mu is None:
-            mu = stationary_distribution(T)
-        EC = n * mu[:, np.newaxis] * T
+            mu = stationary_distribution(T, check_inputs=False)
+        if sparse.issparse(T):
+            D_mu = sparse.diags(mu, 0)
+            EC = n * D_mu.dot(T)
+        else:
+            EC = n * mu[:, np.newaxis] * T
         return EC
 
 
@@ -158,21 +165,29 @@ def ec_matrix_vector(p0, T, n):
         Expected value for transition counts after N steps.
 
     """
-    if (n <= 0):
-        EC = np.zeros(T.shape)
-        return EC
+    if n <= 0:
+        if sparse.issparse(T):
+            return sparse.coo_matrix(T.shape, dtype=float)
+        else:
+            return np.zeros(T.shape)
     else:
         """Probability vector after (k=0) propagations"""
         p_k = 1.0 * p0
         """Sum of vectors after (k=0) propagations"""
         p_sum = 1.0 * p_k
+        """Transpose T to use sparse dot product"""
+        Tt = T.transpose()
         for k in range(n - 1):
             """Propagate one step p_{k} -> p_{k+1}"""
-            p_k = np.dot(p_k, T)
+            p_k = Tt.dot(p_k)
             """Update sum"""
             p_sum += p_k
         """Expected counts"""
-        EC = p_sum[:, np.newaxis] * T
+        if sparse.issparse(T):
+            D_psum = sparse.diags(p_sum, 0)
+            EC = D_psum.dot(T)
+        else:
+            EC = p_sum[:, np.newaxis] * T
         return EC
 
 
@@ -208,7 +223,7 @@ def ec_geometric_series(p0, T, n):
         Expected value for transition counts after N steps.
 
     """
-    if (n <= 0):
+    if n <= 0:
         EC = np.zeros(T.shape)
         return EC
     else:

deeptime/markov/tools/analysis/dense/__init__.py

@@ -1,6 +1,4 @@
-from . import _committor as committor
 from . import _correlations as correlations
-from . import _expectations as expectations
 from . import _hitting_probability as hitting_probability
 from . import _sensitivity as sensitivity
 from . import _pcca as pcca