diff --git a/cmake/FindUnwind.cmake b/cmake/FindUnwind.cmake
index 8d11fb255b0f43070912002ae61015658fe6a440..35fd9c2a19b8f5d3f7cf5086c4e6a105eb3fda23 100644
--- a/cmake/FindUnwind.cmake
+++ b/cmake/FindUnwind.cmake
@@ -39,6 +39,9 @@ if(NOT UNWIND_FOUND)
set(UNWIND_LIBRARIES ${unwind_library_generic} ${unwind_library_target})
+ include(FindPackageHandleStandardArgs)
+ find_package_handle_standard_args(UNWIND DEFAULT_MSG UNWIND_INCLUDE_DIR UNWIND_LIBRARIES)
+
mark_as_advanced(UNWIND_LIBRARIES UNWIND_INCLUDE_DIR)
unset(unwind_search_dir)
diff --git a/docs/model/formulation.tex b/docs/model/formulation.tex
index 0676ff4f9b03d1e1b3df4aa16963c858a4bb14e4..d3eba7a64c4bdf8ad8bcb0da402f9bde9513da3e 100644
--- a/docs/model/formulation.tex
+++ b/docs/model/formulation.tex
@@ -119,10 +119,29 @@ The finite volume method is a natural choice for the solution of the conservatio
%-------------------------------------------------------------------------------
The integral on the lhs of~\eq{eq:cable_balance} can be approximated by assuming that the average transmembrane potential $V$ in $\Omega_i$ is equal to the potential $V_i$ defined at the centre of the segment:
\begin{equation}
- \int_{\Gamma_i}{c_m \pder{V}{t} } \deriv{v} \approx \sigma_i c_m \pder{V_i}{t},
+ \int_{\Gamma_i}{c_m \pder{V}{t} } \deriv{v} \approx \sigma_i \cmi \pder{V_i}{t},
\label{eq:dvdt}
\end{equation}
-where $\sigma_i$ is the surface area of the membrane potential.
+where $\sigma_i$ is the surface area, and $\cmi$ is the average specific membrane capacitance, respectively of the surface $\Gamma_i$.
+
+Each control volume is composed of \emph{sub control volumes}, which are illustrated as the coloured sub-regions in \fig{fig:segment}.
+\begin{equation*}
+ \Omega_i = \bigcup_{j\in\mathcal{N}_i}{\Omega_i^j}.
+\end{equation*}
+Likewise, the surface $\Gamma_i$ is composed of subsufaces as follows
+\begin{equation*}
+ \Gamma_i = \bigcup_{j\in\mathcal{N}_i}{\Gamma_i^j},
+\end{equation*}
+where $\Gamma_i^j$ is the surface of each of the sub-control volumes in $\Omega_i$.
+Thus, the surface area of the CV as
+\begin{equation*}
+ \sigma_i = \sum_{i\in\mathcal{N}_i}{\sigma_i^j},
+\end{equation*}
+where $\sigma_i^j$ is the area of $\Gamma_i^j$, and the average specific membrane capacitance $\cmi$ is
+\begin{equation*}
+ \cmi = \frac{1}{\sigma_i}\sum_{i\in\mathcal{N}_i}{\sigma_i^j c_m^{i,j}}.
+\end{equation*}
+\todo{This is included as a placeholder, we really need more illustrations to show how CV averages are computed for quantities that vary between sub-control volumes of the same CV.}
%-------------------------------------------------------------------------------
\subsubsection{Intra-cellular flux}
@@ -195,7 +214,7 @@ where $a_{i,\ell}$ and $a_{i,r}$ are the radii of at the left and right end of t
%-------------------------------------------------------------------------------
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}
+ \sigma_i \cmi \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)),
\label{eq:ode}
\end{equation}
@@ -228,56 +247,29 @@ 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}
- & \sigma_i V_i^{k+1} + \sum_{j\in\mathcal{N}_i} {\Delta t \alpha_{ij} (V_i^{k+1}-V_j^{k+1})}
+ & \frac{\sigma_i \cmi}{\Delta t} V_i^{k+1} + \sum_{j\in\mathcal{N}_i} {\alpha_{ij} (V_i^{k+1}-V_j^{k+1})}
\nonumber \\
- = & \sigma_i \left( V_i^k - \frac{\Delta t}{c_m}(i_m^{k} - i_e) \right),
+ = & \frac{\sigma_i \cmi}{\Delta t} V_i^k - \sigma_i(i_m^{k} - i_e),
\label{eq:ode_linsys}
\end{align}
where the value
\begin{equation}
- \alpha_{ij} = \alpha_{ji} = \frac{\sigma_{ij}}{ c_m r_L \Delta x_{ij}}
+ \alpha_{ij} = \alpha_{ji} = \frac{\sigma_{ij}}{ 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}},
+ \left[ \frac{\sigma_i \cmi}{\Delta t} + \sum_{j\in\mathcal{N}_i} {\alpha_{ij}} \right] V_i^{k+1}
+ - \sum_{j\in\mathcal{N}_i} { \alpha_{ij} V_j^{k+1}},
+ \label{eq:rhs_linsys}
\end{equation}
which gives the coefficients for the linear system.
-%-------------------------------------------------------------------------------
-\subsubsection{Example: unbranched uniform cable}
-%-------------------------------------------------------------------------------
-For an unrbanched uniform cable of constant radius $a$, with length $L$ and $n$ compartments, the linear system for internal compartments (i.e. not at the end points of the cable) is simplified by the following observations
-\begin{align}
- \Delta x_{ij} &= \Delta x = \frac{L}{n-1}, \nonumber \\
- \sigma_{ij} &= \pi a^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
-\begin{align}
- \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}
-
-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
-\begin{align}
- (1+2\beta)V_1^{k+1} - 2\beta V_{2}^{k+1}, \quad\quad & \text{left} \nonumber \\
- (1+2\beta)V_n^{k+1} - 2\beta V_{n-1}^{k+1}, \quad\quad & \text{right} \nonumber
-\end{align}
-where we note that the ratio $\alpha_{ij}/\sigma_{i}=2\beta$ because the surface area of the control volumes at the boundary are half those on the interior.
+The capacitance of the cell membrane, $\sigma_i \cmi$, varies between control volumes, while the $\alpha_{i,j}$ term in symmetric (i.e. $\alpha_{i,j}=\alpha_{j,i}$.)
+With this in mind, we can see that when the linear system is written in the form~\eq{eq:ode_linsys}, the matrix is symmetric.
+Furthermore, because $\alpha_{i,j} > 0$, the linear system is diagonally dominant for sufficiently small $\Delta t$.
%-------------------------------------------------------------------------------
\subsubsection{The Soma}
diff --git a/docs/model/images/cable.tex b/docs/model/images/cable.tex
index c11e2b740300923efb8985acf7c6e061efe485a2..c2843898c864c723b40e6d04d7c1f73e76c375a5 100644
--- a/docs/model/images/cable.tex
+++ b/docs/model/images/cable.tex
@@ -29,6 +29,10 @@
%\draw [pil, very thin] (-3.4,0) -- ( 3.5, 0);
%\draw [pil, very thin] (-3.2,-0.2) -- (-3.2, 0.5);
+% left sub CV
+\filldraw[white,fill=green!20] (-2,-0.5) -- (-2,0.5) -- (0,0.65) -- (0,-0.65) -- cycle;
+\filldraw[white,fill=blue!20] (0,-0.65) -- (0,0.65) -- ( 2, 0.8) -- ( 2, -0.8) -- cycle;
+
% left volume
\draw [white!60!black] (-6,-0.5) -- (-6, 0.5);
\draw [white!60!black] (-6,-0.5) -- (-2,-0.5);
@@ -74,6 +78,8 @@
\node [] at ( 2.5, 1.15) {$\Gamma_{{i,i+1}}$};
\node [] at ( -1, 1.15) {$\Gamma_{{e}}$};
\node [] at ( 1.7, -1.2) {$\Omega$};
+\node [] at ( -1, 0) {\textcolor{green!50!black}{$\Omega_{i-1}^j$}};
+\node [] at ( 1, 0) {\textcolor{blue}{$\Omega_{i+1}^j$}};
\end{tikzpicture}
\end{document}
diff --git a/docs/model/report.tex b/docs/model/report.tex
index 0f816e1c25a7aa22b7199df27a0cd3d8216f9dc2..14b82f0344e607dec89783de30cff80a33dc761e 100644
--- a/docs/model/report.tex
+++ b/docs/model/report.tex
@@ -71,6 +71,7 @@
\newcommand{\pder}[2]{\frac{\partial{#1}}{\partial{#2}}}
\newcommand{\dder}[2]{\frac{\deriv{#1}}{\deriv{#2}}}
\newcommand{\vv}[1]{\bm{#1}\xspace}
+\newcommand{\cmi}{c_{m,i}}
\newcommand{\unit}[1]{\left[{#1}\right]}
\newcommand{\txtunit}[1]{$\left[{#1}\right]$}
diff --git a/docs/model/symbols.tex b/docs/model/symbols.tex
index baa371ad8efe42d483d2f35325b6d04ef730ca43..7f2bf3a352fe1527aa8647a9bc148bb53bffaddc 100644
--- a/docs/model/symbols.tex
+++ b/docs/model/symbols.tex
@@ -2,18 +2,22 @@
%\subsubsection{Balancing Units}
%-------------------------------------------------------------------------------
Ensuring that units are balanced and correct requires care.
-Take the description of the nonlinear system of ODEs that arises from the finite volume discretisation
+Take the system of linear equations that arises from the finite volume discretisation, specified in \eq{eq:ode_linsys} and \eq{eq:rhs_linsys}
\begin{equation}
\label{eq:linsys_FV}
- V_i^{k+1} + \sum_{j\in\mathcal{N}_i} {\frac{\Delta t \alpha_{ij}}{\sigma_i} (V_i^{k+1}-V_j^{k+1})}
- = V_i^k - \frac{\Delta t}{c_m}(i_m^{k} - i_e).
+ \left[
+ \frac{\sigma_i \cmi}{\Delta t} + \sum_{j\in\mathcal{N}_i} {\alpha_{ij}}
+ \right]
+ V_i^{k+1} - \sum_{j\in\mathcal{N}_i} { \alpha_{ij} V_j^{k+1}}
+ =
+ \frac{\sigma_i \cmi}{\Delta t} V_i^k - \sigma_i(i_m^{k} - i_e).
\end{equation}
+
The choice of units for a parameter, e.g. $\mu m^2$ or $m^2$ for the area $\sigma_{ij}$, introduces a constant of proportionality wherever it is used ($10^{-12}$ in the case of $\mu m^2 \rightarrow m^2$).
Wherever terms are added in \eq{eq:linsys_FV} the units must be checked, and constants of proportionality balanced.
First, appropriate units for each of the parameters and variables are chosen in~\tbl{tbl:units}.
We try to use the same units as NEURON, except for the specific membrane capacitance $c_m$, for which $F\cdot m^{-2}$ is used in place of $nF\cdot mm^{-2}$.
-In \eq{eq:linsys_FV} we choose units of $mV \equiv 10^{-3}V$ for each term because of the $V_i$ terms on either side of the equation.
\begin{table}[hp!]
\begin{tabular}{lllr}
@@ -35,85 +39,129 @@ In \eq{eq:linsys_FV} we choose units of $mV \equiv 10^{-3}V$ for each term becau
\label{tbl:units}
\end{table}
-%------------------------------------------
-\subsubsection{current terms}
-%------------------------------------------
-Membrane current is calculated as follows $i_m = \overline{g}(E-V)$, with units
+\subsection{Left Hand Side}
+First, we calculate the units for the term $\frac{\sigma_i \cmi}{\Delta t}$, as follows
\begin{align}
- \unit{ i_m } &= \unit{ \overline{g} } \unit{ V } \nonumber \\
- &= 10^{4} \cdot A\cdot V^{-1}\cdot m^{-2} \cdot 10^{-3} \cdot V \nonumber \\
- &= 10 \cdot A \cdot m^{-2}. \label{eq:im_unit}
+ \left[ \frac{\sigma_i \cmi}{\Delta t} \right]
+ &=
+ \frac{10^{-12}\cdot m^2 \cdot s\cdot A\cdot V^{-1}\cdot m^{-2} } {10^{-3}\cdot s}
+ \nonumber \\
+ &=
+ 10^{-9} \cdot A \cdot V^{-1}. \label{eq:units_lhs_diag}
\end{align}
-The point process currents are calculated as point sources which must be turned into current densities as follows $i_m = g_s(E-V)/\sigma_i$.
-The units for the synaptic conductance $g_s$ are $\mu S$, so the units are calculated as follows
+The units of $\alpha_{i,j}$ are:
\begin{align}
- \unit{ i_m } &= \unit{ g_s } \unit{ V } \unit{\sigma_i}^-1 \nonumber \\
- &= 10^{-6} \cdot A\cdot V^{-1} \cdot 10^{-3} \cdot 10^{12} \cdot m^{-2} \nonumber \\
- &= 10^{3} \cdot A \cdot m^{-2}, \label{eq:ims_unit}
+ \left[ \alpha_{i,j} \right]
+ &=
+ \left[ \frac{\sigma_{ij}}{ r_L \Delta x_{ij}} \right] \nonumber \\
+ &=
+ \frac{10^{-12}\cdot m^2}{ 10^{-2} \cdot A^{-1}\cdot V\cdot m \cdot 10^{-6}\cdot m}
+ \nonumber \\
+ &=
+ 10^{-4}\cdot A \cdot V^{-1}, \label{eq:units_lhs_lu}
\end{align}
-which must be scaled by $10^{2}$ to match that of of the density channels in \eq{eq:im_unit}.
+which can be scaled by $10^{5}$ to have the same scale as in \eq{eq:units_lhs_diag}.
+Thus, the RHS with scaling for units is:
+\begin{equation}
+ \label{eq:linsys_LHS_scaled}
+ \left[
+ \frac{\sigma_i \cmi}{\Delta t} + \sum_{j\in\mathcal{N}_i} {10^5\cdot\alpha_{ij}}
+ \right]
+ V_i^{k+1} - \sum_{j\in\mathcal{N}_i} { 10^5\cdot\alpha_{ij} V_j^{k+1}}.
+\end{equation}
+The implementation folds the $10^5$ factor into the $\alpha_{ij}$ terms when they are used to calculate the invariant component of the matrix during the initialization phase.
-The injected current $I_e$ has units $nA$, which has to be expressed in terms of current per unit area $i_e=I_e / \sigma_i$ with units
+After this scaling, the units of the LHS are
\begin{align}
- \unit{ i_e } &= \unit{ I_e } \unit{ \sigma_i }^{-1} \nonumber \\
- &= 10^{-9}\cdot A \cdot 10^{12} \cdot m^{-2} \nonumber \\
- &= 10^{3} \cdot A \cdot m ^{-2}, \label{eq:ie_unit}
+ \text{units on LHS}
+ &= (10^{-9} \cdot A \cdot V^{-1})( 10^{-3} \cdot V) \nonumber \\
+ &= 10^{-12} \cdot A. \label{eq:balanced_units}
\end{align}
-which must be scaled by $10^2$ to match $i_m$ in \eq{eq:im_unit}.
-The units for the flux coefficent can be calculated as follows:
-\begin{align}
- \unit{ \frac{\Delta t}{c_m} } &= 10^{-3} \cdot s \cdot s^{-1}\cdot A^{-1}\cdot V\cdot m^2 \nonumber \\
- &= 10^{-3} \cdot A^{-1} \cdot V\cdot m^2. \label{eq:dtcm_unit}
-\end{align}
-From \eq{eq:im_unit} and \eq{eq:dtcm_unit} that the units of the full current term are
-\begin{align}
- \unit{ \frac{\Delta t}{c_m}\left(i_m - i_e\right) }
- &= 10^{-3} \cdot A^{-1} \cdot V\cdot m^2 \cdot 10 \cdot A \cdot m^{-2} \nonumber \\
- &= 10^{-2} \cdot V,
-\end{align}
-which must be scaled by $10$ to match the units of $mV\equiv10^{-3}V$.
-%------------------------------------------
-\subsubsection{flux terms}
-%------------------------------------------
-The coefficients in the linear system have the units
+\subsection{Right Hand Side}
+The first term on the RHS has the same scaling factor as the LHS, so does not need to be changed.
+
+Density channels and point processes describe the membrane current differently in NMODL;
+as current densities ($10\cdot A\cdot m^{-2}$) and currents ($10^{-9}\cdot A$) respectively.
+The current term can be written as follows:
\begin{equation}
- \unit{ \frac{\Delta t\alpha_{ij}}{\sigma_i} }
- =
- \unit{ \frac{\Delta t \sigma_{ij} } {c_m r_L \Delta x_{ij} \sigma_i} }
- =
- \unit{ \frac{\Delta t } {c_m r_L \Delta x_{ij} } },
+ \sigma_i (i_m - i_e) \equiv \sigma_i \bar{i}_m + I_m - I_e,
\end{equation}
-where we we simplify by noting that $\unit{\sigma_{ij}}=\unit{\sigma_i}$.
-The units of the term $c_m r_L$ on the denominator are calculated as follows
-\begin{align}
- \unit{c_m r_L}
- &= s \cdot A \cdot V^{-1} \cdot m^{-2} \cdot 10^{-2} \cdot A^{-1} \cdot V \cdot m \nonumber \\
- &= 10^{-2} \cdot s \cdot m^{-1},
-\end{align}
-so the units of the denominator are
+where $\bar{i}_m$ is the current density contribution from ion channels, $I_m$ is the current contribution from synapses, and $I_e$ is current contribution from electrodes.
+
+The units of the current density as calculated via NMODL are
+\begin{equation}
+ \label{eq:im_unit}
+ \unit{ \bar{i}_m } = \unit{ \overline{g} } \unit{ V }
+ = 10 \cdot A \cdot m^{-2},
+\end{equation}
+so the units of the current from density channels are
+\begin{equation}
+ \unit{\sigma_i \bar{i}_m} = 10^{-12}\cdot{m}^2 \cdot 10 \cdot A \cdot m^{-2} = 10^{-11}\cdot A.
+\end{equation}
+Hence, the $\sigma_i \bar{i}_m$ term must be scaled by 10 to match units.
+
+Likewise the units of synapse and electrode current
+\begin{equation}
+ \label{eq:Im_unit}
+ \unit{ I_e } = \unit{ I_m } = \unit{ g_s } \unit{ V }
+ = 10^{-9}\cdot A,
+\end{equation}
+which must be scaled by $10^3$ to match units.
+
+The properly scaled RHS is
+\begin{equation}
+ \label{eq:linsys_RHS_scaled}
+ \frac{\sigma_i \cmi}{\Delta t} V_i^k -
+ (10\cdot\sigma_i \bar{i}_m + 10^3(I_m - I_e)).
+\end{equation}
+
+\subsection{Putting It Together}
+Hey ho, let's go: from \eq{eq:linsys_LHS_scaled} and \eq{eq:linsys_RHS_scaled} the full scaled linear system is
\begin{align}
- \unit{c_m r_L \Delta x_{ij}}
- &= 10^{-2} \cdot s \cdot m^{-1} \cdot 10^{-6} \cdot m \nonumber \\
- &= 10^{-8} \cdot s,
+ &
+ \left[
+ \frac{\sigma_i \cmi}{\Delta t} + \sum_{j\in\mathcal{N}_i} {10^5\cdot\alpha_{ij}}
+ \right]
+ V_i^{k+1} - \sum_{j\in\mathcal{N}_i} { 10^5\cdot\alpha_{ij} V_j^{k+1}} \nonumber \\
+ & =
+ \frac{\sigma_i \cmi}{\Delta t} V_i^k -
+ (10\cdot\sigma_i \bar{i}_m + 10^3(I_m - I_e)).
\end{align}
-and hence
+This can be expressed more generally in terms of weights
\begin{align}
- \unit{\frac{\Delta t } {c_m r_L \Delta x_{ij} }}
- &= 10^{8} \cdot s^{-1} \cdot 10^{-3} \cdot s \nonumber \\
- &= 10^{5}.
+ &
+ \left[
+ g_i + \sum_{j\in\mathcal{N}_i} {g_{ij}}
+ \right]
+ V_i^{k+1} - \sum_{j\in\mathcal{N}_i} { g_{ij} V_j^{k+1}} \nonumber \\
+ & =
+ g_i V_i^k -
+ (w_i^d \bar{i}_m + w_i^p(I_m - I_e)),
\end{align}
+which can be expressed more compactly as
+\begin{equation}
+ Gv=i,
+\end{equation}
+where $G\in\mathbb{R}^{n\times n}$ is the conductance matrix, and $v, i \in \mathbb{R}^{n}$ are voltage and current vectors respectively.
-So, the terms with $\alpha_{ij}$ must be scaled by $10^5$ to match the units of $mV$.
-%------------------------------------------
-\subsubsection{discretization with scaling}
-%------------------------------------------
-Here is something that I wish the NEURON documentation had provided:
-\begin{align}
-& V_i^{k+1} + \sum_{j\in\mathcal{N}_i} {10^5 \cdot \frac{\Delta t \alpha_{ij}}{\sigma_i} (V_i^{k+1}-V_j^{k+1})} \nonumber \\
-& = V_i^k - 10\cdot \frac{\Delta t}{c_m}(i_m^{k} - 10^2\cdot I_e/\sigma_i).
-\end{align}
+In NestMC the weights are chosen such that the conductance has units $\mu S$, voltage has units $mV$ and current has units $nA$.
+
+ \begin{center}
+
+ \begin{tabular}{llll}
+ \hline
+ weight & value & units & SI \\
+ \hline
+ $g_i$ & $10^{-3}\frac{\sigma_i \cmi}{\Delta t}$ & $\mu S$ & $10^{-6} \cdot A\cdot V^{-1}$ \\
+ $g_{ij}$ & $10^2\alpha_{ij}$ & $\mu S$ & $10^{-6} \cdot A\cdot V^{-1}$ \\
+ $w_i^d$ & $10^{-2}\cdot\sigma_i$ & $10^2\mu m^{-2}$ & $10^{-10}m^2$ \\
+ $w_i^p$ & $1$ & $1$ & $1$ \\
+ \hline
+ \end{tabular}
+
+ \end{center}
%------------------------------------------
\subsection{Supplementary Unit Information}
%------------------------------------------
diff --git a/modcc/cprinter.cpp b/modcc/cprinter.cpp
index a34a8c5432d9ffa99c836cf785892aed790a68d2..beb97a673d28f2b6c24bbb250dd579f945dfd9cb 100644
--- a/modcc/cprinter.cpp
+++ b/modcc/cprinter.cpp
@@ -80,14 +80,14 @@ CPrinter::CPrinter(Module &m, bool o)
text_.add_line("};");
text_.add_line(tname + " ion_" + ion.name + ";");
}
- text_.add_line();
//////////////////////////////////////////////
// constructor
//////////////////////////////////////////////
int num_vars = array_variables.size();
- text_.add_line(class_name + "(view vec_v, view vec_i, const_iview node_index)");
- text_.add_line(": base(vec_v, vec_i, node_index)");
+ text_.add_line();
+ text_.add_line(class_name + "(view vec_v, view vec_i, array&& weights, iarray&& node_index)");
+ text_.add_line(": base(vec_v, vec_i, std::move(node_index))");
text_.add_line("{");
text_.increase_indentation();
text_.add_gutter() << "size_type num_fields = " << num_vars << ";";
@@ -124,8 +124,16 @@ CPrinter::CPrinter(Module &m, bool o)
}
text_.end_line();
}
-
text_.add_line();
+
+ // copy in the weights if this is a density mechanism
+ if (m.kind() == moduleKind::density) {
+ text_.add_line("// add the user-supplied weights for converting from current density");
+ text_.add_line("// to per-compartment current in nA");
+ text_.add_line("memory::copy(weights, weights_(0, size()));");
+ text_.add_line();
+ }
+
text_.add_line("// set initial values for variables and parameters");
for(auto const& var : array_variables) {
double val = var->value();
@@ -349,7 +357,6 @@ CPrinter::CPrinter(Module &m, bool o)
text_.add_line();
text_.add_line("using base::vec_v_;");
text_.add_line("using base::vec_i_;");
- text_.add_line("using base::vec_area_;");
text_.add_line("using base::node_index_;");
text_.add_line();
diff --git a/modcc/cudaprinter.cpp b/modcc/cudaprinter.cpp
index 9d60a3f6630528ce4bdc52d4503ce05316a14c05..c6c33ead52ef2131d3efc6582a603026ee5f3868 100644
--- a/modcc/cudaprinter.cpp
+++ b/modcc/cudaprinter.cpp
@@ -82,9 +82,6 @@ CUDAPrinter::CUDAPrinter(Module &m, bool o)
param_pack.push_back("vec_v_.data()");
param_pack.push_back("vec_i_.data()");
- text_.add_line("T* vec_area;");
- param_pack.push_back("vec_area_.data()");
-
text_.add_line("// node index information");
text_.add_line("I* ni;");
text_.add_line("unsigned long n_;");
@@ -191,7 +188,7 @@ CUDAPrinter::CUDAPrinter(Module &m, bool o)
int num_vars = array_variables.size();
text_.add_line();
- text_.add_line(class_name + "(view vec_v, view vec_i, iarray&& node_index) :");
+ text_.add_line(class_name + "(view vec_v, view vec_i, array&& weights, iarray&& node_index):");
text_.add_line(" base(vec_v, vec_i, std::move(node_index))");
text_.add_line("{");
text_.increase_indentation();
@@ -223,6 +220,7 @@ CUDAPrinter::CUDAPrinter(Module &m, bool o)
array_variables[i]->name() + " = data_("
+ std::to_string(i) + "*field_size, " + std::to_string(i+1) + "*field_size);");
}
+ text_.add_line();
for(auto const& var : array_variables) {
double val = var->value();
@@ -232,6 +230,15 @@ CUDAPrinter::CUDAPrinter(Module &m, bool o)
text_.add_line("memory::fill(" + var->name() + ", " + std::to_string(val) + ");");
}
}
+ text_.add_line();
+
+ // copy in the weights if this is a density mechanism
+ if (m.kind() == moduleKind::density) {
+ text_.add_line("// add the user-supplied weights for converting from current density");
+ text_.add_line("// to per-compartment current in nA");
+ text_.add_line("memory::copy(weights, weights_(0, size()));");
+ text_.add_line();
+ }
text_.decrease_indentation();
text_.add_line("}");
@@ -480,7 +487,6 @@ CUDAPrinter::CUDAPrinter(Module &m, bool o)
text_.add_line("using base::vec_v_;");
text_.add_line("using base::vec_i_;");
- text_.add_line("using base::vec_area_;");
text_.add_line("using base::node_index_;");
text_.add_line();
text_.add_line("param_pack_type param_pack_;");
diff --git a/modcc/module.cpp b/modcc/module.cpp
index 3e8fe3d77da7c2e2be24493cc04c00393f111f03..be76e87c8d19bc7cfe652bf05915963cbe92e12f 100644
--- a/modcc/module.cpp
+++ b/modcc/module.cpp
@@ -551,8 +551,8 @@ bool Module::semantic() {
has_current_update = true;
}
}
- if(has_current_update && kind()==moduleKind::point) {
- block.emplace_back(Parser("current_ = 100. * current_ / area_").parse_line_expression());
+ if(has_current_update && kind()==moduleKind::density) {
+ block.emplace_back(Parser("current_ = weights_ * current_").parse_line_expression());
}
auto v = make_unique<ConstantFolderVisitor>();
@@ -594,6 +594,11 @@ void Module::add_variables_to_symbols() {
create_variable("t", rangeKind::scalar, accessKind::read);
create_variable("dt", rangeKind::scalar, accessKind::read);
+ // density mechanisms use a vector of weights from current densities to
+ // units of nA
+ if (kind()==moduleKind::density) {
+ create_variable("weights_", rangeKind::range, accessKind::read);
+ }
// add indexed variables to the table
auto create_indexed_variable = [this]
@@ -613,8 +618,6 @@ void Module::add_variables_to_symbols() {
accessKind::write, ionKind::none, Location());
create_indexed_variable("v", "vec_v", tok::eq,
accessKind::read, ionKind::none, Location());
- create_indexed_variable("area_", "vec_area", tok::eq,
- accessKind::read, ionKind::none, Location());
// add state variables
for(auto const &var : state_block()) {
diff --git a/src/algorithms.hpp b/src/algorithms.hpp
index 9a3ecdf7f2197e233fd2c848aa86cf124aef3a41..22325d224ce576cdaa32d759d71368845e51c944 100644
--- a/src/algorithms.hpp
+++ b/src/algorithms.hpp
@@ -7,6 +7,7 @@
#include <type_traits>
#include <vector>
+#include <util/compat.hpp>
#include <util/debug.hpp>
#include <util/meta.hpp>
#include <util/range.hpp>
@@ -393,6 +394,28 @@ auto index_into(const Sub& sub, const Super& super)
return util::make_range(begin, end);
}
+/// Binary search, because std::binary_search doesn't return information
+/// about where a match was found.
+template <typename It, typename T>
+It binary_find(It b, It e, const T& value) {
+ auto it = std::lower_bound(b, e, value);
+ return it==e ? e : (*it==value ? it : e);
+}
+
+template <typename Seq, typename T>
+auto binary_find(const Seq& seq, const T& value)
+ -> decltype(binary_find(std::begin(seq), std::end(seq), value))
+{
+ return binary_find(std::begin(seq), compat::end(seq), value);
+}
+
+template <typename Seq, typename T>
+auto binary_find(Seq& seq, const T& value)
+ -> decltype(binary_find(std::begin(seq), std::end(seq), value))
+{
+ return binary_find(std::begin(seq), compat::end(seq), value);
+}
+
} // namespace algorithms
} // namespace mc
} // namespace nest
diff --git a/src/backends/fvm_gpu.hpp b/src/backends/fvm_gpu.hpp
index 0c0af1d491e0e3006f70834bfca1a1f92d031c7a..0718db34a4466d4288a7ccf75efd10bc3335ae11 100644
--- a/src/backends/fvm_gpu.hpp
+++ b/src/backends/fvm_gpu.hpp
@@ -17,7 +17,7 @@ namespace gpu {
template <typename T, typename I>
struct matrix_solve_param_pack {
T* d;
- T* u;
+ const T* u;
T* rhs;
const I* p;
const I* cell_index;
@@ -30,14 +30,13 @@ struct matrix_solve_param_pack {
template <typename T, typename I>
struct matrix_update_param_pack {
T* d;
- T* u;
+ const T* u;
T* rhs;
- const T* sigma;
- const T* alpha_d;
- const T* alpha;
+ const T* invariant_d;
+ const T* cv_capacitance;
+ const T* face_conductance;
const T* voltage;
const T* current;
- const T* cv_capacitance;
I n;
};
@@ -82,29 +81,36 @@ struct backend {
/// Hines matrix assembly interface
struct matrix_assembler {
matrix_update_param_pack<value_type, size_type> params;
- array alpha_d;
+
+ // the invariant part of the matrix diagonal
+ array invariant_d; // [μS]
matrix_assembler() = default;
matrix_assembler(
view d, view u, view rhs, const_iview p,
- const_view sigma, const_view alpha,
- const_view voltage, const_view current, const_view cv_capacitance)
+ const_view cv_capacitance,
+ const_view face_conductance,
+ const_view voltage,
+ const_view current)
{
auto n = d.size();
- host_array alpha_d_tmp(n, 0);
+ host_array invariant_d_tmp(n, 0);
+ // make a copy of the conductance on the host
+ host_array face_conductance_tmp = face_conductance;
for(auto i: util::make_span(1u, n)) {
- alpha_d_tmp[i] += alpha[i];
+ auto gij = face_conductance_tmp[i];
- // add contribution to the diagonal of parent
- alpha_d_tmp[p[i]] += alpha[i];
+ u[i] = -gij;
+ invariant_d_tmp[i] += gij;
+ invariant_d_tmp[p[i]] += gij;
}
- alpha_d = alpha_d_tmp;
+ invariant_d = invariant_d_tmp;
params = {
d.data(), u.data(), rhs.data(),
- sigma.data(), alpha_d.data(), alpha.data(),
- voltage.data(), current.data(), cv_capacitance.data(), size_type(n)};
+ invariant_d.data(), cv_capacitance.data(), face_conductance.data(),
+ voltage.data(), current.data(), size_type(n)};
}
void assemble(value_type dt) {
@@ -151,6 +157,7 @@ struct backend {
static mechanism make_mechanism(
const std::string& name,
view vec_v, view vec_i,
+ const std::vector<value_type>& weights,
const std::vector<size_type>& node_indices)
{
if (!has_mechanism(name)) {
@@ -158,20 +165,20 @@ struct backend {
}
return mech_map_.find(name)->
- second(vec_v, vec_i, memory::make_const_view(node_indices));
+ second(vec_v, vec_i, memory::make_const_view(weights), memory::make_const_view(node_indices));
}
static bool has_mechanism(const std::string& name) { return mech_map_.count(name)>0; }
private:
- using maker_type = mechanism (*)(view, view, iarray&&);
+ using maker_type = mechanism (*)(view, view, array&&, iarray&&);
static std::map<std::string, maker_type> mech_map_;
template <template <typename> class Mech>
- static mechanism maker(view vec_v, view vec_i, iarray&& node_indices) {
+ static mechanism maker(view vec_v, view vec_i, array&& weights, iarray&& node_indices) {
return mechanisms::make_mechanism<Mech<backend>>
- (vec_v, vec_i, std::move(node_indices));
+ (vec_v, vec_i, std::move(weights), std::move(node_indices));
}
};
@@ -219,13 +226,13 @@ __global__
void assemble_matrix(matrix_update_param_pack<T, I> params, T dt) {
auto tid = threadIdx.x + blockDim.x*blockIdx.x;
- T factor_lhs = 1e5*dt;
- T factor_rhs = 10.*dt;
+ T factor = 1e-3/dt;
if(tid < params.n) {
- params.d[tid] = params.sigma[tid] + factor_lhs*params.alpha_d[tid];
- params.u[tid] = -factor_lhs*params.alpha[tid];
- params.rhs[tid] = params.sigma[tid] *
- (params.voltage[tid] - factor_rhs/params.cv_capacitance[tid]*params.current[tid]);
+ auto gi = factor * params.cv_capacitance[tid];
+
+ params.d[tid] = gi + params.invariant_d[tid];
+
+ params.rhs[tid] = gi*params.voltage[tid] - params.current[tid];
}
}
diff --git a/src/backends/fvm_multicore.hpp b/src/backends/fvm_multicore.hpp
index 71a0ba0d06f370ed09f5ce6513fce16223477bc9..8166525ee7d01c89bf77218caf21513a2d2dc422 100644
--- a/src/backends/fvm_multicore.hpp
+++ b/src/backends/fvm_multicore.hpp
@@ -59,53 +59,53 @@ struct backend {
}
}
- // it might be acceptable to have the entire builder defined here
- // because the storage might need to be back end specific
struct matrix_assembler {
- view d;
- view u;
- view rhs;
+ view d; // [μS]
+ view u; // [μS]
+ view rhs; // [nA]
const_iview p;
- const_view sigma;
- const_view alpha;
- const_view voltage;
- const_view current;
- const_view cv_capacitance;
+ const_view cv_capacitance; // [pF]
+ const_view face_conductance; // [μS]
+ const_view voltage; // [mV]
+ const_view current; // [nA]
- array alpha_d;
+ // the invariant part of the matrix diagonal
+ array invariant_d; // [μS]
matrix_assembler() = default;
matrix_assembler(
view d, view u, view rhs, const_iview p,
- const_view sigma, const_view alpha,
- const_view voltage, const_view current, const_view cv_capacitance)
+ const_view cv_capacitance,
+ const_view face_conductance,
+ const_view voltage,
+ const_view current)
:
d{d}, u{u}, rhs{rhs}, p{p},
- sigma{sigma}, alpha{alpha},
- voltage{voltage}, current{current}, cv_capacitance{cv_capacitance}
+ cv_capacitance{cv_capacitance}, face_conductance{face_conductance},
+ voltage{voltage}, current{current}
{
auto n = d.size();
- alpha_d = array(n, 0);
- for(auto i: util::make_span(1u, n)) {
- alpha_d[i] += alpha[i];
+ invariant_d = array(n, 0);
+ for (auto i: util::make_span(1u, n)) {
+ auto gij = face_conductance[i];
- // add contribution to the diagonal of parent
- alpha_d[p[i]] += alpha[i];
+ u[i] = -gij;
+ invariant_d[i] += gij;
+ invariant_d[p[i]] += gij;
}
}
void assemble(value_type dt) {
auto n = d.size();
- value_type factor_lhs = 1e5*dt;
- value_type factor_rhs = 1e1*dt; // units: 10·ms/(F/m^2)·(mA/cm^2) ≡ mV
+ value_type factor = 1e-3/dt;
for (auto i: util::make_span(0u, n)) {
- d[i] = sigma[i] + factor_lhs*alpha_d[i];
- u[i] = -factor_lhs*alpha[i];
- // the RHS of the linear system is
- // cv_area * (V - dt/cm*(im - ie))
- rhs[i] = sigma[i]*(voltage[i] - factor_rhs/cv_capacitance[i]*current[i]);
+ auto gi = factor*cv_capacitance[i];
+
+ d[i] = gi + invariant_d[i];
+
+ rhs[i] = gi*voltage[i] - current[i];
}
}
};
@@ -120,16 +120,19 @@ struct backend {
static mechanism make_mechanism(
const std::string& name,
view vec_v, view vec_i,
+ const std::vector<value_type>& weights,
const std::vector<size_type>& node_indices)
{
if (!has_mechanism(name)) {
throw std::out_of_range("no mechanism in database : " + name);
}
- return mech_map_.find(name)->second(vec_v, vec_i, iarray(node_indices));
+ return mech_map_.find(name)->second(vec_v, vec_i, array(weights), iarray(node_indices));
}
- static bool has_mechanism(const std::string& name) { return mech_map_.count(name)>0; }
+ static bool has_mechanism(const std::string& name) {
+ return mech_map_.count(name)>0;
+ }
static std::string name() {
return "cpu";
@@ -137,13 +140,13 @@ struct backend {
private:
- using maker_type = mechanism (*)(view, view, iarray&&);
+ using maker_type = mechanism (*)(view, view, array&&, iarray&&);
static std::map<std::string, maker_type> mech_map_;
template <template <typename> class Mech>
- static mechanism maker(view vec_v, view vec_i, iarray&& node_indices) {
+ static mechanism maker(view vec_v, view vec_i, array&& weights, iarray&& node_indices) {
return mechanisms::make_mechanism<Mech<backend>>
- (vec_v, vec_i, std::move(node_indices));
+ (vec_v, vec_i, std::move(weights), std::move(node_indices));
}
};
diff --git a/src/cell_group.hpp b/src/cell_group.hpp
index 203423c125cd22ef31f5f7e0e2acaa49f8a26821..606acf586579b2e043ca164d1b60c86e9f7bf78a 100644
--- a/src/cell_group.hpp
+++ b/src/cell_group.hpp
@@ -84,6 +84,10 @@ public:
}
}
+ time_type min_step(time_type dt) {
+ return 0.1*dt;
+ }
+
void advance(time_type tfinal, time_type dt) {
while (cell_.time()<tfinal) {
// take any pending samples
@@ -105,14 +109,22 @@ public:
// look for events in the next time step
time_type tstep = cell_.time()+dt;
tstep = std::min(tstep, tfinal);
-
auto next = events_.pop_if_before(tstep);
- time_type tnext = next ? next->time: tstep;
+
+ // apply events that are due within the smallest allowed time step.
+ while (next && (next->time-cell_.time()) < min_step(dt)) {
+ auto handle = get_target_handle(next->target);
+ cell_.deliver_event(handle, next->weight);
+ next = events_.pop_if_before(tstep);
+ }
// integrate cell state
+ time_type tnext = next ? next->time: tstep;
cell_.advance(tnext - cell_.time());
+
if (!cell_.is_physical_solution()) {
- std::cerr << "warning: solution out of bounds\n";
+ std::cerr << "warning: solution out of bounds for cell "
+ << gid_base_ << " at t " << cell_.time() << " ms\n";
}
PE("events");
@@ -127,13 +139,6 @@ public:
if (next) {
auto handle = get_target_handle(next->target);
cell_.deliver_event(handle, next->weight);
- // apply events that are due within some epsilon of the current
- // time step. This should be a parameter. e.g. with for variable
- // order time stepping, use the minimum possible time step size.
- while(auto e = events_.pop_if_before(cell_.time()+dt/10.)) {
- auto handle = get_target_handle(e->target);
- cell_.deliver_event(handle, e->weight);
- }
}
PL();
}
diff --git a/src/fvm_multicell.hpp b/src/fvm_multicell.hpp
index b60059e4a184758c91726c733956804e76ad9a43..3719648a49ad3de6db1cf2c8fa98815bc999e1ec 100644
--- a/src/fvm_multicell.hpp
+++ b/src/fvm_multicell.hpp
@@ -198,16 +198,16 @@ private:
/// cv_areas_[i] is the surface area of CV i [µm^2]
array 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);
- array face_alpha_; // [µm·m^2/cm/s ≡ 10^5 µm^2/ms]
+ /// CV i and its parent, required when constructing linear system [µS]
+ /// face_conductance_[i] = area_face / (r_L * delta_x);
+ array face_conductance_;
- /// cv_capacitance_[i] is the capacitance of CV i per unit area (i.e. c_m) [F/m^2]
- array cv_capacitance_;
+ /// cv_capacitance_[i] is the capacitance of CV membrane [pF]
+ /// C_m = area*c_m
+ array cv_capacitance_; // units [µm^2*F*m^-2 = pF]
- /// the average current density over the surface of each CV [mA/cm^2]
- /// current_ = i_m - i_e
+ /// the transmembrane current over the surface of each CV [nA]
+ /// I = area*i_m - I_e
array current_;
/// the potential in each CV [mV]
@@ -223,12 +223,42 @@ private:
std::vector<std::pair<const array fvm_multicell::*, size_type>> probes_;
+ /// Compact representation of the control volumes into which a segment is
+ /// decomposed. Used to reconstruct the weights used to convert current
+ /// densities to currents for density channels.
+ struct segment_cv_range {
+ // the contribution to the surface area of the CVs that
+ // are at the beginning and end of the segment
+ std::pair<value_type, value_type> areas;
+
+ // the range of CVs in the segment, excluding the parent CV
+ std::pair<size_type, size_type> segment_cvs;
+
+ // The last CV in the parent segment, which corresponds to the
+ // first CV in this segment.
+ // Set to npos() if there is no parent (i.e. if soma)
+ size_type parent_cv;
+
+ static constexpr size_type npos() {
+ return std::numeric_limits<size_type>::max();
+ }
+
+ // the number of CVs (including the parent)
+ std::size_t size() const {
+ return segment_cvs.second-segment_cvs.first + (parent_cv==npos() ? 0 : 1);
+ }
+
+ bool has_parent() const {
+ return parent_cv != npos();
+ }
+ };
+
// perform area and capacitance calculation on initialization
- void compute_cv_area_unnormalized_capacitance(
+ segment_cv_range compute_cv_area_capacitance(
std::pair<size_type, size_type> comp_ival,
const segment* seg,
const std::vector<size_type>& parent,
- std::vector<value_type>& tmp_face_alpha,
+ std::vector<value_type>& tmp_face_conductance,
std::vector<value_type>& tmp_cv_areas,
std::vector<value_type>& tmp_cv_capacitance
);
@@ -238,11 +268,12 @@ private:
//////////////////////////////// Implementation ////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
template <typename Backend>
-void fvm_multicell<Backend>::compute_cv_area_unnormalized_capacitance(
+typename fvm_multicell<Backend>::segment_cv_range
+fvm_multicell<Backend>::compute_cv_area_capacitance(
std::pair<size_type, size_type> comp_ival,
const segment* seg,
const std::vector<size_type>& parent,
- std::vector<value_type>& tmp_face_alpha,
+ std::vector<value_type>& tmp_face_conductance,
std::vector<value_type>& tmp_cv_areas,
std::vector<value_type>& tmp_cv_capacitance)
{
@@ -251,6 +282,8 @@ void fvm_multicell<Backend>::compute_cv_area_unnormalized_capacitance(
auto ncomp = comp_ival.second-comp_ival.first;
+ segment_cv_range cv_range;
+
if (auto soma = seg->as_soma()) {
// confirm assumption that there is one compartment in soma
if (ncomp!=1) {
@@ -258,9 +291,14 @@ void fvm_multicell<Backend>::compute_cv_area_unnormalized_capacitance(
}
auto i = comp_ival.first;
auto area = math::area_sphere(soma->radius());
+ auto c_m = soma->mechanism("membrane").get("c_m").value;
tmp_cv_areas[i] += area;
- tmp_cv_capacitance[i] += area * soma->mechanism("membrane").get("c_m").value;
+ tmp_cv_capacitance[i] += area*c_m;
+
+ cv_range.segment_cvs = {comp_ival.first, comp_ival.first+1};
+ cv_range.areas = {0.0, area};
+ cv_range.parent_cv = segment_cv_range::npos();
}
else if (auto cable = seg->as_cable()) {
// Loop over each compartment in the cable
@@ -275,9 +313,8 @@ void fvm_multicell<Backend>::compute_cv_area_unnormalized_capacitance(
// the respective control volumes, and the volumes and lengths of
// each half are used to calculate the flux coefficients that
// for the connection between the two control volumes and which
- // (after scaling by inverse capacitance) is stored in
- // `face_alpha[i]`.
- //
+ // is stored in `face_conductance[i]`.
+ //
//
// +------- cv j --------+------- cv i -------+
// | | |
@@ -298,16 +335,22 @@ void fvm_multicell<Backend>::compute_cv_area_unnormalized_capacitance(
auto divs = div_compartments<div_compartment_integrator>(cable, ncomp);
+ // assume that this segment has a parent, which is the case so long
+ // as the soma is the root of all cell trees.
+ cv_range.parent_cv = parent[comp_ival.first];
+ cv_range.segment_cvs = comp_ival;
+ cv_range.areas = {divs(0).left.area, divs(ncomp-1).right.area};
+
for (auto i: util::make_span(comp_ival)) {
const auto& div = divs(i-comp_ival.first);
auto j = parent[i];
// Conductance approximated by weighted harmonic mean of mean
// conductances in each half.
- //
+ //
// Mean conductances:
- // c₠= 1/h₠∫₠A(x)/R dx
- // c₂ = 1/h₂ ∫₂ A(x)/R dx
+ // g₠= 1/h₠∫₠A(x)/R dx
+ // g₂ = 1/h₂ ∫₂ A(x)/R dx
//
// where A(x) is the cross-sectional area, R is the bulk
// resistivity, h is the length of the interval and the
@@ -315,13 +358,18 @@ void fvm_multicell<Backend>::compute_cv_area_unnormalized_capacitance(
// Equivalently, in terms of the semi-compartment volumes
// Vâ‚ and Vâ‚‚:
//
- // câ‚ = 1/R·Vâ‚/hâ‚
- // c₂ = 1/R·V₂/h₂
+ // gâ‚ = 1/R·Vâ‚/hâ‚
+ // g₂ = 1/R·V₂/h₂
//
// Weighted harmonic mean, with h = hâ‚+hâ‚‚:
//
- // c = (hâ‚/h·c₯¹+hâ‚‚/h·c₂¯¹)¯¹
+ // g = (hâ‚/h·g₯¹+hâ‚‚/h·g₂¯¹)¯¹
// = 1/R · hVâ‚Vâ‚‚/(h₂²Vâ‚+h₲Vâ‚‚)
+ //
+ // the following units are used
+ // lengths : μm
+ // areas : μm^2
+ // volumes : μm^3
auto h1 = div.left.length;
auto V1 = div.left.volume;
@@ -330,7 +378,9 @@ void fvm_multicell<Backend>::compute_cv_area_unnormalized_capacitance(
auto h = h1+h2;
auto conductance = 1/r_L*h*V1*V2/(h2*h2*V1+h1*h1*V2);
- tmp_face_alpha[i] = conductance / (c_m * h);
+ // the scaling factor of 10^2 is to convert the quantity
+ // to micro Siemens [μS]
+ tmp_face_conductance[i] = 1e2 * conductance / h;
auto al = div.left.area;
auto ar = div.right.area;
@@ -344,6 +394,8 @@ void fvm_multicell<Backend>::compute_cv_area_unnormalized_capacitance(
else {
throw std::domain_error("FVM lowering encountered unsuported segment type");
}
+
+ return cv_range;
}
template <typename Backend>
@@ -361,6 +413,7 @@ void fvm_multicell<Backend>::initialize(
using util::size;
using util::sort_by;
using util::transform_view;
+ using util::subrange_view;
// count total detectors, targets and probes for validation of handle container sizes
std::size_t detectors_count = 0u;
@@ -383,7 +436,7 @@ void fvm_multicell<Backend>::initialize(
voltage_ = array(ncomp, resting_potential_);
// create maps for mechanism initialization.
- std::map<std::string, std::vector<std::pair<size_type, size_type>>> mech_map;
+ std::map<std::string, std::vector<segment_cv_range>> mech_map;
std::vector<std::vector<cell_lid_type>> syn_mech_map;
std::map<std::string, std::size_t> syn_mech_indices;
@@ -395,10 +448,13 @@ void fvm_multicell<Backend>::initialize(
auto detector_hi = detector_handles.begin();
auto probe_hi = probe_handles.begin();
- // allocate scratch vectors
- std::vector<value_type> tmp_face_alpha(ncomp);
- std::vector<value_type> tmp_cv_areas(ncomp);
- std::vector<value_type> tmp_cv_capacitance(ncomp);
+ // Allocate scratch storage for calculating quantities used to build the
+ // linear system: these will later be copied into target-specific storage
+ // as need be.
+ // Initialize to zero, because the results therin are calculated via accumulation.
+ std::vector<value_type> tmp_face_conductance(ncomp, 0.);
+ std::vector<value_type> tmp_cv_areas(ncomp, 0.);
+ std::vector<value_type> tmp_cv_capacitance(ncomp, 0.);
// Iterate over the input cells and build the indexes etc that descrbe the
// fused cell group. On completion:
@@ -421,7 +477,7 @@ void fvm_multicell<Backend>::initialize(
auto seg_num_compartments =
transform_view(c.segments(), [](const segment_ptr& s) { return s->num_compartments(); });
- auto nseg = seg_num_compartments.size();
+ const auto nseg = seg_num_compartments.size();
std::vector<cell_lid_type> seg_comp_bounds;
auto seg_comp_part =
@@ -431,13 +487,13 @@ void fvm_multicell<Backend>::initialize(
const auto& seg = c.segment(j);
const auto& seg_comp_ival = seg_comp_part[j];
- compute_cv_area_unnormalized_capacitance(
+ auto cv_range = compute_cv_area_capacitance(
seg_comp_ival, seg, group_parent_index,
- tmp_face_alpha, tmp_cv_areas, tmp_cv_capacitance);
+ tmp_face_conductance, tmp_cv_areas, tmp_cv_capacitance);
for (const auto& mech: seg->mechanisms()) {
if (mech.name()!="membrane") {
- mech_map[mech.name()].push_back(seg_comp_ival);
+ mech_map[mech.name()].push_back(cv_range);
}
}
}
@@ -458,8 +514,8 @@ void fvm_multicell<Backend>::initialize(
auto& map_entry = syn_mech_map[syn_mech_index];
- size_type syn_comp = comp_ival.first+find_cv_index(syn.location, graph);
- map_entry.push_back(syn_comp);
+ auto syn_cv = comp_ival.first + find_cv_index(syn.location, graph);
+ map_entry.push_back(syn_cv);
}
// add the stimuli
@@ -501,39 +557,66 @@ void fvm_multicell<Backend>::initialize(
EXPECTS(detectors_size==detectors_count);
EXPECTS(probes_size==probes_count);
- // normalize capacitance across cell
- for (auto i: util::make_span(0, ncomp)) {
- tmp_cv_capacitance[i] /= tmp_cv_areas[i];
- }
-
// store the geometric information in target-specific containers
- face_alpha_ = make_const_view(tmp_face_alpha);
- cv_areas_ = make_const_view(tmp_cv_areas);
- cv_capacitance_ = make_const_view(tmp_cv_capacitance);
+ face_conductance_ = make_const_view(tmp_face_conductance);
+ cv_areas_ = make_const_view(tmp_cv_areas);
+ cv_capacitance_ = make_const_view(tmp_cv_capacitance);
// initalize matrix
matrix_ = matrix_type(group_parent_index, cell_comp_bounds);
matrix_assembler_ = matrix_assembler(
matrix_.d(), matrix_.u(), matrix_.rhs(), matrix_.p(),
- cv_areas_, face_alpha_, voltage_, current_, cv_capacitance_);
+ cv_capacitance_, face_conductance_, voltage_, current_);
// For each density mechanism build the full node index, i.e the list of
// compartments with that mechanism, then build the mechanism instance.
- std::vector<size_type> mech_comp_indices(ncomp);
+ std::vector<size_type> mech_cv_index(ncomp);
+ std::vector<value_type> mech_cv_weight(ncomp);
std::map<std::string, std::vector<size_type>> mech_index_map;
- for (auto& mech: mech_map) {
- mech_comp_indices.clear();
- for (auto comp_ival: mech.second) {
- util::append(mech_comp_indices, make_span(comp_ival));
+ for (auto const& mech: mech_map) {
+ // Clear the pre-allocated storage for mechanism indexes and weights.
+ // Reuse the same vectors each time to have only one malloc and free
+ // outside of the loop for each
+ mech_cv_index.clear();
+ mech_cv_weight.clear();
+
+ const auto& seg_cv_ranges = mech.second;
+ for (auto& rng: seg_cv_ranges) {
+ if (rng.has_parent()) {
+ // locate the parent cv in the partially constructed list of cv indexes
+ auto it = algorithms::binary_find(mech_cv_index, rng.parent_cv);
+ if (it == mech_cv_index.end()) {
+ mech_cv_index.push_back(rng.parent_cv);
+ mech_cv_weight.push_back(0);
+ }
+ auto pos = std::distance(std::begin(mech_cv_index), it);
+
+ // add area contribution to the parent cv for the segment
+ mech_cv_weight[pos] += rng.areas.first;
+ }
+ util::append(mech_cv_index, make_span(rng.segment_cvs));
+ util::append(mech_cv_weight, subrange_view(tmp_cv_areas, rng.segment_cvs));
+
+ // adjust the last CV
+ mech_cv_weight.back() = rng.areas.second;
+
+ EXPECTS(mech_cv_weight.size()==mech_cv_index.size());
+ }
+
+ // Scale the weights to get correct units (see w_i^d in formulation docs)
+ // The units for the density channel weights are [10^2 μm^2 = 10^-10 m^2],
+ // which requires that we scale the areas [μm^2] by 10^-2
+ for (auto& w: mech_cv_weight) {
+ w *= 1e-2;
}
mechanisms_.push_back(
- backend::make_mechanism(mech.first, voltage_, current_, mech_comp_indices)
+ backend::make_mechanism(mech.first, voltage_, current_, mech_cv_weight, mech_cv_index)
);
// save the indices for easy lookup later in initialization
- mech_index_map[mech.first] = mech_comp_indices;
+ mech_index_map[mech.first] = mech_cv_index;
}
// Create point (synapse) mechanisms
@@ -541,22 +624,33 @@ void fvm_multicell<Backend>::initialize(
const auto& mech_name = syni.first;
size_type mech_index = mechanisms_.size();
- auto comp_indices = syn_mech_map[syni.second];
- size_type n_indices = size(comp_indices);
+ auto cv_map = syn_mech_map[syni.second];
+ size_type n_indices = size(cv_map);
// sort indices but keep track of their original order for assigning
// target handles
-
using index_pair = std::pair<cell_lid_type, size_type>;
- auto compartment_index = [](index_pair x) { return x.first; };
+ auto cv_index = [](index_pair x) { return x.first; };
auto target_index = [](index_pair x) { return x.second; };
std::vector<index_pair> permute;
assign_by(permute, make_span(0u, n_indices),
- [&](size_type i) { return index_pair(comp_indices[i], i); });
+ [&](size_type i) { return index_pair(cv_map[i], i); });
+
+ // sort the cv information in order of cv index
+ sort_by(permute, cv_index);
- sort_by(permute, compartment_index);
- assign_by(comp_indices, permute, compartment_index);
+ std::vector<cell_lid_type> cv_indices =
+ assign_from(transform_view(permute, cv_index));
+
+ // Create the mechanism.
+ // An empty weight vector is supplied, because there are no weights applied to point
+ // processes, because their currents are calculated with the target units of [nA]
+ mechanisms_.push_back(
+ backend::make_mechanism(mech_name, voltage_, current_, {}, cv_indices));
+
+ // save the compartment indexes for this synapse type
+ mech_index_map[mech_name] = cv_indices;
// make target handles
std::vector<target_handle> handles(n_indices);
@@ -565,13 +659,6 @@ void fvm_multicell<Backend>::initialize(
}
target_hi = std::copy_n(std::begin(handles), n_indices, target_hi);
targets_count += n_indices;
-
- auto mech = backend::make_mechanism(mech_name, voltage_, current_, comp_indices);
- mech->set_areas(cv_areas_);
- mechanisms_.push_back(std::move(mech));
-
- // save the compartment indexes for this synapse type
- mech_index_map[mech_name] = comp_indices;
}
// confirm user-supplied containers for targets are appropriately sized
@@ -653,18 +740,14 @@ void fvm_multicell<Backend>::advance(double dt) {
PL();
}
- // TODO KERNEL: the stimulus might have to become a "proper" mechanism
- // so that the update kernel is fully implemented on GPU.
-
// add current contributions from stimuli
for (auto& stim : stimuli_) {
auto ie = stim.second.amplitude(t_); // [nA]
auto loc = stim.first;
- // note: current_ in [mA/cm^2], ie in [nA], cv_areas_ in [µm^2].
- // unit scale factor: [nA/µm^2]/[mA/cm^2] = 100
+ // note: current_ and ie have the same units [nA]
if (ie!=0.) {
- current_[loc] = current_[loc] - 100*ie/cv_areas_[loc];
+ current_[loc] = current_[loc] - ie;
}
}
PL();
@@ -672,6 +755,7 @@ void fvm_multicell<Backend>::advance(double dt) {
// solve the linear system
PE("matrix", "setup");
matrix_assembler_.assemble(dt);
+
PL(); PE("solve");
matrix_.solve();
PL();
diff --git a/src/mechanism.hpp b/src/mechanism.hpp
index 5ad78a0937a90933a7d7f59d4ee8c3d12562b884..154253dada02b518cdb21eddf2bc6fbb48025430 100644
--- a/src/mechanism.hpp
+++ b/src/mechanism.hpp
@@ -42,7 +42,9 @@ public:
using ion_type = ion<backend>;
mechanism(view vec_v, view vec_i, iarray&& node_index):
- vec_v_(vec_v), vec_i_(vec_i), node_index_(std::move(node_index))
+ vec_v_(vec_v),
+ vec_i_(vec_i),
+ node_index_(std::move(node_index))
{}
std::size_t size() const {
@@ -63,16 +65,11 @@ public:
virtual bool uses_ion(ionKind) const = 0;
virtual void set_ion(ionKind k, ion_type& i, const std::vector<size_type>& index) = 0;
- void set_areas(view area) {
- vec_area_ = area;
- }
-
virtual mechanismKind kind() const = 0;
view vec_v_;
view vec_i_;
iarray node_index_;
- view vec_area_;
};
template <class Backend>
@@ -82,10 +79,11 @@ template <typename M>
auto make_mechanism(
typename M::view vec_v,
typename M::view vec_i,
+ typename M::array&& weights,
typename M::iarray&& node_indices)
--> decltype(util::make_unique<M>(vec_v, vec_i, std::move(node_indices)))
+-> decltype(util::make_unique<M>(vec_v, vec_i, std::move(weights), std::move(node_indices)))
{
- return util::make_unique<M>(vec_v, vec_i, std::move(node_indices));
+ return util::make_unique<M>(vec_v, vec_i, std::move(weights), std::move(node_indices));
}
} // namespace mechanisms
diff --git a/src/util/rangeutil.hpp b/src/util/rangeutil.hpp
index 7f5bf1ac58bb6c65c42aa1c010e39931872982aa..7dee14b911f9e5bb81d8778d12ec141a0e258ffe 100644
--- a/src/util/rangeutil.hpp
+++ b/src/util/rangeutil.hpp
@@ -51,6 +51,18 @@ subrange_view(Seq& seq, Size bi, Size ei) {
return make_range(b, e);
}
+template <
+ typename Seq,
+ typename Iter = typename sequence_traits<Seq>::iterator,
+ typename Size = typename sequence_traits<Seq>::size_type
+>
+enable_if_t<is_forward_iterator<Iter>::value, range<Iter>>
+subrange_view(Seq& seq, std::pair<Size, Size> index) {
+ Iter b = std::next(std::begin(seq), index.first);
+ Iter e = std::next(b, index.second-index.first);
+ return make_range(b, e);
+}
+
// Append sequence to a container
template <typename Container, typename Seq>
@@ -69,6 +81,30 @@ AssignableContainer& assign(AssignableContainer& c, const Seq& seq) {
return c;
}
+namespace impl {
+ template <typename Seq>
+ struct assign_proxy {
+ assign_proxy(const Seq& seq):
+ ref{seq}
+ {}
+
+ // Convert the sequence to a container of type C.
+ // This requires that C supports construction from a pair of iterators
+ template <typename C>
+ operator C() const {
+ return C(std::begin(ref), std::end(ref));
+ }
+
+ const Seq& ref;
+ };
+}
+
+// Copy-assign sequence to a container
+
+template <typename Seq>
+impl::assign_proxy<Seq> assign_from(const Seq& seq) {
+ return impl::assign_proxy<Seq>(seq);
+}
// Assign sequence to a container with transform `proj`
@@ -233,13 +269,6 @@ Value max_value(const Seq& seq, Compare cmp = Compare{}) {
return m;
}
-template <typename T, typename Seq>
-std::vector<T> make_std_vector(const Seq& seq) {
- auto i = std::begin(seq);
- auto e = std::end(seq);
- return std::vector<T>(i, e);
-}
-
template <typename C, typename Seq>
C make_copy(Seq const& seq) {
return C{std::begin(seq), std::end(seq)};
diff --git a/tests/unit/test_algorithms.cpp b/tests/unit/test_algorithms.cpp
index 3f0830a98650c450fcce348eccb6061703a47a3a..a6b66bfec14d5c91b7ae403e7ed2909c74a20346 100644
--- a/tests/unit/test_algorithms.cpp
+++ b/tests/unit/test_algorithms.cpp
@@ -1,3 +1,4 @@
+#include <iterator>
#include <random>
#include <vector>
@@ -578,3 +579,150 @@ TEST(algorithms, index_into)
EXPECT_EQ(i, *it++);
}
}
+
+TEST(algorithms, binary_find)
+{
+ using nest::mc::algorithms::binary_find;
+
+ // empty containers
+ {
+ std::vector<int> v;
+ EXPECT_TRUE(binary_find(v, 100) == std::end(v));
+ }
+
+ // value not present and greater than all entries
+ {
+ int a[] = {1, 10, 15};
+ EXPECT_TRUE(binary_find(a, 100) == std::end(a));
+
+ std::vector<int> v{1, 10, 15};
+ EXPECT_TRUE(binary_find(v, 100) == std::end(v));
+ }
+
+ // value not present and less than all entries
+ {
+ int a[] = {1, 10, 15};
+ EXPECT_TRUE(binary_find(a, -1) == std::end(a));
+
+ std::vector<int> v{1, 10, 15};
+ EXPECT_TRUE(binary_find(v, -1) == std::end(v));
+ }
+
+ // value not present and inside lower-upper bounds
+ {
+ int a[] = {1, 10, 15};
+ EXPECT_TRUE(binary_find(a, 4) == std::end(a));
+
+ std::vector<int> v{1, 10, 15};
+ EXPECT_TRUE(binary_find(v, 4) == std::end(v));
+ }
+
+ // value is first in range
+ {
+ int a[] = {1, 10, 15};
+ auto ita = binary_find(a, 1);
+ auto found = ita!=std::end(a);
+ EXPECT_TRUE(found);
+ EXPECT_EQ(std::distance(std::begin(a), ita), 0u);
+ if (found) EXPECT_EQ(*ita, 1);
+
+ std::vector<int> v{1, 10, 15};
+ auto itv = binary_find(v, 1);
+ found = itv!=std::end(v);
+ EXPECT_TRUE(found);
+ EXPECT_EQ(std::distance(std::begin(v), itv), 0u);
+ if (found) EXPECT_EQ(*itv, 1);
+ }
+
+ // value is last in range
+ {
+ int a[] = {1, 10, 15};
+ auto ita = binary_find(a, 15);
+ auto found = ita!=std::end(a);
+ EXPECT_TRUE(found);
+ EXPECT_EQ(std::distance(std::begin(a), ita), 2u);
+ if (found) EXPECT_EQ(*ita, 15);
+
+ std::vector<int> v{1, 10, 15};
+ auto itv = binary_find(v, 15);
+ found = itv!=std::end(v);
+ EXPECT_TRUE(found);
+ EXPECT_EQ(std::distance(std::begin(v), itv), 2u);
+ if (found) EXPECT_EQ(*itv, 15);
+ }
+
+ // value is last present and neither first nor last in range
+ {
+ int a[] = {1, 10, 15};
+ auto ita = binary_find(a, 10);
+ auto found = ita!=std::end(a);
+ EXPECT_TRUE(found);
+ EXPECT_EQ(std::distance(std::begin(a), ita), 1u);
+ if (found) EXPECT_EQ(*ita, 10);
+
+ std::vector<int> v{1, 10, 15};
+ auto itv = binary_find(v, 10);
+ found = itv!=std::end(v);
+ EXPECT_TRUE(found);
+ EXPECT_EQ(std::distance(std::begin(v), itv), 1u);
+ if (found) EXPECT_EQ(*itv, 10);
+ }
+
+ // value is last present and neither first nor last in range and range has even size
+ {
+ int a[] = {1, 10, 15, 27};
+ auto ita = binary_find(a, 10);
+ auto found = ita!=std::end(a);
+ EXPECT_TRUE(found);
+ EXPECT_EQ(std::distance(std::begin(a), ita), 1u);
+ if (found) EXPECT_EQ(*ita, 10);
+
+ std::vector<int> v{1, 10, 15, 27};
+ auto itv = binary_find(v, 10);
+ found = itv!=std::end(v);
+ EXPECT_TRUE(found);
+ EXPECT_EQ(std::distance(std::begin(v), itv), 1u);
+ if (found) EXPECT_EQ(*itv, 10);
+ }
+
+ // test for const types
+ // i.e. iterators returned from passing in a const reference to a container
+ // can be compared to a const iterator from the container
+ {
+ std::vector<int> v{1, 10, 15};
+ auto const& vr = v;
+ auto itv = binary_find(vr, 10);
+ auto found = itv!=std::end(vr);
+ EXPECT_TRUE(found);
+ EXPECT_EQ(std::distance(nest::mc::util::cbegin(v), itv), 1u);
+ if (found) EXPECT_EQ(*itv, 10);
+ }
+}
+
+struct int_string {
+ int value;
+
+ friend bool operator<(const int_string& lhs, const std::string& rhs) {
+ return lhs.value<std::stoi(rhs);
+ }
+ friend bool operator<(const std::string& lhs, const int_string& rhs) {
+ return std::stoi(lhs)<rhs.value;
+ }
+ friend bool operator==(const int_string& lhs, const std::string& rhs) {
+ return lhs.value==std::stoi(rhs);
+ }
+ friend bool operator==(const std::string& lhs, const int_string& rhs) {
+ return std::stoi(lhs)==rhs.value;
+ }
+};
+
+TEST(algorithms, binary_find_convert)
+{
+ using nest::mc::algorithms::binary_find;
+
+ std::vector<std::string> values = {"0", "10", "20", "30"};
+ auto it = nest::mc::algorithms::binary_find(values, int_string{20});
+
+ EXPECT_TRUE(it!=values.end());
+ EXPECT_TRUE(std::distance(values.begin(), it)==2u);
+}
diff --git a/tests/unit/test_mechanisms.cpp b/tests/unit/test_mechanisms.cpp
index 58e8a97419f481490bf63dee656128f5592ff62e..fbc0989c66a25098ad816807cf3991e06707cf95 100644
--- a/tests/unit/test_mechanisms.cpp
+++ b/tests/unit/test_mechanisms.cpp
@@ -6,6 +6,7 @@
TEST(mechanisms, helpers) {
using namespace nest::mc;
using size_type = multicore::backend::size_type;
+ using value_type = multicore::backend::value_type;
// verify that the hh and pas channels are available
EXPECT_TRUE(multicore::backend::has_mechanism("hh"));
@@ -13,21 +14,21 @@ TEST(mechanisms, helpers) {
std::vector<size_type> parent_index = {0,0,1,2,3,4,0,6,7,8};
auto node_indices = std::vector<size_type>{0,6,7,8,9};
+ auto weights = std::vector<value_type>(node_indices.size(), 1.0);
auto n = node_indices.size();
multicore::backend::array vec_i(n, 0.);
multicore::backend::array vec_v(n, 0.);
auto mech = multicore::backend::make_mechanism(
- "hh", memory::make_view(vec_v), memory::make_view(vec_i), node_indices);
+ "hh", memory::make_view(vec_v), memory::make_view(vec_i), weights, node_indices);
EXPECT_EQ(mech->name(), "hh");
EXPECT_EQ(mech->size(), 5u);
// check that an out_of_range exception is thrown if an invalid mechanism is requested
ASSERT_THROW(
- multicore::backend::make_mechanism("dachshund", vec_v, vec_i, node_indices),
+ multicore::backend::make_mechanism("dachshund", vec_v, vec_i, weights, node_indices),
std::out_of_range
);
- //0 1 2 3 4 5 6 7 8 9
}
diff --git a/tests/unit/test_range.cpp b/tests/unit/test_range.cpp
index cc73f2750905caaaff1489adb660b1b57dc5590a..d49430ae03be15c510927b1bc387916b4244c962 100644
--- a/tests/unit/test_range.cpp
+++ b/tests/unit/test_range.cpp
@@ -325,6 +325,25 @@ TEST(range, assign) {
EXPECT_EQ("00110", text);
}
+TEST(range, assign_from) {
+ int in[] = {0,1,2};
+
+ {
+ std::vector<int> copy = util::assign_from(in);
+ for (auto i=0u; i<util::size(in); ++i) {
+ EXPECT_EQ(in[i], copy[i]);
+ }
+ }
+
+ {
+ std::vector<int> copy = util::assign_from(
+ util::transform_view(in, [](int i) {return 2*i;}));
+ for (auto i=0u; i<util::size(in); ++i) {
+ EXPECT_EQ(2*in[i], copy[i]);
+ }
+ }
+}
+
TEST(range, sort) {
char cstr[] = "howdy";
diff --git a/tests/unit/test_synapses.cpp b/tests/unit/test_synapses.cpp
index 50deb528ec5f6206b38fe9794beb181c9723cceb..0084ed56b281707ff470a8581c58407840e0588b 100644
--- a/tests/unit/test_synapses.cpp
+++ b/tests/unit/test_synapses.cpp
@@ -46,14 +46,16 @@ TEST(synapses, expsyn_basic_state)
{
using namespace nest::mc;
using size_type = multicore::backend::size_type;
+ using value_type = multicore::backend::value_type;
using synapse_type = mechanisms::expsyn::mechanism_expsyn<multicore::backend>;
auto num_syn = 4;
std::vector<size_type> indexes(num_syn);
+ std::vector<value_type> weights(indexes.size(), 1.0);
synapse_type::array voltage(num_syn, -65.0);
synapse_type::array current(num_syn, 1.0);
- auto mech = mechanisms::make_mechanism<synapse_type>( voltage, current, indexes );
+ auto mech = mechanisms::make_mechanism<synapse_type>(voltage, current, weights, indexes);
auto ptr = dynamic_cast<synapse_type*>(mech.get());
@@ -102,11 +104,13 @@ TEST(synapses, exp2syn_basic_state)
auto num_syn = 4;
using size_type = multicore::backend::size_type;
+ using value_type = multicore::backend::value_type;
std::vector<size_type> indexes(num_syn);
+ std::vector<value_type> weights(indexes.size(), 1.0);
synapse_type::array voltage(num_syn, -65.0);
synapse_type::array current(num_syn, 1.0);
- auto mech = mechanisms::make_mechanism<synapse_type>( voltage, current, indexes );
+ auto mech = mechanisms::make_mechanism<synapse_type>(voltage, current, weights, indexes);
auto ptr = dynamic_cast<synapse_type*>(mech.get());