fix corner case of is_minimal_degree

src/algorithms.hpp

 #pragma once
 
+#include <iostream>
+
 #include <algorithm>
 #include <numeric>
 #include <type_traits>
@@ -85,14 +87,19 @@ namespace algorithms{
             "is_minimal_degree only applies to integral types"
         );
 
+        if(c.size()==0u) {
+            return true;
+        }
+
         using value_type = typename C::value_type;
-        auto i = value_type(0);
-        for(auto v : c) {
-            if(i++<v) {
-                return false;
-            }
+        if(c[0] != value_type(0)) {
+            return false;
         }
-        return true;
+        auto i = value_type(1);
+        auto it = std::find_if(
+            c.begin()+1, c.end(), [&i](value_type v) { return v>=(i++); }
+        );
+        return it==c.end();
     }
 
     template <typename C>

src/fvm.hpp

@@ -149,8 +149,6 @@ fvm_cell<T, I>::fvm_cell(nest::mc::cell const& cell)
                 face_alpha_[i] = area_face  / (c_m * r_L * c.length);
                 cv_capacitance_[i] = c_m;
 
-                std::cout << "radius " << radius_center << ", c_m " << c_m << ", r_L " << r_L << ", dx " << c.length << "\n";
-
                 auto halflen = c.length/2;
 
                 auto al = math::area_frustrum(halflen, left(c.radius), radius_center);

src/mechanism.hpp

+#pragma once
+
+#pragma once
+
+#include <memory>
+#include <string>
+
+#include "matrix.hpp"
+#include "parameter_list.hpp"
+#include "util.hpp"
+
+namespace nest {
+namespace mc {
+
+enum class mechanismKind {point, density};
+
+template <typename T, typename I>
+class mechanism {
+    using value_type  = T;
+    using size_type   = I;
+
+    // define storage types
+    using vector_type = memory::HostVector<value_type>;
+    using view_type   = typename vector_type::view_type;
+    using index_type  = memory::HostVector<size_type>;
+    using index_view  = typename index_type::view_type;
+    //using indexed_view= IndexedView<value_type, size_type, target()>;
+
+    using matrix_type = matrix<value_type, size_type>;
+
+    mechanism(matrix_type *matrix, index_view node_indices)
+    :   matrix_(matrix)
+    ,   node_indices_(node_indices)
+    {}
+
+    std::size_t size() const
+    {
+        return node_indices_.size();
+    }
+
+    virtual void set_params(value_type t_, value_type dt_) = 0;
+    virtual std::string name() const = 0;
+    virtual std::size_t memory() const = 0;
+    virtual void init()     = 0;
+    virtual void state()    = 0;
+    virtual void current()  = 0;
+
+    virtual mechanismKind kind() const = 0;
+
+    matrix_type* matrix_;
+    index_type node_indices_;
+};
+
+template <typename T, typename I>
+using mechanism_ptr = std::unique_ptr<mechanism<T,I>>;
+
+template <typename M>
+mechanism_ptr<typename M::value_type, typename M::size_type>
+make_mechanism(
+    typename M::matrix_type* matrix,
+    typename M::index_view   node_indices
+) {
+    return util::make_unique<M>(matrix, node_indices);
+}
+
+} // namespace nest
+} // namespace mc

src/mechanism_interface.cpp

+#include "mechanism_interface.hpp"
+
+/*  include the mechanisms
+#include <mechanisms/hh.hpp>
+#include <mechanisms/pas.hpp>
+*/
+
+namespace nest {
+namespace mc {
+namespace mechanisms {
+
+std::map<std::string, mechanism_helper<double, int>> mechanism_map;
+
+void setup_mechanism_helpers() {
+    /*  manually insert
+    mechanism_map["hh"]  = mechanisms::hh;
+    mechanism_map["pas"] = mechanisms::pas;
+    */
+}
+
+} // namespace mechanisms
+} // namespace nest
+} // namespace mc
+

src/mechanism_interface.hpp

+#pragma once
+
+#include <map>
+#include <string>
+
+#include "matrix.hpp"
+#include "mechanism.hpp"
+#include "parameter_list.hpp"
+
+namespace nest {
+namespace mc {
+
+///
+template <typename T, typename I>
+struct mechanism_helper {
+    using index_type = memory::HostView<I>;
+    using index_view = typename index_type::view_type;
+    using mechanism_type = mechanism<T, I>;
+    using matrix_type = typename mechanism_type::matrix_type;
+
+    virtual std::string name() const = 0;
+
+    virtual mechanism_type new_mechanism(matrix_type*, index_view) const = 0;
+
+    virtual void set_parameters(mechanism_type&, parameter_list const&) const = 0;
+};
+
+extern std::map<std::string, mechanism_helper<double, int>> mechanism_helpers;
+void setup_mechanism_helpers();
+
+} // namespace nest
+} // namespace mc

src/parameter_list.hpp

@@ -146,8 +146,8 @@ namespace mc {
         membrane_parameters()
         : base("membrane")
         {
-            base::add_parameter({"c_m",   0.01, {0., 1e9}}); // typically 10 nF/mm^2 == 0.01 F/m2
-            base::add_parameter({"r_L", 180.00, {0., 1e9}}); // Ohm.cm
+            base::add_parameter({"r_L",   0.01, {0., 1e9}}); // typically 10 nF/mm^2 == 0.01 F/m2
+            base::add_parameter({"c_m", 180.00, {0., 1e9}}); // Ohm.cm
         }
     };
 

tests/test_algorithms.cpp

@@ -47,11 +47,6 @@ TEST(algorithms, minimal_degree)
         EXPECT_TRUE(nest::mc::algorithms::is_minimal_degree(v));
     }
 
-    {
-        std::vector<int> v = {0, 1, 2, 3, 4};
-        EXPECT_TRUE(nest::mc::algorithms::is_minimal_degree(v));
-    }
-
     {
         std::vector<int> v = {0, 0, 1, 2, 3, 4};
         EXPECT_TRUE(nest::mc::algorithms::is_minimal_degree(v));
@@ -76,6 +71,11 @@ TEST(algorithms, minimal_degree)
         std::vector<int> v = {0, 2};
         EXPECT_FALSE(nest::mc::algorithms::is_minimal_degree(v));
     }
+
+    {
+        std::vector<int> v = {0, 1, 2};
+        EXPECT_FALSE(nest::mc::algorithms::is_minimal_degree(v));
+    }
 }
 
 TEST(algorithms, is_strictly_monotonic_increasing)