From e2348921953615fa52b33c3ddb9f80dea03ba552 Mon Sep 17 00:00:00 2001
From: bcumming <bcumming@cscs.ch>
Date: Tue, 19 Apr 2016 11:50:26 +0200
Subject: [PATCH] fix corner case of is_minimal_degree

---
 src/algorithms.hpp          | 19 +++++++----
 src/fvm.hpp                 |  2 --
 src/mechanism.hpp           | 67 +++++++++++++++++++++++++++++++++++++
 src/mechanism_interface.cpp | 24 +++++++++++++
 src/mechanism_interface.hpp | 32 ++++++++++++++++++
 src/parameter_list.hpp      |  4 +--
 tests/test_algorithms.cpp   | 10 +++---
 7 files changed, 143 insertions(+), 15 deletions(-)
 create mode 100644 src/mechanism.hpp
 create mode 100644 src/mechanism_interface.cpp
 create mode 100644 src/mechanism_interface.hpp

diff --git a/src/algorithms.hpp b/src/algorithms.hpp
index 10ba44d3..bc25b98b 100644
--- a/src/algorithms.hpp
+++ b/src/algorithms.hpp
@@ -1,5 +1,7 @@
 #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>
diff --git a/src/fvm.hpp b/src/fvm.hpp
index 9e84119f..d25d4669 100644
--- a/src/fvm.hpp
+++ b/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);
diff --git a/src/mechanism.hpp b/src/mechanism.hpp
new file mode 100644
index 00000000..34322258
--- /dev/null
+++ b/src/mechanism.hpp
@@ -0,0 +1,67 @@
+#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
diff --git a/src/mechanism_interface.cpp b/src/mechanism_interface.cpp
new file mode 100644
index 00000000..ba24eb06
--- /dev/null
+++ b/src/mechanism_interface.cpp
@@ -0,0 +1,24 @@
+#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
+
diff --git a/src/mechanism_interface.hpp b/src/mechanism_interface.hpp
new file mode 100644
index 00000000..389ef059
--- /dev/null
+++ b/src/mechanism_interface.hpp
@@ -0,0 +1,32 @@
+#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
diff --git a/src/parameter_list.hpp b/src/parameter_list.hpp
index e7dd1bd1..56848466 100644
--- a/src/parameter_list.hpp
+++ b/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
         }
     };
 
diff --git a/tests/test_algorithms.cpp b/tests/test_algorithms.cpp
index 9328b315..cfe058a3 100644
--- a/tests/test_algorithms.cpp
+++ b/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)
-- 
GitLab