diff --git a/.ycm_extra_conf.py b/.ycm_extra_conf.py
index 09a57a749c8438dc4200eeea34b9cb9bc1845af4..b8378fb82ab46999a5a31f6f5ff516f0c6cc5af6 100644
--- a/.ycm_extra_conf.py
+++ b/.ycm_extra_conf.py
@@ -40,7 +40,9 @@ flags = [
'-x',
'c++',
'-I',
- './src',
+ 'src',
+ '-I',
+ '.',
# '-isystem',
# '/Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX10.10.sdk/usr/include/c++/4.2.1',
# '-I',
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/fvm.hpp b/src/fvm.hpp
index db19aaae308840f45445a09be20500bf0138b346..9e84119f526c531f3e80c082ea22d8dc7edd997b 100644
--- a/src/fvm.hpp
+++ b/src/fvm.hpp
@@ -40,11 +40,28 @@ class fvm_cell {
fvm_cell(nest::mc::cell const& cell);
/// build the matrix for a given time step
- void setup_matrx(value_type dt);
-
- matrix_type& matrix()
+ 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_;
+ return matrix_.size();
}
private:
@@ -57,8 +74,21 @@ class fvm_cell {
/// 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_;
+ /// 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_;
+
};
////////////////////////////////////////////////////////////
@@ -67,9 +97,12 @@ 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())
+: 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)}
{
using util::left;
using util::right;
@@ -77,14 +110,6 @@ fvm_cell<T, I>::fvm_cell(nest::mc::cell const& cell)
auto parent_index = matrix_.p();
auto const& segment_index = cell.segment_index();
- // 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()) {
@@ -94,7 +119,9 @@ fvm_cell<T, I>::fvm_cell(nest::mc::cell const& cell)
"FVM lowering encountered soma with non-zero index"
);
}
- cv_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()) {
// loop over each compartment in the cable
@@ -111,13 +138,18 @@ fvm_cell<T, I>::fvm_cell(nest::mc::cell const& cell)
// (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 i = segment_index[seg_idx] + c.index;
auto j = parent_index[i];
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);
+ 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;
@@ -125,6 +157,8 @@ fvm_cell<T, I>::fvm_cell(nest::mc::cell const& cell)
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 {
@@ -132,10 +166,36 @@ fvm_cell<T, I>::fvm_cell(nest::mc::cell const& cell)
}
++seg_idx;
}
+
+ // 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>
-void fvm_cell<T, I>::setup_matrx(T dt)
+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;
@@ -144,6 +204,7 @@ void fvm_cell<T, I>::setup_matrx(T dt)
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
@@ -157,13 +218,9 @@ void fvm_cell<T, I>::setup_matrx(T dt)
// . . .
// l[i] . . d[i]
//
-
- //d(all) = cv_areas_ + dt*(alpha_ + alpha_(p));
- //d[0] = cv_areas_[0];
-
d(all) = cv_areas_;
- for(auto i=1; i<d.size(); ++i) {
- auto a = dt * alpha_[i];
+ for(auto i=1u; i<d.size(); ++i) {
+ auto a = dt * face_alpha_[i];
d[i] += a;
l[i] = -a;
@@ -172,6 +229,12 @@ void fvm_cell<T, I>::setup_matrx(T dt)
// 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/matrix.hpp b/src/matrix.hpp
index 4fbd3346e274c1c712bc99d1bd4ea5c2d583f582..3e5f4442373c24821bc9f6b8499db90a37a87835 100644
--- a/src/matrix.hpp
+++ b/src/matrix.hpp
@@ -57,7 +57,7 @@ class matrix {
}
/// the dimension of the matrix (i.e. the number of rows or colums)
- size_type size() const
+ std::size_t size() const
{
return parent_index_.size();
}
diff --git a/src/parameter_list.cpp b/src/parameter_list.cpp
index 65155100610f20e1d048d80ab732f868236b5690..86c641aa809e6efe1c8f03247cc64f9f1b3a01f0 100644
--- a/src/parameter_list.cpp
+++ b/src/parameter_list.cpp
@@ -85,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
@@ -97,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";
@@ -106,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 74b9b1b3fd396c50f27123c016af1a5ef630810f..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,6 @@ namespace mc {
};
- std::ostream& operator<<(std::ostream& o, parameter_list const& l);
-
///////////////////////////////////////////////////////////////////////////
// predefined parameter sets
///////////////////////////////////////////////////////////////////////////
@@ -159,8 +146,8 @@ namespace mc {
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
+ 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
}
};
@@ -192,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/tests/test_matrix.cpp b/tests/test_matrix.cpp
index 19e6702ac975ea09b01e72d2b9677ebc30368d21..f0271bd40b6ee109489dfd1f1edb781dca49cfa3 100644
--- a/tests/test_matrix.cpp
+++ b/tests/test_matrix.cpp
@@ -16,7 +16,7 @@ TEST(matrix, construct_from_parent_only)
std::vector<int> p = {0,0,1};
matrix_type m{p};
EXPECT_EQ(m.num_cells(), 1);
- EXPECT_EQ(m.size(), 3);
+ EXPECT_EQ(m.size(), 3u);
EXPECT_EQ(p.size(), 3u);
auto mp = m.p();
@@ -31,9 +31,8 @@ TEST(matrix, construct_from_parent_only)
std::vector<int> p = {0,0,1};
matrix_type m{std::move(p)};
EXPECT_EQ(m.num_cells(), 1);
- EXPECT_EQ(m.size(), 3);
+ EXPECT_EQ(m.size(), 3u);
EXPECT_EQ(p.size(), 3u);
- EXPECT_EQ(m.size(), 3);
auto mp = m.p();
EXPECT_EQ(mp[0], 0);
@@ -48,7 +47,7 @@ TEST(matrix, construct_from_parent_only)
std::iota(p.begin()+1, p.end(), 0);
matrix_type m{p};
EXPECT_EQ(m.num_cells(), 1);
- EXPECT_EQ(m.size(), 3);
+ EXPECT_EQ(m.size(), 3u);
EXPECT_EQ(p.size(), 3u);
auto mp = m.p();
@@ -63,7 +62,7 @@ TEST(matrix, construct_from_parent_only)
std::iota(p.begin()+1, p.end(), 0);
matrix_type m{std::move(p)};
EXPECT_EQ(m.num_cells(), 1);
- EXPECT_EQ(m.size(), 3);
+ EXPECT_EQ(m.size(), 3u);
EXPECT_EQ(p.size(), 0u); // 0 implies moved from
auto mp = m.p();
@@ -98,7 +97,7 @@ TEST(matrix, solve)
std::iota(p.begin()+1, p.end(), 0);
matrix_type m{p};
- EXPECT_EQ(m.size(), (int)n);
+ EXPECT_EQ(m.size(), 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 73e7929297f53c2972e52dd3116aecf0affa51ad..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;
@@ -38,11 +88,12 @@ TEST(run, init)
using fvm_cell = fvm::fvm_cell<double, int>;
fvm_cell fvcell(cell);
- EXPECT_EQ(fvcell.matrix().size(), 11);
+ auto& J = fvcell.jacobian();
+ EXPECT_EQ(J.size(), 11u);
- fvcell.setup_matrx(0.01);
+ fvcell.setup_matrix(0.01);
+ std::cout << "areas " << fvcell.cv_areas() << "\n";
- auto& J = fvcell.matrix();
//std::cout << "l" << J.l() << "\n";
//std::cout << "d" << J.d() << "\n";
//std::cout << "u" << J.u() << "\n";