From 5d2a79af25a8de80b8da41e01db733d511073067 Mon Sep 17 00:00:00 2001
From: bcumming <bcumming@cscs.ch>
Date: Thu, 14 Apr 2016 15:36:17 +0200
Subject: [PATCH] matrix construction for FVM method
---
src/algorithms.hpp | 21 +++++--
src/fvm.hpp | 121 ++++++++++++++++++++++++++++++++++-------
src/math.hpp | 7 +++
src/parameter_list.cpp | 8 ++-
src/parameter_list.hpp | 29 ++++++++++
src/segment.hpp | 10 +++-
src/tree.hpp | 2 +-
tests/test_matrix.cpp | 10 ++--
tests/test_run.cpp | 7 +--
tests/test_tree.cpp | 22 ++++----
10 files changed, 187 insertions(+), 50 deletions(-)
diff --git a/src/algorithms.hpp b/src/algorithms.hpp
index 250cade2..10ba44d3 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 bd1339c0..0b4cb1d0 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 ec737fa2..86eccfce 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 a6d401b5..65155100 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 114a18d0..74b9b1b3 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 a73dd921..e339e3ac 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 8a296cec..d517f7e4 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 ba6794f1..19e6702a 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 c4894299..6728cc20 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 097895dd..9bb320d1 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};
--
GitLab