Fixes (#1490)

* raise if n_jobs is non-positive (#1486)

* [plots] deeptime compatibility for plot_flux

* [its] shape fix

* update pybind

* include deeptime as dependency

* fix regspace by using deeptime implementation

* CI: update reference numpy and build distro for pip

Author: clonker

.travis.yml

@@ -17,21 +17,21 @@ env:
     - MACOSX_DEPLOYMENT_TARGET=10.9
     - CONDA_BUILD=YES
   matrix:
-    - CONDA_PY=3.6
-      CONDA_NPY=1.11
-    - CONDA_PY=3.7
-      CONDA_NPY=1.17
+    - CONDA_PY='3.6'
+      CONDA_NPY='1.16'
+    - CONDA_PY='3.7'
+      CONDA_NPY='1.17'
+    - CONDA_PY='3.8'
+      CONDA_NPY='1.17'
 
 matrix:
   exclude:  # test only 3.8 on osx.
-    - env: CONDA_PY=3.5
-      os: osx
     - env: CONDA_PY=3.6
       os: osx
   include:  # test pip install
     - env: CONDA_PY=3.7 CONDA_BUILD=NO python="3.7"
       os: linux
-      dist: xenial # You’ll need to add dist: xenial to use Python 3.7 and higher.
+      dist: bionic # You’ll need to add dist: bionic to use Python 3.7 and higher.
       language: python
       sudo: true
 

devtools/conda-recipe/meta.yaml

@@ -35,9 +35,10 @@ requirements:
     - setuptools
     - pip
     - pybind11
+    - deeptime
 
   run:
-    - bhmm >=0.6.3,<0.7
+    - bhmm >=0.6.3
     - decorator >=4.0.0
     - h5py
     - intel-openmp # [osx]
@@ -52,6 +53,7 @@ requirements:
     - scipy
     - setuptools
     - tqdm
+    - deeptime
 
 test:
   source_files:

ext/pybind11

@@ -37,4 +37,4 @@
 from .kmeans import KmeansClustering
 from .kmeans import MiniBatchKmeansClustering
 from .regspace import RegularSpaceClustering
-from .uniform_time import UniformTimeClustering
+from .uniform_time import UniformTimeClustering

pyemma/coordinates/clustering/__init__.py

@@ -160,8 +160,8 @@ def _transform_array(self, X):
 
         # for performance reasons we pre-center the cluster centers for minRMSD.
         if self.metric == 'minRMSD' and not self._precentered:
-            self.logger.debug("precentering cluster centers for minRMSD.")
-            self._inst.precenter_centers(self.clustercenters)
+            # self.logger.debug("precentering cluster centers for minRMSD.")
+            # self._inst.precenter_centers(self.clustercenters)
             self._precentered = True
 
         dtraj = self._inst.assign(X, self.clustercenters, self.n_jobs)
@@ -279,7 +279,7 @@ def save_dtrajs(self, trajfiles=None, prefix='',
                 output_files.append(name)
         else:
             for i in range(len(self.dtrajs)):
-                if prefix is not '':
+                if prefix != '':
                     name = "%s_%i%s" % (prefix, i, extension)
                 else:
                     name = str(i) + extension

pyemma/coordinates/clustering/interface.py

@@ -31,7 +31,7 @@
 from pyemma.util.exceptions import NotConvergedWarning
 
 import numpy as np
-
+import deeptime as dt
 
 __all__ = ['RegularSpaceClustering']
 
@@ -79,6 +79,10 @@ def __init__(self, dmin, max_centers=1000, metric='euclidean', stride=1, n_jobs=
 
         """
         super(RegularSpaceClustering, self).__init__(metric=metric, n_jobs=n_jobs)
+
+        from ._ext import RMSDMetric
+        dt.clustering.metrics.register("minRMSD", RMSDMetric)
+
         self._converged = False
         self.set_params(dmin=dmin, metric=metric,
                         max_centers=max_centers, stride=stride, skip=skip)
@@ -133,16 +137,17 @@ def _estimate(self, iterable, **kwargs):
         # temporary list to store cluster centers
         clustercenters = []
         used_frames = 0
-        from ._ext import regspace
-        self._inst = regspace.Regspace_f(self.dmin, self.max_centers, self.metric, iterable.ndim)
+        regspace = dt.clustering.RegularSpace(dmin=self.dmin, max_centers=self.max_centers,
+                                              metric=self.metric, n_jobs=self.n_jobs)
+
+        # from ._ext import regspace
         it = iterable.iterator(return_trajindex=False, stride=self.stride,
                                chunk=self.chunksize, skip=self.skip)
         try:
             with it:
                 for X in it:
+                    regspace.partial_fit(X.astype(np.float32, order='C', copy=False), n_jobs=self.n_jobs)
                     used_frames += len(X)
-                    self._inst.cluster(X.astype(np.float32, order='C', copy=False),
-                                       clustercenters, self.n_jobs)
             self._converged = True
         except regspace.MaxCentersReachedException:
             self._converged = False
@@ -156,8 +161,16 @@ def _estimate(self, iterable, **kwargs):
             raise NotConvergedWarning("Used data for centers: %.2f%%" % used_data)
         finally:
             # even if not converged, we store the found centers.
-            new_shape = (len(clustercenters), iterable.ndim)
-            clustercenters = np.array(clustercenters).reshape(new_shape)
+            model = regspace.fetch_model()
+            clustercenters = model.cluster_centers.squeeze().reshape(-1, iterable.ndim)
+            self._inst = dt.clustering.ClusterModel(clustercenters, metric=self.metric)
+            from types import MethodType
+
+            def _assign(self, data, _, n_jobs):
+                out = self.transform(data, n_jobs=n_jobs)
+                return out
+
+            self._inst.assign = MethodType(_assign, self._inst)
             self.update_model_params(clustercenters=clustercenters,
                                      n_clusters=len(clustercenters))
 

pyemma/coordinates/clustering/regspace.py

@@ -2,9 +2,64 @@
 // Created by marscher on 7/17/17.
 //
 
+#include "metric.h"
 #include "regspace.h"
 #include "kmeans.h"
 
+class MaximumMetric : public Metric {
+public:
+
+    double compute_squared_d(const double* xs, const double* ys, std::size_t dim) const override {
+        return _compute(xs, ys, dim);
+    }
+    float compute_squared_f(const float* xs, const float* ys, std::size_t dim) const override {
+        return _compute(xs, ys, dim);
+    }
+private:
+    template<typename T>
+    T _compute(const T* xs, const T* ys, std::size_t dim) const {
+        T result = 0.0;
+        for (size_t i = 0; i < dim; ++i) {
+            auto d = std::abs(xs[i] - ys[i]);
+            if (d > result) {
+                result = d;
+            }
+        }
+        return result*result;
+    }
+};
+
+class RMSDMetric : public Metric {
+public:
+
+    double compute_squared_d(const double* xs, const double* ys, std::size_t dim) const override {
+        std::vector<float> xsCast (xs, xs+dim), ysCast (ys, ys+dim);
+        return _compute(xsCast.data(), ysCast.data(), dim);
+    }
+    float compute_squared_f(const float* xs, const float* ys, std::size_t dim) const override {
+        return _compute(xs, ys, dim);
+    }
+
+private:
+    template<typename T>
+    T _compute(const T* xs, const T* ys, std::size_t dim) const {
+        if (dim % 3 != 0) {
+            throw std::range_error("RMSDMetric is only implemented for input data with a dimension dividable by 3.");
+        }
+        float trace_a, trace_b;
+        auto dim3 = static_cast<const int>(dim / 3);
+        std::vector<float> buffer_b (ys, ys + dim);
+        std::vector<float> buffer_a (xs, xs + dim);
+        inplace_center_and_trace_atom_major(buffer_a.data(), &trace_a, 1, dim3);
+        inplace_center_and_trace_atom_major(buffer_b.data(), &trace_b, 1, dim3);
+
+
+        float msd = msd_atom_major(dim3, dim3, xs,
+                                   buffer_b.data(), trace_a, trace_b, 0, nullptr);
+        return msd;
+    }
+};
+
 using dtype = float;
 
 PYBIND11_MODULE(_ext, m) {
@@ -36,4 +91,7 @@ PYBIND11_MODULE(_ext, m) {
             .def("cluster_loop", &kmeans_f::cluster_loop)
             .def("init_centers_KMpp", &kmeans_f::initCentersKMpp)
             .def("cost_function", &kmeans_f::costFunction);
+    // py::class_<Metric>(m, "Metric");
+    py::object baseMetric = (py::object) py::module_::import("deeptime.clustering._clustering_bindings").attr("Metric");
+    py::class_<RMSDMetric>(m, "RMSDMetric", baseMetric).def(py::init<>());
 }

pyemma/coordinates/clustering/src/clustering_module.cpp

@@ -107,6 +107,7 @@ def dt_model(self, value):
         self._timeunit_model = TimeUnit(self._dt_model)
 
     @property
+    @alias('n_states')
     def nstates(self):
         r"""Returns the number of states.
 

pyemma/msm/models/reactive_flux.py

@@ -544,7 +544,7 @@ def plot_flux(
     """
     from matplotlib import pylab as plt
     F = flux_scale * getattr(flux, attribute_to_plot)
-    c = flux.committor
+    c = flux.forward_committor
     if state_sizes is None:
         state_sizes = flux.stationary_distribution
     plot = NetworkPlot(F, pos=pos, xpos=c, ax=ax)
@@ -554,13 +554,13 @@ def plot_flux(
 
     if isinstance(state_labels, str) and state_labels == 'auto':
         # the first and last element correspond to A and B in ReactiveFlux
-        state_labels = _np.array([str(i) for i in range(flux.nstates)])
+        state_labels = _np.array([str(i) for i in range(flux.n_states)])
         state_labels[_np.array(flux.A)] = "A"
         state_labels[_np.array(flux.B)] = "B"
     elif isinstance(state_labels, (_np.ndarray, list, tuple)):
-        if len(state_labels) != flux.nstates:
+        if len(state_labels) != flux.n_states:
             raise ValueError("length of state_labels({}) has to match length of states({})."
-                             .format(len(state_labels), flux.nstates))
+                             .format(len(state_labels), flux.n_states))
 
     fig = plot.plot_network(
         state_sizes=state_sizes, state_scale=state_scale, state_colors=state_colors,

pyemma/plots/networks.py

@@ -8,6 +8,7 @@ requires = [
  "mdtraj",
  # TODO: include again once https://github.com/pybind/pybind11/issues/1067 is fixed.
  #"pybind11",
+ "deeptime"
  ]