From 6871a9c7b6a9a2dd9613f4f2ced6d586d6901c54 Mon Sep 17 00:00:00 2001 From: bcumming <bcumming@cscs.ch> Date: Thu, 28 Apr 2016 17:51:16 +0200 Subject: [PATCH] big push towards full cell integration --- .ycm_extra_conf.py | 2 + CMakeLists.txt | 1 + data/hh.mod | 78 ++++++++-------- data/passive.mod | 28 ++++++ mechanisms/generate.sh | 4 + mechanisms/mod/hh.mod | 124 +++++++++++++++++++++++++ mechanisms/mod/pas.mod | 28 ++++++ src/CMakeLists.txt | 5 +- src/algorithms.hpp | 53 +++++++++++ src/fvm.hpp | 175 +++++++++++++++++++++++++++++------- src/indexed_view.hpp | 50 +++++++++++ src/ion.hpp | 109 ++++++++++++++++++++++ src/matrix.hpp | 16 +++- src/mechanism.hpp | 31 +++++-- src/mechanism_interface.cpp | 34 +++++-- src/mechanism_interface.hpp | 43 +++++++-- src/parameter_list.hpp | 29 +++++- src/segment.hpp | 25 ++++++ src/util.hpp | 44 +++++++++ tests/CMakeLists.txt | 1 + tests/test_algorithms.cpp | 119 ++++++++++++++++++++++++ tests/test_matrix.cpp | 1 - tests/test_mechanisms.cpp | 34 +++++++ tests/test_run.cpp | 34 +++---- vector | 2 +- 25 files changed, 957 insertions(+), 113 deletions(-) mode change 100755 => 100644 data/hh.mod create mode 100644 data/passive.mod create mode 100755 mechanisms/generate.sh create mode 100644 mechanisms/mod/hh.mod create mode 100644 mechanisms/mod/pas.mod create mode 100644 src/indexed_view.hpp create mode 100644 src/ion.hpp create mode 100644 tests/test_mechanisms.cpp diff --git a/.ycm_extra_conf.py b/.ycm_extra_conf.py index b8378fb8..55d6af4d 100644 --- a/.ycm_extra_conf.py +++ b/.ycm_extra_conf.py @@ -43,6 +43,8 @@ flags = [ 'src', '-I', '.', + '-I', + 'include', # '-isystem', # '/Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX10.10.sdk/usr/include/c++/4.2.1', # '-I', diff --git a/CMakeLists.txt b/CMakeLists.txt index 61278f34..57991c88 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -14,6 +14,7 @@ set(CMAKE_EXPORT_COMPILE_COMMANDS "YES") set(CMAKE_LIBRARY_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/lib) set(CMAKE_ARCHIVE_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/lib) +include_directories(${CMAKE_SOURCE_DIR}/include) include_directories(${CMAKE_SOURCE_DIR}/src) include_directories(${CMAKE_SOURCE_DIR}) diff --git a/data/hh.mod b/data/hh.mod old mode 100755 new mode 100644 index 6020be0c..975d601d --- a/data/hh.mod +++ b/data/hh.mod @@ -1,17 +1,17 @@ TITLE hh.mod squid sodium, potassium, and leak channels -COMMENT - This is the original Hodgkin-Huxley treatment for the set of sodium, - potassium, and leakage channels found in the squid giant axon membrane. - ("A quantitative description of membrane current and its application - conduction and excitation in nerve" J.Physiol. (Lond.) 117:500-544 (1952).) - Membrane voltage is in absolute mV and has been reversed in polarity - from the original HH convention and shifted to reflect a resting potential - of -65 mV. - Remember to set celsius=6.3 (or whatever) in your HOC file. - See squid.hoc for an example of a simulation using this model. - SW Jaslove 6 March, 1992 -ENDCOMMENT +: COMMENT +: This is the original Hodgkin-Huxley treatment for the set of sodium, +: potassium, and leakage channels found in the squid giant axon membrane. +: ("A quantitative description of membrane current and its application +: conduction and excitation in nerve" J.Physiol. (Lond.) 117:500-544 (1952).) +: Membrane voltage is in absolute mV and has been reversed in polarity +: from the original HH convention and shifted to reflect a resting potential +: of -65 mV. +: Remember to set celsius=6.3 (or whatever) in your HOC file. +: See squid.hoc for an example of a simulation using this model. +: SW Jaslove 6 March, 1992 +: ENDCOMMENT UNITS { (mA) = (milliamp) @@ -26,14 +26,14 @@ NEURON { NONSPECIFIC_CURRENT il RANGE gnabar, gkbar, gl, el, gna, gk GLOBAL minf, hinf, ninf, mtau, htau, ntau - THREADSAFE : assigned GLOBALs will be per thread } PARAMETER { - gnabar = .12 (S/cm2) <0,1e9> - gkbar = .036 (S/cm2) <0,1e9> - gl = .0003 (S/cm2) <0,1e9> + gnabar = .12 (S/cm2) : <0,1e9> + gkbar = .036 (S/cm2) : <0,1e9> + gl = .0003 (S/cm2) : <0,1e9> el = -54.3 (mV) + celsius = 37 (degC) } STATE { @@ -42,20 +42,17 @@ STATE { ASSIGNED { v (mV) - celsius (degC) - ena (mV) - ek (mV) gna (S/cm2) gk (S/cm2) - ina (mA/cm2) - ik (mA/cm2) - il (mA/cm2) - minf hinf ninf - mtau (ms) htau (ms) ntau (ms) + minf + hinf + ninf + mtau (ms) + htau (ms) + ntau (ms) } -? currents BREAKPOINT { SOLVE states METHOD cnexp gna = gnabar*m*m*m*h @@ -72,18 +69,15 @@ INITIAL { n = ninf } -DERIVATIVE states -{ +DERIVATIVE states { rates(v) m' = (minf-m)/mtau h' = (hinf-h)/htau n' = (ninf-n)/ntau } -PROCEDURE rates(v(mV)) +PROCEDURE rates(v) :(mV)) { - :Computes rate and other constants at current v. - :Call once from HOC to initialize inf at resting v. LOCAL alpha, beta, sum, q10 q10 = 3^((celsius - 6.3)/10) @@ -103,18 +97,28 @@ PROCEDURE rates(v(mV)) hinf = alpha/sum :"n" potassium activation system - alpha = .01*vtrap(-(v+55),10) + alpha = .01*vtrap(-(v+55),10) beta = .125*exp(-(v+65)/80) sum = alpha + beta ntau = 1/(q10*sum) ninf = alpha/sum } -FUNCTION vtrap(x,y) { :Traps for 0 in denominator of rate eqns. - if (fabs(x/y) < 1e-6) { - vtrap = y*(1 - x/y/2) - }else{ - vtrap = x/(exp(x/y) - 1) - } +FUNCTION vtrap(x,y) { + : Traps for 0 in denominator of rate eqns. + + : Disabled for function inlining... + : there is a very good case for making vtrap a builtin function for the language + : - it is ubiquitous: used in so many models of ion channels + : - platform specific optimizations written in assembly/intrinsics + : required to facilitate vectorization/minimize branch misses + + :if (fabs(x/y) < 1e-6) { + : vtrap = y*(1 - x/y/2) + :}else{ + : vtrap = x/(exp(x/y) - 1) + :} + + vtrap = x/(exp(x/y) - 1) } diff --git a/data/passive.mod b/data/passive.mod new file mode 100644 index 00000000..413bd250 --- /dev/null +++ b/data/passive.mod @@ -0,0 +1,28 @@ +TITLE passive membrane channel + +UNITS { + (mV) = (millivolt) + (mA) = (milliamp) + (S) = (siemens) +} + +NEURON { + SUFFIX pas + NONSPECIFIC_CURRENT i + RANGE g, e +} + +INITIAL {} + +PARAMETER { + g = .001 (S/cm2) :<0,1e9> + e = -70 (mV) +} + +ASSIGNED { + v (mV) +} + +BREAKPOINT { + i = g*(v - e) +} diff --git a/mechanisms/generate.sh b/mechanisms/generate.sh new file mode 100755 index 00000000..9fe83f04 --- /dev/null +++ b/mechanisms/generate.sh @@ -0,0 +1,4 @@ +for mech in pas hh +do + modcc -t cpu -o ../include/mechanisms/$mech.hpp ./mod/$mech.mod +done diff --git a/mechanisms/mod/hh.mod b/mechanisms/mod/hh.mod new file mode 100644 index 00000000..975d601d --- /dev/null +++ b/mechanisms/mod/hh.mod @@ -0,0 +1,124 @@ +TITLE hh.mod squid sodium, potassium, and leak channels + +: COMMENT +: This is the original Hodgkin-Huxley treatment for the set of sodium, +: potassium, and leakage channels found in the squid giant axon membrane. +: ("A quantitative description of membrane current and its application +: conduction and excitation in nerve" J.Physiol. (Lond.) 117:500-544 (1952).) +: Membrane voltage is in absolute mV and has been reversed in polarity +: from the original HH convention and shifted to reflect a resting potential +: of -65 mV. +: Remember to set celsius=6.3 (or whatever) in your HOC file. +: See squid.hoc for an example of a simulation using this model. +: SW Jaslove 6 March, 1992 +: ENDCOMMENT + +UNITS { + (mA) = (milliamp) + (mV) = (millivolt) + (S) = (siemens) +} + +NEURON { + SUFFIX hh + USEION na READ ena WRITE ina + USEION k READ ek WRITE ik + NONSPECIFIC_CURRENT il + RANGE gnabar, gkbar, gl, el, gna, gk + GLOBAL minf, hinf, ninf, mtau, htau, ntau +} + +PARAMETER { + gnabar = .12 (S/cm2) : <0,1e9> + gkbar = .036 (S/cm2) : <0,1e9> + gl = .0003 (S/cm2) : <0,1e9> + el = -54.3 (mV) + celsius = 37 (degC) +} + +STATE { + m h n +} + +ASSIGNED { + v (mV) + + gna (S/cm2) + gk (S/cm2) + minf + hinf + ninf + mtau (ms) + htau (ms) + ntau (ms) +} + +BREAKPOINT { + SOLVE states METHOD cnexp + gna = gnabar*m*m*m*h + ina = gna*(v - ena) + gk = gkbar*n*n*n*n + ik = gk*(v - ek) + il = gl*(v - el) +} + +INITIAL { + rates(v) + m = minf + h = hinf + n = ninf +} + +DERIVATIVE states { + rates(v) + m' = (minf-m)/mtau + h' = (hinf-h)/htau + n' = (ninf-n)/ntau +} + +PROCEDURE rates(v) :(mV)) +{ + LOCAL alpha, beta, sum, q10 + + q10 = 3^((celsius - 6.3)/10) + + :"m" sodium activation system + alpha = .1 * vtrap(-(v+40),10) + beta = 4 * exp(-(v+65)/18) + sum = alpha + beta + mtau = 1/(q10*sum) + minf = alpha/sum + + :"h" sodium inactivation system + alpha = .07 * exp(-(v+65)/20) + beta = 1 / (exp(-(v+35)/10) + 1) + sum = alpha + beta + htau = 1/(q10*sum) + hinf = alpha/sum + + :"n" potassium activation system + alpha = .01*vtrap(-(v+55),10) + beta = .125*exp(-(v+65)/80) + sum = alpha + beta + ntau = 1/(q10*sum) + ninf = alpha/sum +} + +FUNCTION vtrap(x,y) { + : Traps for 0 in denominator of rate eqns. + + : Disabled for function inlining... + : there is a very good case for making vtrap a builtin function for the language + : - it is ubiquitous: used in so many models of ion channels + : - platform specific optimizations written in assembly/intrinsics + : required to facilitate vectorization/minimize branch misses + + :if (fabs(x/y) < 1e-6) { + : vtrap = y*(1 - x/y/2) + :}else{ + : vtrap = x/(exp(x/y) - 1) + :} + + vtrap = x/(exp(x/y) - 1) +} + diff --git a/mechanisms/mod/pas.mod b/mechanisms/mod/pas.mod new file mode 100644 index 00000000..413bd250 --- /dev/null +++ b/mechanisms/mod/pas.mod @@ -0,0 +1,28 @@ +TITLE passive membrane channel + +UNITS { + (mV) = (millivolt) + (mA) = (milliamp) + (S) = (siemens) +} + +NEURON { + SUFFIX pas + NONSPECIFIC_CURRENT i + RANGE g, e +} + +INITIAL {} + +PARAMETER { + g = .001 (S/cm2) :<0,1e9> + e = -70 (mV) +} + +ASSIGNED { + v (mV) +} + +BREAKPOINT { + i = g*(v - e) +} diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 5696287e..8a8286bd 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1,10 +1,11 @@ set(HEADERS swcio.hpp - ) +) set(BASE_SOURCES + cell.cpp + mechanism_interface.cpp parameter_list.cpp swcio.cpp - cell.cpp ) add_library(cellalgo ${BASE_SOURCES} ${HEADERS}) diff --git a/src/algorithms.hpp b/src/algorithms.hpp index bc25b98b..2a549eee 100644 --- a/src/algorithms.hpp +++ b/src/algorithms.hpp @@ -6,6 +6,8 @@ #include <numeric> #include <type_traits> +#include "util.hpp" + /* * Some simple wrappers around stl algorithms to improve readability of code * that uses them. @@ -117,6 +119,57 @@ namespace algorithms{ return true; } + template<typename C> + bool is_sorted(const C& c) + { + return std::is_sorted(c.begin(), c.end()); + } + + template<typename C> + bool is_unique(const C& c) + { + return std::adjacent_find(c.begin(), c.end()) == c.end(); + } + + /// Return and index that maps entries in sub to their corresponding + /// values in super, where sub is a subset of super. + /// + /// Both sets are sorted and have unique entries. + /// Complexity is O(n), where n is size of super + template<typename C> + // C::iterator models forward_iterator + // C::value_type is_integral + C index_into(const C& super, const C& sub) + { + //EXPECTS {s \in super : \forall s \in sub}; + EXPECTS(is_unique(super) && is_unique(sub)); + EXPECTS(is_sorted(super) && is_sorted(sub)); + EXPECTS(sub.size() <= super.size()); + + static_assert( + std::is_integral<typename C::value_type>::value, + "index_into only applies to integral types" + ); + + C out(sub.size()); // out will have one entry for each index in sub + + auto sub_it=sub.begin(); + auto super_it=super.begin(); + auto sub_idx=0u, super_idx = 0u; + + while(sub_it!=sub.end() && super_it!=super.end()) { + if(*sub_it==*super_it) { + out[sub_idx] = super_idx; + ++sub_it; ++sub_idx; + } + ++super_it; ++super_idx; + } + + EXPECTS(sub_idx==sub.size()); + + return out; + } + } // namespace algorithms } // namespace mc } // namespace nest diff --git a/src/fvm.hpp b/src/fvm.hpp index b0eedb41..c6c458c7 100644 --- a/src/fvm.hpp +++ b/src/fvm.hpp @@ -1,13 +1,19 @@ #pragma once #include <algorithm> +#include <map> +#include <string> +#include <vector> +#include <algorithms.hpp> #include <cell.hpp> -#include <segment.hpp> +#include <ion.hpp> #include <math.hpp> #include <matrix.hpp> +#include <mechanism.hpp> +#include <mechanism_interface.hpp> #include <util.hpp> -#include <algorithms.hpp> +#include <segment.hpp> #include <vector/include/Vector.hpp> @@ -24,8 +30,16 @@ class fvm_cell { /// the integral index type using size_type = I; + /// the type used to store matrix information using matrix_type = matrix<value_type, size_type>; + /// mechanism type + using mechanism_type = + nest::mc::mechanisms::mechanism_ptr<value_type, size_type>; + + /// ion species storage + using ion_type = mechanisms::ion<value_type, size_type>; + /// the container used for indexes using index_type = memory::HostVector<size_type>; /// view into index container @@ -44,26 +58,36 @@ class fvm_cell { /// TODO this should be const /// which requires const_view in the vector library - matrix_type& jacobian(); + matrix_type& jacobian() { + return matrix_; + } /// 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(); + vector_view cv_areas() { + return 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(); + /// this is the total capacitance, not per unit area, + /// i.e. equivalent to sigma_i * c_m + vector_view cv_capacitance() { + return cv_capacitance_; + } - std::size_t size() const - { + std::size_t size() const { return matrix_.size(); } + /// return reference to in iterable container of the mechanisms + std::vector<mechanism_type>& mechanisms() { + return mechanisms_; + } + private: /// the linear system for implicit time stepping of cell state @@ -89,11 +113,16 @@ class fvm_cell { /// the potential in mV in each CV vector_type voltage_; + /// the set of mechanisms present in the cell + std::vector<mechanism_type> mechanisms_; + + /// the ion species + std::map<mechanisms::ionKind, ion_type> ions_; }; -//////////////////////////////////////////////////////////// -// Implementation -//////////////////////////////////////////////////////////// +//////////////////////////////////////////////////////////////////////////////// +//////////////////////////////// Implementation //////////////////////////////// +//////////////////////////////////////////////////////////////////////////////// template <typename T, typename I> fvm_cell<T, I>::fvm_cell(nest::mc::cell const& cell) @@ -106,8 +135,16 @@ fvm_cell<T, I>::fvm_cell(nest::mc::cell const& cell) using util::left; using util::right; + // TODO: potential code stink + // matrix_ is a member, but it is not initialized with the other members + // above because it requires the parent_index, which is calculated + // "on the fly" by cell.model(). + // cell.model() is quite expensive, and the information it calculates is + // used elsewhere, so we defer the intialization to inside the constructor + // body. const auto graph = cell.model(); matrix_ = matrix_type(graph.parent_index); + auto parent_index = matrix_.p(); auto const& segment_index = graph.segment_index; @@ -170,27 +207,102 @@ fvm_cell<T, I>::fvm_cell(nest::mc::cell const& cell) 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_; -} + ///////////////////////////////////////////// + // create mechanisms + ///////////////////////////////////////////// -template <typename T, typename I> -typename fvm_cell<T,I>::vector_view -fvm_cell<T,I>::cv_areas() -{ - return cv_areas_; -} + // FIXME : candidate for a private member function -template <typename T, typename I> -typename fvm_cell<T,I>::vector_view -fvm_cell<T,I>::cv_capacitance() -{ - return cv_capacitance_; + // for each mechanism in the cell record the indexes of the segments that + // contain the mechanism + std::map<std::string, std::vector<int>> mech_map; + + for(auto i=0; i<cell.num_segments(); ++i) { + for(const auto& mech : cell.segment(i)->mechanisms()) { + // FIXME : Membrane has to be a proper mechanism, + // because it is exposed via the public interface. + // This if statement is bad + if(mech.name() != "membrane") { + mech_map[mech.name()].push_back(i); + } + } + } + + // Create the mechanism implementations with the state for each mechanism + // instance. + // TODO : this works well for density mechanisms (e.g. ion channels), but + // does it work for point processes (e.g. synapses)? + for(auto& mech : mech_map) { + auto& helper = nest::mc::mechanisms::get_mechanism_helper(mech.first); + + // calculate the number of compartments that contain the mechanism + auto num_comp = 0u; + for(auto seg : mech.second) { + num_comp += segment_index[seg+1] - segment_index[seg]; + } + + // build a vector of the indexes of the compartments that contain + // the mechanism + std::vector<int> compartment_index(num_comp); + auto pos = 0u; + for(auto seg : mech.second) { + auto seg_size = segment_index[seg+1] - segment_index[seg]; + std::iota( + compartment_index.data() + pos, + compartment_index.data() + pos + seg_size, + segment_index[seg] + ); + pos += seg_size; + } + + // instantiate the mechanism + mechanisms_.push_back( + helper->new_mechanism( + &matrix_, + index_view(compartment_index.data(), compartment_index.size()) + ) + ); + } + + ///////////////////////////////////////////// + // build the ion species + // FIXME : this should be a private member function + ///////////////////////////////////////////// + for(auto ion : mechanisms::ion_kinds()) { + auto indexes = + std::make_pair(std::vector<int>(size()), std::vector<int>(size())); + auto ends = + std::make_pair(indexes.first.begin(), indexes.second.begin()); + + // after the loop the range + // [indexes.first.begin(), ends.first) + // will hold the indexes of the compartments that require ion + for(auto& mech : mechanisms_) { + if(mech->uses_ion(ion)) { + ends.second = + std::set_union( + mech->node_index().begin(), mech->node_index().end(), + indexes.first.begin(), ends.first, + indexes.second.begin() + ); + std::swap(indexes.first, indexes.second); + std::swap(ends.first, ends.second); + } + } + + // create the ion state + if(auto n = std::distance(indexes.first.begin(), ends.first)) { + ions_.emplace(ion, index_view(indexes.first.data(), n)); + } + + // join the ion reference in each mechanism into the cell-wide ion state + for(auto& mech : mechanisms_) { + if(mech->uses_ion(ion)) { + mech->set_ion(ion, ions_[ion]); + } + } + } } template <typename T, typename I> @@ -236,6 +348,9 @@ void fvm_cell<T, I>::setup_matrix(T dt) } } + } // namespace fvm } // namespace mc } // namespace nest + + diff --git a/src/indexed_view.hpp b/src/indexed_view.hpp new file mode 100644 index 00000000..0cbc22a9 --- /dev/null +++ b/src/indexed_view.hpp @@ -0,0 +1,50 @@ +#pragma once + +#include <vector/include/Vector.hpp> + +namespace nest { +namespace mc { + +template <typename T, typename I> +struct indexed_view { + using value_type = T; + using size_type = I; + using view_type = typename memory::HostVector<T>::view_type; + using index_view_type = typename memory::HostVector<I>::view_type;; + + view_type view; + index_view_type index; // TODO make this a const view + + indexed_view(view_type v, index_view_type i) + : view(v), + index(i) + {} + + size_type size() const + { + return index.size(); + } + + // TODO + // + // these should either + // - return the internal reference type implementations from the containers + // e.g. the GPU reference that does a copy-on-read from GPU memory + // - or ensure that the result of dereferencing these is properly handled + // i.e. by just returning a value for the const version + + value_type& + operator[] (const size_type i) + { + return view[index[i]]; + } + + value_type const& + operator[] (const size_type i) const + { + return view[index[i]]; + } +}; + +} // namespace mc +} // namespace nest diff --git a/src/ion.hpp b/src/ion.hpp new file mode 100644 index 00000000..2f780616 --- /dev/null +++ b/src/ion.hpp @@ -0,0 +1,109 @@ +#pragma once + +#include <vector/include/Vector.hpp> +#include "indexed_view.hpp" + +namespace nest { +namespace mc { +namespace mechanisms { + +/******************************************************************************* + + ion channels have the following fields : + + --------------------------------------------------- + label Ca Na K name + --------------------------------------------------- + iX ica ina ik current + eX eca ena ek reversal_potential + Xi cai nai ki internal_concentration + Xo cao nao ko external_concentration + gX gca gna gk conductance + --------------------------------------------------- +*******************************************************************************/ + +/// let's enumerate the ion channel types +enum class ionKind {ca, na, k}; + +[[gnu::unused]] static +std::string to_string(ionKind k) +{ + if(k==ionKind::na) return "sodium"; + if(k==ionKind::ca) return "calcium"; + if(k==ionKind::k) return "pottasium"; + return "unkown ion"; +} + +/// and a little helper to iterate over them +[[gnu::unused]] static +std::initializer_list<ionKind> ion_kinds() +{ + return {ionKind::ca, ionKind::na, ionKind::k}; +} + +/// storage for ion channel information in a cell group +template<typename T, typename I> +class ion { +public : + // expose tempalte parameters + using value_type = T; + using size_type = I; + + // define storage types + using vector_type = memory::HostVector<value_type>; + using index_type = memory::HostVector<size_type>; + using vector_view_type = typename vector_type::view_type; + using index_view_type = typename index_type::view_type; + + using indexed_view_type = indexed_view<value_type, size_type>; + + ion() = default; + + ion(index_view_type idx) + : node_index_{idx} + , iX_{idx.size(), std::numeric_limits<value_type>::quiet_NaN()} + , eX_{idx.size(), std::numeric_limits<value_type>::quiet_NaN()} + , Xi_{idx.size(), std::numeric_limits<value_type>::quiet_NaN()} + , Xo_{idx.size(), std::numeric_limits<value_type>::quiet_NaN()} + { } + + std::size_t memory() const { + return 4u*size() * sizeof(value_type) + + size() * sizeof(index_type) + + sizeof(ion); + } + + vector_view_type current() { + return iX_; + } + vector_view_type reversal_potential() { + return eX_; + } + vector_view_type internal_concentration() { + return Xi_; + } + vector_view_type external_concentration() { + return Xo_; + } + + index_type node_index() { + return node_index_; + } + + size_t size() const { + return node_index_.size(); + } + +private : + + index_type node_index_; + vector_type iX_; + vector_type eX_; + vector_type Xi_; + vector_type Xo_; +}; + +} // namespace mechanisms +} // namespace mc +} // namespace nest + diff --git a/src/matrix.hpp b/src/matrix.hpp index 16936fa4..e6d406ee 100644 --- a/src/matrix.hpp +++ b/src/matrix.hpp @@ -67,7 +67,7 @@ class matrix { /// the total memory used to store the matrix std::size_t memory() const { - auto s = 5 * (sizeof(value_type) * size() + sizeof(vector_type)); + auto s = 6 * (sizeof(value_type) * size() + sizeof(vector_type)); s += sizeof(size_type) * (parent_index_.size() + cell_index_.size()) + 2*sizeof(index_type); s += sizeof(matrix); @@ -80,6 +80,11 @@ class matrix { return cell_index_.size() - 1; } + /// FIXME : make modparser use the correct accessors (l,d,u,rhs) instead of these + view_type vec_rhs() { return rhs_; } + view_type vec_d() { return d_; } + view_type vec_v() { return v_; } + /// the vector holding the lower part of the matrix view_type l() { @@ -104,6 +109,12 @@ class matrix { return rhs_; } + /// the vector holding the solution (voltage) + view_type v() + { + return v_; + } + /// the vector holding the parent index index_view p() { @@ -152,6 +163,7 @@ class matrix { d_ = vector_type(n, default_value); u_ = vector_type(n, default_value); rhs_ = vector_type(n, default_value); + v_ = vector_type(n, default_value); } /// the parent indice that describe matrix structure @@ -164,8 +176,10 @@ class matrix { vector_type l_; vector_type d_; vector_type u_; + /// after calling solve, the solution is stored in rhs_ vector_type rhs_; + vector_type v_; // should this be in the matrix? not so sure... }; } // namespace nest diff --git a/src/mechanism.hpp b/src/mechanism.hpp index 34322258..f43abc16 100644 --- a/src/mechanism.hpp +++ b/src/mechanism.hpp @@ -5,17 +5,23 @@ #include <memory> #include <string> +#include "indexed_view.hpp" #include "matrix.hpp" #include "parameter_list.hpp" #include "util.hpp" +#include "ion.hpp" namespace nest { namespace mc { +namespace mechanisms { enum class mechanismKind {point, density}; template <typename T, typename I> class mechanism { + +public: + using value_type = T; using size_type = I; @@ -24,31 +30,39 @@ class mechanism { using view_type = typename vector_type::view_type; using index_type = memory::HostVector<size_type>; using index_view = typename index_type::view_type; - //using indexed_view= IndexedView<value_type, size_type, target()>; + using indexed_view_type = indexed_view<value_type, size_type>; using matrix_type = matrix<value_type, size_type>; + using ion_type = ion<value_type, size_type>; - mechanism(matrix_type *matrix, index_view node_indices) + mechanism(matrix_type *matrix, index_view node_index) : matrix_(matrix) - , node_indices_(node_indices) + , node_index_(node_index) {} std::size_t size() const { - return node_indices_.size(); + return node_index_.size(); + } + + index_view node_index() const + { + return node_index_; } virtual void set_params(value_type t_, value_type dt_) = 0; virtual std::string name() const = 0; virtual std::size_t memory() const = 0; - virtual void init() = 0; - virtual void state() = 0; - virtual void current() = 0; + virtual void nrn_init() = 0; + virtual void nrn_state() = 0; + virtual void nrn_current() = 0; + virtual bool uses_ion(ionKind) const = 0; + virtual void set_ion(ionKind k, ion_type& i) = 0; virtual mechanismKind kind() const = 0; matrix_type* matrix_; - index_type node_indices_; + index_type node_index_; }; template <typename T, typename I> @@ -63,5 +77,6 @@ make_mechanism( return util::make_unique<M>(matrix, node_indices); } +} // namespace mechanisms } // namespace nest } // namespace mc diff --git a/src/mechanism_interface.cpp b/src/mechanism_interface.cpp index ba24eb06..10fad45f 100644 --- a/src/mechanism_interface.cpp +++ b/src/mechanism_interface.cpp @@ -1,21 +1,41 @@ #include "mechanism_interface.hpp" -/* include the mechanisms +// +// include the mechanisms +// + #include <mechanisms/hh.hpp> #include <mechanisms/pas.hpp> -*/ namespace nest { namespace mc { namespace mechanisms { -std::map<std::string, mechanism_helper<double, int>> mechanism_map; +std::map<std::string, mechanism_helper_ptr<value_type, index_type>> mechanism_helpers; void setup_mechanism_helpers() { - /* manually insert - mechanism_map["hh"] = mechanisms::hh; - mechanism_map["pas"] = mechanisms::pas; - */ + mechanism_helpers["pas"] = + make_mechanism_helper< + mechanisms::pas::helper<value_type, index_type> + >(); + + mechanism_helpers["hh"] = + make_mechanism_helper< + mechanisms::hh::helper<value_type, index_type> + >(); +} + +mechanism_helper_ptr<value_type, index_type>& +get_mechanism_helper(const std::string& name) +{ + auto helper = mechanism_helpers.find(name); + if(helper==mechanism_helpers.end()) { + throw std::out_of_range( + nest::mc::util::pprintf("there is no mechanism named \'%\'", name) + ); + } + + return helper->second; } } // namespace mechanisms diff --git a/src/mechanism_interface.hpp b/src/mechanism_interface.hpp index 389ef059..33b6e7f2 100644 --- a/src/mechanism_interface.hpp +++ b/src/mechanism_interface.hpp @@ -9,24 +9,51 @@ namespace nest { namespace mc { +namespace mechanisms { -/// +using value_type = double; +using index_type = int; + +/// helper type for building mechanisms +/// the use of abstract base classes everywhere is a bit ugly template <typename T, typename I> struct mechanism_helper { - using index_type = memory::HostView<I>; + using value_type = T; + using size_type = I; + using index_type = memory::HostVector<I>; using index_view = typename index_type::view_type; - using mechanism_type = mechanism<T, I>; - using matrix_type = typename mechanism_type::matrix_type; + using mechanism_ptr_type = mechanism_ptr<T, I>; + using matrix_type = typename mechanism<T,I>::matrix_type; virtual std::string name() const = 0; - virtual mechanism_type new_mechanism(matrix_type*, index_view) const = 0; + virtual mechanism_ptr<T,I> new_mechanism(matrix_type*, index_view) const = 0; - virtual void set_parameters(mechanism_type&, parameter_list const&) const = 0; + virtual void set_parameters(mechanism_ptr_type&, parameter_list const&) const = 0; }; -extern std::map<std::string, mechanism_helper<double, int>> mechanism_helpers; +template <typename T, typename I> +using mechanism_helper_ptr = + std::unique_ptr<mechanism_helper<T,I>>; + +template <typename M> +mechanism_helper_ptr<typename M::value_type, typename M::size_type> +make_mechanism_helper() +{ + return util::make_unique<M>(); +} + +// for now use a global variable for the map of mechanism helpers +extern std::map< + std::string, + mechanism_helper_ptr<value_type, index_type> +> mechanism_helpers; + void setup_mechanism_helpers(); -} // namespace nest +mechanism_helper_ptr<value_type, index_type>& +get_mechanism_helper(const std::string& name); + +} // namespace mechanisms } // namespace mc +} // namespace nest diff --git a/src/parameter_list.hpp b/src/parameter_list.hpp index 56848466..1080ac49 100644 --- a/src/parameter_list.hpp +++ b/src/parameter_list.hpp @@ -169,13 +169,36 @@ namespace mc { hh_parameters() : base("hh") { - base::add_parameter({"gnabar", 0.12, {0., 1e9}}); - base::add_parameter({"gkbar", 0.036, {0., 1e9}}); - base::add_parameter({"gl", 0.0003,{0., 1e9}}); + base::add_parameter({"gnabar", 0.12, {0, 1e9}}); + base::add_parameter({"gkbar", 0.036, {0, 1e9}}); + base::add_parameter({"gl", 0.0003,{0, 1e9}}); base::add_parameter({"el", -54.3}); } }; + /// parameters for passive channel + class pas_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; + + pas_parameters() + : base("pas") + { + base::add_parameter({"g", 0.001, {0, 1e9}}); + base::add_parameter({"e", -70}); + } + }; + } // namespace mc } // namespace nest diff --git a/src/segment.hpp b/src/segment.hpp index e339e3ac..fd1d8ded 100644 --- a/src/segment.hpp +++ b/src/segment.hpp @@ -113,6 +113,31 @@ class segment { return *it; } + const parameter_list& mechanism(std::string const& n) const + { + auto it = std::find_if( + mechanisms_.begin(), mechanisms_.end(), + [&n](parameter_list const& l){return l.name() == n;} + ); + if(it==mechanisms_.end()) { + throw std::out_of_range( + "attempt to access a parameter that is not defined in a segment" + ); + } + + return *it; + } + + std::vector<parameter_list>& mechanisms() + { + return mechanisms_; + } + + const std::vector<parameter_list>& mechanisms() const + { + return mechanisms_; + } + protected: segmentKind kind_; diff --git a/src/util.hpp b/src/util.hpp index d9239279..3dae43f8 100644 --- a/src/util.hpp +++ b/src/util.hpp @@ -2,6 +2,11 @@ #include "vector/include/Vector.hpp" +#ifdef DEBUG +#define EXPECTS(expression) assert(expression) +#else +#define EXPECTS(expression) +#endif /* using memory::util::red; @@ -80,6 +85,45 @@ namespace util { std::false_type >::type {}; + + // printf with variadic templates for simplified string creation + // mostly used to simplify error string creation + [[gnu::unused]] static + std::string pprintf(const char *s) { + std::string errstring; + while(*s) { + if(*s == '%' && s[1]!='%') { + // instead of throwing an exception, replace with ?? + //throw std::runtime_error("pprintf: the number of arguments did not match the format "); + errstring += "<?>"; + } + else { + errstring += *s; + } + ++s; + } + return errstring; + } + + template <typename T, typename ... Args> + std::string pprintf(const char *s, T value, Args... args) { + std::string errstring; + while(*s) { + if(*s == '%' && s[1]!='%') { + std::stringstream str; + str << value; + errstring += str.str(); + errstring += pprintf(++s, args...); + return errstring; + } + else { + errstring += *s; + ++s; + } + } + return errstring; + } + } // namespace util } // namespace mc } // namespace nest diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index b555fc46..918e00d6 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -17,6 +17,7 @@ set(TEST_SOURCES test_cell.cpp test_compartments.cpp test_matrix.cpp + test_mechanisms.cpp test_point.cpp test_run.cpp test_segment.cpp diff --git a/tests/test_algorithms.cpp b/tests/test_algorithms.cpp index cfe058a3..e6568807 100644 --- a/tests/test_algorithms.cpp +++ b/tests/test_algorithms.cpp @@ -169,3 +169,122 @@ TEST(algorithms, is_positive) ) ); } + +TEST(algorithms, is_unique) +{ + EXPECT_TRUE( + nest::mc::algorithms::is_unique( + std::vector<int>{} + ) + ); + EXPECT_TRUE( + nest::mc::algorithms::is_unique( + std::vector<int>{0} + ) + ); + EXPECT_TRUE( + nest::mc::algorithms::is_unique( + std::vector<int>{0,1,100} + ) + ); + EXPECT_FALSE( + nest::mc::algorithms::is_unique( + std::vector<int>{0,0} + ) + ); + EXPECT_FALSE( + nest::mc::algorithms::is_unique( + std::vector<int>{0,1,2,2,3,4} + ) + ); + EXPECT_FALSE( + nest::mc::algorithms::is_unique( + std::vector<int>{0,1,2,3,4,4} + ) + ); +} + +TEST(algorithms, is_sorted) +{ + EXPECT_TRUE( + nest::mc::algorithms::is_sorted( + std::vector<int>{} + ) + ); + EXPECT_TRUE( + nest::mc::algorithms::is_sorted( + std::vector<int>{100} + ) + ); + EXPECT_TRUE( + nest::mc::algorithms::is_sorted( + std::vector<int>{0,1,2} + ) + ); + EXPECT_TRUE( + nest::mc::algorithms::is_sorted( + std::vector<int>{0,2,100} + ) + ); + EXPECT_TRUE( + nest::mc::algorithms::is_sorted( + std::vector<int>{0,0} + ) + ); + EXPECT_TRUE( + nest::mc::algorithms::is_sorted( + std::vector<int>{0,1,2,2,2,2,3,4,5,5,5} + ) + ); + EXPECT_FALSE( + nest::mc::algorithms::is_sorted( + std::vector<int>{0,1,2,1} + ) + ); + EXPECT_FALSE( + nest::mc::algorithms::is_sorted( + std::vector<int>{1,0} + ) + ); +} + +TEST(algorithms, index_into) +{ + using C = std::vector<int>; + + // by default index_into assumes that the inputs satisfy + // quite a strong set of prerequisites + // + // TODO: test that the EXPECTS() catch bad inputs when DEBUG mode is enabled + // put this in a seperate unit test + auto tests = { + std::make_pair(C{}, C{}), + std::make_pair(C{100}, C{}), + std::make_pair(C{0,1,3,4,6,7,10,11}, C{0,4,6,7,11}), + std::make_pair(C{0,1,3,4,6,7,10,11}, C{0}), + std::make_pair(C{0,1,3,4,6,7,10,11}, C{11}), + std::make_pair(C{0,1,3,4,6,7,10,11}, C{4}), + std::make_pair(C{0,1,3,4,6,7,10,11}, C{0,11}), + std::make_pair(C{0,1,3,4,6,7,10,11}, C{4,11}), + std::make_pair(C{0,1,3,4,6,7,10,11}, C{0,1,3,4,6,7,10,11}) + }; + + auto test_result = [] (const C& super, const C& sub, const C& index) { + if(sub.size()!=index.size()) return false; + for(auto i=0u; i<sub.size(); ++i) { + if(index[i]>=C::value_type(super.size())) return false; + if(super[index[i]]!=sub[i]) return false; + } + return true; + }; + + for(auto& t : tests) { + EXPECT_TRUE( + test_result( + t.first, t.second, + nest::mc::algorithms::index_into(t.first, t.second) + ) + ); + } +} + diff --git a/tests/test_matrix.cpp b/tests/test_matrix.cpp index f0271bd4..73663d68 100644 --- a/tests/test_matrix.cpp +++ b/tests/test_matrix.cpp @@ -116,7 +116,6 @@ TEST(matrix, solve) EXPECT_NEAR(0., std::sqrt(err), 1e-8); } - } } diff --git a/tests/test_mechanisms.cpp b/tests/test_mechanisms.cpp new file mode 100644 index 00000000..f4e5a4d2 --- /dev/null +++ b/tests/test_mechanisms.cpp @@ -0,0 +1,34 @@ +#include "gtest.h" + +#include <mechanism_interface.hpp> +#include <matrix.hpp> + +TEST(mechanisms, helpers) { + nest::mc::mechanisms::setup_mechanism_helpers(); + + EXPECT_EQ(nest::mc::mechanisms::mechanism_helpers.size(), 2u); + + // verify that the hh and pas channels are available + EXPECT_EQ(nest::mc::mechanisms::get_mechanism_helper("hh")->name(), "hh"); + EXPECT_EQ(nest::mc::mechanisms::get_mechanism_helper("pas")->name(), "pas"); + + // check that an out_of_range exception is thrown if an invalid mechanism is + // requested + ASSERT_THROW( + nest::mc::mechanisms::get_mechanism_helper("dachshund"), + std::out_of_range + ); + + //0 1 2 3 4 5 6 7 8 9 + std::vector<int> parent_index = {0,0,1,2,3,4,0,6,7,8}; + memory::HostVector<int> node_indices = std::vector<int>{0,6,7,8,9}; + + nest::mc::matrix<double, int> matrix(parent_index); + + auto& helper = nest::mc::mechanisms::get_mechanism_helper("hh"); + auto mech = helper->new_mechanism(&matrix, node_indices); + + EXPECT_EQ(mech->name(), "hh"); + EXPECT_EQ(mech->size(), 5u); + EXPECT_EQ(mech->matrix_, &matrix); +} diff --git a/tests/test_run.cpp b/tests/test_run.cpp index c4e3b102..55c9c6fc 100644 --- a/tests/test_run.cpp +++ b/tests/test_run.cpp @@ -13,11 +13,18 @@ TEST(run, cable) // 1um radius and 4mm long, all in cm cell.add_cable(0, segmentKind::dendrite, 1e-4, 1e-4, 4e-1); + 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.model().tree.num_segments(), 3u); - std::cout << cell.segment(1)->area() << " is the area\n"; - EXPECT_EQ(cell.model().tree.num_segments(), 2u); + // add passive to all 3 segments in the cell + for(auto& seg :cell.segments()) { + seg->add_mechanism(pas_parameters()); + } cell.soma()->add_mechanism(hh_parameters()); + cell.segment(2)->add_mechanism(hh_parameters()); auto& soma_hh = cell.soma()->mechanism("hh"); @@ -27,23 +34,24 @@ TEST(run, cable) soma_hh.set("el", -54.387); cell.segment(1)->set_compartments(4); + cell.segment(2)->set_compartments(4); using fvm_cell = fvm::fvm_cell<double, int>; fvm_cell fvcell(cell); auto& J = fvcell.jacobian(); - EXPECT_EQ(J.size(), 5u); + EXPECT_EQ(J.size(), 9u); 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; + //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"; + //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.; @@ -91,11 +99,7 @@ TEST(run, init) 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"; + //std::cout << "areas " << fvcell.cv_areas() << "\n"; J.rhs()(memory::all) = 1.; J.rhs()[0] = 10.; diff --git a/vector b/vector index e95bcbf1..8153d030 160000 --- a/vector +++ b/vector @@ -1 +1 @@ -Subproject commit e95bcbf1de6f833b5113a82dcc6d47bb7f41c397 +Subproject commit 8153d030abba20ddd174da9804406139c393a808 -- GitLab