diff --git a/src/algorithms.hpp b/src/algorithms.hpp
index 250cade2905ffe8aa45a9554dab196412d574aff..10ba44d30a7721e2ca58cf7e98bf399d99c05f24 100644
--- a/src/algorithms.hpp
+++ b/src/algorithms.hpp
@@ -35,6 +35,11 @@ namespace algorithms{
template <typename C>
C make_index(C const& c)
{
+ static_assert(
+ std::is_integral<typename C::value_type>::value,
+ "make_index only applies to integral types"
+ );
+
C out(c.size()+1);
out[0] = 0;
std::partial_sum(c.begin(), c.end(), out.begin()+1);
@@ -69,15 +74,21 @@ namespace algorithms{
);
}
- template <typename C>
+ template <
+ typename C,
+ typename = typename std::enable_if<std::is_integral<typename C::value_type>::value>
+ >
bool is_minimal_degree(C const& c)
{
static_assert(
std::is_integral<typename C::value_type>::value,
"is_minimal_degree only applies to integral types"
);
- for(auto i=0; i<c.size(); ++i) {
- if(i<c[i]) {
+
+ using value_type = typename C::value_type;
+ auto i = value_type(0);
+ for(auto v : c) {
+ if(i++<v) {
return false;
}
}
@@ -91,8 +102,8 @@ namespace algorithms{
std::is_integral<typename C::value_type>::value,
"is_positive only applies to integral types"
);
- for(auto i=0; i<c.size(); ++i) {
- if(c[i]<1) {
+ for(auto v : c) {
+ if(v<1) {
return false;
}
}
diff --git a/src/fvm.hpp b/src/fvm.hpp
index bd1339c01c09539e52539c418728fe9952863486..0b4cb1d09a63df689d7aa2c88908e7cff079c05d 100644
--- a/src/fvm.hpp
+++ b/src/fvm.hpp
@@ -29,11 +29,29 @@ class fvm_cell {
/// view into index container
using index_view = typename index_type::view_type;
+ /// the container used for values
+ using vector_type = memory::HostVector<value_type>;
+ /// view into value container
+ using vector_view = typename vector_type::view_type;
+
+ /// constructor
fvm_cell(nest::mc::cell const& cell);
+ /// build the matrix for a given time step
+ void setup_matrx(value_type dt);
+
private:
+ /// the linear system for implicit time stepping of cell state
matrix<value_type, size_type> matrix_;
+
+ /// cv_areas_[i] is the surface area of CV i
+ vector_type cv_areas_;
+
+ /// alpha_[i] is the following value at the CV face between
+ /// CV i and its parent, required when constructing linear system
+ /// alpha_[i] = area_face / (c_m * r_L * delta_x);
+ vector_type alpha_;
};
////////////////////////////////////////////////////////////
@@ -43,42 +61,69 @@ class fvm_cell {
template <typename T, typename I>
fvm_cell<T, I>::fvm_cell(nest::mc::cell const& cell)
: matrix_(cell.parent_index())
+, cv_areas_(matrix_.size())
+, alpha_(matrix_.size())
{
+ using util::left;
+ using util::right;
+
+ auto parent_index = matrix_.p();
auto const& segment_index = cell.segment_index();
- auto num_fv = segment_index.back();
- // storage for the volume and surface area of the finite volumes
- std::vector<value_type> fv_areas(num_fv);
- std::vector<value_type> fv_volumes(num_fv);
+ // Use the membrane parameters for the first segment everywhere
+ // in the cell to start with.
+ // This has to be extended to use compartment/segment specific
+ // membrane properties.
+ auto membrane_params = cell.segments()[0]->mechanism("membrane");
+ auto c_m = membrane_params.get("c_m").value;
+ auto r_L = membrane_params.get("r_L").value;
auto seg_idx = 0;
for(auto const& s : cell.segments()) {
if(auto soma = s->as_soma()) {
- // make the assumption that the soma is at 0
+ // assert the assumption that the soma is at 0
if(seg_idx!=0) {
throw std::domain_error(
"FVM lowering encountered soma with non-zero index"
);
}
- fv_volumes[0] += math::volume_sphere(soma->radius());
- fv_areas[0] += math::area_sphere (soma->radius());
+ cv_areas_[0] += math::area_sphere(soma->radius());
+ // d[0] += cv_areas_[0];
}
else if(auto cable = s->as_cable()) {
- using util::left;
- using util::right;
-
+ // loop over each compartment in the cable
+ // each compartment has the face between two CVs at its centre
+ // the centers of the CVs are the end points of the compartment
+ //
+ // __________________________________
+ // | ........ | .cvleft. | cv |
+ // | ........ L ........ C R
+ // |__________|__________|__________|
+ //
+ // The compartment has end points marked L and R (left and right).
+ // The left compartment is assumed to be closer to the soma
+ // (i.e. it follows the minimal degree ordering)
+ // The face is at the center, marked C.
+ // The full control volume to the left (marked with .)
for(auto c : cable->compartments()) {
- auto rhs = segment_index[seg_idx] + c.index;
- auto lhs = matrix_.p()[rhs];
+ auto i = segment_index[seg_idx] + c.index;
+ auto j = parent_index[i];
- auto rad_C = math::mean(left(c.radius), right(c.radius));
- auto len = c.length/2;
+ auto radius_center = math::mean(c.radius);
+ auto area_face = math::area_circle( radius_center );
+ alpha_[i] = area_face / (c_m * r_L * c.length);
- fv_volumes[lhs] += math::volume_frustrum(len, left(c.radius), rad_C);
- fv_areas[lhs] += math::area_frustrum (len, left(c.radius), rad_C);
+ auto halflen = c.length/2;
- fv_volumes[rhs] += math::volume_frustrum(len, right(c.radius), rad_C);
- fv_areas[rhs] += math::area_frustrum (len, right(c.radius), rad_C);
+ auto al = math::area_frustrum(halflen, left(c.radius), radius_center);
+ auto ar = math::area_frustrum(halflen, right(c.radius), radius_center);
+ cv_areas_[j] += al;
+ cv_areas_[i] += ar;
+
+ // d[j] += al;
+ // d[i] += ar;
+ // l[i] -= alpha_ij;
+ // u[i] -= alpha_ij;
}
}
else {
@@ -88,11 +133,45 @@ fvm_cell<T, I>::fvm_cell(nest::mc::cell const& cell)
}
++seg_idx;
}
+}
- for(auto i=0; i<num_fv; ++i) {
- //auto area = fv_areas[i];
+template <typename T, typename I>
+void fvm_cell<T, I>::setup_matrx(T dt)
+{
+ using memory::all;
+
+ // convenience accesors to matrix storage
+ auto l = matrix_.l();
+ auto d = matrix_.d();
+ auto u = matrix_.u();
+ auto p = matrix_.p();
+
+ // The matrix has the following layout in memory
+ // where j is the parent index of i, i.e. i<j
+ //
+ // d[i] is the diagonal entry at a_ii
+ // u[i] is the upper triangle entry at a_ji
+ // l[i] is the lower triangle entry at a_ij
+ //
+ // d[j] . . u[i]
+ // . . .
+ // . . .
+ // l[i] . . d[i]
+ //
+ d(all) = 0;
+ for(auto i : d.range()) {
+ auto a = dt * alpha_[i];
+
+ // add area of CV and contribution from face with parent CV to diagonal
+ d[i] += cv_areas_[i] + a;
+
+ // add contribution to the diagonal of parent
+ d[p[i]] += a;
+
+ // note the symmetry
+ l[i] = -a;
+ u[i] = -a;
}
-
}
} // namespace fvm
diff --git a/src/math.hpp b/src/math.hpp
index ec737fa2e9dd3c2b0ce49f6708da3e235c1172ba..86eccfce67e7f53fe6645b15808c201c88634e19 100644
--- a/src/math.hpp
+++ b/src/math.hpp
@@ -1,5 +1,6 @@
#pragma once
+#include <utility>
#include <cmath>
namespace nest {
@@ -17,6 +18,12 @@ namespace math {
return (a+b) / T(2);
}
+ template <typename T>
+ T constexpr mean(std::pair<T,T> const& p)
+ {
+ return (p.first+p.second) / T(2);
+ }
+
template <typename T>
T constexpr square(T a)
{
diff --git a/src/parameter_list.cpp b/src/parameter_list.cpp
index a6d401b588c6eb8f543dd5f005201410cbaab827..65155100610f20e1d048d80ab732f868236b5690 100644
--- a/src/parameter_list.cpp
+++ b/src/parameter_list.cpp
@@ -45,7 +45,13 @@ namespace mc {
parameter& parameter_list::get(std::string const& n)
{
- return *find_by_name(n);
+ auto it = find_by_name(n);
+ if(it==parameters_.end()) {
+ throw std::domain_error(
+ "parameter list does not contain parameter"
+ );
+ }
+ return *it;
}
std::string const& parameter_list::name() const {
diff --git a/src/parameter_list.hpp b/src/parameter_list.hpp
index 114a18d07198aaa2b8230206e517ff4817bda445..74b9b1b3fd396c50f27123c016af1a5ef630810f 100644
--- a/src/parameter_list.hpp
+++ b/src/parameter_list.hpp
@@ -136,6 +136,35 @@ namespace mc {
std::ostream& operator<<(std::ostream& o, parameter_list const& l);
+ ///////////////////////////////////////////////////////////////////////////
+ // predefined parameter sets
+ ///////////////////////////////////////////////////////////////////////////
+
+ /// default set of parameters for the cell membrane that are added to every
+ /// segment when it is created.
+ class membrane_parameters
+ : public parameter_list
+ {
+ public:
+
+ using base = parameter_list;
+
+ using base::value_type;
+
+ using base::set;
+ using base::get;
+ using base::parameters;
+ using base::has_parameter;
+
+ membrane_parameters()
+ : base("membrane")
+ {
+ 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
+ }
+ };
+
+ /// parameters for the classic Hodgkin & Huxley model (1952)
class hh_parameters
: public parameter_list
{
diff --git a/src/segment.hpp b/src/segment.hpp
index a73dd9219cba674f90fb03d4e11ea44e0ab29dfd..e339e3ac4bc974477a624fe698d11bf8eb521f26 100644
--- a/src/segment.hpp
+++ b/src/segment.hpp
@@ -171,13 +171,13 @@ class soma_segment : public segment
: radius_{r}
{
kind_ = segmentKind::soma;
+ mechanisms_.push_back(membrane_parameters());
}
soma_segment(value_type r, point_type const &c)
: soma_segment(r)
{
center_ = c;
- kind_ = segmentKind::soma;
}
value_type volume() const override
@@ -245,6 +245,9 @@ class cable_segment : public segment
radii_ = std::move(r);
lengths_ = std::move(lens);
+
+ // add default membrane parameters
+ mechanisms_.push_back(membrane_parameters());
}
cable_segment(
@@ -253,7 +256,7 @@ class cable_segment : public segment
value_type r2,
value_type len
)
- : cable_segment{k, {r1, r2}, {len}}
+ : cable_segment{k, std::vector<value_type>{r1, r2}, std::vector<value_type>{len}}
{ }
// constructor that lets the user describe the cable as a
@@ -269,6 +272,9 @@ class cable_segment : public segment
radii_ = std::move(r);
locations_ = std::move(p);
update_lengths();
+
+ // add default membrane parameters
+ mechanisms_.push_back(membrane_parameters());
}
// helper that lets user specify with the radius and location
diff --git a/src/tree.hpp b/src/tree.hpp
index 8a296cece8c44f750b2a6798780449231e7cf878..d517f7e4570fc9278bab22ff974a3baaa964ea50 100644
--- a/src/tree.hpp
+++ b/src/tree.hpp
@@ -372,7 +372,7 @@ std::vector<int> make_parent_index(tree const& t, C const& counts)
auto num_compartments = index.back();
std::vector<int> parent_index(num_compartments);
auto pos = 0;
- for(auto i : range(0, t.num_nodes())) {
+ for(int i : range(0, t.num_nodes())) {
// get the parent of this segment
// taking care for the case where the root node has -1 as its parent
auto parent = t.parent(i);
diff --git a/tests/test_matrix.cpp b/tests/test_matrix.cpp
index ba6794f1bd4ab095cc542bee340bb8a4e5ff15cc..19e6702ac975ea09b01e72d2b9677ebc30368d21 100644
--- a/tests/test_matrix.cpp
+++ b/tests/test_matrix.cpp
@@ -17,7 +17,7 @@ TEST(matrix, construct_from_parent_only)
matrix_type m{p};
EXPECT_EQ(m.num_cells(), 1);
EXPECT_EQ(m.size(), 3);
- EXPECT_EQ(p.size(), 3);
+ EXPECT_EQ(p.size(), 3u);
auto mp = m.p();
EXPECT_EQ(mp[0], 0);
@@ -32,7 +32,7 @@ TEST(matrix, construct_from_parent_only)
matrix_type m{std::move(p)};
EXPECT_EQ(m.num_cells(), 1);
EXPECT_EQ(m.size(), 3);
- EXPECT_EQ(p.size(), 3);
+ EXPECT_EQ(p.size(), 3u);
EXPECT_EQ(m.size(), 3);
auto mp = m.p();
@@ -49,7 +49,7 @@ TEST(matrix, construct_from_parent_only)
matrix_type m{p};
EXPECT_EQ(m.num_cells(), 1);
EXPECT_EQ(m.size(), 3);
- EXPECT_EQ(p.size(), 3);
+ EXPECT_EQ(p.size(), 3u);
auto mp = m.p();
EXPECT_EQ(mp[0], 0);
@@ -64,7 +64,7 @@ TEST(matrix, construct_from_parent_only)
matrix_type m{std::move(p)};
EXPECT_EQ(m.num_cells(), 1);
EXPECT_EQ(m.size(), 3);
- EXPECT_EQ(p.size(), 0); // 0 implies moved from
+ EXPECT_EQ(p.size(), 0u); // 0 implies moved from
auto mp = m.p();
EXPECT_EQ(mp[0], 0);
@@ -98,7 +98,7 @@ TEST(matrix, solve)
std::iota(p.begin()+1, p.end(), 0);
matrix_type m{p};
- EXPECT_EQ(m.size(), n);
+ EXPECT_EQ(m.size(), (int)n);
EXPECT_EQ(m.num_cells(), 1);
m.d()(memory::all) = 2;
diff --git a/tests/test_run.cpp b/tests/test_run.cpp
index c489429926c39d8ff424f46bb3d6e24dc2a4b591..6728cc20889792fe9e67883a4ab64777234bcd09 100644
--- a/tests/test_run.cpp
+++ b/tests/test_run.cpp
@@ -39,7 +39,6 @@ TEST(run, init)
using fvm_cell = fvm::fvm_cell<double, int>;
fvm_cell fvcell(cell);
-
}
// test out the parameter infrastructure
@@ -54,16 +53,16 @@ TEST(run, parameters)
// add_parameter() returns a bool that indicates whether
// it was able to successfull add the parameter
EXPECT_EQ(list.add_parameter(std::move(p)), true);
- EXPECT_EQ(list.num_parameters(), 1u);
+ EXPECT_EQ(list.num_parameters(), 1);
// test in place construction of a parameter
EXPECT_EQ(list.add_parameter({"b", -3.0}), true);
- EXPECT_EQ(list.num_parameters(), 2u);
+ EXPECT_EQ(list.num_parameters(), 2);
// check that adding a parameter that already exists returns false
// and does not increase the number of parameters
EXPECT_EQ(list.add_parameter({"b", -3.0}), false);
- EXPECT_EQ(list.num_parameters(), 2u);
+ EXPECT_EQ(list.num_parameters(), 2);
auto &parms = list.parameters();
EXPECT_EQ(parms[0].name, "a");
diff --git a/tests/test_tree.cpp b/tests/test_tree.cpp
index 097895dd7cb962ba8cdd3db1bcea09b461285112..9bb320d1e4cf5f1b0e8b1f5e320384f32935d728 100644
--- a/tests/test_tree.cpp
+++ b/tests/test_tree.cpp
@@ -267,10 +267,10 @@ TEST(tree, make_parent_index)
std::vector<int> counts = {5};
nest::mc::tree t(parent_index);
auto new_parent_index = make_parent_index(t, counts);
- EXPECT_EQ(new_parent_index.size(), counts[0]);
+ EXPECT_EQ(new_parent_index.size(), (unsigned)counts[0]);
EXPECT_EQ(new_parent_index[0], 0);
- for(auto i=1; i<new_parent_index.size(); ++i) {
- EXPECT_EQ(new_parent_index[i], i-1);
+ for(auto i=1u; i<new_parent_index.size(); ++i) {
+ EXPECT_EQ((unsigned)new_parent_index[i], i-1);
}
}
// some trees with single compartment per segment
@@ -281,13 +281,13 @@ TEST(tree, make_parent_index)
// 1
std::vector<int>{0,0},
// 0
- // / \
+ // / \.
// 1 2
std::vector<int>{0,0,0},
// 0
- // / \
+ // / \.
// 1 4
- // / \ |\
+ // / \ |\.
// 2 3 5 6
std::vector<int>{0,0,0,1,1,2,2}
};
@@ -301,15 +301,15 @@ TEST(tree, make_parent_index)
// a tree with multiple compartments per segment
//
// 0
- // / \
+ // / \.
// 1 8
- // / \
+ // / \.
// 2 9
- // /
+ // /.
// 3
- // / \
+ // / \.
// 4 6
- // / \
+ // / \.
// 5 7
{
std::vector<int> parent_index = {0,0,1,2,3,4,3,6,0,8};