diff --git a/arbor/CMakeLists.txt b/arbor/CMakeLists.txt index 33d410b1208bee626f8d166b77975845cadd79b0..e54d2d5089268661911881efcee00ec720d49c34 100644 --- a/arbor/CMakeLists.txt +++ b/arbor/CMakeLists.txt @@ -50,6 +50,7 @@ set(arbor_sources if(ARB_WITH_CUDA) list(APPEND arbor_sources backends/gpu/mechanism.cpp + backends/gpu/mechanism.cu backends/gpu/shared_state.cpp backends/gpu/stimulus.cpp backends/gpu/stimulus.cu diff --git a/arbor/backends/gpu/mechanism.cpp b/arbor/backends/gpu/mechanism.cpp index c7c629c4dff237dcb51b26a8bda751ebef2e8e84..d73e46e9cfeb60371cce8f4e3f12be6aa1336f03 100644 --- a/arbor/backends/gpu/mechanism.cpp +++ b/arbor/backends/gpu/mechanism.cpp @@ -56,6 +56,7 @@ void mechanism::instantiate(unsigned id, backend::shared_state& shared, const layout& pos_data) { + mult_in_place_ = !pos_data.multiplicity.empty(); mechanism_id_ = id; width_ = pos_data.cv.size(); @@ -122,9 +123,10 @@ void mechanism::instantiate(unsigned id, } // Allocate and initialize index vectors, viz. node_index_ and any ion indices. - // (First sub-array of indices_ is used for node_index_.) + // (First sub-array of indices_ is used for node_index_, last sub-array used for multiplicity_ if it is not empty) - indices_ = iarray((1+num_ions_)*width_padded_); + size_type num_elements = mult_in_place_ ? 2 + num_ions_ : 1 + num_ions_; + indices_ = iarray((num_elements)*width_padded_); memory::copy(make_const_view(pos_data.cv), device_view(indices_.data(), width_)); pp->node_index_ = indices_.data(); @@ -148,6 +150,11 @@ void mechanism::instantiate(unsigned id, ion_index_ptr = index_start; memory::copy(make_const_view(mech_ion_index), device_view(index_start, width_)); } + + if (mult_in_place_) { + memory::copy(make_const_view(pos_data.multiplicity), device_view(indices_.data() + width_padded_, width_)); + pp->multiplicity_ = indices_.data() + (num_ions_ + 1)*width_padded_; + } } void mechanism::set_parameter(const std::string& key, const std::vector<fvm_value_type>& values) { @@ -178,5 +185,20 @@ void mechanism::set_global(const std::string& key, fvm_value_type value) { } } +void multiply_in_place(fvm_value_type* s, const fvm_index_type* p, int n); + +void mechanism::initialize() { + nrn_init(); + mechanism_ppack_base* pp = ppack_ptr(); + auto states = state_table(); + + if(mult_in_place_) { + for (auto& state: states) { + multiply_in_place(*state.second, pp->multiplicity_, pp->width_); + } + } +} + + } // namespace multicore } // namespace arb diff --git a/arbor/backends/gpu/mechanism.cu b/arbor/backends/gpu/mechanism.cu new file mode 100644 index 0000000000000000000000000000000000000000..4fa1f998a5f90b41d96302759c1a52306332acd1 --- /dev/null +++ b/arbor/backends/gpu/mechanism.cu @@ -0,0 +1,28 @@ +#include <iostream> +#include <backends/event.hpp> +#include <backends/multi_event_stream_state.hpp> +#include <backends/gpu/cuda_common.hpp> +#include <backends/gpu/math_cu.hpp> +#include <backends/gpu/mechanism_ppack_base.hpp> +#include <backends/gpu/reduce_by_key.hpp> + +namespace arb { +namespace gpu { + +__global__ +void multiply_in_place_(fvm_value_type* s, const fvm_index_type* p, int n) { + int tid_ = threadIdx.x + blockDim.x*blockIdx.x; + if (tid_<n) { + s[tid_] *= p[tid_]; + } +} + +void multiply_in_place(fvm_value_type* s, const fvm_index_type* p, int n) { + unsigned block_dim = 128; + unsigned grid_dim = gpu::impl::block_count(n, block_dim); + + multiply_in_place_<<<grid_dim, block_dim>>>(s, p, n); +} + +} // namespace gpu +} // namespace arb diff --git a/arbor/backends/gpu/mechanism.hpp b/arbor/backends/gpu/mechanism.hpp index d9da7aa12358fb2fd09ff71ee2c3c9407271adf7..98e70ce382b99d94942dcc9cc705d92dd2add6e5 100644 --- a/arbor/backends/gpu/mechanism.hpp +++ b/arbor/backends/gpu/mechanism.hpp @@ -59,6 +59,8 @@ public: void set_global(const std::string& key, fvm_value_type value) override; + void initialize() override; + protected: size_type width_ = 0; // Instance width (number of CVs/sites) size_type num_ions_ = 0; @@ -78,6 +80,7 @@ protected: iarray indices_; array data_; + bool mult_in_place_; // Generated mechanism field, global and ion table lookup types. // First component is name, second is pointer to corresponing member in @@ -87,6 +90,9 @@ protected: using global_table_entry = std::pair<const char*, value_type*>; using mechanism_global_table = std::vector<global_table_entry>; + using state_table_entry = std::pair<const char*, value_type**>; + using mechanism_state_table = std::vector<state_table_entry>; + using field_table_entry = std::pair<const char*, value_type**>; using mechanism_field_table = std::vector<field_table_entry>; @@ -99,6 +105,8 @@ protected: using ion_index_entry = std::pair<ionKind, const index_type**>; using mechanism_ion_index_table = std::vector<ion_index_entry>; + virtual void nrn_init() = 0; + // Generated mechanisms must implement the following methods, together with // fingerprint(), clone(), kind(), nrn_init(), nrn_state(), nrn_current() // and deliver_events() (if required) from arb::mechanism. @@ -109,6 +117,7 @@ protected: virtual mechanism_field_table field_table() { return {}; } virtual mechanism_field_default_table field_default_table() { return {}; } virtual mechanism_global_table global_table() { return {}; } + virtual mechanism_state_table state_table() { return {}; } virtual mechanism_ion_state_table ion_state_table() { return {}; } virtual mechanism_ion_index_table ion_index_table() { return {}; } diff --git a/arbor/backends/gpu/mechanism_ppack_base.hpp b/arbor/backends/gpu/mechanism_ppack_base.hpp index fb256579dbacd3325094a014c98c1e85f754d643..16518c095fe373ded2042c5fd9557ee98e832eb0 100644 --- a/arbor/backends/gpu/mechanism_ppack_base.hpp +++ b/arbor/backends/gpu/mechanism_ppack_base.hpp @@ -38,6 +38,8 @@ struct mechanism_ppack_base { const value_type* temperature_degC_; const index_type* node_index_; + const index_type* multiplicity_; + const value_type* weight_; }; diff --git a/arbor/backends/multicore/mechanism.cpp b/arbor/backends/multicore/mechanism.cpp index 5dc5c61878aadb47907acb473b31391d80e0cdaa..97aa457e3bb40686983c49f1c65632f93da83a26 100644 --- a/arbor/backends/multicore/mechanism.cpp +++ b/arbor/backends/multicore/mechanism.cpp @@ -64,6 +64,7 @@ void copy_extend(const Source& source, Dest&& dest, const Fill& fill) { void mechanism::instantiate(unsigned id, backend::shared_state& shared, const layout& pos_data) { using util::make_range; + mult_in_place_ = !pos_data.multiplicity.empty(); util::padded_allocator<> pad(shared.alignment); mechanism_id_ = id; width_ = pos_data.cv.size(); @@ -131,10 +132,16 @@ void mechanism::instantiate(unsigned id, backend::shared_state& shared, const la // * For indices in the padded tail of ion index maps, set index to last valid ion index. node_index_ = iarray(width_padded_, pad); + copy_extend(pos_data.cv, node_index_, pos_data.cv.back()); copy_extend(pos_data.weight, make_range(data_.data(), data_.data()+width_padded_), 0); index_constraints_ = make_constraint_partition(node_index_, width_, simd_width); + if (mult_in_place_) { + multiplicity_ = iarray(width_padded_, pad); + copy_extend(pos_data.multiplicity, multiplicity_, 1); + } + for (auto i: ion_index_table()) { util::optional<ion_state&> oion = value_by_key(shared.ion_data, i.first); if (!oion) { @@ -150,7 +157,6 @@ void mechanism::instantiate(unsigned id, backend::shared_state& shared, const la arb_assert(compatible_index_constraints(node_index_, ion_index, simd_width)); } - } void mechanism::set_parameter(const std::string& key, const std::vector<fvm_value_type>& values) { @@ -183,5 +189,19 @@ void mechanism::set_global(const std::string& key, fvm_value_type value) { } } +void mechanism::initialize() { + nrn_init(); + + auto states = state_table(); + + if(mult_in_place_) { + for (auto& state: states) { + for (std::size_t j = 0; j < width_; ++j) { + (*state.second)[j] *= multiplicity_[j]; + } + } + } +} + } // namespace multicore } // namespace arb diff --git a/arbor/backends/multicore/mechanism.hpp b/arbor/backends/multicore/mechanism.hpp index af08a9318ee8af6e1d6f650d0c12bab3833b89f8..21493e947e37796bff1138aac3dc5fd6f2c3d02a 100644 --- a/arbor/backends/multicore/mechanism.hpp +++ b/arbor/backends/multicore/mechanism.hpp @@ -55,6 +55,7 @@ public: } void instantiate(fvm_size_type id, backend::shared_state& shared, const layout& w) override; + void initialize() override; void deliver_events() override { // Delegate to derived class, passing in event queue state. @@ -84,6 +85,8 @@ protected: // Per-mechanism index and weight data, excepting ion indices. iarray node_index_; + iarray multiplicity_; + bool mult_in_place_; constraint_partition index_constraints_; const value_type* weight_; // Points within data_ after instantiation. @@ -99,6 +102,9 @@ protected: using global_table_entry = std::pair<const char*, value_type*>; using mechanism_global_table = std::vector<global_table_entry>; + using state_table_entry = std::pair<const char*, value_type**>; + using mechanism_state_table = std::vector<state_table_entry>; + using field_table_entry = std::pair<const char*, value_type**>; using mechanism_field_table = std::vector<field_table_entry>; @@ -111,6 +117,8 @@ protected: using ion_index_entry = std::pair<ionKind, iarray*>; using mechanism_ion_index_table = std::vector<ion_index_entry>; + virtual void nrn_init() = 0; + // Generated mechanisms must implement the following methods, together with // fingerprint(), clone(), kind(), nrn_init(), nrn_state(), nrn_current() // and deliver_events() (if required) from arb::mechanism. @@ -121,6 +129,7 @@ protected: virtual mechanism_field_table field_table() { return {}; } virtual mechanism_field_default_table field_default_table() { return {}; } virtual mechanism_global_table global_table() { return {}; } + virtual mechanism_state_table state_table() { return {}; } virtual mechanism_ion_state_table ion_state_table() { return {}; } virtual mechanism_ion_index_table ion_index_table() { return {}; } diff --git a/arbor/fvm_layout.cpp b/arbor/fvm_layout.cpp index 24c69019e21350dbda9c80891d0346e2ee7a8868..7e43669cd0336c7bdb5e50cba64be1ce0fd114a5 100644 --- a/arbor/fvm_layout.cpp +++ b/arbor/fvm_layout.cpp @@ -46,6 +46,14 @@ namespace { segment_index = algorithms::make_index(counts); } }; + + struct cv_param { + fvm_size_type cv; + std::vector<fvm_value_type> params; + fvm_size_type target; + }; + + ARB_DEFINE_LEXICOGRAPHIC_ORDERING(cv_param,(a.cv,a.params,a.target),(b.cv,b.params,b.target)) } // namespace // Cable segment discretization @@ -249,7 +257,7 @@ fvm_discretization fvm_discretize(const std::vector<mc_cell>& cells) { // IIb. Density mechanism CVs, parameter values; ion channel default concentration contributions. // IIc. Point mechanism CVs, parameter values, and targets. -fvm_mechanism_data fvm_build_mechanism_data(const mechanism_catalogue& catalogue, const std::vector<mc_cell>& cells, const fvm_discretization& D) { +fvm_mechanism_data fvm_build_mechanism_data(const mechanism_catalogue& catalogue, const std::vector<mc_cell>& cells, const fvm_discretization& D, bool coalesce_syanpses) { using util::assign; using util::sort_by; using util::optional; @@ -271,7 +279,7 @@ fvm_mechanism_data fvm_build_mechanism_data(const mechanism_catalogue& catalogue // 3. Pointer to mechanism metadata from catalogue. struct density_mech_data { - std::vector<std::pair<size_type, const mechanism_desc*>> segments; // + std::vector<std::pair<size_type, const mechanism_desc*>> segments; // string_set paramset; const mechanism_info* info = nullptr; }; @@ -562,19 +570,69 @@ fvm_mechanism_data fvm_build_mechanism_data(const mechanism_catalogue& catalogue fvm_mechanism_config& config = mechdata.mechanisms[name]; config.kind = mechanismKind::point; - assign(config.cv, transform_view(cv_order, [&](size_type j) { return points[j].cv; })); - assign(config.target, transform_view(cv_order, [&](size_type j) { return points[j].target_index; })); + // Generate config.cv: contains cv of group of synapses that can be coalesced into one instance + // Generate config.param_values: contains parameters of group of synapses that can be coalesced into one instance + // Generate multiplicity: contains number of synapses in each coalesced group of synapses + // Generate target: contains the synapse target number + if(catalogue[name].linear && coalesce_syanpses) { + // cv_param_vec used to lexicographically sort the cv, parameters and target, which are stored in that order + std::vector<cv_param> cv_param_vec(cv_order.size()); + + for (unsigned i = 0; i < cv_order.size(); ++i) { + auto loc = cv_order[i]; + std::vector<value_type> p(nparam); + for (auto pidx: make_span(0, nparam)) { + value_type pdefault = param_default[pidx]; + const std::string& pname = param_name[pidx]; + p[pidx] = value_by_key(points[loc].desc->values(), pname).value_or(pdefault); + } + cv_param_vec[i] = cv_param{points[loc].cv, p, points[loc].target_index}; + } - config.param_values.resize(nparam); - for (auto pidx: make_span(0, nparam)) { - value_type pdefault = param_default[pidx]; - const std::string& pname = param_name[pidx]; + std::sort(cv_param_vec.begin(), cv_param_vec.end()); - config.param_values[pidx].first = pname; + auto identical_synapse = [](const cv_param& i, const cv_param& j) { + return i.cv==j.cv && i.params==j.params; + }; - auto& values = config.param_values[pidx].second; - assign(values, transform_view(cv_order, - [&](size_type j) { return value_by_key(points[j].desc->values(), pname).value_or(pdefault); })); + config.param_values.resize(nparam); + for (auto pidx: make_span(0, nparam)) { + config.param_values[pidx].first = param_name[pidx]; + } + config.target.reserve(cv_param_vec.size()); + + for (auto i = cv_param_vec.begin(), j = i; i!=cv_param_vec.end(); i = j) { + ++j; + while (j!=cv_param_vec.end() && identical_synapse(*i, *j)) ++j; + + auto mergeable = util::make_range(i, j); + + config.cv.push_back((*i).cv); + for (auto pidx: make_span(0, nparam)) { + config.param_values[pidx].second.push_back((*i).params[pidx]); + } + config.multiplicity.push_back(mergeable.size()); + + for (auto e: mergeable) { + config.target.push_back(e.target); + } + } + } + else { + assign(config.cv, transform_view(cv_order, [&](size_type j) { return points[j].cv; })); + assign(config.target, transform_view(cv_order, [&](size_type j) { return points[j].target_index; })); + + config.param_values.resize(nparam); + for (auto pidx: make_span(0, nparam)) { + value_type pdefault = param_default[pidx]; + const std::string& pname = param_name[pidx]; + + config.param_values[pidx].first = pname; + + auto& values = config.param_values[pidx].second; + assign(values, transform_view(cv_order, + [&](size_type j) { return value_by_key(points[j].desc->values(), pname).value_or(pdefault); })); + } } } diff --git a/arbor/fvm_layout.hpp b/arbor/fvm_layout.hpp index 4ab94980feed388102cd5f807ae02cbefdceccff..abf7389e2dc1fb3c33f00406ae9ab29fa4dadd35 100644 --- a/arbor/fvm_layout.hpp +++ b/arbor/fvm_layout.hpp @@ -102,6 +102,9 @@ struct fvm_mechanism_config { // duplicates for point mechanisms. std::vector<index_type> cv; + // Coalesced synapse multiplier + std::vector<index_type> multiplicity; + // Normalized area contribution in corresponding CV (density mechanisms only). std::vector<value_type> norm_area; @@ -137,6 +140,6 @@ struct fvm_mechanism_data { std::size_t ntarget = 0; }; -fvm_mechanism_data fvm_build_mechanism_data(const mechanism_catalogue& catalogue, const std::vector<mc_cell>& cells, const fvm_discretization& D); +fvm_mechanism_data fvm_build_mechanism_data(const mechanism_catalogue& catalogue, const std::vector<mc_cell>& cells, const fvm_discretization& D, bool coalesce = true); } // namespace arb diff --git a/arbor/fvm_lowered_cell_impl.hpp b/arbor/fvm_lowered_cell_impl.hpp index 70c132637cbaab4e7d60257fd970c549910f82f8..15d95a44944f7224c80ce4769b3bf17754236a98 100644 --- a/arbor/fvm_lowered_cell_impl.hpp +++ b/arbor/fvm_lowered_cell_impl.hpp @@ -157,7 +157,7 @@ void fvm_lowered_cell_impl<Backend>::reset() { set_tmin(0); for (auto& m: mechanisms_) { - m->nrn_init(); + m->initialize(); } update_ion_state(); @@ -387,7 +387,7 @@ void fvm_lowered_cell_impl<B>::initialize( // Discretize mechanism data. - fvm_mechanism_data mech_data = fvm_build_mechanism_data(*catalogue, cells, D); + fvm_mechanism_data mech_data = fvm_build_mechanism_data(*catalogue, cells, D, global_props.coalesce_synapses); // Discritize and build gap junction info @@ -421,8 +421,12 @@ void fvm_lowered_cell_impl<B>::initialize( mechanism::layout layout; layout.cv = config.cv; + layout.multiplicity = config.multiplicity; layout.weight.resize(layout.cv.size()); + std::vector<fvm_index_type> multiplicity_divs; + auto multiplicity_part = util::make_partition(multiplicity_divs, layout.multiplicity); + // Mechanism weights are F·α where α ∈ [0, 1] is the proportional // contribution in the CV, and F is the scaling factor required // to convert from the mechanism current contribution units to A/m². @@ -431,13 +435,16 @@ void fvm_lowered_cell_impl<B>::initialize( // Point mechanism contributions are in [nA]; CV area A in [µm^2]. // F = 1/A * [nA/µm²] / [A/m²] = 1000/A. - for (auto i: count_along(layout.cv)) { + for (auto i: count_along(config.cv)) { auto cv = layout.cv[i]; layout.weight[i] = 1000/D.cv_area[cv]; // (builtin stimulus, for example, has no targets) + if (!config.target.empty()) { - target_handles[config.target[i]] = target_handle(mech_id, i, cv_to_intdom[cv]); + for (auto j: make_span(multiplicity_part[i])) { + target_handles[config.target[j]] = target_handle(mech_id, i, cv_to_intdom[cv]); + } } } } diff --git a/arbor/include/arbor/mc_cell.hpp b/arbor/include/arbor/mc_cell.hpp index 8ac2adddec06133acfa01d68d5ce6e4e2fbb788b..82aa8cb58e806e1111619ef256b6c2543363f42c 100644 --- a/arbor/include/arbor/mc_cell.hpp +++ b/arbor/include/arbor/mc_cell.hpp @@ -90,6 +90,7 @@ struct mc_cell_global_properties { double temperature_K = constant::hh_squid_temp; // [K] double init_membrane_potential_mV = -65; // [mV] + bool coalesce_synapses = true; }; /// high-level abstract representation of a cell and its segments diff --git a/arbor/include/arbor/mechanism.hpp b/arbor/include/arbor/mechanism.hpp index ff28b43df1f2716f620bde858adeb0df571635cc..8e17ad66bc51c1b71f3dfd053857da7c6012c2c3 100644 --- a/arbor/include/arbor/mechanism.hpp +++ b/arbor/include/arbor/mechanism.hpp @@ -29,6 +29,8 @@ public: struct layout { std::vector<fvm_index_type> cv; // Maps in-instance index to CV index. std::vector<fvm_value_type> weight; // Maps in-instance index to compartment contribution. + std::vector<fvm_index_type> multiplicity; // Number of logical point processes at in-instance index; + // if empty point processes are not coalesced, all multipliers are 1 }; // Return fingerprint of mechanism dynamics source description for validation/replication. @@ -61,7 +63,7 @@ public: virtual void set_parameter(const std::string& key, const std::vector<fvm_value_type>& values) = 0; // Simulation interfaces: - virtual void nrn_init() = 0; + virtual void initialize() = 0; virtual void nrn_state() = 0; virtual void nrn_current() = 0; virtual void deliver_events() {}; diff --git a/arbor/include/arbor/mechinfo.hpp b/arbor/include/arbor/mechinfo.hpp index 41a3c631f6987e1c28a2fe2ebd051aa24a55dbed..dce3c70dc24bd5b732f2ef251fa8417db9f49fdd 100644 --- a/arbor/include/arbor/mechinfo.hpp +++ b/arbor/include/arbor/mechinfo.hpp @@ -60,6 +60,8 @@ struct mechanism_info { std::unordered_map<ionKind, ion_dependency> ions; mechanism_fingerprint fingerprint; + + bool linear = false; }; } // namespace arb diff --git a/cmake/CompilerOptions.cmake b/cmake/CompilerOptions.cmake index 4717ec0aee086fd4fd435126f2c941c56070aafb..03dadabc5c7a6de851d5e5248914bf4a1e64004d 100644 --- a/cmake/CompilerOptions.cmake +++ b/cmake/CompilerOptions.cmake @@ -31,6 +31,14 @@ set(CXXOPT_WALL $<IF:$<CXX_COMPILER_ID:Clang>,-Wno-potentially-evaluated-expression,> $<IF:$<CXX_COMPILER_ID:AppleClang>,-Wno-potentially-evaluated-expression,> + # Clang: + # + # * Disable 'unused-function' warning: this will flag the unused + # functions generated by `ARB_DEFINE_LEXICOGRAPHIC_ORDERING` + + $<IF:$<CXX_COMPILER_ID:Clang>,-Wno-unused-function,> + $<IF:$<CXX_COMPILER_ID:AppleClang>,-Wno-unused-function,> + # * Clang erroneously warns that T is an 'unused type alias' # in code like this: # struct X { diff --git a/modcc/module.cpp b/modcc/module.cpp index 81312784081f2bf80f15b78d191eb8b8db502966..93e5ac7c414e595cb97ba88643b6b6cdc974c404 100644 --- a/modcc/module.cpp +++ b/modcc/module.cpp @@ -141,6 +141,8 @@ bool Module::semantic() { // symbol table. // Returns false if a symbol name clashes with the name of a symbol that // is already in the symbol table. + bool linear = true; + std::vector<std::string> state_vars; auto move_symbols = [this] (std::vector<symbol_ptr>& symbol_list) { for(auto& symbol: symbol_list) { bool is_found = (symbols_.find(symbol->name()) != symbols_.end()); @@ -231,6 +233,7 @@ bool Module::semantic() { auto state = sym.second->is_variable(); if (!state || !state->is_state()) continue; + state_vars.push_back(state->name()); auto shadowed = state->shadows(); if (!shadowed) continue; @@ -334,6 +337,13 @@ bool Module::semantic() { } else { deriv->body()->accept(solver.get()); + for (auto& s: deriv->body()->statements()) { + if(s->is_assignment() && !state_vars.empty()) { + linear_test_result r = linear_test(s->is_assignment()->rhs(), state_vars); + linear &= r.is_linear; + linear &= r.is_homogeneous; + } + } } if (auto solve_block = solver->as_block(false)) { @@ -382,6 +392,15 @@ bool Module::semantic() { //.......................................................... NrnCurrentRewriter nrn_current_rewriter; breakpoint->accept(&nrn_current_rewriter); + + for (auto& s: breakpoint->body()->statements()) { + if(s->is_assignment() && !state_vars.empty()) { + linear_test_result r = linear_test(s->is_assignment()->rhs(), state_vars); + linear &= r.is_linear; + linear &= r.is_homogeneous; + } + } + auto nrn_current_block = nrn_current_rewriter.as_block(); if (!nrn_current_block) { append_errors(nrn_current_rewriter.errors()); @@ -400,6 +419,32 @@ bool Module::semantic() { return false; } + if (has_symbol("net_receive", symbolKind::procedure)) { + auto net_rec_api = make_empty_api_method("net_rec_api", "net_receive"); + if (net_rec_api.second) { + for (auto &s: net_rec_api.second->body()->statements()) { + if (s->is_assignment()) { + for (const auto &id: state_vars) { + auto coef = symbolic_pdiff(s->is_assignment()->rhs(), id); + if(coef->is_number()) { + if (!s->is_assignment()->lhs()->is_identifier()) { + error(pprintf("Left hand side of assignment is not an identifier")); + return false; + } + linear &= s->is_assignment()->lhs()->is_identifier()->name() == id ? + coef->is_number()->value() == 1 : + coef->is_number()->value() == 0; + } + else { + linear = false; + } + } + } + } + } + } + linear_ = linear; + // Perform semantic analysis and inlining on function and procedure bodies // in order to inline calls inside the newly crated API methods. semantic_func_proc(); diff --git a/modcc/module.hpp b/modcc/module.hpp index 8b650870efc28d6b8feadc96ca30b9341cf73842..dbe8c1ffcb6008bd0bcd26b48b103070fbece99e 100644 --- a/modcc/module.hpp +++ b/modcc/module.hpp @@ -103,6 +103,7 @@ public: return find_ion(k) != neuron_block().ions.end(); }; + bool is_linear() const { return linear_; } private: moduleKind kind_; @@ -116,6 +117,7 @@ private: UnitsBlock units_block_; ParameterBlock parameter_block_; AssignedBlock assigned_block_; + bool linear_; // AST storage. std::vector<symbol_ptr> callables_; diff --git a/modcc/printer/cprinter.cpp b/modcc/printer/cprinter.cpp index 19eda4d04dd649fa15dfa0856a7dd52e8b389814..1c1ce229350cefccddefc94de01ea3b92c945c95 100644 --- a/modcc/printer/cprinter.cpp +++ b/modcc/printer/cprinter.cpp @@ -258,6 +258,19 @@ std::string emit_cpp_source(const Module& module_, const printer_options& opt) { } out << popindent << "\n};" << popindent << "\n}\n"; + out << + "mechanism_state_table state_table() override {\n" << indent << + "return {" << indent; + + sep.reset(); + for (const auto& array: vars.arrays) { + auto memb = array->name(); + if(array->is_state()) { + out << sep << "{" << quote(memb) << ", &" << memb << "}"; + } + } + out << popindent << "\n};" << popindent << "\n}\n"; + } if (!ion_deps.empty()) { diff --git a/modcc/printer/cudaprinter.cpp b/modcc/printer/cudaprinter.cpp index 1013d5a80cd6efc59d84133e672c677505283c61..a5304524a2444abf0ef24379eaf883eace420c52 100644 --- a/modcc/printer/cudaprinter.cpp +++ b/modcc/printer/cudaprinter.cpp @@ -160,6 +160,20 @@ std::string emit_cuda_cpp_source(const Module& module_, const printer_options& o } out << popindent << "\n};" << popindent << "\n}\n"; + out << + "mechanism_state_table state_table() override {\n" << indent << + "return {" << indent; + + sep.reset(); + for (const auto& array: vars.arrays) { + auto memb = array->name(); + if(array->is_state()) { + out << sep << "{" << quote(memb) << ", &pp_." << memb << "}"; + } + } + out << popindent << "\n};" << popindent << "\n}\n"; + + } if (!ion_deps.empty()) { diff --git a/modcc/printer/infoprinter.cpp b/modcc/printer/infoprinter.cpp index 0d64f56a360708c586f2b041496efd4ddcd125fa..7c69918e834fb5b84b36eedbefe7481e493f4379 100644 --- a/modcc/printer/infoprinter.cpp +++ b/modcc/printer/infoprinter.cpp @@ -116,7 +116,8 @@ std::string build_info_header(const Module& m, const printer_options& opt) { std::string fingerprint = "<placeholder>"; out << popindent << "\n" "},\n" - "// fingerprint\n" << quote(fingerprint) << "\n" + "// fingerprint\n" << quote(fingerprint) << ",\n" + "// linear, homogeneous mechanism\n" << m.is_linear() << "\n" << popindent << "};\n" "\n" diff --git a/test/unit-modcc/test.mod b/test/unit-modcc/mod_files/test0.mod similarity index 100% rename from test/unit-modcc/test.mod rename to test/unit-modcc/mod_files/test0.mod diff --git a/test/unit-modcc/mod_files/test1.mod b/test/unit-modcc/mod_files/test1.mod new file mode 100644 index 0000000000000000000000000000000000000000..ce62502aa6920c0ae3615c657d6acdd116ee9ea9 --- /dev/null +++ b/test/unit-modcc/mod_files/test1.mod @@ -0,0 +1,38 @@ +NEURON { +POINT_PROCESS expsyn +RANGE tau, e +NONSPECIFIC_CURRENT i +} + +UNITS { +(mV) = (millivolt) +} + +PARAMETER { +tau = 2.0 (ms) : the default for Neuron is 0.1 +e = 0 (mV) +} + +ASSIGNED {} + +STATE { +g +} + +INITIAL { +g=0 +} + +BREAKPOINT { +SOLVE state METHOD cnexp +i = g*(v - e) +} + +DERIVATIVE state { +g' = -g/tau +} + +NET_RECEIVE(weight) { +g = g + weight +} + diff --git a/test/unit-modcc/mod_files/test2.mod b/test/unit-modcc/mod_files/test2.mod new file mode 100644 index 0000000000000000000000000000000000000000..a66a8b7174982901ae9189ebeb608d0640b004bf --- /dev/null +++ b/test/unit-modcc/mod_files/test2.mod @@ -0,0 +1,25 @@ +UNITS { +(mV) = (millivolt) +(S) = (siemens) +} + +NEURON { +SUFFIX pas +NONSPECIFIC_CURRENT i +RANGE g, e +} + +INITIAL {} + +PARAMETER { +g = .001 (S/cm2) +e = -65 (mV) : we use -65 for the ball and stick model, instead of Neuron default of -70 +} + +ASSIGNED { +v (mV) +} + +BREAKPOINT { +i = g*(v - e) +} diff --git a/test/unit-modcc/mod_files/test3.mod b/test/unit-modcc/mod_files/test3.mod new file mode 100644 index 0000000000000000000000000000000000000000..e43f6ce39ec5c64cd7d047888af59fef001b1567 --- /dev/null +++ b/test/unit-modcc/mod_files/test3.mod @@ -0,0 +1,38 @@ +NEURON { +POINT_PROCESS expsyn +RANGE tau, e +NONSPECIFIC_CURRENT i +} + +UNITS { +(mV) = (millivolt) +} + +PARAMETER { +tau = 2.0 (ms) : the default for Neuron is 0.1 +e = 0 (mV) +} + +ASSIGNED {} + +STATE { +g +} + +INITIAL { +g=0 +} + +BREAKPOINT { +SOLVE state METHOD cnexp +i = g*(v - e) +} + +DERIVATIVE state { +g' = -g/tau + tau +} + +NET_RECEIVE(weight) { +g = g + weight +} + diff --git a/test/unit-modcc/mod_files/test4.mod b/test/unit-modcc/mod_files/test4.mod new file mode 100644 index 0000000000000000000000000000000000000000..68b1f4a9c9ea35c169d662318cb4cc4aeb15599e --- /dev/null +++ b/test/unit-modcc/mod_files/test4.mod @@ -0,0 +1,38 @@ +NEURON { +POINT_PROCESS expsyn +RANGE tau, e +NONSPECIFIC_CURRENT i +} + +UNITS { +(mV) = (millivolt) +} + +PARAMETER { +tau = 2.0 (ms) : the default for Neuron is 0.1 +e = 0 (mV) +} + +ASSIGNED {} + +STATE { +g +} + +INITIAL { +g=0 +} + +BREAKPOINT { +SOLVE state METHOD cnexp +i = g*(v - e) +} + +DERIVATIVE state { +g' = -g/tau + tau +} + +NET_RECEIVE(weight) { +g = g * weight +} + diff --git a/test/unit-modcc/mod_files/test5.mod b/test/unit-modcc/mod_files/test5.mod new file mode 100644 index 0000000000000000000000000000000000000000..fc0574ffbeec3427ab38d019ca9e84119c083472 --- /dev/null +++ b/test/unit-modcc/mod_files/test5.mod @@ -0,0 +1,38 @@ +NEURON { +POINT_PROCESS expsyn +RANGE tau, e +NONSPECIFIC_CURRENT i +} + +UNITS { +(mV) = (millivolt) +} + +PARAMETER { +tau = 2.0 (ms) : the default for Neuron is 0.1 +e = 0 (mV) +} + +ASSIGNED {} + +STATE { +g +} + +INITIAL { +g=0 +} + +BREAKPOINT { +SOLVE state METHOD cnexp +i = g*(v - e) + 2 +} + +DERIVATIVE state { +g' = -g/tau +} + +NET_RECEIVE(weight) { +g = g + weight +} + diff --git a/test/unit-modcc/test_module.cpp b/test/unit-modcc/test_module.cpp index e82f0058350397ea9fe7ee8cc3bd3be9e368b8ce..e71964525f99d5b607df91c4b3780a7af032358f 100644 --- a/test/unit-modcc/test_module.cpp +++ b/test/unit-modcc/test_module.cpp @@ -3,7 +3,7 @@ #include "module.hpp" TEST(Module, open) { - Module m(io::read_all(DATADIR "/test.mod"), "test.mod"); + Module m(io::read_all(DATADIR "/mod_files/test0.mod"), "test0.mod"); if(!m.buffer().size()) { std::cout << "skipping Module.open test because unable to open input file" << std::endl; return; @@ -15,3 +15,31 @@ TEST(Module, open) { EXPECT_NE(t.type, tok::reserved); } } + +TEST(Module, linear_mechanisms) { + for(int i = 1; i < 6; i++) + { + auto file_name = "test"+std::to_string(i)+".mod"; + + Module m(io::read_all(DATADIR "/mod_files/" + file_name), file_name); + if (!m.buffer().size()) { + std::cout << "skipping Module.open test because unable to open input file" << std::endl; + return; + } + + Parser p(m, false); + if (!p.parse()) { + std::cout << "problem with parsing input file" << std::endl; + return; + } + + m.semantic(); + + if(i < 3) { + EXPECT_TRUE(m.is_linear()); + } + else { + EXPECT_FALSE(m.is_linear()); + } + } +} diff --git a/test/unit-modcc/test_parser.cpp b/test/unit-modcc/test_parser.cpp index 6e83fc28f60c9ec9a6ea5efe2c490cdaade173ba..65a77c3a926ed5fbd7d36f6b94f0d2e0210c802f 100644 --- a/test/unit-modcc/test_parser.cpp +++ b/test/unit-modcc/test_parser.cpp @@ -71,7 +71,7 @@ template <typename RetUniqPtr> } TEST(Parser, full_file) { - Module m(io::read_all(DATADIR "/test.mod"), "test.mod"); + Module m(io::read_all(DATADIR "/mod_files/test0.mod"), "test0.mod"); if (m.buffer().size()==0) { std::cout << "skipping Parser.full_file test because unable to open input file" << std::endl; return; diff --git a/test/unit/test_fvm_layout.cpp b/test/unit/test_fvm_layout.cpp index d81e37457c03359ad5cd0137bd81eeb7aa3c74c3..a3d6039e13325652953985c0292932abbbe04948 100644 --- a/test/unit/test_fvm_layout.cpp +++ b/test/unit/test_fvm_layout.cpp @@ -312,9 +312,10 @@ TEST(fvm_layout, mech_index) { EXPECT_TRUE(testing::seq_almost_eq<double>(norm_area, hh_config.norm_area)); // Three expsyn synapses, two 0.4 along segment 1, and one 0.4 along segment 5. + // These two synapses can be coalesced into 1 synapse // 0.4 along => second (non-parent) CV for segment. - EXPECT_EQ(ivec({2, 2, 15}), expsyn_config.cv); + EXPECT_EQ(ivec({2, 15}), expsyn_config.cv); // One exp2syn synapse, 0.4 along segment 4. @@ -331,6 +332,151 @@ TEST(fvm_layout, mech_index) { EXPECT_EQ(ivec({0,5}), M.ions.at(ionKind::k).cv); } +TEST(fvm_layout, coalescing_synapses) { + using ivec = std::vector<fvm_index_type>; + using fvec = std::vector<fvm_value_type>; + + auto syn_desc = [&](const char* name, double val0, double val1) { + mechanism_desc m(name); + m.set("e", val0); + m.set("tau", val1); + return m; + }; + + auto syn_desc_2 = [&](const char* name, double val0, double val1) { + mechanism_desc m(name); + m.set("e", val0); + m.set("tau1", val1); + return m; + }; + + { + mc_cell cell = make_cell_ball_and_stick(); + + // Add synapses of two varieties. + cell.add_synapse({1, 0.3}, "expsyn"); + cell.add_synapse({1, 0.5}, "expsyn"); + cell.add_synapse({1, 0.7}, "expsyn"); + cell.add_synapse({1, 0.9}, "expsyn"); + + fvm_discretization D = fvm_discretize({cell}); + fvm_mechanism_data M = fvm_build_mechanism_data(global_default_catalogue(), {cell}, D); + + auto &expsyn_config = M.mechanisms.at("expsyn"); + EXPECT_EQ(ivec({1, 2, 3, 4}), expsyn_config.cv); + EXPECT_EQ(ivec({1, 1, 1, 1}), expsyn_config.multiplicity); + } + { + mc_cell cell = make_cell_ball_and_stick(); + + // Add synapses of two varieties. + cell.add_synapse({1, 0.3}, "expsyn"); + cell.add_synapse({1, 0.5}, "exp2syn"); + cell.add_synapse({1, 0.7}, "expsyn"); + cell.add_synapse({1, 0.9}, "exp2syn"); + + fvm_discretization D = fvm_discretize({cell}); + fvm_mechanism_data M = fvm_build_mechanism_data(global_default_catalogue(), {cell}, D); + + auto &expsyn_config = M.mechanisms.at("expsyn"); + EXPECT_EQ(ivec({1, 3}), expsyn_config.cv); + EXPECT_EQ(ivec({1, 1}), expsyn_config.multiplicity); + + auto &exp2syn_config = M.mechanisms.at("exp2syn"); + EXPECT_EQ(ivec({2, 4}), exp2syn_config.cv); + EXPECT_EQ(ivec({1, 1}), exp2syn_config.multiplicity); + } + { + mc_cell cell = make_cell_ball_and_stick(); + + // Add synapses of two varieties. + cell.add_synapse({1, 0.3}, "expsyn"); + cell.add_synapse({1, 0.3}, "expsyn"); + cell.add_synapse({1, 0.7}, "expsyn"); + cell.add_synapse({1, 0.7}, "expsyn"); + + fvm_discretization D = fvm_discretize({cell}); + fvm_mechanism_data M = fvm_build_mechanism_data(global_default_catalogue(), {cell}, D); + + auto &expsyn_config = M.mechanisms.at("expsyn"); + EXPECT_EQ(ivec({1, 3}), expsyn_config.cv); + EXPECT_EQ(ivec({2, 2}), expsyn_config.multiplicity); + } + { + mc_cell cell = make_cell_ball_and_stick(); + + // Add synapses of two varieties. + cell.add_synapse({1, 0.3}, syn_desc("expsyn", 0, 0.2)); + cell.add_synapse({1, 0.3}, syn_desc("expsyn", 0, 0.2)); + cell.add_synapse({1, 0.3}, syn_desc("expsyn", 0.1, 0.2)); + cell.add_synapse({1, 0.7}, syn_desc("expsyn", 0.1, 0.2)); + + fvm_discretization D = fvm_discretize({cell}); + fvm_mechanism_data M = fvm_build_mechanism_data(global_default_catalogue(), {cell}, D); + + auto &expsyn_config = M.mechanisms.at("expsyn"); + EXPECT_EQ(ivec({1, 1, 3}), expsyn_config.cv); + EXPECT_EQ(ivec({2, 1, 1}), expsyn_config.multiplicity); + EXPECT_EQ(fvec({0, 0.1, 0.1}), expsyn_config.param_values[0].second); + EXPECT_EQ(fvec({0.2, 0.2, 0.2}), expsyn_config.param_values[1].second); + } + { + mc_cell cell = make_cell_ball_and_stick(); + + // Add synapses of two varieties. + cell.add_synapse({1, 0.7}, syn_desc("expsyn", 0, 3)); + cell.add_synapse({1, 0.7}, syn_desc("expsyn", 1, 3)); + cell.add_synapse({1, 0.7}, syn_desc("expsyn", 0, 3)); + cell.add_synapse({1, 0.7}, syn_desc("expsyn", 1, 3)); + cell.add_synapse({1, 0.3}, syn_desc("expsyn", 0, 2)); + cell.add_synapse({1, 0.3}, syn_desc("expsyn", 1, 2)); + cell.add_synapse({1, 0.3}, syn_desc("expsyn", 0, 2)); + cell.add_synapse({1, 0.3}, syn_desc("expsyn", 1, 2)); + + fvm_discretization D = fvm_discretize({cell}); + fvm_mechanism_data M = fvm_build_mechanism_data(global_default_catalogue(), {cell}, D); + + auto &expsyn_config = M.mechanisms.at("expsyn"); + EXPECT_EQ(ivec({1, 1, 3, 3}), expsyn_config.cv); + EXPECT_EQ(ivec({4, 6, 5, 7, 0, 2, 1, 3}), expsyn_config.target); + EXPECT_EQ(ivec({2, 2, 2, 2}), expsyn_config.multiplicity); + EXPECT_EQ(fvec({0, 1, 0, 1}), expsyn_config.param_values[0].second); + EXPECT_EQ(fvec({2, 2, 3, 3}), expsyn_config.param_values[1].second); + } + { + mc_cell cell = make_cell_ball_and_stick(); + + // Add synapses of two varieties. + cell.add_synapse({1, 0.3}, syn_desc("expsyn", 1, 2)); + cell.add_synapse({1, 0.3}, syn_desc_2("exp2syn", 4, 1)); + cell.add_synapse({1, 0.3}, syn_desc("expsyn", 1, 2)); + cell.add_synapse({1, 0.3}, syn_desc("expsyn", 5, 1)); + cell.add_synapse({1, 0.3}, syn_desc_2("exp2syn", 1, 3)); + cell.add_synapse({1, 0.3}, syn_desc("expsyn", 1, 2)); + cell.add_synapse({1, 0.7}, syn_desc_2("exp2syn", 2, 2)); + cell.add_synapse({1, 0.7}, syn_desc_2("exp2syn", 2, 1)); + cell.add_synapse({1, 0.7}, syn_desc_2("exp2syn", 2, 1)); + cell.add_synapse({1, 0.7}, syn_desc_2("exp2syn", 2, 2)); + + fvm_discretization D = fvm_discretize({cell}); + fvm_mechanism_data M = fvm_build_mechanism_data(global_default_catalogue(), {cell}, D); + + auto &expsyn_config = M.mechanisms.at("expsyn"); + EXPECT_EQ(ivec({1, 1}), expsyn_config.cv); + EXPECT_EQ(ivec({0, 2, 5, 3}), expsyn_config.target); + EXPECT_EQ(ivec({3, 1}), expsyn_config.multiplicity); + EXPECT_EQ(fvec({1, 5}), expsyn_config.param_values[0].second); + EXPECT_EQ(fvec({2, 1}), expsyn_config.param_values[1].second); + + auto &exp2syn_config = M.mechanisms.at("exp2syn"); + EXPECT_EQ(ivec({1, 1, 3, 3}), exp2syn_config.cv); + EXPECT_EQ(ivec({4, 1, 7, 8, 6, 9}), exp2syn_config.target); + EXPECT_EQ(ivec({1, 1, 2, 2}), exp2syn_config.multiplicity); + EXPECT_EQ(fvec({1, 4, 2, 2}), exp2syn_config.param_values[0].second); + EXPECT_EQ(fvec({3, 1, 1, 2}), exp2syn_config.param_values[1].second); + } +} + TEST(fvm_layout, synapse_targets) { std::vector<mc_cell> cells = two_cell_system(); diff --git a/test/unit/test_mech_temperature.cpp b/test/unit/test_mech_temperature.cpp index e4e10226bf40b489f1f07a6bca40cedb5cae8260..3fcb4caaabcea667f33ca4f761f323668858ebaa 100644 --- a/test/unit/test_mech_temperature.cpp +++ b/test/unit/test_mech_temperature.cpp @@ -44,7 +44,7 @@ void run_celsius_test() { // expect 0 value in state 'c' after init: - celsius_test->nrn_init(); + celsius_test->initialize(); std::vector<fvm_value_type> expected_c_values(ncv, 0.); EXPECT_EQ(expected_c_values, mechanism_field(celsius_test.get(), "c")); @@ -62,7 +62,7 @@ void run_celsius_test() { temperature_C = temperature_K-273.15; shared_state->reset(-65., temperature_K); - celsius_test->nrn_init(); + celsius_test->initialize(); celsius_test->nrn_state(); expected_c_values.assign(ncv, temperature_C); diff --git a/test/unit/test_mechcat.cpp b/test/unit/test_mechcat.cpp index c19418f50e39c5bb531eb2a03ad905785098592b..44df5f86616958064a42bb2387c10dd47f924e35 100644 --- a/test/unit/test_mechcat.cpp +++ b/test/unit/test_mechcat.cpp @@ -69,7 +69,7 @@ struct common_impl: concrete_mechanism<B> { overrides_[key] = v; } - void nrn_init() override {} + void initialize() override {} void nrn_state() override {} void nrn_current() override {} void deliver_events() override {} @@ -286,7 +286,7 @@ TEST(mechcat, instantiate) { // write its specialized global variables to shared state, but we do in // these tests for testing purposes. - mechanism::layout layout = {{0u, 1u, 2u}, {1., 2., 1.}}; + mechanism::layout layout = {{0u, 1u, 2u}, {1., 2., 1.}, {1u, 1u, 1u}}; bar_backend::shared_state bar_state; auto cat = build_fake_catalogue(); diff --git a/test/unit/test_synapses.cpp b/test/unit/test_synapses.cpp index 9219edbd3f8fd2941837b3099a4510304bd1783a..c433099fa67d0d2c25084bc1e4658ff5ed5642e6 100644 --- a/test/unit/test_synapses.cpp +++ b/test/unit/test_synapses.cpp @@ -95,10 +95,11 @@ TEST(synapses, syn_basic_state) { state.set_dt(); std::vector<index_type> syn_cv(num_syn, 0); + std::vector<index_type> syn_mult(num_syn, 1); std::vector<value_type> syn_weight(num_syn, 1.0); - expsyn->instantiate(0, state, {syn_cv, syn_weight}); - exp2syn->instantiate(1, state, {syn_cv, syn_weight}); + expsyn->instantiate(0, state, {syn_cv, syn_weight, syn_mult}); + exp2syn->instantiate(1, state, {syn_cv, syn_weight, syn_mult}); // Parameters initialized to default values? @@ -130,10 +131,10 @@ TEST(synapses, syn_basic_state) { // Initialize state then check g, A, B have been set to zero. - expsyn->nrn_init(); + expsyn->initialize(); EXPECT_TRUE(all_equal_to(mechanism_field(expsyn, "g"), 0.)); - exp2syn->nrn_init(); + exp2syn->initialize(); EXPECT_TRUE(all_equal_to(mechanism_field(exp2syn, "A"), 0.)); EXPECT_TRUE(all_equal_to(mechanism_field(exp2syn, "B"), 0.));