From 7b0b69f3eae977a7ac0e4d2ae32b825456887255 Mon Sep 17 00:00:00 2001
From: noraabiakar <nora.abiakar@gmail.com>
Date: Mon, 18 Feb 2019 16:44:11 +0100
Subject: [PATCH] Integrate over integration-domains instead of cells in the
lowered/back-end implementation (#688)
Fixes #683.
* Distinguish between per-cell matrix solver and per-'integration domain' state in the FVM back-end, adding new index for FVM back-ends mapping CVs to integration domain index.
* Move to one event queue per integration domain, avoiding extra logic for synchronisation of times across gap-junction connected cells. This requires a merging of event lanes in `mc_cell_group` when constructing the staged events in the multi-event queue.
* Add synapses to gap junction example.
---
arbor/backends/event.hpp | 8 +-
arbor/backends/gpu/fvm.hpp | 2 +-
arbor/backends/gpu/matrix_assemble.cu | 20 +-
arbor/backends/gpu/matrix_fine.cu | 10 +-
arbor/backends/gpu/matrix_fine.hpp | 3 +-
arbor/backends/gpu/matrix_state_fine.hpp | 13 +-
arbor/backends/gpu/matrix_state_flat.hpp | 21 +-
.../backends/gpu/matrix_state_interleaved.hpp | 19 +-
arbor/backends/gpu/mechanism.cpp | 2 +-
arbor/backends/gpu/shared_state.cpp | 40 ++--
arbor/backends/gpu/shared_state.cu | 39 +--
arbor/backends/gpu/shared_state.hpp | 21 +-
arbor/backends/gpu/threshold_watcher.cu | 8 +-
arbor/backends/gpu/threshold_watcher.hpp | 10 +-
arbor/backends/multicore/fvm.hpp | 2 +-
arbor/backends/multicore/matrix_state.hpp | 14 +-
arbor/backends/multicore/mechanism.cpp | 2 +-
arbor/backends/multicore/shared_state.cpp | 59 ++---
arbor/backends/multicore/shared_state.hpp | 21 +-
.../backends/multicore/threshold_watcher.hpp | 8 +-
arbor/fvm_layout.cpp | 1 +
arbor/fvm_lowered_cell.hpp | 2 +-
arbor/fvm_lowered_cell_impl.hpp | 81 ++++++-
arbor/matrix.hpp | 9 +-
arbor/mc_cell_group.cpp | 96 +++-----
arbor/mc_cell_group.hpp | 18 +-
example/gap_junctions/gap_junctions.cpp | 76 ++++--
example/gap_junctions/parameters.hpp | 22 +-
example/gap_junctions/readme.md | 42 ++++
test/unit/CMakeLists.txt | 1 -
test/unit/test_fvm_lowered.cpp | 226 ++++++++++++++++--
test/unit/test_matrix.cpp | 12 +-
test/unit/test_matrix.cu | 8 +-
test/unit/test_matrix_cpuvsgpu.cpp | 6 +-
test/unit/test_mc_cell_group.cpp | 149 ------------
test/unit/test_mech_temperature.cpp | 5 +-
test/unit/test_probe.cpp | 3 +-
test/unit/test_supercell.cpp | 71 ------
test/unit/test_synapses.cpp | 4 +-
39 files changed, 602 insertions(+), 552 deletions(-)
delete mode 100644 test/unit/test_supercell.cpp
diff --git a/arbor/backends/event.hpp b/arbor/backends/event.hpp
index deba95ee..cdf4b5cc 100644
--- a/arbor/backends/event.hpp
+++ b/arbor/backends/event.hpp
@@ -13,11 +13,11 @@ namespace arb {
struct target_handle {
cell_local_size_type mech_id; // mechanism type identifier (per cell group).
cell_local_size_type mech_index; // instance of the mechanism
- cell_size_type cell_index; // which cell (acts as index into e.g. vec_t)
+ cell_size_type intdom_index; // which integration domain (acts as index into e.g. vec_t)
target_handle() {}
- target_handle(cell_local_size_type mech_id, cell_local_size_type mech_index, cell_size_type cell_index):
- mech_id(mech_id), mech_index(mech_index), cell_index(cell_index) {}
+ target_handle(cell_local_size_type mech_id, cell_local_size_type mech_index, cell_size_type intdom_index):
+ mech_id(mech_id), mech_index(mech_index), intdom_index(intdom_index) {}
};
struct deliverable_event {
@@ -33,7 +33,7 @@ struct deliverable_event {
// Stream index accessor function for multi_event_stream:
inline cell_size_type event_index(const deliverable_event& ev) {
- return ev.handle.cell_index;
+ return ev.handle.intdom_index;
}
// Subset of event information required for mechanism delivery.
diff --git a/arbor/backends/gpu/fvm.hpp b/arbor/backends/gpu/fvm.hpp
index cc82bed0..f19a71bd 100644
--- a/arbor/backends/gpu/fvm.hpp
+++ b/arbor/backends/gpu/fvm.hpp
@@ -64,7 +64,7 @@ struct backend {
const execution_context& context)
{
return threshold_watcher(
- state.cv_to_cell.data(),
+ state.cv_to_intdom.data(),
state.time.data(),
state.time_to.data(),
state.voltage.data(),
diff --git a/arbor/backends/gpu/matrix_assemble.cu b/arbor/backends/gpu/matrix_assemble.cu
index fde63e1a..bcc589d2 100644
--- a/arbor/backends/gpu/matrix_assemble.cu
+++ b/arbor/backends/gpu/matrix_assemble.cu
@@ -24,14 +24,15 @@ void assemble_matrix_flat(
const T* cv_capacitance,
const T* area,
const I* cv_to_cell,
- const T* dt_cell,
+ const T* dt_intdom,
+ const I* cell_to_intdom,
unsigned n)
{
const unsigned tid = threadIdx.x + blockDim.x*blockIdx.x;
if (tid<n) {
auto cid = cv_to_cell[tid];
- auto dt = dt_cell[cid];
+ auto dt = dt_intdom[cell_to_intdom[cid]];
// Note: dt==0 case is expected only at the end of a mindelay/2
// integration period, and consequently divergence is unlikely
@@ -73,7 +74,8 @@ void assemble_matrix_interleaved(
const I* sizes,
const I* starts,
const I* matrix_to_cell,
- const T* dt_cell,
+ const T* dt_intdom,
+ const I* cell_to_intdom,
unsigned padded_size, unsigned num_mtx)
{
static_assert(BlockWidth*LoadWidth==Threads,
@@ -107,7 +109,7 @@ void assemble_matrix_interleaved(
if (permuted_cid<num_mtx) {
auto cid = matrix_to_cell[permuted_cid];
- dt = dt_cell[cid];
+ dt = dt_intdom[cell_to_intdom[cid]];
// The 1e-3 is a constant of proportionality required to ensure that the
// conductance (gi) values have units μS (micro-Siemens).
@@ -155,7 +157,8 @@ void assemble_matrix_flat(
const fvm_value_type* cv_capacitance,
const fvm_value_type* area,
const fvm_index_type* cv_to_cell,
- const fvm_value_type* dt_cell,
+ const fvm_value_type* dt_intdom,
+ const fvm_index_type* cell_to_intdom,
unsigned n)
{
constexpr unsigned block_dim = 128;
@@ -165,7 +168,7 @@ void assemble_matrix_flat(
<fvm_value_type, fvm_index_type>
<<<grid_dim, block_dim>>>
(d, rhs, invariant_d, voltage, current, cv_capacitance,
- area, cv_to_cell, dt_cell, n);
+ area, cv_to_cell, dt_intdom, cell_to_intdom, n);
}
//template <typename T, typename I, unsigned BlockWidth, unsigned LoadWidth, unsigned Threads>
@@ -180,7 +183,8 @@ void assemble_matrix_interleaved(
const fvm_index_type* sizes,
const fvm_index_type* starts,
const fvm_index_type* matrix_to_cell,
- const fvm_value_type* dt_cell,
+ const fvm_value_type* dt_intdom,
+ const fvm_index_type* cell_to_intdom,
unsigned padded_size, unsigned num_mtx)
{
constexpr unsigned bd = impl::matrices_per_block();
@@ -195,7 +199,7 @@ void assemble_matrix_interleaved(
<<<grid_dim, block_dim>>>
(d, rhs, invariant_d, voltage, current, cv_capacitance, area,
sizes, starts, matrix_to_cell,
- dt_cell, padded_size, num_mtx);
+ dt_intdom, cell_to_intdom, padded_size, num_mtx);
}
} // namespace gpu
diff --git a/arbor/backends/gpu/matrix_fine.cu b/arbor/backends/gpu/matrix_fine.cu
index 764b39d5..f3ad5372 100644
--- a/arbor/backends/gpu/matrix_fine.cu
+++ b/arbor/backends/gpu/matrix_fine.cu
@@ -53,7 +53,8 @@ void assemble_matrix_fine(
const T* cv_capacitance,
const T* area,
const I* cv_to_cell,
- const T* dt_cell,
+ const T* dt_intdom,
+ const I* cell_to_intdom,
const I* perm,
unsigned n)
{
@@ -61,7 +62,7 @@ void assemble_matrix_fine(
if (tid<n) {
auto cid = cv_to_cell[tid];
- auto dt = dt_cell[cid];
+ auto dt = dt_intdom[cell_to_intdom[cid]];
if (dt>0) {
// The 1e-3 is a constant of proportionality required to ensure that the
@@ -264,7 +265,8 @@ void assemble_matrix_fine(
const fvm_value_type* cv_capacitance,
const fvm_value_type* area,
const fvm_index_type* cv_to_cell,
- const fvm_value_type* dt_cell,
+ const fvm_value_type* dt_intdom,
+ const fvm_index_type* cell_to_intdom,
const fvm_index_type* perm,
unsigned n)
{
@@ -273,7 +275,7 @@ void assemble_matrix_fine(
kernels::assemble_matrix_fine<<<num_blocks, block_dim>>>(
d, rhs, invariant_d, voltage, current, cv_capacitance, area,
- cv_to_cell, dt_cell,
+ cv_to_cell, dt_intdom, cell_to_intdom,
perm, n);
}
diff --git a/arbor/backends/gpu/matrix_fine.hpp b/arbor/backends/gpu/matrix_fine.hpp
index 547b80e2..896bb28b 100644
--- a/arbor/backends/gpu/matrix_fine.hpp
+++ b/arbor/backends/gpu/matrix_fine.hpp
@@ -58,7 +58,8 @@ void assemble_matrix_fine(
const fvm_value_type* cv_capacitance,
const fvm_value_type* area,
const fvm_index_type* cv_to_cell,
- const fvm_value_type* dt_cell,
+ const fvm_value_type* dt_intdom,
+ const fvm_index_type* cell_to_intdom,
const fvm_index_type* perm,
unsigned n);
diff --git a/arbor/backends/gpu/matrix_state_fine.hpp b/arbor/backends/gpu/matrix_state_fine.hpp
index c457fe24..f5c44cbe 100644
--- a/arbor/backends/gpu/matrix_state_fine.hpp
+++ b/arbor/backends/gpu/matrix_state_fine.hpp
@@ -104,6 +104,8 @@ public:
// for storing the solution in unpacked format
array solution_;
+ iarray cell_to_intdom;
+
// the maximum nuber of branches in each level per block
unsigned max_branches_per_level;
@@ -139,7 +141,8 @@ public:
const std::vector<size_type>& cell_cv_divs,
const std::vector<value_type>& cap,
const std::vector<value_type>& face_conductance,
- const std::vector<value_type>& area)
+ const std::vector<value_type>& area,
+ const std::vector<size_type>& cell_intdom)
{
using util::make_span;
constexpr unsigned npos = unsigned(-1);
@@ -421,6 +424,7 @@ public:
// to be stored in flat format
cv_capacitance = memory::make_const_view(cap);
invariant_d = memory::make_const_view(invariant_d_tmp);
+ cell_to_intdom = memory::make_const_view(cell_intdom);
// calculte the cv -> cell mappings
std::vector<size_type> cv_to_cell_tmp(matrix_size);
@@ -434,10 +438,10 @@ public:
// Assemble the matrix
// Afterwards the diagonal and RHS will have been set given dt, voltage and current
- // dt_cell [ms] (per cell)
+ // dt_intdom [ms] (per cell)
// voltage [mV]
// current [nA]
- void assemble(const_view dt_cell, const_view voltage, const_view current) {
+ void assemble(const_view dt_intdom, const_view voltage, const_view current) {
assemble_matrix_fine(
d.data(),
rhs.data(),
@@ -447,7 +451,8 @@ public:
cv_capacitance.data(),
cv_area.data(),
cv_to_cell.data(),
- dt_cell.data(),
+ dt_intdom.data(),
+ cell_to_intdom.data(),
perm.data(),
size());
}
diff --git a/arbor/backends/gpu/matrix_state_flat.hpp b/arbor/backends/gpu/matrix_state_flat.hpp
index e4ddc666..6f4cbd01 100644
--- a/arbor/backends/gpu/matrix_state_flat.hpp
+++ b/arbor/backends/gpu/matrix_state_flat.hpp
@@ -30,7 +30,8 @@ void assemble_matrix_flat(
const fvm_value_type* cv_capacitance,
const fvm_value_type* cv_area,
const fvm_index_type* cv_to_cell,
- const fvm_value_type* dt_cell,
+ const fvm_value_type* dt_intdom,
+ const fvm_index_type* cell_to_intdom,
unsigned n);
/// matrix state
@@ -57,6 +58,8 @@ struct matrix_state_flat {
array face_conductance; // [μS]
array cv_area; // [μm^2]
+ iarray cell_to_intdom;
+
// the invariant part of the matrix diagonal
array invariant_d; // [μS]
@@ -66,7 +69,8 @@ struct matrix_state_flat {
const std::vector<index_type>& cell_cv_divs,
const std::vector<value_type>& cv_cap,
const std::vector<value_type>& face_cond,
- const std::vector<value_type>& area):
+ const std::vector<value_type>& area,
+ const std::vector<index_type>& cell_to_intdom):
parent_index(memory::make_const_view(p)),
cell_cv_divs(memory::make_const_view(cell_cv_divs)),
cv_to_cell(p.size()),
@@ -74,7 +78,8 @@ struct matrix_state_flat {
u(p.size()),
rhs(p.size()),
cv_capacitance(memory::make_const_view(cv_cap)),
- cv_area(memory::make_const_view(area))
+ cv_area(memory::make_const_view(area)),
+ cell_to_intdom(memory::make_const_view(cell_to_intdom))
{
arb_assert(cv_cap.size() == size());
arb_assert(face_cond.size() == size());
@@ -115,15 +120,15 @@ struct matrix_state_flat {
// Assemble the matrix
// Afterwards the diagonal and RHS will have been set given dt, voltage and current.
- // dt_cell [ms] (per cell)
- // voltage [mV]
- // current [nA]
- void assemble(const_view dt_cell, const_view voltage, const_view current) {
+ // dt_intdom [ms] (per integration domain)
+ // voltage [mV]
+ // current [nA]
+ void assemble(const_view dt_intdom, const_view voltage, const_view current) {
// perform assembly on the gpu
assemble_matrix_flat(
d.data(), rhs.data(), invariant_d.data(), voltage.data(),
current.data(), cv_capacitance.data(), cv_area.data(),
- cv_to_cell.data(), dt_cell.data(), size());
+ cv_to_cell.data(), dt_intdom.data(), cell_to_intdom.data(), size());
}
void solve() {
diff --git a/arbor/backends/gpu/matrix_state_interleaved.hpp b/arbor/backends/gpu/matrix_state_interleaved.hpp
index 2ac1989e..864b77eb 100644
--- a/arbor/backends/gpu/matrix_state_interleaved.hpp
+++ b/arbor/backends/gpu/matrix_state_interleaved.hpp
@@ -28,7 +28,8 @@ void assemble_matrix_interleaved(
const fvm_index_type* sizes,
const fvm_index_type* starts,
const fvm_index_type* matrix_to_cell,
- const fvm_value_type* dt_cell,
+ const fvm_value_type* dt_intdom,
+ const fvm_index_type* cell_to_intdom,
unsigned padded_size, unsigned num_mtx);
// host side wrapper for interleaved matrix solver kernel
@@ -122,6 +123,8 @@ struct matrix_state_interleaved {
// the invariant part of the matrix diagonal
array invariant_d; // [μS]
+ iarray cell_to_intdom;
+
// the length of a vector required to store values for one
// matrix with padding
unsigned padded_size;
@@ -144,7 +147,8 @@ struct matrix_state_interleaved {
const std::vector<index_type>& cell_cv_divs,
const std::vector<value_type>& cv_cap,
const std::vector<value_type>& face_cond,
- const std::vector<value_type>& area)
+ const std::vector<value_type>& area,
+ const std::vector<index_type>& cell_intdom)
{
arb_assert(cv_cap.size() == p.size());
arb_assert(face_cond.size() == p.size());
@@ -248,6 +252,7 @@ struct matrix_state_interleaved {
matrix_sizes = memory::make_const_view(sizes_p);
matrix_index = memory::make_const_view(cell_to_cv_p);
matrix_to_cell_index = memory::make_const_view(perm);
+ cell_to_intdom = memory::make_const_view(cell_intdom);
// Allocate space for storing the un-interleaved solution.
solution_ = array(p.size());
@@ -259,16 +264,16 @@ struct matrix_state_interleaved {
// Assemble the matrix
// Afterwards the diagonal and RHS will have been set given dt, voltage and current.
- // dt_cell [ms] (per cell)
- // voltage [mV] (per compartment)
- // current density [A.m^-2] (per compartment)
- void assemble(const_view dt_cell, const_view voltage, const_view current) {
+ // dt_intdom [ms] (per integration domain)
+ // voltage [mV] (per compartment)
+ // current density [A.m^-2] (per compartment)
+ void assemble(const_view dt_intdom, const_view voltage, const_view current) {
assemble_matrix_interleaved
(d.data(), rhs.data(), invariant_d.data(),
voltage.data(), current.data(), cv_capacitance.data(), cv_area.data(),
matrix_sizes.data(), matrix_index.data(),
matrix_to_cell_index.data(),
- dt_cell.data(), padded_matrix_size(), num_matrices());
+ dt_intdom.data(), cell_to_intdom.data(), padded_matrix_size(), num_matrices());
}
diff --git a/arbor/backends/gpu/mechanism.cpp b/arbor/backends/gpu/mechanism.cpp
index 9a70b562..c7c629c4 100644
--- a/arbor/backends/gpu/mechanism.cpp
+++ b/arbor/backends/gpu/mechanism.cpp
@@ -68,7 +68,7 @@ void mechanism::instantiate(unsigned id,
pp->width_ = width_;
- pp->vec_ci_ = shared.cv_to_cell.data();
+ pp->vec_ci_ = shared.cv_to_intdom.data();
pp->vec_t_ = shared.time.data();
pp->vec_t_to_ = shared.time_to.data();
pp->vec_dt_ = shared.dt_cv.data();
diff --git a/arbor/backends/gpu/shared_state.cpp b/arbor/backends/gpu/shared_state.cpp
index 5760b02c..fe499023 100644
--- a/arbor/backends/gpu/shared_state.cpp
+++ b/arbor/backends/gpu/shared_state.cpp
@@ -35,11 +35,9 @@ void update_time_to_impl(
std::size_t n, fvm_value_type* time_to, const fvm_value_type* time,
fvm_value_type dt, fvm_value_type tmax);
-void sync_time_to_impl(std::size_t n, fvm_value_type* time_to, const fvm_index_type* time_deps);
-
void set_dt_impl(
- fvm_size_type ncell, fvm_size_type ncomp, fvm_value_type* dt_cell, fvm_value_type* dt_comp,
- const fvm_value_type* time_to, const fvm_value_type* time, const fvm_index_type* cv_to_cell);
+ fvm_size_type nintdom, fvm_size_type ncomp, fvm_value_type* dt_intdom, fvm_value_type* dt_comp,
+ const fvm_value_type* time_to, const fvm_value_type* time, const fvm_index_type* cv_to_intdom);
void add_gj_current_impl(
fvm_size_type n_gj, const fvm_gap_junction* gj, const fvm_value_type* v, fvm_value_type* i);
@@ -112,26 +110,24 @@ void ion_state::zero_current() {
// Shared state methods:
shared_state::shared_state(
- fvm_size_type n_cell,
- const std::vector<fvm_index_type>& cv_to_cell_vec,
- const std::vector<fvm_index_type>& time_dep_vec,
+ fvm_size_type n_intdom,
+ const std::vector<fvm_index_type>& cv_to_intdom_vec,
const std::vector<fvm_gap_junction>& gj_vec,
unsigned // alignment parameter ignored.
):
- n_cell(n_cell),
- n_cv(cv_to_cell_vec.size()),
+ n_intdom(n_intdom),
+ n_cv(cv_to_intdom_vec.size()),
n_gj(gj_vec.size()),
- cv_to_cell(make_const_view(cv_to_cell_vec)),
- time_dep(make_const_view(time_dep_vec)),
+ cv_to_intdom(make_const_view(cv_to_intdom_vec)),
gap_junctions(make_const_view(gj_vec)),
- time(n_cell),
- time_to(n_cell),
- dt_cell(n_cell),
+ time(n_intdom),
+ time_to(n_intdom),
+ dt_intdom(n_intdom),
dt_cv(n_cv),
voltage(n_cv),
current_density(n_cv),
temperature_degC(1),
- deliverable_events(n_cell)
+ deliverable_events(n_intdom)
{}
void shared_state::add_ion(
@@ -177,15 +173,11 @@ void shared_state::ions_nernst_reversal_potential(fvm_value_type temperature_K)
}
void shared_state::update_time_to(fvm_value_type dt_step, fvm_value_type tmax) {
- update_time_to_impl(n_cell, time_to.data(), time.data(), dt_step, tmax);
-}
-
-void shared_state::sync_time_to() {
- sync_time_to_impl(n_cell, time_to.data(), time_dep.data());
+ update_time_to_impl(n_intdom, time_to.data(), time.data(), dt_step, tmax);
}
void shared_state::set_dt() {
- set_dt_impl(n_cell, n_cv, dt_cell.data(), dt_cv.data(), time_to.data(), time.data(), cv_to_cell.data());
+ set_dt_impl(n_intdom, n_cv, dt_intdom.data(), dt_cv.data(), time_to.data(), time.data(), cv_to_intdom.data());
}
void shared_state::add_gj_current() {
@@ -193,7 +185,7 @@ void shared_state::add_gj_current() {
}
std::pair<fvm_value_type, fvm_value_type> shared_state::time_bounds() const {
- return minmax_value_impl(n_cell, time.data());
+ return minmax_value_impl(n_intdom, time.data());
}
std::pair<fvm_value_type, fvm_value_type> shared_state::voltage_bounds() const {
@@ -206,10 +198,10 @@ void shared_state::take_samples(const sample_event_stream::state& s, array& samp
// Debug interface
std::ostream& operator<<(std::ostream& o, shared_state& s) {
- o << " cv_to_cell " << s.cv_to_cell << "\n";
+ o << " cv_to_intdom " << s.cv_to_intdom << "\n";
o << " time " << s.time << "\n";
o << " time_to " << s.time_to << "\n";
- o << " dt_cell " << s.dt_cell << "\n";
+ o << " dt_intdom " << s.dt_intdom << "\n";
o << " dt_cv " << s.dt_cv << "\n";
o << " voltage " << s.voltage << "\n";
o << " current " << s.current_density << "\n";
diff --git a/arbor/backends/gpu/shared_state.cu b/arbor/backends/gpu/shared_state.cu
index 806418b9..a5901cf5 100644
--- a/arbor/backends/gpu/shared_state.cu
+++ b/arbor/backends/gpu/shared_state.cu
@@ -34,24 +34,6 @@ void init_concentration_impl(unsigned n, T* Xi, T* Xo, const T* weight_Xi, const
}
}
-template <typename T, typename I>
-__global__ void sync_time_to_impl(unsigned n, T* time_to, const I* time_deps) {
- unsigned i = threadIdx.x+blockIdx.x*blockDim.x;
- if (i<n) {
- if (time_deps[i] > 0) {
- auto min_t = time_to[i];
- for (int j = 1; j < time_deps[i]; j++) {
- if (time_to[i+j] < min_t) {
- min_t = time_to[i+j];
- }
- }
- for (int j = 0; j < time_deps[i]; j++) {
- time_to[i+j] = min_t;
- }
- }
- }
-}
-
template <typename T>
__global__ void update_time_to_impl(unsigned n, T* time_to, const T* time, T dt, T tmax) {
unsigned i = threadIdx.x+blockIdx.x*blockDim.x;
@@ -131,15 +113,6 @@ void init_concentration_impl(
kernel::init_concentration_impl<<<nblock, block_dim>>>(n, Xi, Xo, weight_Xi, weight_Xo, c_int, c_ext);
}
-void sync_time_to_impl(std::size_t n, fvm_value_type* time_to, const fvm_index_type* time_deps)
-{
- if (!n) return;
-
- constexpr int block_dim = 128;
- int nblock = block_count(n, block_dim);
- kernel::sync_time_to_impl<<<nblock, block_dim>>>(n, time_to, time_deps);
-}
-
void update_time_to_impl(
std::size_t n, fvm_value_type* time_to, const fvm_value_type* time,
fvm_value_type dt, fvm_value_type tmax)
@@ -152,17 +125,17 @@ void update_time_to_impl(
}
void set_dt_impl(
- fvm_size_type ncell, fvm_size_type ncomp, fvm_value_type* dt_cell, fvm_value_type* dt_comp,
- const fvm_value_type* time_to, const fvm_value_type* time, const fvm_index_type* cv_to_cell)
+ fvm_size_type nintdom, fvm_size_type ncomp, fvm_value_type* dt_intdom, fvm_value_type* dt_comp,
+ const fvm_value_type* time_to, const fvm_value_type* time, const fvm_index_type* cv_to_intdom)
{
- if (!ncell || !ncomp) return;
+ if (!nintdom || !ncomp) return;
constexpr int block_dim = 128;
- int nblock = block_count(ncell, block_dim);
- kernel::vec_minus<<<nblock, block_dim>>>(ncell, dt_cell, time_to, time);
+ int nblock = block_count(nintdom, block_dim);
+ kernel::vec_minus<<<nblock, block_dim>>>(nintdom, dt_intdom, time_to, time);
nblock = block_count(ncomp, block_dim);
- kernel::gather<<<nblock, block_dim>>>(ncomp, dt_comp, dt_cell, cv_to_cell);
+ kernel::gather<<<nblock, block_dim>>>(ncomp, dt_comp, dt_intdom, cv_to_intdom);
}
void add_gj_current_impl(
diff --git a/arbor/backends/gpu/shared_state.hpp b/arbor/backends/gpu/shared_state.hpp
index c66ac3f1..778d82b2 100644
--- a/arbor/backends/gpu/shared_state.hpp
+++ b/arbor/backends/gpu/shared_state.hpp
@@ -65,16 +65,15 @@ struct ion_state {
};
struct shared_state {
- fvm_size_type n_cell = 0; // Number of distinct cells (integration domains).
+ fvm_size_type n_intdom = 0; // Number of distinct integration domains.
fvm_size_type n_cv = 0; // Total number of CVs.
fvm_size_type n_gj = 0; // Total number of GJs.
- iarray cv_to_cell; // Maps CV index to cell index.
- iarray time_dep; // Provides information about supercells
+ iarray cv_to_intdom; // Maps CV index to intdom index.
gjarray gap_junctions; // Stores gap_junction info.
- array time; // Maps cell index to integration start time [ms].
- array time_to; // Maps cell index to integration stop time [ms].
- array dt_cell; // Maps cell index to (stop time) - (start time) [ms].
+ array time; // Maps intdom index to integration start time [ms].
+ array time_to; // Maps intdom index to integration stop time [ms].
+ array dt_intdom; // Maps intdom index to (stop time) - (start time) [ms].
array dt_cv; // Maps CV index to dt [ms].
array voltage; // Maps CV index to membrane voltage [mV].
array current_density; // Maps CV index to current density [A/m²].
@@ -87,9 +86,8 @@ struct shared_state {
shared_state() = default;
shared_state(
- fvm_size_type n_cell,
- const std::vector<fvm_index_type>& cv_to_cell_vec,
- const std::vector<fvm_index_type>& time_dep_vec,
+ fvm_size_type n_intdom,
+ const std::vector<fvm_index_type>& cv_to_intdom_vec,
const std::vector<fvm_gap_junction>& gj_vec,
unsigned align
);
@@ -109,10 +107,7 @@ struct shared_state {
// Set time_to to earliest of time+dt_step and tmax.
void update_time_to(fvm_value_type dt_step, fvm_value_type tmax);
- // Synchrnize the time_to for supercells.
- void sync_time_to();
-
- // Set the per-cell and per-compartment dt from time_to - time.
+ // Set the per-intdom and per-compartment dt from time_to - time.
void set_dt();
// Update gap_junction state
diff --git a/arbor/backends/gpu/threshold_watcher.cu b/arbor/backends/gpu/threshold_watcher.cu
index c18022ee..6377f009 100644
--- a/arbor/backends/gpu/threshold_watcher.cu
+++ b/arbor/backends/gpu/threshold_watcher.cu
@@ -30,7 +30,7 @@ inline T lerp(T a, T b, T u) {
__global__
void test_thresholds_impl(
int size,
- const fvm_index_type* cv_to_cell, const fvm_value_type* t_after, const fvm_value_type* t_before,
+ const fvm_index_type* cv_to_intdom, const fvm_value_type* t_after, const fvm_value_type* t_before,
stack_storage<threshold_crossing>& stack,
fvm_index_type* is_crossed, fvm_value_type* prev_values,
const fvm_index_type* cv_index, const fvm_value_type* values, const fvm_value_type* thresholds)
@@ -43,7 +43,7 @@ void test_thresholds_impl(
if (i<size) {
// Test for threshold crossing
const auto cv = cv_index[i];
- const auto cell = cv_to_cell[cv];
+ const auto cell = cv_to_intdom[cv];
const auto v_prev = prev_values[i];
const auto v = values[cv];
const auto thresh = thresholds[i];
@@ -86,7 +86,7 @@ extern void reset_crossed_impl(
void test_thresholds_impl(
int size,
- const fvm_index_type* cv_to_cell, const fvm_value_type* t_after, const fvm_value_type* t_before,
+ const fvm_index_type* cv_to_intdom, const fvm_value_type* t_after, const fvm_value_type* t_before,
stack_storage<threshold_crossing>& stack,
fvm_index_type* is_crossed, fvm_value_type* prev_values,
const fvm_index_type* cv_index, const fvm_value_type* values, const fvm_value_type* thresholds)
@@ -95,7 +95,7 @@ void test_thresholds_impl(
constexpr int block_dim = 128;
const int grid_dim = impl::block_count(size, block_dim);
kernel::test_thresholds_impl<<<grid_dim, block_dim>>>(
- size, cv_to_cell, t_after, t_before, stack, is_crossed, prev_values, cv_index, values, thresholds);
+ size, cv_to_intdom, t_after, t_before, stack, is_crossed, prev_values, cv_index, values, thresholds);
}
}
diff --git a/arbor/backends/gpu/threshold_watcher.hpp b/arbor/backends/gpu/threshold_watcher.hpp
index e0668a9f..7b1ad17b 100644
--- a/arbor/backends/gpu/threshold_watcher.hpp
+++ b/arbor/backends/gpu/threshold_watcher.hpp
@@ -22,7 +22,7 @@ namespace gpu {
void test_thresholds_impl(
int size,
- const fvm_index_type* cv_to_cell, const fvm_value_type* t_after, const fvm_value_type* t_before,
+ const fvm_index_type* cv_to_intdom, const fvm_value_type* t_after, const fvm_value_type* t_before,
stack_storage<threshold_crossing>& stack,
fvm_index_type* is_crossed, fvm_value_type* prev_values,
const fvm_index_type* cv_index, const fvm_value_type* values, const fvm_value_type* thresholds);
@@ -44,7 +44,7 @@ public:
threshold_watcher(const execution_context& ctx): stack_(ctx.gpu) {}
threshold_watcher(
- const fvm_index_type* cv_to_cell,
+ const fvm_index_type* cv_to_intdom,
const fvm_value_type* t_before,
const fvm_value_type* t_after,
const fvm_value_type* values,
@@ -52,7 +52,7 @@ public:
const std::vector<fvm_value_type>& thresholds,
const execution_context& ctx
):
- cv_to_cell_(cv_to_cell),
+ cv_to_intdom_(cv_to_intdom),
t_before_(t_before),
t_after_(t_after),
values_(values),
@@ -109,7 +109,7 @@ public:
if (size()>0) {
test_thresholds_impl(
(int)size(),
- cv_to_cell_, t_after_, t_before_,
+ cv_to_intdom_, t_after_, t_before_,
stack_.storage(),
is_crossed_.data(), v_prev_.data(),
cv_index_.data(), values_, thresholds_.data());
@@ -129,7 +129,7 @@ public:
private:
/// Non-owning pointers to gpu-side cv-to-cell map, per-cell time data,
/// and the values for to test against thresholds.
- const fvm_index_type* cv_to_cell_ = nullptr;
+ const fvm_index_type* cv_to_intdom_ = nullptr;
const fvm_value_type* t_before_ = nullptr;
const fvm_value_type* t_after_ = nullptr;
const fvm_value_type* values_ = nullptr;
diff --git a/arbor/backends/multicore/fvm.hpp b/arbor/backends/multicore/fvm.hpp
index 1cbae3ef..1f9c019d 100644
--- a/arbor/backends/multicore/fvm.hpp
+++ b/arbor/backends/multicore/fvm.hpp
@@ -51,7 +51,7 @@ struct backend {
const execution_context& context)
{
return threshold_watcher(
- state.cv_to_cell.data(),
+ state.cv_to_intdom.data(),
state.time.data(),
state.time_to.data(),
state.voltage.data(),
diff --git a/arbor/backends/multicore/matrix_state.hpp b/arbor/backends/multicore/matrix_state.hpp
index a4d318aa..9e54b8eb 100644
--- a/arbor/backends/multicore/matrix_state.hpp
+++ b/arbor/backends/multicore/matrix_state.hpp
@@ -29,6 +29,8 @@ public:
array face_conductance; // [μS]
array cv_area; // [μm^2]
+ iarray cell_to_intdom;
+
// the invariant part of the matrix diagonal
array invariant_d; // [μS]
@@ -38,13 +40,15 @@ public:
const std::vector<index_type>& cell_cv_divs,
const std::vector<value_type>& cap,
const std::vector<value_type>& cond,
- const std::vector<value_type>& area):
+ const std::vector<value_type>& area,
+ const std::vector<index_type>& cell_to_intdom):
parent_index(p.begin(), p.end()),
cell_cv_divs(cell_cv_divs.begin(), cell_cv_divs.end()),
d(size(), 0), u(size(), 0), rhs(size()),
cv_capacitance(cap.begin(), cap.end()),
face_conductance(cond.begin(), cond.end()),
- cv_area(area.begin(), area.end())
+ cv_area(area.begin(), area.end()),
+ cell_to_intdom(cell_to_intdom.begin(), cell_to_intdom.end())
{
arb_assert(cap.size() == size());
arb_assert(cond.size() == size());
@@ -70,16 +74,16 @@ public:
// Assemble the matrix
// Afterwards the diagonal and RHS will have been set given dt, voltage and current.
- // dt_cell [ms] (per cell)
+ // dt_intdom [ms] (per integration domain)
// voltage [mV] (per compartment)
// current density [A.m^-2] (per compartment)
- void assemble(const_view dt_cell, const_view voltage, const_view current) {
+ void assemble(const_view dt_intdom, const_view voltage, const_view current) {
auto cell_cv_part = util::partition_view(cell_cv_divs);
const index_type ncells = cell_cv_part.size();
// loop over submatrices
for (auto m: util::make_span(0, ncells)) {
- auto dt = dt_cell[m];
+ auto dt = dt_intdom[cell_to_intdom[m]];
if (dt>0) {
value_type factor = 1e-3/dt;
diff --git a/arbor/backends/multicore/mechanism.cpp b/arbor/backends/multicore/mechanism.cpp
index 63bd3449..5dc5c618 100644
--- a/arbor/backends/multicore/mechanism.cpp
+++ b/arbor/backends/multicore/mechanism.cpp
@@ -70,7 +70,7 @@ void mechanism::instantiate(unsigned id, backend::shared_state& shared, const la
// Assign non-owning views onto shared state:
- vec_ci_ = shared.cv_to_cell.data();
+ vec_ci_ = shared.cv_to_intdom.data();
vec_t_ = shared.time.data();
vec_t_to_ = shared.time_to.data();
vec_dt_ = shared.dt_cv.data();
diff --git a/arbor/backends/multicore/shared_state.cpp b/arbor/backends/multicore/shared_state.cpp
index 3bb73441..6313d429 100644
--- a/arbor/backends/multicore/shared_state.cpp
+++ b/arbor/backends/multicore/shared_state.cpp
@@ -115,35 +115,31 @@ void ion_state::zero_current() {
// shared_state methods:
shared_state::shared_state(
- fvm_size_type n_cell,
- const std::vector<fvm_index_type>& cv_to_cell_vec,
- const std::vector<fvm_index_type>& time_dep_vec,
+ fvm_size_type n_intdom,
+ const std::vector<fvm_index_type>& cv_to_intdom_vec,
const std::vector<fvm_gap_junction>& gj_vec,
unsigned align
):
alignment(min_alignment(align)),
alloc(alignment),
- n_cell(n_cell),
- n_cv(cv_to_cell_vec.size()),
+ n_intdom(n_intdom),
+ n_cv(cv_to_intdom_vec.size()),
n_gj(gj_vec.size()),
- cv_to_cell(math::round_up(n_cv, alignment), pad(alignment)),
- time_dep(n_cell),
+ cv_to_intdom(math::round_up(n_cv, alignment), pad(alignment)),
gap_junctions(math::round_up(n_gj, alignment), pad(alignment)),
- time(n_cell, pad(alignment)),
- time_to(n_cell, pad(alignment)),
- dt_cell(n_cell, pad(alignment)),
+ time(n_intdom, pad(alignment)),
+ time_to(n_intdom, pad(alignment)),
+ dt_intdom(n_intdom, pad(alignment)),
dt_cv(n_cv, pad(alignment)),
voltage(n_cv, pad(alignment)),
current_density(n_cv, pad(alignment)),
temperature_degC(NAN),
- deliverable_events(n_cell)
+ deliverable_events(n_intdom)
{
- std::copy(time_dep_vec.begin(), time_dep_vec.end(), time_dep.begin());
-
- // For indices in the padded tail of cv_to_cell, set index to last valid cell index.
+ // For indices in the padded tail of cv_to_intdom, set index to last valid intdom index.
if (n_cv>0) {
- std::copy(cv_to_cell_vec.begin(), cv_to_cell_vec.end(), cv_to_cell.begin());
- std::fill(cv_to_cell.begin() + n_cv, cv_to_cell.end(), cv_to_cell_vec.back());
+ std::copy(cv_to_intdom_vec.begin(), cv_to_intdom_vec.end(), cv_to_intdom.begin());
+ std::fill(cv_to_intdom.begin() + n_cv, cv_to_intdom.end(), cv_to_intdom_vec.back());
}
if (n_gj>0) {
std::copy(gj_vec.begin(), gj_vec.end(), gap_junctions.begin());
@@ -192,24 +188,9 @@ void shared_state::ions_nernst_reversal_potential(fvm_value_type temperature_K)
i.second.nernst(temperature_K);
}
}
-void shared_state::sync_time_to() {
- for (fvm_size_type i = 0; i<n_cell; i++) {
- if (!time_dep[i]) continue;
-
- fvm_value_type min_t = time_to[i];
- for (int j = 1; j < time_dep[i]; j++) {
- if (time_to[i+j] < min_t) {
- min_t = time_to[i+j];
- }
- }
- for (int j = 0; j < time_dep[i]; j++) {
- time_to[i+j] = min_t;
- }
- }
-}
void shared_state::update_time_to(fvm_value_type dt_step, fvm_value_type tmax) {
- for (fvm_size_type i = 0; i<n_cell; i+=simd_width) {
+ for (fvm_size_type i = 0; i<n_intdom; i+=simd_width) {
simd_value_type t(time.data()+i);
t = min(t+dt_step, simd_value_type(tmax));
t.copy_to(time_to.data()+i);
@@ -217,18 +198,18 @@ void shared_state::update_time_to(fvm_value_type dt_step, fvm_value_type tmax) {
}
void shared_state::set_dt() {
- for (fvm_size_type j = 0; j<n_cell; j+=simd_width) {
+ for (fvm_size_type j = 0; j<n_intdom; j+=simd_width) {
simd_value_type t(time.data()+j);
simd_value_type t_to(time_to.data()+j);
auto dt = t_to-t;
- dt.copy_to(dt_cell.data()+j);
+ dt.copy_to(dt_intdom.data()+j);
}
for (fvm_size_type i = 0; i<n_cv; i+=simd_width) {
- simd_index_type cell_idx(cv_to_cell.data()+i);
+ simd_index_type intdom_idx(cv_to_intdom.data()+i);
- simd_value_type dt(simd::indirect(dt_cell.data(), cell_idx));
+ simd_value_type dt(simd::indirect(dt_intdom.data(), intdom_idx));
dt.copy_to(dt_cv.data()+i);
}
}
@@ -272,12 +253,12 @@ void shared_state::take_samples(
std::ostream& operator<<(std::ostream& out, const shared_state& s) {
using io::csv;
- out << "n_cell " << s.n_cell << "\n";
+ out << "n_intdom " << s.n_intdom << "\n";
out << "n_cv " << s.n_cv << "\n";
- out << "cv_to_cell " << csv(s.cv_to_cell) << "\n";
+ out << "cv_to_intdom " << csv(s.cv_to_intdom) << "\n";
out << "time " << csv(s.time) << "\n";
out << "time_to " << csv(s.time_to) << "\n";
- out << "dt_cell " << csv(s.dt_cell) << "\n";
+ out << "dt_intdom " << csv(s.dt_intdom) << "\n";
out << "dt_cv " << csv(s.dt_cv) << "\n";
out << "voltage " << csv(s.voltage) << "\n";
out << "current " << csv(s.current_density) << "\n";
diff --git a/arbor/backends/multicore/shared_state.hpp b/arbor/backends/multicore/shared_state.hpp
index 40be9eb2..fc2c0e7f 100644
--- a/arbor/backends/multicore/shared_state.hpp
+++ b/arbor/backends/multicore/shared_state.hpp
@@ -84,16 +84,15 @@ struct shared_state {
unsigned alignment = 1; // Alignment and padding multiple.
util::padded_allocator<> alloc; // Allocator with corresponging alignment/padding.
- fvm_size_type n_cell = 0; // Number of distinct cells (integration domains).
+ fvm_size_type n_intdom = 0; // Number of integration domains.
fvm_size_type n_cv = 0; // Total number of CVs.
fvm_size_type n_gj = 0; // Total number of GJs.
- iarray cv_to_cell; // Maps CV index to cell index.
- iarray time_dep; // Provides information about supercells
+ iarray cv_to_intdom; // Maps CV index to integration domain index.
gjarray gap_junctions; // Stores gap_junction info.
- array time; // Maps cell index to integration start time [ms].
- array time_to; // Maps cell index to integration stop time [ms].
- array dt_cell; // Maps cell index to (stop time) - (start time) [ms].
+ array time; // Maps intdom index to integration start time [ms].
+ array time_to; // Maps intdom index to integration stop time [ms].
+ array dt_intdom; // Maps index to (stop time) - (start time) [ms].
array dt_cv; // Maps CV index to dt [ms].
array voltage; // Maps CV index to membrane voltage [mV].
array current_density; // Maps CV index to current density [A/m²].
@@ -106,9 +105,8 @@ struct shared_state {
shared_state() = default;
shared_state(
- fvm_size_type n_cell,
- const std::vector<fvm_index_type>& cv_to_cell_vec,
- const std::vector<fvm_index_type>& time_dep_vec,
+ fvm_size_type n_intdom,
+ const std::vector<fvm_index_type>& cv_to_intdom_vec,
const std::vector<fvm_gap_junction>& gj_vec,
unsigned align
);
@@ -128,10 +126,7 @@ struct shared_state {
// Set time_to to earliest of time+dt_step and tmax.
void update_time_to(fvm_value_type dt_step, fvm_value_type tmax);
- // Synchrnize the time_to for supercells.
- void sync_time_to();
-
- // Set the per-cell and per-compartment dt from time_to - time.
+ // Set the per-integration domain and per-compartment dt from time_to - time.
void set_dt();
// Update gap_junction state
diff --git a/arbor/backends/multicore/threshold_watcher.hpp b/arbor/backends/multicore/threshold_watcher.hpp
index 68b1d61d..e356f486 100644
--- a/arbor/backends/multicore/threshold_watcher.hpp
+++ b/arbor/backends/multicore/threshold_watcher.hpp
@@ -18,7 +18,7 @@ public:
threshold_watcher(const execution_context& ctx) {}
threshold_watcher(
- const fvm_index_type* cv_to_cell,
+ const fvm_index_type* cv_to_intdom,
const fvm_value_type* t_before,
const fvm_value_type* t_after,
const fvm_value_type* values,
@@ -26,7 +26,7 @@ public:
const std::vector<fvm_value_type>& thresholds,
const execution_context& context
):
- cv_to_cell_(cv_to_cell),
+ cv_to_intdom_(cv_to_intdom),
t_before_(t_before),
t_after_(t_after),
values_(values),
@@ -66,7 +66,7 @@ public:
void test() {
for (fvm_size_type i = 0; i<n_cv_; ++i) {
auto cv = cv_index_[i];
- auto cell = cv_to_cell_[cv];
+ auto cell = cv_to_intdom_[cv];
auto v_prev = v_prev_[i];
auto v = values_[cv];
auto thresh = thresholds_[i];
@@ -104,7 +104,7 @@ public:
private:
/// Non-owning pointers to cv-to-cell map, per-cell time data,
/// and the values for to test against thresholds.
- const fvm_index_type* cv_to_cell_ = nullptr;
+ const fvm_index_type* cv_to_intdom_ = nullptr;
const fvm_value_type* t_before_ = nullptr;
const fvm_value_type* t_after_ = nullptr;
const fvm_value_type* values_ = nullptr;
diff --git a/arbor/fvm_layout.cpp b/arbor/fvm_layout.cpp
index 9d2d1607..24c69019 100644
--- a/arbor/fvm_layout.cpp
+++ b/arbor/fvm_layout.cpp
@@ -111,6 +111,7 @@ namespace {
//
fvm_discretization fvm_discretize(const std::vector<mc_cell>& cells) {
+
using value_type = fvm_value_type;
using index_type = fvm_index_type;
using size_type = fvm_size_type;
diff --git a/arbor/fvm_lowered_cell.hpp b/arbor/fvm_lowered_cell.hpp
index e49c610b..5f77ca1e 100644
--- a/arbor/fvm_lowered_cell.hpp
+++ b/arbor/fvm_lowered_cell.hpp
@@ -28,8 +28,8 @@ struct fvm_lowered_cell {
virtual void initialize(
const std::vector<cell_gid_type>& gids,
- const std::vector<int>& deps,
const recipe& rec,
+ std::vector<fvm_index_type>& cell_to_intdom,
std::vector<target_handle>& target_handles,
probe_association_map<probe_handle>& probe_map) = 0;
diff --git a/arbor/fvm_lowered_cell_impl.hpp b/arbor/fvm_lowered_cell_impl.hpp
index 31edf2bf..70c13263 100644
--- a/arbor/fvm_lowered_cell_impl.hpp
+++ b/arbor/fvm_lowered_cell_impl.hpp
@@ -10,9 +10,11 @@
#include <cmath>
#include <iterator>
#include <memory>
+#include <queue>
#include <stdexcept>
#include <utility>
#include <vector>
+#include <unordered_set>
#include <arbor/assert.hpp>
#include <arbor/common_types.hpp>
@@ -50,8 +52,8 @@ public:
void initialize(
const std::vector<cell_gid_type>& gids,
- const std::vector<int>& deps,
const recipe& rec,
+ std::vector<fvm_index_type>& cell_to_intdom,
std::vector<target_handle>& target_handles,
probe_association_map<probe_handle>& probe_map) override;
@@ -67,6 +69,13 @@ public:
const recipe& rec,
const fvm_discretization& D);
+ // Generates indom index for every gid, guarantees that gids belonging to the same supercell are in the same intdom
+ // Fills cell_to_intdom map; returns number of intdoms
+ fvm_size_type fvm_intdom(
+ const recipe& rec,
+ const std::vector<cell_gid_type>& gids,
+ std::vector<fvm_index_type>& cell_to_intdom);
+
value_type time() const override { return tmin_; }
//Exposed for testing purposes
@@ -215,7 +224,6 @@ fvm_integration_result fvm_lowered_cell_impl<Backend>::integrate(
state_->update_time_to(dt_max, tfinal);
state_->deliverable_events.event_time_if_before(state_->time_to);
- state_->sync_time_to();
state_->set_dt();
PL();
@@ -230,7 +238,7 @@ fvm_integration_result fvm_lowered_cell_impl<Backend>::integrate(
// Integrate voltage by matrix solve.
PE(advance_integrate_matrix_build);
- matrix_.assemble(state_->dt_cell, state_->voltage, state_->current_density);
+ matrix_.assemble(state_->dt_intdom, state_->voltage, state_->current_density);
PL();
PE(advance_integrate_matrix_solve);
matrix_.solve();
@@ -312,8 +320,8 @@ void fvm_lowered_cell_impl<B>::assert_voltage_bounded(fvm_value_type bound) {
template <typename B>
void fvm_lowered_cell_impl<B>::initialize(
const std::vector<cell_gid_type>& gids,
- const std::vector<int>& deps,
const recipe& rec,
+ std::vector<fvm_index_type>& cell_to_intdom,
std::vector<target_handle>& target_handles,
probe_association_map<probe_handle>& probe_map)
{
@@ -363,12 +371,19 @@ void fvm_lowered_cell_impl<B>::initialize(
check_voltage_mV = global_props.membrane_voltage_limit_mV;
+ auto num_intdoms = fvm_intdom(rec, gids, cell_to_intdom);
+
// Discretize cells, build matrix.
fvm_discretization D = fvm_discretize(cells);
+
+ std::vector<index_type> cv_to_intdom(D.ncomp);
+ std::transform(D.cv_to_cell.begin(), D.cv_to_cell.end(), cv_to_intdom.begin(),
+ [&cell_to_intdom](index_type i){ return cell_to_intdom[i]; });
+
arb_assert(D.ncell == ncell);
- matrix_ = matrix<backend>(D.parent_cv, D.cell_cv_bounds, D.cv_capacitance, D.face_conductance, D.cv_area);
- sample_events_ = sample_event_stream(ncell);
+ matrix_ = matrix<backend>(D.parent_cv, D.cell_cv_bounds, D.cv_capacitance, D.face_conductance, D.cv_area, cell_to_intdom);
+ sample_events_ = sample_event_stream(num_intdoms);
// Discretize mechanism data.
@@ -385,7 +400,7 @@ void fvm_lowered_cell_impl<B>::initialize(
util::transform_view(keys(mech_data.mechanisms),
[&](const std::string& name) { return mech_instance(name)->data_alignment(); }));
- state_ = std::make_unique<shared_state>(ncell, D.cv_to_cell, deps, gj_vector, data_alignment? data_alignment: 1u);
+ state_ = std::make_unique<shared_state>(num_intdoms, cv_to_intdom, gj_vector, data_alignment? data_alignment: 1u);
// Instantiate mechanisms and ions.
@@ -422,7 +437,7 @@ void fvm_lowered_cell_impl<B>::initialize(
// (builtin stimulus, for example, has no targets)
if (!config.target.empty()) {
- target_handles[config.target[i]] = target_handle(mech_id, i, D.cv_to_cell[cv]);
+ target_handles[config.target[i]] = target_handle(mech_id, i, cv_to_intdom[cv]);
}
}
}
@@ -531,4 +546,54 @@ std::vector<fvm_gap_junction> fvm_lowered_cell_impl<B>::fvm_gap_junctions(
return v;
}
+template <typename B>
+fvm_size_type fvm_lowered_cell_impl<B>::fvm_intdom(
+ const recipe& rec,
+ const std::vector<cell_gid_type>& gids,
+ std::vector<fvm_index_type>& cell_to_intdom) {
+
+ cell_to_intdom.resize(gids.size());
+
+ std::unordered_map<cell_gid_type, cell_size_type> gid_to_loc;
+ for (auto i: util::count_along(gids)) {
+ gid_to_loc[gids[i]] = i;
+ }
+
+ std::unordered_set<cell_gid_type> visited;
+ std::queue<cell_gid_type> intdomq;
+ cell_size_type intdom_id = 0;
+
+ for (auto gid: gids) {
+ if (visited.count(gid)) continue;
+ visited.insert(gid);
+
+ intdomq.push(gid);
+ while (!intdomq.empty()) {
+ auto g = intdomq.front();
+ intdomq.pop();
+
+ cell_to_intdom[gid_to_loc[g]] = intdom_id;
+
+ for (auto gj: rec.gap_junctions_on(g)) {
+ cell_gid_type peer =
+ gj.local.gid==g? gj.peer.gid:
+ gj.peer.gid==g? gj.local.gid:
+ throw bad_cell_description(cell_kind::cable1d_neuron, g);
+
+ if (!gid_to_loc.count(peer)) {
+ throw gj_unsupported_domain_decomposition(g, peer);
+ }
+
+ if (!visited.count(peer)) {
+ visited.insert(peer);
+ intdomq.push(peer);
+ }
+ }
+ }
+ intdom_id++;
+ }
+
+ return intdom_id;
+}
+
} // namespace arb
diff --git a/arbor/matrix.hpp b/arbor/matrix.hpp
index bab4c6ef..94958c00 100644
--- a/arbor/matrix.hpp
+++ b/arbor/matrix.hpp
@@ -35,10 +35,12 @@ public:
const std::vector<index_type>& ci,
const std::vector<value_type>& cv_capacitance,
const std::vector<value_type>& face_conductance,
- const std::vector<value_type>& cv_area):
+ const std::vector<value_type>& cv_area,
+ const std::vector<index_type>& cell_to_intdom):
parent_index_(pi.begin(), pi.end()),
cell_index_(ci.begin(), ci.end()),
- state_(pi, ci, cv_capacitance, face_conductance, cv_area)
+ cell_to_intdom_(cell_to_intdom.begin(), cell_to_intdom.end()),
+ state_(pi, ci, cv_capacitance, face_conductance, cv_area, cell_to_intdom)
{
arb_assert(cell_index_[num_cells()] == index_type(parent_index_.size()));
}
@@ -81,6 +83,9 @@ private:
/// indexes that point to the start of each cell in the index
iarray cell_index_;
+ /// indexes that point to the index of the integration domain
+ iarray cell_to_intdom_;
+
public:
// Provide via public interface to make testing much easier.
// If you modify this directly without knowing what you are doing,
diff --git a/arbor/mc_cell_group.cpp b/arbor/mc_cell_group.cpp
index ad6673ff..ec8371ab 100644
--- a/arbor/mc_cell_group.cpp
+++ b/arbor/mc_cell_group.cpp
@@ -24,10 +24,12 @@
namespace arb {
+ARB_DEFINE_LEXICOGRAPHIC_ORDERING(arb::target_handle,(a.mech_id,a.mech_index,a.intdom_index),(b.mech_id,b.mech_index,b.intdom_index))
+ARB_DEFINE_LEXICOGRAPHIC_ORDERING(arb::deliverable_event,(a.time,a.handle,a.weight),(b.time,b.handle,b.weight))
+
mc_cell_group::mc_cell_group(const std::vector<cell_gid_type>& gids, const recipe& rec, fvm_lowered_cell_ptr lowered):
- lowered_(std::move(lowered))
+ gids_(gids), lowered_(std::move(lowered))
{
- generate_deps_gids(rec, gids);
// Default to no binning of events
set_binning_policy(binning_kind::none, 0);
@@ -47,7 +49,7 @@ mc_cell_group::mc_cell_group(const std::vector<cell_gid_type>& gids, const recip
target_handles_.reserve(n_targets);
// Construct cell implementation, retrieving handles and maps.
- lowered_->initialize(gids_, deps_, rec, target_handles_, probe_map_);
+ lowered_->initialize(gids_, rec, cell_to_intdom_, target_handles_, probe_map_);
// Create a list of the global identifiers for the spike sources
for (auto source_gid: gids_) {
@@ -58,58 +60,6 @@ mc_cell_group::mc_cell_group(const std::vector<cell_gid_type>& gids, const recip
spike_sources_.shrink_to_fit();
}
-// Fills gids_ and deps_: gids_ are sorted such that members of the same supercell are consecutive
-void mc_cell_group::generate_deps_gids(const recipe& rec, std::vector<cell_gid_type> gids) {
-
- std::unordered_map<cell_gid_type, cell_size_type> gid_to_loc;
- for (auto i: util::count_along(gids)) {
- gid_to_loc[gids[i]] = i;
- }
-
- deps_.reserve(gids_.size());
- std::unordered_set<cell_gid_type> visited;
- std::queue<cell_gid_type> scq;
-
- for (auto gid: gids) {
- if (visited.count(gid)) continue;
- visited.insert(gid);
-
- cell_size_type sc_size = 0;
- scq.push(gid);
- while (!scq.empty()) {
- auto g = scq.front();
- scq.pop();
-
- gids_.push_back(g);
- ++sc_size;
-
- for (auto gj: rec.gap_junctions_on(g)) {
- cell_gid_type peer =
- gj.local.gid==g? gj.peer.gid:
- gj.peer.gid==g? gj.local.gid:
- throw bad_cell_description(cell_kind::cable1d_neuron, g);
-
- if (!gid_to_loc.count(peer)) {
- throw gj_unsupported_domain_decomposition(g, peer);
- }
-
- if (!visited.count(peer)) {
- visited.insert(peer);
- scq.push(peer);
- }
- }
- }
-
- deps_.push_back(sc_size>1? sc_size: 0);
- deps_.insert(deps_.end(), sc_size-1, 0);
- }
-
- perm_gids_.reserve(gids_.size());
- for (auto gid: gids_) {
- perm_gids_.push_back(gid_to_loc[gid]);
- }
-}
-
void mc_cell_group::reset() {
spikes_.clear();
@@ -135,21 +85,51 @@ void mc_cell_group::advance(epoch ep, time_type dt, const event_lane_subrange& e
PE(advance_eventsetup);
staged_events_.clear();
+
+ fvm_index_type ev_begin = 0, ev_mid = 0, ev_end = 0;
// skip event binning if empty lanes are passed
if (event_lanes.size()) {
- for (auto lid: util::count_along(gids_)) {
- auto& lane = event_lanes[perm_gids_[lid]];
+
+ std::vector<cell_size_type> idx_sorted_by_intdom(cell_to_intdom_.size());
+ std::iota(idx_sorted_by_intdom.begin(), idx_sorted_by_intdom.end(), 0);
+ util::sort_by(idx_sorted_by_intdom, [&](cell_size_type i) { return cell_to_intdom_[i]; });
+
+ /// Event merging on integration domain could benefit from the use of the logic from `tree_merge_events`
+ fvm_index_type prev_intdom = -1;
+ for (auto i: util::count_along(gids_)) {
+ unsigned count_staged = 0;
+
+ auto lid = idx_sorted_by_intdom[i];
+ auto& lane = event_lanes[lid];
+ auto curr_intdom = cell_to_intdom_[lid];
+
for (auto e: lane) {
if (e.time>=ep.tfinal) break;
e.time = binners_[lid].bin(e.time, tstart);
auto h = target_handles_[target_handle_divisions_[lid]+e.target.index];
auto ev = deliverable_event(e.time, h, e.weight);
staged_events_.push_back(ev);
+ count_staged++;
}
+
+ ev_end += count_staged;
+
+ if (curr_intdom != prev_intdom) {
+ ev_begin = ev_end - count_staged;
+ prev_intdom = curr_intdom;
+ }
+ else {
+ std::inplace_merge(staged_events_.begin() + ev_begin,
+ staged_events_.begin() + ev_mid,
+ staged_events_.begin() + ev_end);
+ }
+
+ ev_mid = ev_end;
}
}
PL();
+
// Create sample events and delivery information.
//
// For each (schedule, sampler, probe set) in the sampler association
@@ -196,7 +176,7 @@ void mc_cell_group::advance(epoch ep, time_type dt, const event_lane_subrange& e
call_info.push_back({sa.sampler, pid, p.tag, n_samples, n_samples+n_times});
for (auto t: sample_times) {
- sample_event ev{t, cell_index, {p.handle, n_samples++}};
+ sample_event ev{t, (cell_gid_type)cell_to_intdom_[cell_index], {p.handle, n_samples++}};
sample_events.push_back(ev);
}
}
diff --git a/arbor/mc_cell_group.hpp b/arbor/mc_cell_group.hpp
index 4efd99d0..f22d3b22 100644
--- a/arbor/mc_cell_group.hpp
+++ b/arbor/mc_cell_group.hpp
@@ -57,26 +57,12 @@ public:
void remove_all_samplers() override;
- void generate_deps_gids(const recipe& rec, std::vector<cell_gid_type>);
-
- std::vector<cell_gid_type> get_gids() {
- return gids_;
- }
-
- std::vector<int> get_dependencies() {
- return deps_;
- }
-
private:
// List of the gids of the cells in the group.
std::vector<cell_gid_type> gids_;
- // Permutation of the gids of the cells in the group.
- // perm_gids_[i] is the index of gids_[i] in the original gids vector passed to the ctor
- std::vector<cell_size_type> perm_gids_;
-
- // List of the dependencies of the cells in the group.
- std::vector<int> deps_;
+ // Map from gid to integration domain id
+ std::vector<fvm_index_type> cell_to_intdom_;
// Hash table for converting gid to local index
std::unordered_map<cell_gid_type, cell_gid_type> gid_index_map_;
diff --git a/example/gap_junctions/gap_junctions.cpp b/example/gap_junctions/gap_junctions.cpp
index 94a697f7..e37bb671 100644
--- a/example/gap_junctions/gap_junctions.cpp
+++ b/example/gap_junctions/gap_junctions.cpp
@@ -21,6 +21,9 @@
#include <arbor/recipe.hpp>
#include <arbor/version.hpp>
+#include <arborenv/concurrency.hpp>
+#include <arborenv/gpu_env.hpp>
+
#include <sup/ioutil.hpp>
#include <sup/json_meter.hpp>
@@ -43,18 +46,18 @@ using arb::cell_probe_address;
void write_trace_json(const std::vector<arb::trace_data<double>>& trace, unsigned rank);
// Generate a cell.
-arb::mc_cell gj_cell(double delay, double duration);
+arb::mc_cell gj_cell(cell_gid_type gid, unsigned ncells, double stim_duration);
class gj_recipe: public arb::recipe {
public:
gj_recipe(gap_params params): params_(params) {}
cell_size_type num_cells() const override {
- return params_.cells_per_ring * params_.num_rings;
+ return params_.n_cells_per_cable * params_.n_cables;
}
arb::util::unique_any get_cell_description(cell_gid_type gid) const override {
- return gj_cell(params_.delay*gid, params_.duration);
+ return gj_cell(gid, params_.n_cells_per_cable, params_.stim_duration);
}
cell_kind get_cell_kind(cell_gid_type gid) const override {
@@ -66,9 +69,16 @@ public:
return 1;
}
- // The cell has one target synapse, which will be connected to cell gid-1.
+ // The cell has one target synapse, which will be connected to a cell in another cable.
cell_size_type num_targets(cell_gid_type gid) const override {
- return 0;
+ return 1;
+ }
+
+ std::vector<arb::cell_connection> connections_on(cell_gid_type gid) const override {
+ if(gid % params_.n_cells_per_cable || (int)gid - 1 < 0) {
+ return{};
+ }
+ return {arb::cell_connection({gid - 1, 0}, {gid, 0}, params_.event_weight, params_.event_min_delay)};
}
// There is one probe (for measuring voltage at the soma) on the cell.
@@ -94,14 +104,22 @@ public:
std::vector<arb::gap_junction_connection> gap_junctions_on(cell_gid_type gid) const override{
std::vector<arb::gap_junction_connection> conns;
- cell_gid_type ring_start_gid = (gid/params_.cells_per_ring) * params_.cells_per_ring;
- cell_gid_type next_cell = (gid + 1) % params_.cells_per_ring + ring_start_gid;
- cell_gid_type prev_cell = (gid - 1 + params_.cells_per_ring) % params_.cells_per_ring + ring_start_gid;
+ int cable_begin = (gid/params_.n_cells_per_cable) * params_.n_cells_per_cable;
+ int cable_end = cable_begin + params_.n_cells_per_cable;
+
+ int next_cell = gid + 1;
+ int prev_cell = gid - 1;
+
+ // Soma is connected to the prev cell's dendrite
+ // Dendrite is connected to the next cell's soma
+ // Gap junction conductance in μS
- // Our soma is connected to the next cell's dendrite
- // Our dendrite is connected to the prev cell's soma
- conns.push_back(arb::gap_junction_connection({next_cell, 1}, {gid, 0}, 0.0037)); // 1 is the id of the dendrite junction
- conns.push_back(arb::gap_junction_connection({prev_cell, 0}, {gid, 1}, 0.0037)); // 0 is the id of the soma junction
+ if (next_cell < cable_end) {
+ conns.push_back(arb::gap_junction_connection({(cell_gid_type)next_cell, 0}, {gid, 1}, 0.015));
+ }
+ if (prev_cell >= cable_begin) {
+ conns.push_back(arb::gap_junction_connection({(cell_gid_type)prev_cell, 1}, {gid, 0}, 0.015));
+ }
return conns;
}
@@ -157,16 +175,26 @@ int main(int argc, char** argv) {
try {
bool root = true;
+ arb::proc_allocation resources;
+ if (auto nt = arbenv::get_env_num_threads()) {
+ resources.num_threads = nt;
+ }
+ else {
+ resources.num_threads = arbenv::thread_concurrency();
+ }
+
#ifdef ARB_MPI_ENABLED
arbenv::with_mpi guard(argc, argv, false);
- auto context = arb::make_context(arb::proc_allocation(), MPI_COMM_WORLD);
+ resources.gpu_id = arbenv::find_private_gpu(MPI_COMM_WORLD);
+ auto context = arb::make_context(resources, MPI_COMM_WORLD);
{
int rank;
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
root = rank==0;
}
#else
- auto context = arb::make_context();
+ resources.gpu_id = arbenv::default_gpu();
+ auto context = arb::make_context(resources);
#endif
#ifdef ARB_PROFILE_ENABLED
@@ -225,7 +253,7 @@ int main(int argc, char** argv) {
std::cout << "running simulation" << std::endl;
// Run the simulation for 100 ms, with time steps of 0.025 ms.
- sim.run(params.duration, 0.025);
+ sim.run(params.sim_duration, 0.025);
meters.checkpoint("model-run", context);
@@ -234,7 +262,7 @@ int main(int argc, char** argv) {
// Write spikes to file
if (root) {
std::cout << "\n" << ns << " spikes generated at rate of "
- << params.duration/ns << " ms between spikes\n";
+ << params.sim_duration/ns << " ms between spikes\n";
std::ofstream fid("spikes.gdf");
if (!fid.good()) {
std::cerr << "Warning: unable to open file spikes.gdf for spike output\n";
@@ -275,6 +303,7 @@ void write_trace_json(const std::vector<arb::trace_data<double>>& trace, unsigne
json["name"] = "gj demo: cell " + std::to_string(i);
json["units"] = "mV";
json["cell"] = std::to_string(i);
+ json["group"] = std::to_string(rank);
json["probe"] = "0";
auto &jt = json["data"]["time"];
@@ -290,7 +319,7 @@ void write_trace_json(const std::vector<arb::trace_data<double>>& trace, unsigne
}
}
-arb::mc_cell gj_cell(double delay, double duration) {
+arb::mc_cell gj_cell(cell_gid_type gid, unsigned ncell, double stim_duration) {
arb::mc_cell cell;
arb::mechanism_desc nax("nax");
@@ -327,11 +356,16 @@ arb::mc_cell gj_cell(double delay, double duration) {
cell.add_detector({0,0}, 10);
- cell.add_gap_junction({0, 1}); //ggap in μS
- cell.add_gap_junction({1, 1}); //ggap in μS
+ cell.add_gap_junction({0, 1});
+ cell.add_gap_junction({1, 1});
+
+ if (!gid) {
+ arb::i_clamp stim(0, stim_duration, 0.4);
+ cell.add_stimulus({0, 0.5}, stim);
+ }
- arb::i_clamp stim(delay, duration, 0.2);
- cell.add_stimulus({1, 0.95}, stim);
+ // Add a synapse to the mid point of the first dendrite.
+ cell.add_synapse({1, 0.5}, "expsyn");
return cell;
}
diff --git a/example/gap_junctions/parameters.hpp b/example/gap_junctions/parameters.hpp
index 2c2e7487..f7dc5f05 100644
--- a/example/gap_junctions/parameters.hpp
+++ b/example/gap_junctions/parameters.hpp
@@ -13,11 +13,13 @@ struct gap_params {
gap_params() = default;
std::string name = "default";
- unsigned cells_per_ring = 100;
- unsigned num_rings = 10;
- double duration = 100;
- double delay = 0.5;
- bool print_all = false;
+ unsigned n_cables = 3;
+ unsigned n_cells_per_cable = 5;
+ double stim_duration = 30;
+ double event_min_delay = 10;
+ double event_weight = 0.05;
+ double sim_duration = 100;
+ bool print_all = true;
};
gap_params read_options(int argc, char** argv) {
@@ -44,10 +46,12 @@ gap_params read_options(int argc, char** argv) {
json << f;
param_from_json(params.name, "name", json);
- param_from_json(params.cells_per_ring, "num-cells", json);
- param_from_json(params.num_rings, "num-rings", json);
- param_from_json(params.duration, "duration", json);
- param_from_json(params.delay, "delay", json);
+ param_from_json(params.n_cables, "n-cables", json);
+ param_from_json(params.n_cells_per_cable, "n-cells-per-cable", json);
+ param_from_json(params.stim_duration, "stim-duration", json);
+ param_from_json(params.event_min_delay, "event-min-delay", json);
+ param_from_json(params.event_weight, "event-weight", json);
+ param_from_json(params.sim_duration, "sim-duration", json);
param_from_json(params.print_all, "print-all", json);
if (!json.empty()) {
diff --git a/example/gap_junctions/readme.md b/example/gap_junctions/readme.md
index 1d32a192..a338020c 100644
--- a/example/gap_junctions/readme.md
+++ b/example/gap_junctions/readme.md
@@ -1,3 +1,45 @@
# Gap Junctions Example
A miniapp that demonstrates how to describe how to build a network with gap junctions.
+
+##Structure:
+Cells are structured into groups that are inter-connected by gap junctions; Groups are
+connected by synapses. The first cell of the first group (top left in diagram) has a
+current stimulus.
+
+
+```
+c --gj-- c --gj-- c --gj-- c --gj-- c
+ |
+ syn
+ |
+c --gj-- c --gj-- c --gj-- c --gj-- c
+|
+syn
+|
+c --gj-- c --gj-- c --gj-- c --gj-- c
+```
+
+
+## Tunable parameters
+* _n_cables_: number of groups of cells connected by gap junctions.
+* _n_cells_per_cable_: number of cells in a group.
+* _stim_duration_: duration that the stimulus on the first cell is on.
+* _event_min_delay_: minimum delay of the network.
+* _event_weight_: weight of an event.
+* _sim_duration_: duration of the simulation.
+* _print_all_: print the voltages of all cells in nerwork.
+
+An example parameter file is:
+```
+{
+ "name": "small test",
+ "n-cables": 3,
+ "n-cells-per-cable": 5,
+ "stim-duration": 30,
+ "event-min-delay": 10,
+ "event-weight": 0.05,
+ "sim-duration": 100,
+ "print-all": false
+}
+```
\ No newline at end of file
diff --git a/test/unit/CMakeLists.txt b/test/unit/CMakeLists.txt
index 1148c002..977372b9 100644
--- a/test/unit/CMakeLists.txt
+++ b/test/unit/CMakeLists.txt
@@ -103,7 +103,6 @@ set(unit_sources
test_spike_store.cpp
test_stats.cpp
test_strprintf.cpp
- test_supercell.cpp
test_swcio.cpp
test_synapses.cpp
test_thread.cpp
diff --git a/test/unit/test_fvm_lowered.cpp b/test/unit/test_fvm_lowered.cpp
index 84e2e226..8e2c2f25 100644
--- a/test/unit/test_fvm_lowered.cpp
+++ b/test/unit/test_fvm_lowered.cpp
@@ -78,6 +78,128 @@ ACCESS_BIND(\
using namespace arb;
+class gap_recipe_0: public recipe {
+public:
+ gap_recipe_0() {}
+
+ cell_size_type num_cells() const override {
+ return size_;
+ }
+
+ arb::util::unique_any get_cell_description(cell_gid_type gid) const override {
+ mc_cell c;
+ c.add_soma(20);
+ c.add_gap_junction({0, 1});
+ return {std::move(c)};
+ }
+
+ cell_kind get_cell_kind(cell_gid_type gid) const override {
+ return cell_kind::cable1d_neuron;
+ }
+ std::vector<gap_junction_connection> gap_junctions_on(cell_gid_type gid) const override {
+ switch (gid) {
+ case 0 :
+ return {gap_junction_connection({5, 0}, {0, 0}, 0.1)};
+ case 2 :
+ return {
+ gap_junction_connection({3, 0}, {2, 0}, 0.1),
+ };
+ case 3 :
+ return {
+ gap_junction_connection({7, 0}, {3, 0}, 0.1),
+ gap_junction_connection({3, 0}, {2, 0}, 0.1)
+ };
+ case 5 :
+ return {gap_junction_connection({5, 0}, {0, 0}, 0.1)};
+ case 7 :
+ return {
+ gap_junction_connection({3, 0}, {7, 0}, 0.1),
+ };
+ default :
+ return {};
+ }
+ }
+
+private:
+ cell_size_type size_ = 12;
+};
+
+class gap_recipe_1: public recipe {
+public:
+ gap_recipe_1() {}
+
+ cell_size_type num_cells() const override {
+ return size_;
+ }
+
+ arb::util::unique_any get_cell_description(cell_gid_type) const override {
+ mc_cell c;
+ c.add_soma(20);
+ return {std::move(c)};
+ }
+
+ cell_kind get_cell_kind(cell_gid_type gid) const override {
+ return cell_kind::cable1d_neuron;
+ }
+
+private:
+ cell_size_type size_ = 12;
+};
+
+class gap_recipe_2: public recipe {
+public:
+ gap_recipe_2() {}
+
+ cell_size_type num_cells() const override {
+ return size_;
+ }
+
+ arb::util::unique_any get_cell_description(cell_gid_type) const override {
+ mc_cell c;
+ c.add_soma(20);
+ c.add_gap_junction({0,1});
+ return {std::move(c)};
+ }
+
+ cell_kind get_cell_kind(cell_gid_type gid) const override {
+ return cell_kind::cable1d_neuron;
+ }
+ std::vector<gap_junction_connection> gap_junctions_on(cell_gid_type gid) const override {
+ switch (gid) {
+ case 0 :
+ return {
+ gap_junction_connection({2, 0}, {0, 0}, 0.1),
+ gap_junction_connection({3, 0}, {0, 0}, 0.1),
+ gap_junction_connection({5, 0}, {0, 0}, 0.1)
+ };
+ case 2 :
+ return {
+ gap_junction_connection({0, 0}, {2, 0}, 0.1),
+ gap_junction_connection({3, 0}, {2, 0}, 0.1),
+ gap_junction_connection({5, 0}, {2, 0}, 0.1)
+ };
+ case 3 :
+ return {
+ gap_junction_connection({0, 0}, {3, 0}, 0.1),
+ gap_junction_connection({2, 0}, {3, 0}, 0.1),
+ gap_junction_connection({5, 0}, {3, 0}, 0.1)
+ };
+ case 5 :
+ return {
+ gap_junction_connection({2, 0}, {5, 0}, 0.1),
+ gap_junction_connection({3, 0}, {5, 0}, 0.1),
+ gap_junction_connection({0, 0}, {5, 0}, 0.1)
+ };
+ default :
+ return {};
+ }
+ }
+
+private:
+ cell_size_type size_ = 12;
+};
+
+
TEST(fvm_lowered, matrix_init)
{
execution_context context;
@@ -92,10 +214,11 @@ TEST(fvm_lowered, matrix_init)
cell.segment(1)->set_compartments(10);
std::vector<target_handle> targets;
+ std::vector<fvm_index_type> cell_to_intdom;
probe_association_map<probe_handle> probe_map;
fvm_cell fvcell(context);
- fvcell.initialize({0}, {0}, cable1d_recipe(cell), targets, probe_map);
+ fvcell.initialize({0}, cable1d_recipe(cell), cell_to_intdom, targets, probe_map);
auto& J = fvcell.*private_matrix_ptr;
EXPECT_EQ(J.size(), 11u);
@@ -137,10 +260,11 @@ TEST(fvm_lowered, target_handles) {
cells[1].add_detector({0, 0}, 3.3);
std::vector<target_handle> targets;
+ std::vector<fvm_index_type> cell_to_intdom;
probe_association_map<probe_handle> probe_map;
fvm_cell fvcell(context);
- fvcell.initialize({0, 1}, {0, 0}, cable1d_recipe(cells), targets, probe_map);
+ fvcell.initialize({0, 1}, cable1d_recipe(cells), cell_to_intdom, targets, probe_map);
mechanism* expsyn = find_mechanism(fvcell, "expsyn");
ASSERT_TRUE(expsyn);
@@ -154,19 +278,19 @@ TEST(fvm_lowered, target_handles) {
EXPECT_EQ(expsyn_id, targets[0].mech_id);
EXPECT_EQ(1u, targets[0].mech_index);
- EXPECT_EQ(0u, targets[0].cell_index);
+ EXPECT_EQ(0u, targets[0].intdom_index);
EXPECT_EQ(expsyn_id, targets[1].mech_id);
EXPECT_EQ(0u, targets[1].mech_index);
- EXPECT_EQ(0u, targets[1].cell_index);
+ EXPECT_EQ(0u, targets[1].intdom_index);
EXPECT_EQ(exp2syn_id, targets[2].mech_id);
EXPECT_EQ(0u, targets[2].mech_index);
- EXPECT_EQ(1u, targets[2].cell_index);
+ EXPECT_EQ(1u, targets[2].intdom_index);
EXPECT_EQ(expsyn_id, targets[3].mech_id);
EXPECT_EQ(2u, targets[3].mech_index);
- EXPECT_EQ(1u, targets[3].cell_index);
+ EXPECT_EQ(1u, targets[3].intdom_index);
}
TEST(fvm_lowered, stimulus) {
@@ -196,6 +320,7 @@ TEST(fvm_lowered, stimulus) {
// The implementation of the stimulus is tested by creating a lowered cell, then
// testing that the correct currents are injected at the correct control volumes
// as during the stimulus windows.
+ std::vector<fvm_index_type> cell_to_intdom(cells.size(), 0);
fvm_discretization D = fvm_discretize(cells);
const auto& A = D.cv_area;
@@ -204,7 +329,7 @@ TEST(fvm_lowered, stimulus) {
probe_association_map<probe_handle> probe_map;
fvm_cell fvcell(context);
- fvcell.initialize({0}, {0}, cable1d_recipe(cells), targets, probe_map);
+ fvcell.initialize({0}, cable1d_recipe(cells), cell_to_intdom, targets, probe_map);
mechanism* stim = find_mechanism(fvcell, "_builtin_stimulus");
ASSERT_TRUE(stim);
@@ -294,11 +419,12 @@ TEST(fvm_lowered, derived_mechs) {
// Test initialization and global parameter values.
std::vector<target_handle> targets;
+ std::vector<fvm_index_type> cell_to_intdom;
probe_association_map<probe_handle> probe_map;
execution_context context;
fvm_cell fvcell(context);
- fvcell.initialize({0, 1, 2}, {0, 0, 0}, rec, targets, probe_map);
+ fvcell.initialize({0, 1, 2}, rec, cell_to_intdom, targets, probe_map);
// Both mechanisms will have the same internal name, "test_kin1".
@@ -401,10 +527,11 @@ TEST(fvm_lowered, weighted_write_ion) {
c.segments()[3] ->add_mechanism("test_ca");
std::vector<target_handle> targets;
+ std::vector<fvm_index_type> cell_to_intdom;
probe_association_map<probe_handle> probe_map;
fvm_cell fvcell(context);
- fvcell.initialize({0}, {0}, cable1d_recipe(c), targets, probe_map);
+ fvcell.initialize({0}, cable1d_recipe(c), cell_to_intdom, targets, probe_map);
auto& state = *(fvcell.*private_state_ptr).get();
auto& ion = state.ion_data.at(ionKind::ca);
@@ -591,9 +718,13 @@ TEST(fvm_lowered, gj_coords_complex) {
cells.push_back(std::move(c1));
cells.push_back(std::move(c2));
- fvm_discretization D = fvm_discretize(cells);
+ std::vector<fvm_index_type> cell_to_intdom;
+
std::vector<cell_gid_type> gids = {0, 1, 2};
+ fvcell.fvm_intdom(rec, gids, cell_to_intdom);
+ fvm_discretization D = fvm_discretize(cells);
+
auto GJ = fvcell.fvm_gap_junctions(cells, gids, rec, D);
EXPECT_EQ(10u, GJ.size());
@@ -661,18 +792,22 @@ TEST(fvm_lowered, cell_group_gj) {
c.add_soma(2.1);
if (i % 2 == 0) {
c.add_gap_junction({0, 1});
- if (i < 10) {
- cell_group0.push_back(std::move(c));
- }
- else {
- cell_group1.push_back(std::move(c));
- }
}
-
+ if (i < 10) {
+ cell_group0.push_back(std::move(c));
+ }
+ else {
+ cell_group1.push_back(std::move(c));
+ }
}
- std::vector<cell_gid_type> gids_cg0 = { 0, 2, 4, 6, 8};
- std::vector<cell_gid_type> gids_cg1 = {10,12,14,16,18};
+ std::vector<cell_gid_type> gids_cg0 = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9};
+ std::vector<cell_gid_type> gids_cg1 = {10, 11, 12, 13, 14, 15, 16, 17, 18, 19};
+
+ std::vector<fvm_index_type> cell_to_intdom0, cell_to_intdom1;
+
+ auto num_dom0 = fvcell.fvm_intdom(rec, gids_cg0, cell_to_intdom0);
+ auto num_dom1 = fvcell.fvm_intdom(rec, gids_cg1, cell_to_intdom1);
fvm_discretization D0 = fvm_discretize(cell_group0);
fvm_discretization D1 = fvm_discretize(cell_group1);
@@ -683,10 +818,61 @@ TEST(fvm_lowered, cell_group_gj) {
EXPECT_EQ(10u, GJ0.size());
EXPECT_EQ(10u, GJ1.size());
- std::vector<pair> expected_loc = {{0, 1}, {0, 4}, {1, 2}, {1, 0}, {2, 3} ,{2, 1}, {3, 4}, {3, 2}, {4, 0}, {4, 3}};
+ std::vector<pair> expected_loc = {{0, 2}, {0, 8}, {2, 4}, {2, 0}, {4, 6} ,{4, 2}, {6, 8}, {6, 4}, {8, 0}, {8, 6}};
for (unsigned i = 0; i < GJ0.size(); i++) {
EXPECT_EQ(expected_loc[i], GJ0[i].loc);
EXPECT_EQ(expected_loc[i], GJ1[i].loc);
}
+
+ std::vector<fvm_index_type> expected_doms= {0u, 1u, 0u, 2u, 0u, 3u, 0u, 4u, 0u, 5u};
+ EXPECT_EQ(6u, num_dom0);
+ EXPECT_EQ(6u, num_dom1);
+
+ EXPECT_EQ(expected_doms, cell_to_intdom0);
+ EXPECT_EQ(expected_doms, cell_to_intdom1);
+
}
+
+TEST(fvm_lowered, integration_domains) {
+ {
+ execution_context context;
+ fvm_cell fvcell(context);
+
+ std::vector<cell_gid_type> gids = {11u, 5u, 2u, 3u, 0u, 8u, 7u};
+ std::vector<fvm_index_type> cell_to_intdom;
+
+ auto num_dom = fvcell.fvm_intdom(gap_recipe_0(), gids, cell_to_intdom);
+ std::vector<fvm_index_type> expected_doms= {0u, 1u, 2u, 2u, 1u, 3u, 2u};
+
+ EXPECT_EQ(4u, num_dom);
+ EXPECT_EQ(expected_doms, cell_to_intdom);
+ }
+ {
+ execution_context context;
+ fvm_cell fvcell(context);
+
+ std::vector<cell_gid_type> gids = {11u, 5u, 2u, 3u, 0u, 8u, 7u};
+ std::vector<fvm_index_type> cell_to_intdom;
+
+ auto num_dom = fvcell.fvm_intdom(gap_recipe_1(), gids, cell_to_intdom);
+ std::vector<fvm_index_type> expected_doms= {0u, 1u, 2u, 3u, 4u, 5u, 6u};
+
+ EXPECT_EQ(7u, num_dom);
+ EXPECT_EQ(expected_doms, cell_to_intdom);
+ }
+ {
+ execution_context context;
+ fvm_cell fvcell(context);
+
+ std::vector<cell_gid_type> gids = {5u, 2u, 3u, 0u};
+ std::vector<fvm_index_type> cell_to_intdom;
+
+ auto num_dom = fvcell.fvm_intdom(gap_recipe_2(), gids, cell_to_intdom);
+ std::vector<fvm_index_type> expected_doms= {0u, 0u, 0u, 0u};
+
+ EXPECT_EQ(1u, num_dom);
+ EXPECT_EQ(expected_doms, cell_to_intdom);
+ }
+}
+
diff --git a/test/unit/test_matrix.cpp b/test/unit/test_matrix.cpp
index c977f88b..eb72dacf 100644
--- a/test/unit/test_matrix.cpp
+++ b/test/unit/test_matrix.cpp
@@ -23,7 +23,7 @@ using vvec = std::vector<value_type>;
TEST(matrix, construct_from_parent_only)
{
std::vector<index_type> p = {0,0,1};
- matrix_type m(p, {0, 3}, vvec(3), vvec(3), vvec(3));
+ matrix_type m(p, {0, 3}, vvec(3), vvec(3), vvec(3), {0});
EXPECT_EQ(m.num_cells(), 1u);
EXPECT_EQ(m.size(), 3u);
EXPECT_EQ(p.size(), 3u);
@@ -41,7 +41,7 @@ TEST(matrix, solve_host)
// trivial case : 1x1 matrix
{
- matrix_type m({0}, {0,1}, vvec(1), vvec(1), vvec(1));
+ matrix_type m({0}, {0,1}, vvec(1), vvec(1), vvec(1), {0});
auto& state = m.state_;
fill(state.d, 2);
fill(state.u, -1);
@@ -57,7 +57,7 @@ TEST(matrix, solve_host)
for(auto n : make_span(2, 1001)) {
auto p = std::vector<index_type>(n);
std::iota(p.begin()+1, p.end(), 0);
- matrix_type m(p, {0, n}, vvec(n), vvec(n), vvec(n));
+ matrix_type m(p, {0, n}, vvec(n), vvec(n), vvec(n), {0});
EXPECT_EQ(m.size(), (unsigned)n);
EXPECT_EQ(m.num_cells(), 1u);
@@ -94,7 +94,8 @@ TEST(matrix, zero_diagonal)
// Three matrices, sizes 3, 3 and 2, with no branching.
std::vector<index_type> p = {0, 0, 1, 3, 3, 5, 5};
std::vector<index_type> c = {0, 3, 5, 7};
- matrix_type m(p, c, vvec(7), vvec(7), vvec(7));
+ std::vector<index_type> i = {0, 1, 2};
+ matrix_type m(p, c, vvec(7), vvec(7), vvec(7), i);
EXPECT_EQ(7u, m.size());
EXPECT_EQ(3u, m.num_cells());
@@ -129,6 +130,7 @@ TEST(matrix, zero_diagonal_assembled)
// Three matrices, sizes 3, 3 and 2, with no branching.
std::vector<index_type> p = {0, 0, 1, 3, 3, 5, 5};
std::vector<index_type> c = {0, 3, 5, 7};
+ std::vector<index_type> s = {0, 1, 2};
// Face conductances.
vvec g = {0, 1, 1, 0, 1, 0, 2};
@@ -152,7 +154,7 @@ TEST(matrix, zero_diagonal_assembled)
// Expected solution:
// x = [ 4 5 6 7 8 9 10]
- matrix_type m(p, c, Cm, g, area);
+ matrix_type m(p, c, Cm, g, area, s);
m.assemble(dt, v, i);
m.solve();
diff --git a/test/unit/test_matrix.cu b/test/unit/test_matrix.cu
index b3424ad2..8d01632d 100644
--- a/test/unit/test_matrix.cu
+++ b/test/unit/test_matrix.cu
@@ -256,6 +256,7 @@ TEST(matrix, backends)
std::vector<I> p;
std::vector<I> cell_cv_divs;
+ std::vector<I> cell_to_intdom;
for (auto m=0; m<num_mtx; ++m) {
auto &p_ref = p_base[m%2];
auto first = p.size();
@@ -263,6 +264,7 @@ TEST(matrix, backends)
p.push_back(i + first);
}
cell_cv_divs.push_back(first);
+ cell_to_intdom.push_back(m);
}
cell_cv_divs.push_back(p.size());
@@ -286,9 +288,9 @@ TEST(matrix, backends)
std::generate(i.begin(), i.end(), [&](){return dist(gen);});
// Make the reference matrix and the gpu matrix
- auto flat = state_flat(p, cell_cv_divs, Cm, g, area); // flat
- auto intl = state_intl(p, cell_cv_divs, Cm, g, area); // interleaved
- auto fine = state_fine(p, cell_cv_divs, Cm, g, area); // interleaved
+ auto flat = state_flat(p, cell_cv_divs, Cm, g, area, cell_to_intdom); // flat
+ auto intl = state_intl(p, cell_cv_divs, Cm, g, area, cell_to_intdom); // interleaved
+ auto fine = state_fine(p, cell_cv_divs, Cm, g, area, cell_to_intdom); // interleaved
// Set the integration times for the cells to be between 0.01 and 0.02 ms.
std::vector<T> dt(num_mtx, 0);
diff --git a/test/unit/test_matrix_cpuvsgpu.cpp b/test/unit/test_matrix_cpuvsgpu.cpp
index 861c83dd..9d90e062 100644
--- a/test/unit/test_matrix_cpuvsgpu.cpp
+++ b/test/unit/test_matrix_cpuvsgpu.cpp
@@ -86,6 +86,7 @@ TEST(matrix, assemble)
std::vector<I> p;
std::vector<I> cell_index;
+ std::vector<I> intdom_index;
for (auto m=0; m<num_mtx; ++m) {
auto &p_ref = p_base[m%p_base.size()];
auto first = p.size();
@@ -93,6 +94,7 @@ TEST(matrix, assemble)
p.push_back(i + first);
}
cell_index.push_back(first);
+ intdom_index.push_back(m);
}
cell_index.push_back(p.size());
@@ -112,8 +114,8 @@ TEST(matrix, assemble)
std::vector<T> area(group_size, 1e3);
// Make the reference matrix and the gpu matrix
- auto m_mc = mc_state( p, cell_index, Cm, g, area); // on host
- auto m_gpu = gpu_state(p, cell_index, Cm, g, area); // on gpu
+ auto m_mc = mc_state( p, cell_index, Cm, g, area, intdom_index); // on host
+ auto m_gpu = gpu_state(p, cell_index, Cm, g, area, intdom_index); // on gpu
// Set the integration times for the cells to be between 0.1 and 0.2 ms.
std::vector<T> dt(num_mtx);
diff --git a/test/unit/test_mc_cell_group.cpp b/test/unit/test_mc_cell_group.cpp
index 70077f32..11e33300 100644
--- a/test/unit/test_mc_cell_group.cpp
+++ b/test/unit/test_mc_cell_group.cpp
@@ -28,128 +28,6 @@ namespace {
return c;
}
-
- class gap_recipe_0: public recipe {
- public:
- gap_recipe_0() {}
-
- cell_size_type num_cells() const override {
- return size_;
- }
-
- arb::util::unique_any get_cell_description(cell_gid_type gid) const override {
- mc_cell c;
- c.add_soma(20);
- c.add_gap_junction({0, 1});
- return {std::move(c)};
- }
-
- cell_kind get_cell_kind(cell_gid_type gid) const override {
- return cell_kind::cable1d_neuron;
- }
- std::vector<gap_junction_connection> gap_junctions_on(cell_gid_type gid) const override {
- switch (gid) {
- case 0 :
- return {gap_junction_connection({5, 0}, {0, 0}, 0.1)};
- case 2 :
- return {
- gap_junction_connection({3, 0}, {2, 0}, 0.1),
- };
- case 3 :
- return {
- gap_junction_connection({7, 0}, {3, 0}, 0.1),
- gap_junction_connection({3, 0}, {2, 0}, 0.1)
- };
- case 5 :
- return {gap_junction_connection({5, 0}, {0, 0}, 0.1)};
- case 7 :
- return {
- gap_junction_connection({3, 0}, {7, 0}, 0.1),
- };
- default :
- return {};
- }
- }
-
- private:
- cell_size_type size_ = 12;
- };
-
- class gap_recipe_1: public recipe {
- public:
- gap_recipe_1() {}
-
- cell_size_type num_cells() const override {
- return size_;
- }
-
- arb::util::unique_any get_cell_description(cell_gid_type) const override {
- mc_cell c;
- c.add_soma(20);
- return {std::move(c)};
- }
-
- cell_kind get_cell_kind(cell_gid_type gid) const override {
- return cell_kind::cable1d_neuron;
- }
-
- private:
- cell_size_type size_ = 12;
- };
-
- class gap_recipe_2: public recipe {
- public:
- gap_recipe_2() {}
-
- cell_size_type num_cells() const override {
- return size_;
- }
-
- arb::util::unique_any get_cell_description(cell_gid_type) const override {
- mc_cell c;
- c.add_soma(20);
- c.add_gap_junction({0,1});
- return {std::move(c)};
- }
-
- cell_kind get_cell_kind(cell_gid_type gid) const override {
- return cell_kind::cable1d_neuron;
- }
- std::vector<gap_junction_connection> gap_junctions_on(cell_gid_type gid) const override {
- switch (gid) {
- case 0 :
- return {
- gap_junction_connection({2, 0}, {0, 0}, 0.1),
- gap_junction_connection({3, 0}, {0, 0}, 0.1),
- gap_junction_connection({5, 0}, {0, 0}, 0.1)
- };
- case 2 :
- return {
- gap_junction_connection({0, 0}, {2, 0}, 0.1),
- gap_junction_connection({3, 0}, {2, 0}, 0.1),
- gap_junction_connection({5, 0}, {2, 0}, 0.1)
- };
- case 3 :
- return {
- gap_junction_connection({0, 0}, {3, 0}, 0.1),
- gap_junction_connection({2, 0}, {3, 0}, 0.1),
- gap_junction_connection({5, 0}, {3, 0}, 0.1)
- };
- case 5 :
- return {
- gap_junction_connection({2, 0}, {5, 0}, 0.1),
- gap_junction_connection({3, 0}, {5, 0}, 0.1),
- gap_junction_connection({0, 0}, {5, 0}, 0.1)
- };
- default :
- return {};
- }
- }
-
- private:
- cell_size_type size_ = 12;
- };
-
}
ACCESS_BIND(
@@ -207,30 +85,3 @@ TEST(mc_cell_group, sources) {
}
}
-TEST(mc_cell_group, generated_gids_deps_) {
- {
- std::vector<cell_gid_type> gids = {11u, 5u, 2u, 3u, 0u, 8u, 7u};
- mc_cell_group group{gids, gap_recipe_0(), lowered_cell()};
-
- std::vector<cell_gid_type> expected_gids = {11u, 5u, 0u, 2u, 3u, 7u, 8u};
- std::vector<int> expected_deps = {0, 2, 0, 3, 0, 0, 0};
- EXPECT_EQ(expected_gids, group.get_gids());
- EXPECT_EQ(expected_deps, group.get_dependencies());
- }
- {
- std::vector<cell_gid_type> gids = {11u, 5u, 2u, 3u, 0u, 8u, 7u};
- mc_cell_group group{gids, gap_recipe_1(), lowered_cell()};
-
- std::vector<int> expected_deps = {0, 0, 0, 0, 0, 0, 0};
- EXPECT_EQ(gids, group.get_gids());
- EXPECT_EQ(expected_deps, group.get_dependencies());
- }
- {
- std::vector<cell_gid_type> gids = {5u, 2u, 3u, 0u};
- mc_cell_group group{gids, gap_recipe_2(), lowered_cell()};
-
- std::vector<int> expected_deps = {4, 0, 0, 0};
- EXPECT_EQ(gids, group.get_gids());
- EXPECT_EQ(expected_deps, group.get_dependencies());
- }
-}
diff --git a/test/unit/test_mech_temperature.cpp b/test/unit/test_mech_temperature.cpp
index 81a11394..e4e10226 100644
--- a/test/unit/test_mech_temperature.cpp
+++ b/test/unit/test_mech_temperature.cpp
@@ -22,13 +22,12 @@ void run_celsius_test() {
fvm_size_type ncell = 1;
fvm_size_type ncv = 3;
- std::vector<fvm_index_type> cv_to_cell(ncv, 0);
+ std::vector<fvm_index_type> cv_to_intdom(ncv, 0);
std::vector<fvm_gap_junction> gj = {};
- std::vector<int> deps = {0};
auto celsius_test = cat.instance<backend>("celsius_test");
auto shared_state = std::make_unique<typename backend::shared_state>(
- ncell, cv_to_cell, deps, gj, celsius_test->data_alignment());
+ ncell, cv_to_intdom, gj, celsius_test->data_alignment());
mechanism::layout layout;
layout.weight.assign(ncv, 1.);
diff --git a/test/unit/test_probe.cpp b/test/unit/test_probe.cpp
index 2e52784a..6f449cb1 100644
--- a/test/unit/test_probe.cpp
+++ b/test/unit/test_probe.cpp
@@ -37,10 +37,11 @@ TEST(probe, fvm_lowered_cell) {
rec.add_probe(0, 30, cell_probe_address{loc2, cell_probe_address::membrane_current});
std::vector<target_handle> targets;
+ std::vector<fvm_index_type> cell_to_intdom;
probe_association_map<probe_handle> probe_map;
fvm_cell lcell(context);
- lcell.initialize({0}, {0}, rec, targets, probe_map);
+ lcell.initialize({0}, rec, cell_to_intdom, targets, probe_map);
EXPECT_EQ(3u, rec.num_probes(0));
EXPECT_EQ(3u, probe_map.size());
diff --git a/test/unit/test_supercell.cpp b/test/unit/test_supercell.cpp
deleted file mode 100644
index f4bebb2c..00000000
--- a/test/unit/test_supercell.cpp
+++ /dev/null
@@ -1,71 +0,0 @@
-#include "../gtest.h"
-
-#include <cmath>
-#include <tuple>
-#include <vector>
-
-#include <arbor/constants.hpp>
-#include <arbor/mechcat.hpp>
-#include <arbor/util/optional.hpp>
-#include <arbor/mc_cell.hpp>
-
-#include "backends/multicore/fvm.hpp"
-#include "backends/multicore/mechanism.hpp"
-#include "util/maputil.hpp"
-#include "util/range.hpp"
-
-#include "common.hpp"
-#include "mech_private_field_access.hpp"
-
-using namespace arb;
-
-using backend = ::arb::multicore::backend;
-using shared_state = backend::shared_state;
-using value_type = backend::value_type;
-using size_type = backend::size_type;
-
-// Access to more mechanism protected data:
-
-ACCESS_BIND(const value_type* multicore::mechanism::*, vec_v_ptr, &multicore::mechanism::vec_v_)
-ACCESS_BIND(value_type* multicore::mechanism::*, vec_i_ptr, &multicore::mechanism::vec_i_)
-
-TEST(supercell, sync_time_to) {
- using value_type = multicore::backend::value_type;
- using index_type = multicore::backend::index_type;
-
- int num_cell = 10;
-
- std::vector<fvm_gap_junction> gj = {};
- std::vector<index_type> deps = {4, 0, 0, 0, 3, 0, 0, 2, 0, 0};
-
- shared_state state(num_cell, std::vector<index_type>(num_cell, 0), deps, gj, 1u);
-
- state.time_to = {0.3, 0.1, 0.2, 0.4, 0.5, 0.6, 0.1, 0.1, 0.6, 0.9};
- state.sync_time_to();
-
- std::vector<value_type> expected = {0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.9};
-
- for (unsigned i = 0; i < state.time_to.size(); i++) {
- EXPECT_EQ(expected[i], state.time_to[i]);
- }
-
- state.time_dep = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
- state.time_to = {0.3, 0.1, 0.2, 0.4, 0.5, 0.6, 0.1, 0.1, 0.6, 0.9};
- expected = {0.3, 0.1, 0.2, 0.4, 0.5, 0.6, 0.1, 0.1, 0.6, 0.9};
- state.sync_time_to();
-
- for (unsigned i = 0; i < state.time_to.size(); i++) {
- EXPECT_EQ(expected[i], state.time_to[i]);
- }
-
- state.time_dep = {10, 0, 0, 0, 0, 0, 0, 0, 0, 0};
- state.time_to = {0.3, 0.1, 0.2, 0.4, 0.5, 0.6, 0.1, 0.1, 0.6, 0.9};
- expected = {0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1};
- state.sync_time_to();
-
- for (unsigned i = 0; i < state.time_to.size(); i++) {
- EXPECT_EQ(expected[i], state.time_to[i]);
- }
-
-}
-
diff --git a/test/unit/test_synapses.cpp b/test/unit/test_synapses.cpp
index 5210a7d4..9219edbd 100644
--- a/test/unit/test_synapses.cpp
+++ b/test/unit/test_synapses.cpp
@@ -77,7 +77,7 @@ TEST(synapses, syn_basic_state) {
int num_syn = 4;
int num_comp = 4;
- int num_cell = 1;
+ int num_intdom = 1;
auto expsyn = unique_cast<multicore::mechanism>(global_default_catalogue().instance<backend>("expsyn"));
ASSERT_TRUE(expsyn);
@@ -87,7 +87,7 @@ TEST(synapses, syn_basic_state) {
std::vector<fvm_gap_junction> gj = {};
auto align = std::max(expsyn->data_alignment(), exp2syn->data_alignment());
- shared_state state(num_cell, std::vector<index_type>(num_comp, 0), std::vector<index_type>(num_cell, 0), gj, align);
+ shared_state state(num_intdom, std::vector<index_type>(num_comp, 0), gj, align);
state.reset(-65., constant::hh_squid_temp);
fill(state.current_density, 1.0);
--
GitLab