diff --git a/.ycm_extra_conf.py b/.ycm_extra_conf.py
index 156103100b8885036931a2cbb148adc8371a05a3..b8378fb82ab46999a5a31f6f5ff516f0c6cc5af6 100644
--- a/.ycm_extra_conf.py
+++ b/.ycm_extra_conf.py
@@ -40,13 +40,11 @@ flags = [
'-x',
'c++',
'-I',
- '.',
- '-isystem',
- '/Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX10.10.sdk/usr/include/c++/4.2.1',
- '-isystem',
- '/Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX10.10.sdk/usr/include/c++/4.2.1/bits',
- '-I',
'src',
+ '-I',
+ '.',
+# '-isystem',
+# '/Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX10.10.sdk/usr/include/c++/4.2.1',
# '-I',
# '/usr/include/c++/4.9.2',
# '-isystem',
diff --git a/docs/formulation.tex b/docs/formulation.tex
index 88e8e0a55f9a5bff24865cdc316c8f49ae814d78..0676ff4f9b03d1e1b3df4aa16963c858a4bb14e4 100644
--- a/docs/formulation.tex
+++ b/docs/formulation.tex
@@ -193,7 +193,7 @@ where $a_{i,\ell}$ and $a_{i,r}$ are the radii of at the left and right end of t
%-------------------------------------------------------------------------------
\subsubsection{Putting it all together}
%-------------------------------------------------------------------------------
-By substituting the volume averaging of the temporal derivative in~\eq{eq:dvdt} approximations for the flux over the surfaces in~\eq{eq:q_ij} and~\eq{eq:cv_volume} respectively into the conservation equation~\eq{eq:cable_balance} we get the following ODE defined for each node in the cell
+By substituting the volume averaging of the temporal derivative in~\eq{eq:dvdt} approximations for the flux over the surfaces in~\eq{eq:J_ij} and~\eq{eq:cv_volume} respectively into the conservation equation~\eq{eq:cable_balance} we get the following ODE defined for each node in the cell
\begin{equation}
\sigma_i c_m \dder{V_i}{t}
= -\sum_{j\in\mathcal{N}_i} {\frac{\sigma_{i,j}}{r_L \Delta x_{i,j}} (V_i-V_j)} - \sigma_i\cdot(i_m(V_i) - i_e(x_i)),
@@ -228,18 +228,25 @@ The current $i_m$ is often a nonlinear function of voltage, so if it was formula
The equations can be rearranged to have all unknown voltage values on the lhs, and values that can be calculated directly on the rhs:
\begin{align}
- & V_i^{k+1} + \sum_{j\in\mathcal{N}_i} {\frac{\alpha_{ij}}{\sigma_i} (V_i^{k+1}-V_j^{k+1})}
+ & \sigma_i V_i^{k+1} + \sum_{j\in\mathcal{N}_i} {\Delta t \alpha_{ij} (V_i^{k+1}-V_j^{k+1})}
\nonumber \\
- = & V_i^k - \frac{\Delta t}{c_m}(i_m^{k} - i_e),
+ = & \sigma_i \left( V_i^k - \frac{\Delta t}{c_m}(i_m^{k} - i_e) \right),
\label{eq:ode_linsys}
\end{align}
where the value
\begin{equation}
- \alpha_{ij} = \alpha_{ji} = \frac{\Delta t \sigma_{ij}}{ c_m r_L \Delta x_{ij}}
+ \alpha_{ij} = \alpha_{ji} = \frac{\sigma_{ij}}{ c_m r_L \Delta x_{ij}}
\label{eq:alpha_linsys}
\end{equation}
is a constant that can be computed for each interface between adjacent compartments during set up.
+The left hand side of \eq{eq:ode_linsys} can be rearranged
+\begin{equation}
+ \left[ \sigma_i + \sum_{j\in\mathcal{N}_i} {\Delta t \alpha_{ij}} \right] V_i^{k+1}
+ - \sum_{j\in\mathcal{N}_i} { \Delta t \alpha_{ij} V_j^{k+1}},
+\end{equation}
+which gives the coefficients for the linear system.
+
%-------------------------------------------------------------------------------
\subsubsection{Example: unbranched uniform cable}
%-------------------------------------------------------------------------------
@@ -247,18 +254,22 @@ For an unrbanched uniform cable of constant radius $a$, with length $L$ and $n$
\begin{align}
\Delta x_{ij} &= \Delta x = \frac{L}{n-1}, \nonumber \\
\sigma_{ij} &= \pi a^2, \nonumber \\
- \sigma_{i} &= 2 \pi a \Delta x, \nonumber \\
- \alpha_{ij} &= \frac{\pi a^2\Delta t}{c_m r_L\Delta x}, \nonumber \\
- \frac{\alpha_{ij}}{\sigma_i}
- &= \frac{a\Delta t}{2c_m r_L\Delta x^2}. \nonumber
+ \sigma_{i} &= \sigma_s = 2 \pi a \Delta x, \nonumber \\
+ \alpha_{ij} &= \alpha = \frac{\pi a^2}{c_m r_L\Delta x}, \nonumber
\end{align}
-With these simplifications, the lhs of the linear system is
+With these simplifications, the LHS of the linear system is
\begin{align}
- & V_i^{k+1} + \beta (V_i^{k+1}-V_{i+1}^{k+1}) + \beta (V_i^{k+1}-V_{i-1}^{k+1})
- \nonumber \\
- = & (1+2\beta)V_i^{k+1} - \beta V_{i+1}^{k+1} - \beta V_{i-1}^{k+1}.
+ \left[\sigma_s + 2\Delta t\alpha\right] & V_{i}^{k+1}
+ - \Delta t \alpha V_{i+1}^{k+1}
+ - \Delta t \alpha V_{i-1}^{k+1}
+ \nonumber \\
+ \left[\sigma_s/2 + \Delta t\alpha\right] & V_{1}^{k+1}
+ - \Delta t \alpha V_{2}^{k+1}
+ \nonumber \\
+ \left[\sigma_s/2 + \Delta t\alpha\right] & V_{n}^{k+1}
+ - \Delta t \alpha V_{n-1}^{k+1}
+ \nonumber
\end{align}
-where $\beta=\frac{a\Delta t}{2c_m r_L\Delta x^2}$.
The end points of the cable, i.e. the compartments for $x_1$ and $x_n$, have to be handled differently.
If we assume that a no-flux boundary condition, i.e. $\vv{J}\cdot\vv{n}=0$, is imposed at the end of the cable, the lhs of the linear system are
diff --git a/src/algorithms.hpp b/src/algorithms.hpp
index 2521091545deefa26c7c7e3acf9abebfd23eef02..797f471bc51afac58f650d728ea57cfec445e4ab 100644
--- a/src/algorithms.hpp
+++ b/src/algorithms.hpp
@@ -36,6 +36,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);
@@ -70,15 +75,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;
}
}
@@ -92,8 +103,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 e72bd5e82cf0098c22fc5d92d20c61d1a76e0ba9..9e84119f526c531f3e80c082ea22d8dc7edd997b 100644
--- a/src/fvm.hpp
+++ b/src/fvm.hpp
@@ -2,14 +2,14 @@
#include <algorithm>
-#include "cell.hpp"
-#include "segment.hpp"
+#include <cell.hpp>
+#include <segment.hpp>
+#include <math.hpp>
+#include <matrix.hpp>
+#include <util.hpp>
+#include <algorithms.hpp>
-#include "math.hpp"
-#include "util.hpp"
-#include "algorithms.hpp"
-
-#include "../vector/include/Vector.hpp"
+#include <vector/include/Vector.hpp>
namespace nest {
namespace mc {
@@ -24,16 +24,71 @@ class fvm_cell {
/// the integral index type
using size_type = I;
+ using matrix_type = matrix<value_type, size_type>;
+
/// the container used for indexes
using index_type = memory::HostVector<size_type>;
/// 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_matrix(value_type dt);
+
+ /// TODO this should be const
+ /// which requires const_view in the vector library
+ matrix_type& jacobian();
+
+ /// TODO this should be const
+ /// return list of CV areas in :
+ /// um^2
+ /// 1e-6.mm^2
+ /// 1e-8.cm^2
+ vector_view cv_areas();
+
+ /// TODO this should be const
+ /// return the capacitance of each CV surface
+ /// note that this is the total capacitance, not per unit area
+ /// which is equivalent to sigma_i * c_m
+ vector_view cv_capacitance();
+
+ std::size_t size() const
+ {
+ return matrix_.size();
+ }
+
private:
- index_type parent_index_;
+ /// the linear system for implicit time stepping of cell state
+ matrix_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
+ /// face_alpha_[i] = area_face / (c_m * r_L * delta_x);
+ vector_type face_alpha_;
+
+ /// cv_capacitance_[i] is the capacitance of CV i per unit area (i.e. c_m)
+ vector_type cv_capacitance_;
+
+ /// the average current over the surface of each CV
+ /// current_ = i_m - i_e
+ /// so the total current over the surface of CV i is
+ /// current_[i] * cv_areas_
+ vector_type current_;
+
+ /// the potential in mV in each CV
+ vector_type voltage_;
+
};
////////////////////////////////////////////////////////////
@@ -42,55 +97,144 @@ 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_ {size(), T(0)}
+, face_alpha_ {size(), T(0)}
+, cv_capacitance_{size(), T(0)}
+, current_ {size(), T(0)}
+, voltage_ {size(), T(0)}
{
- parent_index_ = cell.parent_index();
- auto const& segment_index = cell.segment_index();
- auto num_fv = segment_index.back();
+ using util::left;
+ using util::right;
- // 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);
+ auto parent_index = matrix_.p();
+ auto const& segment_index = cell.segment_index();
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());
+ auto area = math::area_sphere(soma->radius());
+ cv_areas_[0] += area;
+ cv_capacitance_[0] += area * soma->mechanism("membrane").get("c_m").value;
}
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 .)
+ auto c_m = cable->mechanism("membrane").get("c_m").value;
+ auto r_L = cable->mechanism("membrane").get("r_L").value;
for(auto c : cable->compartments()) {
- auto rhs = segment_index[seg_idx] + c.index;
- auto lhs = parent_index_[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 );
+ face_alpha_[i] = area_face / (c_m * r_L * c.length);
+ cv_capacitance_[i] = c_m;
- fv_volumes[lhs] += math::volume_frustrum(len, left(c.radius), rad_C);
- fv_areas[lhs] += math::area_frustrum (len, left(c.radius), rad_C);
+ std::cout << "radius " << radius_center << ", c_m " << c_m << ", r_L " << r_L << ", dx " << c.length << "\n";
- 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 halflen = c.length/2;
+
+ 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;
+ cv_capacitance_[j] += al * c_m;
+ cv_capacitance_[i] += ar * c_m;
}
}
else {
- throw std::domain_error(
- "FVM lowering encountered unsuported segment type"
- );
+ throw std::domain_error("FVM lowering encountered unsuported segment type");
}
++seg_idx;
}
- //std::cout << "areas " << algorithms::sum(fv_volumes) << ", " << cell.volume() << "\n";
- //std::cout << "volumes " << algorithms::sum(fv_areas) << ", " << cell.area() << "\n";
+ // normalize the capacitance by cv_area
+ for(auto i=0u; i<size(); ++i) {
+ cv_capacitance_[i] /= cv_areas_[i];
+ }
+}
+
+template <typename T, typename I>
+typename fvm_cell<T,I>::matrix_type&
+fvm_cell<T,I>::jacobian()
+{
+ return matrix_;
+}
+
+template <typename T, typename I>
+typename fvm_cell<T,I>::vector_view
+fvm_cell<T,I>::cv_areas()
+{
+ return cv_areas_;
+}
+
+template <typename T, typename I>
+typename fvm_cell<T,I>::vector_view
+fvm_cell<T,I>::cv_capacitance()
+{
+ return cv_capacitance_;
+}
+
+template <typename T, typename I>
+void fvm_cell<T, I>::setup_matrix(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();
+ auto rhs = matrix_.rhs();
+
+ // 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) = cv_areas_;
+ for(auto i=1u; i<d.size(); ++i) {
+ auto a = dt * face_alpha_[i];
+
+ d[i] += a;
+ l[i] = -a;
+ u[i] = -a;
+
+ // add contribution to the diagonal of parent
+ d[p[i]] += a;
+ }
+
+ // the RHS of the linear system is
+ // sigma_i * (V[i] - dt/cm*(im - ie))
+ for(auto i=0u; i<d.size(); ++i) {
+ rhs[i] = cv_areas_[i] * (voltage_[i] - dt/cv_capacitance_[i]*current_[i]);
+ }
}
} // 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/matrix.hpp b/src/matrix.hpp
index 1b1daad2fc17af2d8a3ad87ced5fc03435deca96..3e5f4442373c24821bc9f6b8499db90a37a87835 100644
--- a/src/matrix.hpp
+++ b/src/matrix.hpp
@@ -1,6 +1,10 @@
#pragma once
-#include "../vector/include/Vector.hpp"
+#include <type_traits>
+
+#include "util.hpp"
+
+#include <vector/include/Vector.hpp>
namespace nest {
namespace mc {
@@ -20,67 +24,146 @@ class matrix {
using index_type = memory::HostVector<size_type>;
using index_view = typename index_type::view_type;
+ /// construct matrix for one or more cells, with combined parent index
+ /// and a cell index
+ template <
+ typename LHS, typename RHS,
+ typename = typename
+ std::enable_if<
+ util::is_container<LHS>::value &&
+ util::is_container<RHS>::value
+ >
+ >
+ matrix(LHS&& pi, RHS&& ci)
+ : parent_index_(std::forward<LHS>(pi))
+ , cell_index_(std::forward<RHS>(ci))
+ {
+ setup();
+ }
+
+ /// construct matrix for a single cell described by a parent index
+ template <
+ typename IDX,
+ typename = typename
+ std::enable_if< util::is_container<IDX>::value >
+ >
+ matrix(IDX&& pi)
+ : parent_index_(std::forward<IDX>(pi))
+ , cell_index_(2)
+ {
+ cell_index_[0] = 0;
+ cell_index_[1] = size();
+ setup();
+ }
/// the dimension of the matrix (i.e. the number of rows or colums)
- size_type size() const
+ std::size_t size() const
{
- return parent_indices_.size();
+ return parent_index_.size();
}
/// the total memory used to store the matrix
std::size_t memory() const
{
- auto s = sizeof(value_type) * data_.size();
- s += sizeof(size_type) * parent_indices_.size();
- return s + sizeof(matrix);
+ auto s = 5 * (sizeof(value_type) * size() + sizeof(vector_type));
+ s += sizeof(size_type) * (parent_index_.size() + cell_index_.size())
+ + 2*sizeof(index_type);
+ s += sizeof(matrix);
+ return s;
}
/// the number of cell matrices that have been packed together
size_type num_cells() const
{
- return cell_index_.size();
+ return cell_index_.size() - 1;
}
- private:
+ /// the vector holding the lower part of the matrix
+ view_type l()
+ {
+ return l_;
+ }
+
+ /// the vector holding the diagonal of the matrix
+ view_type d()
+ {
+ return d_;
+ }
+
+ /// the vector holding the upper part of the matrix
+ view_type u()
+ {
+ return u_;
+ }
+
+ /// the vector holding the right hand side of the linear equation system
+ view_type rhs()
+ {
+ return rhs_;
+ }
+
+ /// the vector holding the parent index
+ index_view p()
+ {
+ return parent_index_;
+ }
- /// calculate the number of values per sub-array plus padding
- size_type reservation() const
+ /// solve the linear system
+ /// upon completion the solution is stored in the RHS storage
+ /// and can be accessed via rhs()
+ void solve()
{
- constexpr auto alignment = vector_type::coordinator_type::alignment();
- auto const padding = memory::impl::get_padding<value_type>(alignment, size());
- return size() + padding;
+ index_view const& p = parent_index_;
+ auto const ncells = num_cells();
+
+ // loop over submatrices
+ for(auto m=0; m<ncells; ++m) {
+ auto first = cell_index_[m];
+ auto last = cell_index_[m+1];
+
+ // backward sweep
+ for(auto i=last-1; i>first; --i) {
+ auto factor = l_[i] / d_[i];
+ d_[p[i]] -= factor * u_[i];
+ rhs_[p[i]] -= factor * rhs_[i];
+ }
+ rhs_[first] /= d_[first];
+
+ // forward sweep
+ for(auto i=first+1; i<last; ++i) {
+ rhs_[i] -= u_[i] * rhs_[p[i]];
+ rhs_[i] /= d_[i];
+ }
+ }
}
- /// subdivide the memory in data_ into the sub-arrays that store the matrix
- /// diagonals and rhs/solution vectors
+ private:
+
+ /// allocate memory for storing matrix and right hand side vector
void setup()
{
const auto n = size();
- const auto n_alloc = reservation();
-
- // point sub-vectors into storage
- a_ = data_( 0, n);
- d_ = data_( n_alloc, n_alloc + n);
- b_ = data_(2*n_alloc, 2*n_alloc + n);
- rhs_ = data_(3*n_alloc, 3*n_alloc + n);
- v_ = data_(4*n_alloc, 4*n_alloc + n);
- }
+ constexpr auto default_value
+ = std::numeric_limits<value_type>::quiet_NaN();
- // the parent indices are stored as an index array
- index_type parent_indices_;
- index_type cell_index_;
+ l_ = vector_type(n, default_value);
+ d_ = vector_type(n, default_value);
+ u_ = vector_type(n, default_value);
+ rhs_ = vector_type(n, default_value);
+ }
- size_type num_cells_;
+ /// the parent indice that describe matrix structure
+ index_type parent_index_;
- // all the data fields point into a single block of storage
- vector_type data_;
+ /// indexes that point to the start of each cell in the index
+ index_type cell_index_;
- // the individual data files that point into block storage
- view_type a_;
- view_type d_;
- view_type b_;
- view_type rhs_;
- view_type v_;
+ /// storage for lower, diagonal, upper and rhs
+ vector_type l_;
+ vector_type d_;
+ vector_type u_;
+ /// after calling solve, the solution is stored in rhs_
+ vector_type rhs_;
};
} // namespace nest
diff --git a/src/parameter_list.cpp b/src/parameter_list.cpp
index a6d401b588c6eb8f543dd5f005201410cbaab827..86c641aa809e6efe1c8f03247cc64f9f1b3a01f0 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 {
@@ -79,8 +85,7 @@ namespace mc {
} // namespace mc
} // namespace nest
-/*
-static std::ostream&
+std::ostream&
operator<<(std::ostream& o, nest::mc::parameter const& p)
{
return o
@@ -91,7 +96,7 @@ operator<<(std::ostream& o, nest::mc::parameter const& p)
<< ")";
}
-static std::ostream&
+std::ostream&
operator<<(std::ostream& o, nest::mc::parameter_list const& l)
{
o << "parameters \"" << l.name() << "\" :\n";
@@ -100,4 +105,3 @@ operator<<(std::ostream& o, nest::mc::parameter_list const& l)
}
return o;
}
-*/
diff --git a/src/parameter_list.hpp b/src/parameter_list.hpp
index 114a18d07198aaa2b8230206e517ff4817bda445..e7dd1bd1cb1da945d48d8aa11fca8c07674a030f 100644
--- a/src/parameter_list.hpp
+++ b/src/parameter_list.hpp
@@ -54,15 +54,6 @@ namespace mc {
value_type max;
};
- template <typename T>
- std::ostream& operator<<(std::ostream& o, value_range<T> const& r)
- {
- return
- o << "[ "
- << (r.has_lower_bound() ? std::to_string(r.min) : "-inf") << ", "
- << (r.has_upper_bound() ? std::to_string(r.max) : "inf") << "]";
- }
-
struct parameter {
using value_type = double;
using range_type = value_range<value_type>;
@@ -91,8 +82,6 @@ namespace mc {
range_type range;
};
- std::ostream& operator<<(std::ostream& o, parameter const& p);
-
// Use a dumb container class for now
// might have to use a more sophisticated interface in the future if need be
class parameter_list {
@@ -134,8 +123,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({"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
+ }
+ };
+
+ /// parameters for the classic Hodgkin & Huxley model (1952)
class hh_parameters
: public parameter_list
{
@@ -163,3 +179,22 @@ namespace mc {
} // namespace mc
} // namespace nest
+template <typename T>
+std::ostream& operator<<(std::ostream& o, nest::mc::value_range<T> const& r)
+{
+ o << "[";
+ if(r.has_lower_bound())
+ o << r.min;
+ else
+ o<< "-inf";
+ o << ", ";
+ if(r.has_upper_bound())
+ o << r.max;
+ else
+ o<< "inf";
+ return o << "]";
+}
+
+std::ostream& operator<<(std::ostream& o, nest::mc::parameter const& p);
+std::ostream& operator<<(std::ostream& o, nest::mc::parameter_list const& l);
+
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/src/util.hpp b/src/util.hpp
index cb68c2d1f7aa05c5ad4da7a37a5a1abae11c19dc..d92392799cf98da27ee1ca62788686a308c02c1d 100644
--- a/src/util.hpp
+++ b/src/util.hpp
@@ -65,6 +65,21 @@ namespace util {
return p.second;
}
+ template <typename T>
+ struct is_std_vector : std::false_type {};
+
+ template <typename T, typename C>
+ struct is_std_vector<std::vector<T,C>> : std::true_type {};
+
+ template <typename T>
+ struct is_container :
+ std::conditional<
+ is_std_vector<typename std::decay<T>::type>::value ||
+ memory::is_array<T>::value,
+ std::true_type,
+ std::false_type
+ >::type
+ {};
} // namespace util
} // namespace mc
} // namespace nest
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 15d3ec160e9f9a701564c802ee188ebd8845b48e..b555fc465f3e3292ae2fa48f11dfad17bbe0e6e5 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -13,10 +13,11 @@ set(TEST_SOURCES
gtest-all.cpp
# unit tests
+ test_algorithms.cpp
test_cell.cpp
test_compartments.cpp
+ test_matrix.cpp
test_point.cpp
- test_algorithms.cpp
test_run.cpp
test_segment.cpp
test_stimulus.cpp
diff --git a/tests/test_algorithms.cpp b/tests/test_algorithms.cpp
index 3e59db7e5557701f48ebac0d508efe6559b8a9ba..91e43c8b86717ce536110caa7394b9b8b6d74d0d 100644
--- a/tests/test_algorithms.cpp
+++ b/tests/test_algorithms.cpp
@@ -2,8 +2,8 @@
#include <vector>
-#include "../src/algorithms.hpp"
-#include "../src/util.hpp"
+#include <algorithms.hpp>
+#include <util.hpp>
TEST(algorithms, sum)
{
diff --git a/tests/test_matrix.cpp b/tests/test_matrix.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..f0271bd40b6ee109489dfd1f1edb781dca49cfa3
--- /dev/null
+++ b/tests/test_matrix.cpp
@@ -0,0 +1,122 @@
+#include <numeric>
+
+#include "gtest.h"
+
+#include <matrix.hpp>
+#include <math.hpp>
+#include <vector/include/Vector.hpp>
+
+TEST(matrix, construct_from_parent_only)
+{
+ using matrix_type = nest::mc::matrix<double, int>;
+
+ // pass parent index as a std::vector by reference
+ // the input should not be moved from
+ {
+ std::vector<int> p = {0,0,1};
+ matrix_type m{p};
+ EXPECT_EQ(m.num_cells(), 1);
+ EXPECT_EQ(m.size(), 3u);
+ EXPECT_EQ(p.size(), 3u);
+
+ auto mp = m.p();
+ EXPECT_EQ(mp[0], 0);
+ EXPECT_EQ(mp[1], 0);
+ EXPECT_EQ(mp[2], 1);
+ }
+
+ // pass parent index as a std::vector by rvalue reference
+ // the input should not be moved from
+ {
+ std::vector<int> p = {0,0,1};
+ matrix_type m{std::move(p)};
+ EXPECT_EQ(m.num_cells(), 1);
+ EXPECT_EQ(m.size(), 3u);
+ EXPECT_EQ(p.size(), 3u);
+
+ auto mp = m.p();
+ EXPECT_EQ(mp[0], 0);
+ EXPECT_EQ(mp[1], 0);
+ EXPECT_EQ(mp[2], 1);
+ }
+
+ // pass parent index as a HostVector by reference
+ // expect that the input is not moved from
+ {
+ memory::HostVector<int> p(3, 0);
+ std::iota(p.begin()+1, p.end(), 0);
+ matrix_type m{p};
+ EXPECT_EQ(m.num_cells(), 1);
+ EXPECT_EQ(m.size(), 3u);
+ EXPECT_EQ(p.size(), 3u);
+
+ auto mp = m.p();
+ EXPECT_EQ(mp[0], 0);
+ EXPECT_EQ(mp[1], 0);
+ EXPECT_EQ(mp[2], 1);
+ }
+ // pass parent index as a HostVector by rvalue reference
+ // expect that the input is moved from
+ {
+ memory::HostVector<int> p(3, 0);
+ std::iota(p.begin()+1, p.end(), 0);
+ matrix_type m{std::move(p)};
+ EXPECT_EQ(m.num_cells(), 1);
+ EXPECT_EQ(m.size(), 3u);
+ EXPECT_EQ(p.size(), 0u); // 0 implies moved from
+
+ auto mp = m.p();
+ EXPECT_EQ(mp[0], 0);
+ EXPECT_EQ(mp[1], 0);
+ EXPECT_EQ(mp[2], 1);
+ }
+}
+
+TEST(matrix, solve)
+{
+ using matrix_type = nest::mc::matrix<double, int>;
+
+ // trivial case : 1x1 matrix
+ {
+ matrix_type m{std::vector<int>{0}};
+
+ m.d()(memory::all) = 2;
+ m.l()(memory::all) = -1;
+ m.u()(memory::all) = -1;
+ m.rhs()(memory::all) = 1;
+
+ m.solve();
+
+ EXPECT_EQ(m.rhs()[0], 0.5);
+ }
+ // matrices in the range of 2x2 to 1000x1000
+ {
+ using namespace nest::mc::math;
+ for(auto n : memory::Range(2,1001)) {
+ auto p = std::vector<int>(n);
+ std::iota(p.begin()+1, p.end(), 0);
+ matrix_type m{p};
+
+ EXPECT_EQ(m.size(), n);
+ EXPECT_EQ(m.num_cells(), 1);
+
+ m.d()(memory::all) = 2;
+ m.l()(memory::all) = -1;
+ m.u()(memory::all) = -1;
+ m.rhs()(memory::all) = 1;
+
+ m.solve();
+
+ auto x = m.rhs();
+ auto err = square(std::fabs(2.*x[0] - x[1] - 1.));
+ for(auto i : memory::Range(1,n-1)) {
+ err += square(std::fabs(2.*x[i] - x[i-1] - x[i+1] - 1.));
+ }
+ err += square(std::fabs(2.*x[n-1] - x[n-2] - 1.));
+
+ EXPECT_NEAR(0., std::sqrt(err), 1e-8);
+ }
+
+ }
+}
+
diff --git a/tests/test_run.cpp b/tests/test_run.cpp
index c489429926c39d8ff424f46bb3d6e24dc2a4b591..712fca4346e2a1780a9c1c3dca080381eeb8218d 100644
--- a/tests/test_run.cpp
+++ b/tests/test_run.cpp
@@ -3,6 +3,56 @@
#include "../src/cell.hpp"
#include "../src/fvm.hpp"
+TEST(run, cable)
+{
+ using namespace nest::mc;
+
+ nest::mc::cell cell;
+
+ cell.add_soma(6e-4); // 6um in cm
+
+ // 1um radius and 4mm long, all in cm
+ cell.add_cable(0, segmentKind::dendrite, 1e-4, 1e-4, 4e-1);
+
+ std::cout << cell.segment(1)->area() << " is the area\n";
+ EXPECT_EQ(cell.tree().num_segments(), 2u);
+
+ cell.soma()->add_mechanism(hh_parameters());
+
+ auto& soma_hh = cell.soma()->mechanism("hh");
+
+ soma_hh.set("gnabar", 0.12);
+ soma_hh.set("gkbar", 0.036);
+ soma_hh.set("gl", 0.0003);
+ soma_hh.set("el", -54.387);
+
+ cell.segment(1)->set_compartments(4);
+
+ using fvm_cell = fvm::fvm_cell<double, int>;
+ fvm_cell fvcell(cell);
+ auto& J = fvcell.jacobian();
+ EXPECT_EQ(J.size(), 5u);
+
+ fvcell.setup_matrix(0.02);
+ EXPECT_EQ(fvcell.cv_areas().size(), J.size());
+
+ auto& cable_parms = cell.segment(1)->mechanism("membrane");
+ std::cout << soma_hh << std::endl;
+ std::cout << cable_parms << std::endl;
+
+ std::cout << "l " << J.l() << "\n";
+ std::cout << "d " << J.d() << "\n";
+ std::cout << "u " << J.u() << "\n";
+ std::cout << "p " << J.p() << "\n";
+
+ J.rhs()(memory::all) = 1.;
+ J.rhs()[0] = 10.;
+
+ J.solve();
+
+ //std::cout << "x" << J.rhs() << "\n";
+}
+
TEST(run, init)
{
using namespace nest::mc;
@@ -34,12 +84,26 @@ TEST(run, init)
EXPECT_EQ(cell.soma()->mechanism("hh").get("gl").value, 0.0003);
EXPECT_EQ(cell.soma()->mechanism("hh").get("el").value, -54.3);
-
- cell.segment(1)->set_compartments(200);
+ cell.segment(1)->set_compartments(10);
using fvm_cell = fvm::fvm_cell<double, int>;
fvm_cell fvcell(cell);
+ auto& J = fvcell.jacobian();
+ EXPECT_EQ(J.size(), 11u);
+
+ fvcell.setup_matrix(0.01);
+ std::cout << "areas " << fvcell.cv_areas() << "\n";
+
+ //std::cout << "l" << J.l() << "\n";
+ //std::cout << "d" << J.d() << "\n";
+ //std::cout << "u" << J.u() << "\n";
+
+ J.rhs()(memory::all) = 1.;
+ J.rhs()[0] = 10.;
+
+ J.solve();
+ //std::cout << "x" << J.rhs() << "\n";
}
// test out the parameter infrastructure
@@ -54,16 +118,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};
diff --git a/vector b/vector
index c2fc6318f6d65ae9d7aad60309f3b53d590a00ec..e95bcbf1de6f833b5113a82dcc6d47bb7f41c397 160000
--- a/vector
+++ b/vector
@@ -1 +1 @@
-Subproject commit c2fc6318f6d65ae9d7aad60309f3b53d590a00ec
+Subproject commit e95bcbf1de6f833b5113a82dcc6d47bb7f41c397