diff --git a/arbor/communication/communicator.cpp b/arbor/communication/communicator.cpp
index 9fa4a0d76424d9389ffaec3bb40b1db9d02df86a..6f816bf73854d876bfc5b5b7dec7fbad7bd8110e 100644
--- a/arbor/communication/communicator.cpp
+++ b/arbor/communication/communicator.cpp
@@ -26,15 +26,23 @@ communicator::communicator(const recipe& rec,
const domain_decomposition& dom_dec,
const label_resolution_map& source_resolution_map,
const label_resolution_map& target_resolution_map,
- execution_context& ctx)
-{
- distributed_ = ctx.distributed;
- thread_pool_ = ctx.thread_pool;
+ execution_context& ctx): num_total_cells_{rec.num_cells()},
+ num_local_cells_{dom_dec.num_local_cells()},
+ num_local_groups_{dom_dec.num_groups()},
+ num_domains_{(cell_size_type) ctx.distributed->size()},
+ distributed_{ctx.distributed},
+ thread_pool_{ctx.thread_pool} {
+ update_connections(rec, dom_dec, source_resolution_map, target_resolution_map);
+}
- num_domains_ = distributed_->size();
- num_local_groups_ = dom_dec.num_groups();
- num_local_cells_ = dom_dec.num_local_cells();
- auto num_total_cells = rec.num_cells();
+void communicator::update_connections(const connectivity& rec,
+ const domain_decomposition& dom_dec,
+ const label_resolution_map& source_resolution_map,
+ const label_resolution_map& target_resolution_map) {
+ // Forget all lingering information
+ connections_.clear();
+ connection_part_.clear();
+ index_divisions_.clear();
// For caching information about each cell
struct gid_info {
@@ -81,10 +89,11 @@ communicator::communicator(const recipe& rec,
for (const auto& cell: gid_infos) {
for (auto c: cell.conns) {
- if (c.source.gid >= num_total_cells) {
- throw arb::bad_connection_source_gid(cell.gid, c.source.gid, num_total_cells);
+ auto sgid = c.source.gid;
+ if (sgid >= num_total_cells_) {
+ throw arb::bad_connection_source_gid(cell.gid, sgid, num_total_cells_);
}
- const auto src = dom_dec.gid_domain(c.source.gid);
+ const auto src = dom_dec.gid_domain(sgid);
src_domains.push_back(src);
src_counts[src]++;
}
@@ -95,17 +104,18 @@ communicator::communicator(const recipe& rec,
// the connections as partitioned by the domain of their source gid.
connections_.resize(n_cons);
util::make_partition(connection_part_, src_counts);
- auto offsets = connection_part_;
- std::size_t pos = 0;
+ auto offsets = connection_part_; // Copy, as we use this as the list of current target indices to write into
+ auto src_domain = src_domains.begin();
auto target_resolver = resolver(&target_resolution_map);
for (const auto& cell: gid_infos) {
+ auto index = cell.index_on_domain;
auto source_resolver = resolver(&source_resolution_map);
for (const auto& c: cell.conns) {
- const auto i = offsets[src_domains[pos]]++;
auto src_lid = source_resolver.resolve(c.source);
auto tgt_lid = target_resolver.resolve({cell.gid, c.dest});
- connections_[i] = {{c.source.gid, src_lid}, tgt_lid, c.weight, c.delay, cell.index_on_domain};
- ++pos;
+ auto offset = offsets[*src_domain]++;
+ ++src_domain;
+ connections_[offset] = {{c.source.gid, src_lid}, tgt_lid, c.weight, c.delay, index};
}
}
@@ -119,9 +129,7 @@ communicator::communicator(const recipe& rec,
// This is num_domains_ independent sorts, so it can be parallelized trivially.
const auto& cp = connection_part_;
threading::parallel_for::apply(0, num_domains_, thread_pool_.get(),
- [&](cell_size_type i) {
- util::sort(util::subrange_view(connections_, cp[i], cp[i+1]));
- });
+ [&](cell_size_type i) { util::sort(util::subrange_view(connections_, cp[i], cp[i+1])); });
}
std::pair<cell_size_type, cell_size_type> communicator::group_queue_range(cell_size_type i) {
@@ -130,11 +138,9 @@ std::pair<cell_size_type, cell_size_type> communicator::group_queue_range(cell_s
}
time_type communicator::min_delay() {
- auto local_min = std::numeric_limits<time_type>::max();
- for (auto& con : connections_) {
- local_min = std::min(local_min, con.delay());
- }
-
+ auto local_min = std::accumulate(connections_.begin(), connections_.end(),
+ std::numeric_limits<time_type>::max(),
+ [](auto&& acc, auto&& el) { return std::min(acc, time_type(el.delay)); });
return distributed_->min(local_min);
}
@@ -153,29 +159,21 @@ gathered_vector<spike> communicator::exchange(std::vector<spike> local_spikes) {
return global_spikes;
}
-void communicator::make_event_queues(
- const gathered_vector<spike>& global_spikes,
- std::vector<pse_vector>& queues)
-{
+void communicator::make_event_queues(const gathered_vector<spike>& global_spikes,
+ std::vector<pse_vector>& queues) {
arb_assert(queues.size()==num_local_cells_);
-
- using util::subrange_view;
- using util::make_span;
- using util::make_range;
-
+ // Predicate for partitioning
+ struct spike_pred {
+ bool operator()(const spike& spk, const cell_member_type& src) { return spk.source < src; }
+ bool operator()(const cell_member_type& src, const spike& spk) { return src < spk.source; }
+ };
const auto& sp = global_spikes.partition();
const auto& cp = connection_part_;
- for (auto dom: make_span(num_domains_)) {
- auto cons = subrange_view(connections_, cp[dom], cp[dom+1]);
- auto spks = subrange_view(global_spikes.values(), sp[dom], sp[dom+1]);
-
- struct spike_pred {
- bool operator()(const spike& spk, const cell_member_type& src)
- {return spk.source<src;}
- bool operator()(const cell_member_type& src, const spike& spk)
- {return src<spk.source;}
- };
-
+ for (auto dom: util::make_span(num_domains_)) {
+ auto cons = util::subrange_view(connections_, cp[dom], cp[dom+1]);
+ auto spks = util::subrange_view(global_spikes.values(), sp[dom], sp[dom+1]);
+ auto sp = spks.begin(), se = spks.end();
+ auto cn = cons.begin(), ce = cons.end();
// We have a choice of whether to walk spikes or connections:
// i.e., we can iterate over the spikes, and for each spike search
// the for connections that have the same source; or alternatively
@@ -186,27 +184,19 @@ void communicator::make_event_queues(
// complexity of order max(S log(C), C log(S)), where S is the
// number of spikes, and C is the number of connections.
if (cons.size()<spks.size()) {
- auto sp = spks.begin();
- auto cn = cons.begin();
- while (cn!=cons.end() && sp!=spks.end()) {
- auto sources = std::equal_range(sp, spks.end(), cn->source(), spike_pred());
- for (auto s: make_range(sources)) {
- queues[cn->index_on_domain()].push_back(cn->make_event(s));
- }
-
+ while (cn!=ce && sp!=se) {
+ auto& queue = queues[cn->index_on_domain];
+ auto src = cn->source;
+ auto sources = std::equal_range(sp, se, src, spike_pred());
+ for (auto s: util::make_range(sources)) queue.push_back(cn->make_event(s));
sp = sources.first;
++cn;
}
}
else {
- auto cn = cons.begin();
- auto sp = spks.begin();
- while (cn!=cons.end() && sp!=spks.end()) {
- auto targets = std::equal_range(cn, cons.end(), sp->source);
- for (auto c: make_range(targets)) {
- queues[c.index_on_domain()].push_back(c.make_event(*sp));
- }
-
+ while (cn!=ce && sp!=se) {
+ auto targets = std::equal_range(cn, ce, sp->source);
+ for (auto c: util::make_range(targets)) queues[c.index_on_domain].push_back(c.make_event(*sp));
cn = targets.first;
++sp;
}
@@ -218,6 +208,10 @@ std::uint64_t communicator::num_spikes() const {
return num_spikes_;
}
+void communicator::set_num_spikes(std::uint64_t n) {
+ num_spikes_ = n;
+}
+
cell_size_type communicator::num_local_cells() const {
return num_local_cells_;
}
diff --git a/arbor/communication/communicator.hpp b/arbor/communication/communicator.hpp
index 61872c86d911bc31d6f61f717a848b79677b72aa..7a4becd3b19295a174b929922353112ef8d1b7c8 100644
--- a/arbor/communication/communicator.hpp
+++ b/arbor/communication/communicator.hpp
@@ -62,6 +62,7 @@ public:
/// Returns the total number of global spikes over the duration of the simulation
std::uint64_t num_spikes() const;
+ void set_num_spikes(std::uint64_t n);
cell_size_type num_local_cells() const;
@@ -69,10 +70,16 @@ public:
void reset();
+ void update_connections(const connectivity& rec,
+ const domain_decomposition& dom_dec,
+ const label_resolution_map& source_resolution_map,
+ const label_resolution_map& target_resolution_map);
+
private:
- cell_size_type num_local_cells_;
- cell_size_type num_local_groups_;
- cell_size_type num_domains_;
+ cell_size_type num_total_cells_ = 0;
+ cell_size_type num_local_cells_ = 0;
+ cell_size_type num_local_groups_ = 0;
+ cell_size_type num_domains_ = 0;
std::vector<connection> connections_;
std::vector<cell_size_type> connection_part_;
std::vector<cell_size_type> index_divisions_;
diff --git a/arbor/connection.hpp b/arbor/connection.hpp
index 0779fda646d2daab5dab89bc662ccf58165f83cd..7fab0162bf35ffa87652b941181b06ef2d2fe3fa 100644
--- a/arbor/connection.hpp
+++ b/arbor/connection.hpp
@@ -15,52 +15,33 @@ public:
float w,
float d,
cell_gid_type didx=cell_gid_type(-1)):
- source_(src),
- destination_(dest),
- weight_(w),
- delay_(d),
- index_on_domain_(didx)
+ source(src),
+ destination(dest),
+ weight(w),
+ delay(d),
+ index_on_domain(didx)
{}
- float weight() const { return weight_; }
- time_type delay() const { return delay_; }
+ spike_event make_event(const spike& s) { return { destination, s.time + delay, weight}; }
- cell_member_type source() const { return source_; }
- cell_lid_type destination() const { return destination_; }
- cell_size_type index_on_domain() const { return index_on_domain_; }
-
- spike_event make_event(const spike& s) {
- return {destination_, s.time + delay_, weight_};
- }
-
-private:
- cell_member_type source_;
- cell_lid_type destination_;
- float weight_;
- float delay_;
- cell_size_type index_on_domain_;
+ cell_member_type source;
+ cell_lid_type destination;
+ float weight;
+ float delay;
+ cell_size_type index_on_domain;
};
// connections are sorted by source id
// these operators make for easy interopability with STL algorithms
-
-static inline bool operator<(const connection& lhs, const connection& rhs) {
- return lhs.source() < rhs.source();
-}
-
-static inline bool operator<(const connection& lhs, cell_member_type rhs) {
- return lhs.source() < rhs;
-}
-
-static inline bool operator<(cell_member_type lhs, const connection& rhs) {
- return lhs < rhs.source();
-}
+static inline bool operator<(const connection& lhs, const connection& rhs) { return lhs.source < rhs.source; }
+static inline bool operator<(const connection& lhs, cell_member_type rhs) { return lhs.source < rhs; }
+static inline bool operator<(cell_member_type lhs, const connection& rhs) { return lhs < rhs.source; }
} // namespace arb
static inline std::ostream& operator<<(std::ostream& o, arb::connection const& con) {
- return o << "con [" << con.source() << " -> " << con.destination()
- << " : weight " << con.weight()
- << ", delay " << con.delay()
- << ", index " << con.index_on_domain() << "]";
+ return o << "con [" << con.source << " -> " << con.destination
+ << " : weight " << con.weight
+ << ", delay " << con.delay
+ << ", index " << con.index_on_domain << "]";
}
diff --git a/arbor/domain_decomposition.cpp b/arbor/domain_decomposition.cpp
index de443ab3d76f50881e16dda57dfe2418e6817f87..aa22082120a42c53a1359d8e1ea11a7c94fef2a7 100644
--- a/arbor/domain_decomposition.cpp
+++ b/arbor/domain_decomposition.cpp
@@ -15,7 +15,7 @@
namespace arb {
domain_decomposition::domain_decomposition(
const recipe& rec,
- const context& ctx,
+ context ctx,
const std::vector<group_description>& groups)
{
struct partition_gid_domain {
diff --git a/arbor/execution_context.cpp b/arbor/execution_context.cpp
index 0437024b7b759d336a097ed79773570bc9122c69..4b555f1c29158ea1d155d4affdde98061347b91f 100644
--- a/arbor/execution_context.cpp
+++ b/arbor/execution_context.cpp
@@ -14,10 +14,6 @@
namespace arb {
-void execution_context_deleter::operator()(execution_context* p) const {
- delete p;
-}
-
execution_context::execution_context(const proc_allocation& resources):
distributed(make_local_context()),
thread_pool(std::make_shared<threading::task_system>(resources.num_threads)),
@@ -26,7 +22,7 @@ execution_context::execution_context(const proc_allocation& resources):
{}
ARB_ARBOR_API context make_context(const proc_allocation& p) {
- return context(new execution_context(p));
+ return std::make_shared<execution_context>(p);
}
#ifdef ARB_HAVE_MPI
@@ -40,7 +36,7 @@ execution_context::execution_context(const proc_allocation& resources, MPI_Comm
template <>
ARB_ARBOR_API context make_context<MPI_Comm>(const proc_allocation& p, MPI_Comm comm) {
- return context(new execution_context(p, comm));
+ return std::make_shared<execution_context>(p, comm);
}
#endif
template <>
@@ -55,30 +51,30 @@ execution_context::execution_context(
template <>
ARB_ARBOR_API context make_context(const proc_allocation& p, dry_run_info d) {
- return context(new execution_context(p, d));
+ return std::make_shared<execution_context>(p, d);
}
-ARB_ARBOR_API std::string distribution_type(const context& ctx) {
+ARB_ARBOR_API std::string distribution_type(context ctx) {
return ctx->distributed->name();
}
-ARB_ARBOR_API bool has_gpu(const context& ctx) {
+ARB_ARBOR_API bool has_gpu(context ctx) {
return ctx->gpu->has_gpu();
}
-ARB_ARBOR_API unsigned num_threads(const context& ctx) {
+ARB_ARBOR_API unsigned num_threads(context ctx) {
return ctx->thread_pool->get_num_threads();
}
-ARB_ARBOR_API unsigned num_ranks(const context& ctx) {
+ARB_ARBOR_API unsigned num_ranks(context ctx) {
return ctx->distributed->size();
}
-ARB_ARBOR_API unsigned rank(const context& ctx) {
+ARB_ARBOR_API unsigned rank(context ctx) {
return ctx->distributed->id();
}
-ARB_ARBOR_API bool has_mpi(const context& ctx) {
+ARB_ARBOR_API bool has_mpi(context ctx) {
return ctx->distributed->name() == "MPI";
}
diff --git a/arbor/include/arbor/context.hpp b/arbor/include/arbor/context.hpp
index 75ba004bc5c5a92a86a1c58aaddbf2d0d2459bdb..f4fadc389d40c97835d9511f7d389f51f122eaf4 100644
--- a/arbor/include/arbor/context.hpp
+++ b/arbor/include/arbor/context.hpp
@@ -47,14 +47,8 @@ struct proc_allocation {
struct execution_context;
// arb::context is an opaque handle for the execution context for use
-// in the public API, implemented as a unique pointer.
-//
-// As execution_context is an incomplete type, an explicit deleter must be
-// provided.
-struct ARB_ARBOR_API execution_context_deleter {
- void operator()(execution_context*) const;
-};
-using context = std::unique_ptr<execution_context, execution_context_deleter>;
+// in the public API, implemented as a shared pointer.
+using context = std::shared_ptr<execution_context>;
// Helpers for creating contexts. These are implemented in the back end.
@@ -68,11 +62,11 @@ ARB_ARBOR_API context make_context(const proc_allocation& resources, Comm comm);
// Queries for properties of execution resources in a context.
-ARB_ARBOR_API std::string distribution_type(const context&);
-ARB_ARBOR_API bool has_gpu(const context&);
-ARB_ARBOR_API unsigned num_threads(const context&);
-ARB_ARBOR_API bool has_mpi(const context&);
-ARB_ARBOR_API unsigned num_ranks(const context&);
-ARB_ARBOR_API unsigned rank(const context&);
+ARB_ARBOR_API std::string distribution_type(context);
+ARB_ARBOR_API bool has_gpu(context);
+ARB_ARBOR_API unsigned num_threads(context);
+ARB_ARBOR_API bool has_mpi(context);
+ARB_ARBOR_API unsigned num_ranks(context);
+ARB_ARBOR_API unsigned rank(context);
}
diff --git a/arbor/include/arbor/domain_decomposition.hpp b/arbor/include/arbor/domain_decomposition.hpp
index 0de6383ae7bf77341008386413eb890c4715a6d3..21accd76c89410be3c3fbd624640df18b9c494ea 100644
--- a/arbor/include/arbor/domain_decomposition.hpp
+++ b/arbor/include/arbor/domain_decomposition.hpp
@@ -36,7 +36,10 @@ struct group_description {
class ARB_ARBOR_API domain_decomposition {
public:
domain_decomposition() = delete;
- domain_decomposition(const recipe& rec, const context& ctx, const std::vector<group_description>& groups);
+ domain_decomposition(const recipe& rec, context ctx, const std::vector<group_description>& groups);
+
+ domain_decomposition(const domain_decomposition&) = default;
+ domain_decomposition& operator=(const domain_decomposition&) = default;
int gid_domain(cell_gid_type gid) const;
int num_domains() const;
diff --git a/arbor/include/arbor/load_balance.hpp b/arbor/include/arbor/load_balance.hpp
index 1c07398095294f2ff0924ff2a5eed203d4a8283f..31614403a40665e4024cce42466b48c1424ae64e 100644
--- a/arbor/include/arbor/load_balance.hpp
+++ b/arbor/include/arbor/load_balance.hpp
@@ -22,6 +22,6 @@ using partition_hint_map = std::unordered_map<cell_kind, partition_hint>;
ARB_ARBOR_API domain_decomposition partition_load_balance(
const recipe& rec,
- const context& ctx,
+ context ctx,
const partition_hint_map& hint_map = {});
} // namespace arb
diff --git a/arbor/include/arbor/profile/meter_manager.hpp b/arbor/include/arbor/profile/meter_manager.hpp
index 985e43c6e80c831fdfeacd99f3162a3a7aee200a..1b57c44f75bf2f14a7fa752eaceab0f20173f157 100644
--- a/arbor/include/arbor/profile/meter_manager.hpp
+++ b/arbor/include/arbor/profile/meter_manager.hpp
@@ -26,7 +26,7 @@ struct measurement {
std::string name;
std::string units;
std::vector<std::vector<double>> measurements;
- measurement(std::string, std::string, const std::vector<double>&, const context&);
+ measurement(std::string, std::string, const std::vector<double>&, context);
};
class ARB_ARBOR_API meter_manager {
@@ -41,8 +41,8 @@ private:
public:
meter_manager();
- void start(const context& ctx);
- void checkpoint(std::string name, const context& ctx);
+ void start(context ctx);
+ void checkpoint(std::string name, context ctx);
const std::vector<std::unique_ptr<meter>>& meters() const;
const std::vector<std::string>& checkpoint_names() const;
@@ -59,7 +59,7 @@ struct meter_report {
std::vector<std::string> hosts;
};
-ARB_ARBOR_API meter_report make_meter_report(const meter_manager& manager, const context& ctx);
+ARB_ARBOR_API meter_report make_meter_report(const meter_manager& manager, context ctx);
ARB_ARBOR_API std::ostream& operator<<(std::ostream& o, const meter_report& report);
} // namespace profile
diff --git a/arbor/include/arbor/profile/profiler.hpp b/arbor/include/arbor/profile/profiler.hpp
index 8d7d02c7b59cde44ed2a350d8231e4c3cf029f2d..fd8dd980c59fa9d4e5c9988b549b50f9a12a171a 100644
--- a/arbor/include/arbor/profile/profiler.hpp
+++ b/arbor/include/arbor/profile/profiler.hpp
@@ -36,7 +36,7 @@ struct profile {
// TODO: remove declaration and update the docs
void profiler_clear();
-ARB_ARBOR_API void profiler_initialize(context& ctx);
+ARB_ARBOR_API void profiler_initialize(context ctx);
ARB_ARBOR_API void profiler_enter(std::size_t region_id);
ARB_ARBOR_API void profiler_leave();
diff --git a/arbor/include/arbor/recipe.hpp b/arbor/include/arbor/recipe.hpp
index ed5e621850c9f9de547056278568741db21db822..bf17317e42d531e31f41fee83e23a56672225c99 100644
--- a/arbor/include/arbor/recipe.hpp
+++ b/arbor/include/arbor/recipe.hpp
@@ -61,28 +61,43 @@ struct gap_junction_connection {
peer(std::move(peer)), local(std::move(local)), weight(g) {}
};
-class ARB_ARBOR_API recipe {
-public:
- virtual cell_size_type num_cells() const = 0;
-
- // Cell description type will be specific to cell kind of cell with given gid.
- virtual util::unique_any get_cell_description(cell_gid_type gid) const = 0;
- virtual cell_kind get_cell_kind(cell_gid_type) const = 0;
-
- virtual std::vector<event_generator> event_generators(cell_gid_type) const {
+struct ARB_ARBOR_API has_gap_junctions {
+ virtual std::vector<gap_junction_connection> gap_junctions_on(cell_gid_type) const {
return {};
}
+};
+
+struct ARB_ARBOR_API has_synapses {
virtual std::vector<cell_connection> connections_on(cell_gid_type) const {
return {};
}
- virtual std::vector<gap_junction_connection> gap_junctions_on(cell_gid_type) const {
+};
+
+struct ARB_ARBOR_API has_probes {
+ virtual std::vector<probe_info> get_probes(cell_gid_type gid) const {
return {};
}
+};
- virtual std::vector<probe_info> get_probes(cell_gid_type gid) const {
+struct ARB_ARBOR_API has_generators {
+ virtual std::vector<event_generator> event_generators(cell_gid_type) const {
return {};
}
+};
+// Toppings allow updating a simulation
+struct ARB_ARBOR_API connectivity: public has_synapses, has_generators {
+ virtual ~connectivity() {}
+};
+
+// Recipes allow building a simulation by lazy queries
+struct ARB_ARBOR_API recipe: public has_gap_junctions, has_probes, connectivity {
+ // number of cells to build
+ virtual cell_size_type num_cells() const = 0;
+ // Cell description type will be specific to cell kind of cell with given gid.
+ virtual util::unique_any get_cell_description(cell_gid_type gid) const = 0;
+ // Query cell kind per gid
+ virtual cell_kind get_cell_kind(cell_gid_type) const = 0;
// Global property type will be specific to given cell kind.
virtual std::any get_global_properties(cell_kind) const { return std::any{}; };
diff --git a/arbor/include/arbor/simulation.hpp b/arbor/include/arbor/simulation.hpp
index 205c8f940f0edf775bf6afa646669021060e35a5..a96e0cd7a2277a350fe29d8e44c2af377faa9bae 100644
--- a/arbor/include/arbor/simulation.hpp
+++ b/arbor/include/arbor/simulation.hpp
@@ -28,11 +28,13 @@ class simulation_state;
class ARB_ARBOR_API simulation {
public:
- simulation(const recipe& rec, const context& ctx, const domain_decomposition& decomp);
+ simulation(const recipe& rec, context ctx, const domain_decomposition& decomp);
simulation(const recipe& rec,
- const context& ctx=make_context(),
- std::function<domain_decomposition(const recipe&, const context&)> balancer=[](auto& r, auto& c) { return partition_load_balance(r, c); }): simulation(rec, ctx, balancer(rec, ctx)) {}
+ context ctx=make_context(),
+ std::function<domain_decomposition(const recipe&, context)> balancer=[](auto& r, auto c) { return partition_load_balance(r, c); }): simulation(rec, ctx, balancer(rec, ctx)) {}
+
+ void update(const connectivity& rec);
void reset();
diff --git a/arbor/include/arbor/spike.hpp b/arbor/include/arbor/spike.hpp
index 69b1558ec959db60beb1ef6a0230a7c1730156ea..f804dfb5781c67d8bd34860acb224a99289d9ace 100644
--- a/arbor/include/arbor/spike.hpp
+++ b/arbor/include/arbor/spike.hpp
@@ -17,7 +17,7 @@ struct basic_spike {
basic_spike() = default;
basic_spike(id_type s, time_type t):
- source(s), time(t)
+ source(std::move(s)), time(t)
{}
friend bool operator==(const basic_spike& l, const basic_spike& r) {
diff --git a/arbor/label_resolution.hpp b/arbor/label_resolution.hpp
index c015524e4a1f80bc609597902d89d69e2e3976d6..dbfe2014a9ce7517512cd46020e69d9840cc9b5c 100644
--- a/arbor/label_resolution.hpp
+++ b/arbor/label_resolution.hpp
@@ -77,7 +77,7 @@ public:
lid_hopefully at(unsigned idx) const;
};
- label_resolution_map() = delete;
+ label_resolution_map() = default;
explicit label_resolution_map(const cell_labels_and_gids&);
const range_set& at(const cell_gid_type& gid, const cell_tag_type& tag) const;
diff --git a/arbor/partition_load_balance.cpp b/arbor/partition_load_balance.cpp
index 71467aabff1fc7e7cbb7a471a522f2b7046a1b74..8718c743d111eb83506ab8876321e4a804314bbc 100644
--- a/arbor/partition_load_balance.cpp
+++ b/arbor/partition_load_balance.cpp
@@ -22,7 +22,7 @@ namespace arb {
ARB_ARBOR_API domain_decomposition partition_load_balance(
const recipe& rec,
- const context& ctx,
+ context ctx,
const partition_hint_map& hint_map)
{
using util::make_span;
diff --git a/arbor/profile/meter_manager.cpp b/arbor/profile/meter_manager.cpp
index c2540c84267c9e08154113a515cf8638649091e6..2886c1c21dcb44e629e8259920a20b0a9db31eea 100644
--- a/arbor/profile/meter_manager.cpp
+++ b/arbor/profile/meter_manager.cpp
@@ -24,7 +24,7 @@ double mean(const C& c) {
measurement::measurement(std::string n, std::string u,
const std::vector<double>& readings,
- const context& ctx):
+ context ctx):
name(std::move(n)), units(std::move(u))
{
auto dist = ctx->distributed;
@@ -54,7 +54,7 @@ meter_manager::meter_manager() {
}
};
-void meter_manager::start(const context& ctx) {
+void meter_manager::start(context ctx) {
arb_assert(!started_);
started_ = true;
@@ -70,7 +70,7 @@ void meter_manager::start(const context& ctx) {
start_time_ = timer_type::tic();
};
-void meter_manager::checkpoint(std::string name, const context& ctx) {
+void meter_manager::checkpoint(std::string name, context ctx) {
arb_assert(started_);
// Record the time taken on this domain since the last checkpoint
@@ -101,7 +101,7 @@ const std::vector<double>& meter_manager::times() const {
// Build a report of meters, for use at the end of a simulation
// for output to file or analysis.
-ARB_ARBOR_API meter_report make_meter_report(const meter_manager& manager, const context& ctx) {
+ARB_ARBOR_API meter_report make_meter_report(const meter_manager& manager, context ctx) {
meter_report report;
// Add the times to the meter outputs
diff --git a/arbor/profile/profiler.cpp b/arbor/profile/profiler.cpp
index c5e0b641dfcb59ad2a9fe6896ff6e156d613b15e..be5b2aff19fd53e9b3978d91b2803af2e1047ce9 100644
--- a/arbor/profile/profiler.cpp
+++ b/arbor/profile/profiler.cpp
@@ -357,7 +357,7 @@ ARB_ARBOR_API void profiler_enter(region_id_type region_id) {
profiler::get_global_profiler().enter(region_id);
}
-ARB_ARBOR_API void profiler_initialize(context& ctx) {
+ARB_ARBOR_API void profiler_initialize(context ctx) {
profiler::get_global_profiler().initialize(ctx->thread_pool);
}
diff --git a/arbor/simulation.cpp b/arbor/simulation.cpp
index 562db162668b0a4f4a55e448058c4a73b9f5e501..a822ae616b1f798c76cbcc1402ec326009859021 100644
--- a/arbor/simulation.cpp
+++ b/arbor/simulation.cpp
@@ -90,7 +90,9 @@ ARB_ARBOR_API void merge_cell_events(
class simulation_state {
public:
- simulation_state(const recipe& rec, const domain_decomposition& decomp, execution_context ctx);
+ simulation_state(const recipe& rec, const domain_decomposition& decomp, context ctx);
+
+ void update(const connectivity& rec);
void reset();
@@ -116,6 +118,8 @@ public:
spike_export_function global_export_callback_;
spike_export_function local_export_callback_;
epoch_function epoch_callback_;
+ label_resolution_map source_resolution_map_;
+ label_resolution_map target_resolution_map_;
private:
// Record last computed epoch (integration interval).
@@ -137,6 +141,8 @@ private:
std::unordered_map<cell_gid_type, gid_local_info> gid_to_local_;
communicator communicator_;
+ context ctx_;
+ domain_decomposition ddc_;
task_system_handle task_system_;
@@ -183,11 +189,13 @@ private:
simulation_state::simulation_state(
const recipe& rec,
const domain_decomposition& decomp,
- execution_context ctx
+ context ctx
):
- task_system_(ctx.thread_pool),
- local_spikes_({thread_private_spike_store(ctx.thread_pool), thread_private_spike_store(ctx.thread_pool)})
-{
+ ctx_{ctx},
+ ddc_{decomp},
+ task_system_(ctx->thread_pool),
+ local_spikes_({thread_private_spike_store(ctx->thread_pool),
+ thread_private_spike_store(ctx->thread_pool)}) {
// Generate the cell groups in parallel, with one task per cell group.
auto num_groups = decomp.num_groups();
cell_groups_.resize(num_groups);
@@ -197,7 +205,7 @@ simulation_state::simulation_state(
[&](cell_group_ptr& group, int i) {
const auto& group_info = decomp.group(i);
cell_label_range sources, targets;
- auto factory = cell_kind_implementation(group_info.kind, group_info.backend, ctx);
+ auto factory = cell_kind_implementation(group_info.kind, group_info.backend, *ctx_);
group = factory(group_info.gids, rec, sources, targets);
cg_sources[i] = cell_labels_and_gids(std::move(sources), group_info.gids);
@@ -209,44 +217,44 @@ simulation_state::simulation_state(
local_sources.append(cg_sources.at(i));
local_targets.append(cg_targets.at(i));
}
- auto global_sources = ctx.distributed->gather_cell_labels_and_gids(local_sources);
-
- auto source_resolution_map = label_resolution_map(std::move(global_sources));
- auto target_resolution_map = label_resolution_map(std::move(local_targets));
+ auto global_sources = ctx->distributed->gather_cell_labels_and_gids(local_sources);
- communicator_ = arb::communicator(rec, decomp, source_resolution_map, target_resolution_map, ctx);
-
- const auto num_local_cells = communicator_.num_local_cells();
+ source_resolution_map_ = label_resolution_map(std::move(global_sources));
+ target_resolution_map_ = label_resolution_map(std::move(local_targets));
+ communicator_ = communicator(rec, ddc_, source_resolution_map_, target_resolution_map_, *ctx_);
+ update(rec);
+ epoch_.reset();
+}
+void simulation_state::update(const connectivity& rec) {
+ communicator_.update_connections(rec, ddc_, source_resolution_map_, target_resolution_map_);
// Use half minimum delay of the network for max integration interval.
t_interval_ = communicator_.min_delay()/2;
+ const auto num_local_cells = communicator_.num_local_cells();
// Initialize empty buffers for pending events for each local cell
pending_events_.resize(num_local_cells);
-
+ // Forget old generators, if present
+ event_generators_.clear();
event_generators_.resize(num_local_cells);
cell_size_type lidx = 0;
cell_size_type grpidx = 0;
-
- auto target_resolution_map_ptr = std::make_shared<label_resolution_map>(std::move(target_resolution_map));
- for (const auto& group_info: decomp.groups()) {
+ auto target_resolution_map_ptr = std::make_shared<label_resolution_map>(target_resolution_map_);
+ for (const auto& group_info: ddc_.groups()) {
for (auto gid: group_info.gids) {
// Store mapping of gid to local cell index.
- gid_to_local_[gid] = gid_local_info{lidx, grpidx};
-
- // Resolve event_generator targets.
- // Each event generator gets their own resolver state.
+ gid_to_local_[gid] = {lidx, grpidx};
+ // Resolve event_generator targets; each event generator gets their own resolver state.
auto event_gens = rec.event_generators(gid);
for (auto& g: event_gens) {
- g.resolve_label([target_resolution_map_ptr, event_resolver=resolver(target_resolution_map_ptr.get()), gid]
- (const cell_local_label_type& label) mutable {
+ g.resolve_label([target_resolution_map_ptr,
+ event_resolver=resolver(target_resolution_map_ptr.get()),
+ gid] (const cell_local_label_type& label) mutable {
return event_resolver.resolve({gid, label});
});
}
-
// Set up the event generators for cell gid.
event_generators_[lidx] = event_gens;
-
++lidx;
}
++grpidx;
@@ -257,10 +265,9 @@ simulation_state::simulation_state(
// the following epoch. In each buffer there is one lane for each local cell.
event_lanes_[0].resize(num_local_cells);
event_lanes_[1].resize(num_local_cells);
-
- epoch_.reset();
}
+
void simulation_state::reset() {
epoch_ = epoch();
@@ -510,16 +517,18 @@ void simulation_state::inject_events(const cse_vector& events) {
simulation::simulation(
const recipe& rec,
- const context& ctx,
+ context ctx,
const domain_decomposition& decomp)
{
- impl_.reset(new simulation_state(rec, decomp, *ctx));
+ impl_.reset(new simulation_state(rec, decomp, ctx));
}
void simulation::reset() {
impl_->reset();
}
+void simulation::update(const connectivity& rec) { impl_->update(rec); }
+
time_type simulation::run(time_type tfinal, time_type dt) {
if (dt <= 0.0) {
throw domain_error("Finite time-step must be supplied.");
diff --git a/doc/concepts/domdec.rst b/doc/concepts/domdec.rst
index 835152f2a9c332448185e9beff48c6e3e4037b3e..f4dac781a639972c421f640d563838be4eadc66a 100644
--- a/doc/concepts/domdec.rst
+++ b/doc/concepts/domdec.rst
@@ -3,30 +3,64 @@
Domain decomposition
====================
-An Arbor simulation requires a :ref:`modelrecipe`, a :ref:`(hardware) context <modelhardware>`, and a domain decomposition. The Recipe contains the neuroscientific model, the hardware context describes the computational resources you are going to execute the simulation on, and the domain decomposition describes how Arbor will use the hardware. Since the context and domain decomposition may seem closely related at first, it might be instructive to see how recipes are used by Arbor:
+An Arbor simulation requires a :ref:`modelrecipe`, a :ref:`(hardware) context
+<modelhardware>`, and a domain decomposition. The Recipe contains the
+neuroscientific model, the hardware context describes the computational
+resources you are going to execute the simulation on, and the domain
+decomposition describes how Arbor will use the hardware. Since the context and
+domain decomposition may seem closely related at first, it might be instructive
+to see how recipes are used by Arbor:
.. raw:: html
:file: domdec-diag-1.html
-A *domain decomposition* describes the distribution of the model over the available computational resources.
-The description partitions the cells in the model as follows:
+A *domain decomposition* describes the distribution of the model over the
+available computational resources. The description partitions the cells in the
+model as follows:
-* group the cells into cell groups of the same :ref:`kind <modelcellkind>` of cell;
+* group the cells into cell groups of the same :ref:`kind <modelcellkind>` of
+ cell;
* assign each cell group to either a CPU core or GPU on a specific MPI rank.
-The number of cells in each cell group depends on different factors, including the type of the cell, and whether the
-cell group will run on a CPU core or the GPU. The domain decomposition is solely responsible for describing the distribution
-of cells across cell groups and domains.
+The number of cells in each cell group depends on different factors, including
+the type of the cell, and whether the cell group will run on a CPU core or the
+GPU. The domain decomposition is solely responsible for describing the
+distribution of cells across cell groups and domains.
-The domain decomposition can be built manually by the modeler, or an automatic load balancer can be used.
+The domain decomposition can be built manually by the modeler, or an automatic
+load balancer can be used.
+We define some terms as used in the context of connectivity
+
+.. glossary::
+
+ connection
+ Tuple of ``(source, target, weight, delay)`` describing an
+ axon/synapse connection as travelling time (`delay`) and attenuation
+ (`weight`) between two sites `source = (gid, spike_detector)` and `target
+ = (gid, synapse)` where `spike_detector` and `synapse` are string labels.
+
+ cell_group
+ List of same-kinded cells that share some information. Must not be split
+ across domains.
+
+ domain
+ Produced by a `load_balancer`, a list of all `cell_groups`
+ located on the same hardware. A `communicator` deals with the full set of
+ cells of one `domain`.
+
+ domain_decomposition
+ List of domains; distributed across MPI processes.
Load balancers
--------------
-A *load balancer* generates the domain decomposition using the model :ref:`recipe <modelrecipe>` and a description of
-the available computational resources on which the model will run described by an :ref:`execution context <modelcontext>`.
-Currently Arbor provides one load balancer and more will be added over time.
+A *load balancer* generates the domain decomposition using the model
+:ref:`recipe <modelrecipe>` and a description of the available computational
+resources on which the model will run described by an :ref:`execution context
+<modelcontext>`. Currently Arbor provides one automatic load balancer and a
+method for manually listing out the cells per MPI task. More approaches might be
+added over time.
API
---
diff --git a/doc/concepts/interconnectivity.rst b/doc/concepts/interconnectivity.rst
index 4fee5f696332593f582e29b1d02f442c351850db..8c56064801cb1bfbe9467165a2c6d685a6220f10 100644
--- a/doc/concepts/interconnectivity.rst
+++ b/doc/concepts/interconnectivity.rst
@@ -12,9 +12,67 @@ Connections only capture the propagation delay and attenuation associated with s
connectivity: the biophysical modelling of the chemical synapses themselves is the
responsibility of the target cell model.
-Connection sites and gap junction sites are defined on locations on cells as part of the
-:ref:`cell description <modelcelldesc>`.
-A recipe lets you define which sites are connected to which.
+Connection sites and gap junction sites are defined on locations on cells as
+part of the :ref:`cell description <modelcelldesc>`.
+These sites as such are not connected yet, however the :ref:`recipe <modelrecipe>`
+exposes a number of callbacks to form connections and gap junctions between sites.
+The recipe callbacks are interrogated during simulation creation.
+
+In addition, simulations may update their connectivity by building a new
+connection table outside calls to `run`, for example
+
+.. code-block:: python
+
+ rec = recipe()
+ dec = arb.domain_decomposition(rec, ctx)
+ sim = arb.simulation(rec, ctx, dec)
+
+ # run simulation for 0.25ms with the basic connectivity
+ sim.run(0.25, 0.025)
+
+ # extend the recipe to more connections
+ rec.add_connections()
+ # use `connections_on` to build a new connection table
+ sim.update_connections(rec)
+
+ # run simulation for 0.25ms with the extended connectivity
+ sim.run(0.5, 0.025)
+
+This will completely replace the old table, previous connections to be retained
+must be explicitly included in the updated callback. This can also be used to
+update connection weights and delays. Note, however, that there is currently no
+way to introduce new sites to the simulation, nor any changes to gap junctions.
+
+The ``update_connections`` method accepts either a full ``recipe`` (but will
+**only** use the ``connections_on`` and ``events_generators`` callbacks) or a
+``connectivity``, which is a reduced recipe exposing only the relevant callbacks.
+Currently ``connectivity`` is only available in C++; Python users have to pass a
+full recipe.
+
+.. warn::
+
+ The semantics of connection updates are subtle and might produce surprising
+ results if handled carelessly. In particular, spikes in-flight over a
+ connection will *always* be delivered, even if the connection has been
+ deleted before the time of delivery has passed (`= t_emitted +
+ connection_delay`). As Arbor's connection model joins processes on the axon,
+ the synaptic cleft, and the receiving synapse into a simple pair `(weight,
+ delay)` it is unclear 'where' the action potential is located at the time of
+ deletion relative to the locus of disconnection. Thus, it was decided to
+ deliver spike events regardless. This is will not cause issues when the
+ transition is slow and smooth, ie weights decays over time towards a small
+ value and then the connection is removed. However, drastic and/or frequent
+ changes across busy synapses might cause unexpected behaviour.
+
+.. note::
+
+ Arbor uses a lazily constructed network (from the ``recipe`` callbacks) for
+ good reason; storing the full connectivity (for all ``gids``) in the
+ ``recipe`` can lead to probitivively large memory footprints. Keep this in
+ mind when designing your connectivity and heed the consequences of doing I/O
+ in these callbacks. This is doubly important when using models with dynamic
+ connectivity where the temptation to store all connections is even larger and
+ each call to ``update`` will re-evaluate the corresponding callbacks.
.. _modelconnections:
@@ -35,6 +93,14 @@ A recipe lets you define which sites are connected to which.
a :gen:`global_label`; a target identified using a :gen:`local_label` (:gen:`gid` of target is
the argument of the recipe method); a connection delay and a connection weight.
+ .. code-block:: python
+
+ def connections_on(self, gid):
+ if gid + 1 < self.num_cells():
+ return [arbor.connection((gid + 1, "spike-source"), "synapse", weight, delay)]
+ else:
+ return []
+
spike
action potential
Spikes travel over :term:`connections <connection>`. In a synapse, they generate an event.
@@ -71,14 +137,34 @@ A recipe lets you define which sites are connected to which.
Similarly to `Connections`, Gap Junctions in Arbor are defined in two steps:
- 1. Create labeled **gap junction sites** on two separate cells as part of their
- :ref:`cell descriptions <modelcelldesc>` in the :ref:`recipe <modelrecipe>`.
- Each labeled group of gap junctions may contain multiple items on possibly multiple locations
- on the cell.
- 2. Declare the Gap Junction connections in the recipe *on the local cell*: from a peer **gap junction site**
- identified using a :gen:`global_label`; to a local **gap junction site** identified using a
- :gen:`local_label` (:gen:`gid` of the site is implicitly known); and a unit-less connection weight.
- Two of these connections are needed, on each of the peer and local cells.
+ 1. Create labeled **gap junction sites** on two separate cells as part of
+ their :ref:`cell descriptions <modelcelldesc>` in the :ref:`recipe
+ <modelrecipe>`. Each labeled group of gap junctions may contain multiple
+ items on possibly multiple locations on the cell.
+ 2. Declare the Gap Junction connections in the recipe *on the local cell*:
+ from a peer **gap junction site** identified using a
+ :gen:`global_label`; to a local **gap junction site** identified using
+ a :gen:`local_label` (:gen:`gid` of the site is implicitly known); and
+ a unit-less connection weight. Two of these connections are needed, on
+ each of the peer and local cells. The callback `gap_junctions_on`
+ returns a list of these items, eg
+
+ .. code-block:: python
+
+ def gap_junctions_on(self, gid):
+ n = self.num_cells
+ if gid + 1 < n and gid > 0:
+ return [arbor.gap_junction_connection((gid + 1, "gj"), "gj", weight),
+ arbor.gap_junction_connection((gid - 1, "gj"), "gj", weight),]
+ elif gid + 1 < n:
+ return [arbor.gap_junction_connection((gid + 1, "gj"), "gj", weight),]
+ if gid > 0:
+ return [arbor.gap_junction_connection((gid - 1, "gj"), "gj", weight),]
+ else:
+ return []
+
+ Note that gap junction connections are symmetrical and thus the above
+ example generates two connections, one incoming and one outgoing.
.. Note::
Only cable cells support gap junctions as of now.
diff --git a/doc/dev/communication.rst b/doc/dev/communication.rst
new file mode 100644
index 0000000000000000000000000000000000000000..81661345337bec61e45ac0b45680c47a53ee02f5
--- /dev/null
+++ b/doc/dev/communication.rst
@@ -0,0 +1,300 @@
+.. _communication:
+
+Communication in Arbor
+======================
+
+Communication between cells (and thus cell groups) is facilitated using discrete
+events which we call `spikes` in analogy to the underlying neurobiological
+process. Spikes are the only way to communicate between different cell kinds, eg
+cable cells and integrate-and-fire. In accordance to current theory, spikes are
+reduced to a simple time value. The connection is abstracted into a weight and a
+delay; modelling all axonal processes. While this may seem crude, it is a well
+supported model and commonly used in neuronal simulations.
+
+Connections are formed between sources (cable cell: spike detectors) and targets
+(cable cell: synapses). During runtime, events from all sources are concatenated
+on all MPI ranks using ``Allgatherv`` and targets are responsible for selecting
+events they have subscribed to. This is optimised for by sorting events locally
+(by source) and relying on the process layout to convert this into a *globally*
+sorted array.
+
+Exchange of Spikes
+==================
+
+We will start by discussing the exchange of spikes, an extended summary of the
+full communicator class can be found in the next section. Communication of
+spikes is facilitated through the ``communicator::exchange`` method. which is
+
+.. code-block:: c++
+
+ // aka cell_member_type, lexicographic order provided
+ struct id_type {
+ gid_type gid;
+ lid_type index;
+ };
+
+ struct spike {
+ id_type source;
+ double time;
+ };
+
+ using spike_vec = std::vector<spike>;
+
+ // From gathered_vector
+ template <typename T>
+ struct gathered_vector {
+ std::vector<T> values_;
+ std::vector<int> partition_;
+ };
+
+ using g_spike_vec = gathered_vector<spike>;
+
+ g_spike_vec communicator::exchange(spike_vec local_spikes) { // Take by value, since we modify (sort) anyhow.
+ // sort the spikes in ascending order of source gid
+ util::sort_by(local_spikes, [](spike s){ return s.source; });
+ // global all-to-all to gather a local copy of the global spike list on each node.
+ auto global_spikes = distributed_->gather_spikes(local_spikes);
+ num_spikes_ += global_spikes.size();
+ return global_spikes;
+ }
+
+ // From mpi_context
+ g_spike_vec gather_spikes(const spike_ve& local_spikes) const {
+ return mpi::gather_all_with_partition(local_spikes, comm_);
+ }
+
+ // From mpi
+ g_spike_vec gather_all_with_partition(const spike_vec& values, MPI_Comm comm) {
+ using traits = mpi_traits<spike>;
+ // Collect all sender's counts, scale to bytes
+ std::vector<int> counts = gather_all(int(values.size()), comm);
+ for (auto& c: counts) c *= traits::count();
+ // Scan-left to get offsets
+ std::vector<int> displs;
+ util::make_partition(displs, counts);
+ // Allocate recv buffer and exchange
+ spike_vec buffer(displs.back()/traits::count());
+ MPI_Allgatherv(// send buffer count to send MPI datatype of spike
+ values.data(), counts[rank(comm)], traits::mpi_type(),
+ // recv buffer count of each sender offset of senders MPI datatype of spike
+ buffer.data(), counts.data(), displs.data(), traits::mpi_type(),
+ comm);
+ // Scale back to sizeof(spike)
+ for (auto& d: displs) d /= traits::count();
+ return {std::move(buffer), std::move(displs)};
+ }
+
+Note that these snippets have been simplified and shortened from the
+actual code, they are intended for education only.
+
+After ``exchange`` is done, each process has received an object like
+this
+
+.. code-block:: c++
+
+ {
+ spikes: { {s0, t0}, ...},
+ offsets: { 0, ... }
+ }
+
+Now ``spikes`` is the array of all spikes during the last *epoch* where
+each sub-array of spikes is sorted, ie between ``offsets[ix]`` and
+``offest[ix+1]``. The ``offsets`` array has a length of MPI task count
+and its ``i``'th element gives the position of the first spike sent by
+task ``i``.
+
+Distribution of Events to Targets
+=================================
+
+Having received the generated spikes, the concatenated data is converted
+into events on each local cell group. This is done asynchronously with
+computation of the next cell state. In ``simulation.cpp`` we find
+
+.. code-block:: c++
+
+ auto exchange = [this](epoch prev) {
+ // Collate locally generated spikes.
+ auto all_local_spikes = local_spikes(prev.id).gather();
+ // Gather generated spikes across all ranks.
+ auto global_spikes = communicator_.exchange(all_local_spikes);
+ // Append events formed from global spikes to per-cell pending event queues.
+ communicator_.make_event_queues(global_spikes, pending_events_);
+ };
+
+which uses this
+
+.. code-block:: c++
+
+ // Check each global spike in turn to see it generates local events.
+ // If so, make the events and insert them into the appropriate event list.
+ //
+ // Takes reference to a vector of event lists as an argument, with one list
+ // for each local cell group. On completion, the events in each list are
+ // all events that must be delivered to targets in that cell group as a
+ // result of the global spike exchange, plus any events that were already
+ // in the list.
+ void make_event_queues(const gathered_vector<spike>& global_spikes,
+ std::vector<pse_vector>& queues) {
+ // Predicate for partitioning
+ struct spike_pred {
+ bool operator()(const spike& spk, const cell_member_type& src) { return spk.source < src; }
+ bool operator()(const cell_member_type& src, const spike& spk) { return src < spk.source; }
+ };
+
+ const auto& sp = global_spikes.partition();
+ for (auto dom: util::make_span(num_domains_)) {
+ // Fetch connections and spikes per integration domain
+ auto cons = util::subrange_view(connections_, connection_part_[dom], connection_part_[dom + 1]);
+ auto spks = util::subrange_view(global_spikes.values(), sp[dom], sp[dom + 1]);
+ auto sp = spks.begin(), se = spks.end();
+ auto cn = cons.begin(), ce = cons.end();
+ // We have a choice of whether to walk spikes or connections:
+ // i.e., we can iterate over the spikes, and for each spike search
+ // the for connections that have the same source; or alternatively
+ // for each connection, we can search the list of spikes for spikes
+ // with the same source.
+ //
+ // We iterate over whichever set is the smallest, which has
+ // complexity of order max(S log(C), C log(S)), where S is the
+ // number of spikes, and C is the number of connections.
+ if (cons.size() < spks.size()) {
+ while (cn != ce && sp != ce) {
+ auto src = cn->source(); // Source for connection
+ auto cix = cn->index_on_domain(); // Queue index for connection
+ // Given a source src split the range [sp, spks.end) into a pair sources=[l, h]
+ // st *l is the last element not smaller than src
+ // and *h is the first element greater than src.
+ // 'Greater' and 'smaller' are defined via the predicate above.
+ // The range [sp, spks.end) must be (partially) ordered wrt the predicate.
+ auto sources = std::equal_range(sp, se, src, spike_pred());
+ // Consequently, the range returned is the range of equal spike sources,
+ // we pick out ours and add all of them to the appropriate queue.
+ for (auto s: util::make_range(sources)) queues[cix].push_back(cn->make_event(s));
+ // now, move to next
+ sp = sources.first;
+ ++cn;
+ }
+ }
+ else {
+ while (cn != ce && sp != se) {
+ auto targets = std::equal_range(cn, ce, sp->source);
+ for (auto c: util::make_range(targets)) queues[c.index_on_domain()].push_back(c.make_event(*sp));
+ cn = targets.first;
+ ++sp;
+ }
+ }
+ }
+ }
+
+After ``make_event_queues`` there is one queue per cell and each queue
+is filled with a time ordered list of events for that cell. We now need
+to understand the actual connection table stored in
+
+.. code-block:: c++
+
+ struct connection {
+ spike_event make_event(const spike& s) {
+ return { destination_, s.time + delay_, weight_};
+ }
+
+ cell_member_type source;
+ cell_lid_type destination;
+ float weight;
+ float delay;
+ cell_size_type index_on_domain;
+ };
+
+ struct communicator {
+ // [...]
+ cell_size_type num_domains_;
+ std::vector<connection> connections_;
+ std::vector<cell_size_type> connection_part_;
+ // [...]
+ };
+
+The ``connections`` vector is a list of connections partitioned by the
+domain (as in domain decomposition) of their source's ``gid``, while
+``connection_part`` stores the partioning indices.
+
+Building the Connection Table
+=============================
+
+The table of connections on the local rank is built during the construction of
+the ``communicator`` object
+
+.. code-block:: c++
+
+ communicator::communicator(const recipe& rec,
+ const domain_decomposition& dom_dec,
+ const label_resolution_map& source_resolution_map,
+ const label_resolution_map& target_resolution_map,
+ execution_context& ctx);
+
+After that process,
+
+.. code-block:: c++
+
+ struct communicator {
+ // ...
+ std::vector<connection> connections_;
+ std::vector<cell_size_type> connection_part_;
+ };
+
+will contain all connections in ``connections_`` partitioned by the
+domain of the source's ``gid`` in ``dom_dec``. Beginnings of the
+respective partitions are pointed to by the indices in
+``connection_part_``.
+
+The algorithm for building is slightly obscured by caching and the use
+of labels and resolving them via ``target_/source_resolution_map`` to
+local ids on the respective source and target cells.
+
+.. note::
+
+ The ``label_resolution_map`` class is used to translate from labels at the
+ user facing API layers to Arbor's internal mappings in the vein of
+ ``(cell_gid, item_offset)``, where ``item_offset`` is an automatically
+ assigned integer ID. Textual labels are created by calls to ``place``
+ as in this example
+
+ .. code-block:: c++
+
+ auto d = arb::decor{};
+ d.place("..."_ls, arb::synapse{"..."}, "synapse-label");
+
+The construction is performed in-place
+
+.. code-block:: c++
+
+ // Allocate space for our connections
+ connections_.resize(n_cons);
+ // We have pre-computed `src_counts`, connection_part_ will now hold the starting indices
+ // of each `domain`.
+ util::make_partition(connection_part_, src_counts);
+ // Copy, as we use this as the list of the currently available next free target slots in
+ // `connections_`
+ auto offsets = connection_part_;
+ auto target_resolver = resolver(&target_resolution_map);
+ auto src_domain = src_domains.begin();
+ for (const auto& cell: gid_infos) {
+ auto index = cell.index_on_domain;
+ auto source_resolver = resolver(&source_resolution_map);
+ for (const auto& c: cell.conns) {
+ // Compute index representation of labels
+ auto src_lid = source_resolver.resolve(c.source);
+ auto tgt_lid = target_resolver.resolve({cell.gid, c.dest});
+ // Get offset of current source and bump to next free slot
+ auto offset = offsets[*src_domain]++;
+ // Write connection info into slot
+ connections_[offset] = {{c.source.gid, src_lid}, tgt_lid, c.weight, c.delay, index};
+ // Next source domain
+ ++src_domain;
+ }
+ }
+ // Now
+ // * all slots in `connections_` are filled.
+ // * `offsets` points at the ends of each partition.
+
+
+Next, each *partition* is sorted independently according to their
+source's ``gid``.
diff --git a/doc/dev/index.rst b/doc/dev/index.rst
index 0174b15f05839b48ee98987f8523c78e16dcf015..9de6651c00cceb6a33f90d01afaa6c1bf7dd6952 100644
--- a/doc/dev/index.rst
+++ b/doc/dev/index.rst
@@ -16,6 +16,7 @@ Here we document internal components of Arbor. These pages can be useful if you'
cable_cell
cell_groups
+ communication
debug
matrix_solver
simd_api
diff --git a/doc/python/simulation.rst b/doc/python/simulation.rst
index 533041d7dd8a0d45651196ac5fd171248bc6424a..818c73633d6599d752446833e1246f628d872231 100644
--- a/doc/python/simulation.rst
+++ b/doc/python/simulation.rst
@@ -63,12 +63,22 @@ over the local and distributed hardware resources (see :ref:`pydomdec`). Then, t
**Constructor:**
- .. function:: simulation(recipe, domain_decomposition, context)
+ .. function:: simulation(recipe, context, domain_decomposition)
- Initialize the model described by an :py:class:`arbor.recipe`, with cells and network distributed according to :py:class:`arbor.domain_decomposition`, and computational resources described by :py:class:`arbor.context`.
+ Initialize the model described by an :py:class:`arbor.recipe`, with
+ cells and network distributed according to
+ :py:class:`arbor.domain_decomposition`, and computational resources
+ described by :py:class:`arbor.context`.
**Updating Model State:**
+ .. function:: update_connections(recipe)
+
+ Rebuild the connection table as described by
+ :py:class:`arbor.recipe::connections_on` The recipe must differ **only**
+ in the return value of its :py:func:`connections_on` when compared to
+ the original recipe used to construct the simulation object.
+
.. function:: reset()
Reset the state of the simulation to its initial state.
diff --git a/example/CMakeLists.txt b/example/CMakeLists.txt
index 752737ec64ef94096dd29339f797f928d18cdb19..f13932683922fb197373aa82b92f613ffb8e98d9 100644
--- a/example/CMakeLists.txt
+++ b/example/CMakeLists.txt
@@ -11,5 +11,6 @@ add_subdirectory(ring)
add_subdirectory(gap_junctions)
add_subdirectory(single)
add_subdirectory(probe-demo)
+add_subdirectory(plasticity)
add_subdirectory(lfp)
add_subdirectory(diffusion)
diff --git a/example/brunel/brunel.cpp b/example/brunel/brunel.cpp
index 351e656fe03f4979ea3bc1403fb29d4c6fd0bcff..e1e19b3b8399cd6a50297ce4e05d9b2085fd8f80 100644
--- a/example/brunel/brunel.cpp
+++ b/example/brunel/brunel.cpp
@@ -70,7 +70,7 @@ std::ostream& operator<<(std::ostream& o, const cl_options& opt);
std::optional<cl_options> read_options(int argc, char** argv);
-void banner(const context& ctx);
+void banner(context ctx);
// Samples m unique values in interval [start, end) - gid.
// We exclude gid because we don't want self-loops.
@@ -249,7 +249,7 @@ int main(int argc, char** argv) {
simulation sim(recipe,
context,
- [&hints](auto& r, auto& c) { return partition_load_balance(r, c, hints); });
+ [&hints](auto& r, auto c) { return partition_load_balance(r, c, hints); });
// Set up spike recording.
std::vector<arb::spike> recorded_spikes;
@@ -296,7 +296,7 @@ int main(int argc, char** argv) {
return 0;
}
-void banner(const context& ctx) {
+void banner(context ctx) {
std::cout << "==========================================\n";
std::cout << " Brunel model miniapp\n";
std::cout << " - distributed : " << arb::num_ranks(ctx)
diff --git a/example/plasticity/CMakeLists.txt b/example/plasticity/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..904ea0943a4b09562d5a4d0423fccc5b95fdb72c
--- /dev/null
+++ b/example/plasticity/CMakeLists.txt
@@ -0,0 +1,4 @@
+add_executable(plasticity EXCLUDE_FROM_ALL plasticity.cpp)
+add_dependencies(examples plasticity)
+
+target_link_libraries(plasticity PRIVATE arbor arborio arborenv)
diff --git a/example/plasticity/branch_cell.hpp b/example/plasticity/branch_cell.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..527f26012c4a47dac2b285d25c66347f18054e55
--- /dev/null
+++ b/example/plasticity/branch_cell.hpp
@@ -0,0 +1,138 @@
+#pragma once
+
+#include <array>
+#include <random>
+
+#include <nlohmann/json.hpp>
+
+#include <arborio/label_parse.hpp>
+
+#include <arbor/cable_cell.hpp>
+#include <arbor/cable_cell_param.hpp>
+#include <arbor/common_types.hpp>
+#include <arbor/morph/segment_tree.hpp>
+
+#include <string>
+#include <sup/json_params.hpp>
+
+using namespace arborio::literals;
+
+// Parameters used to generate the random cell morphologies.
+struct cell_parameters {
+ cell_parameters() = default;
+
+ // Maximum number of levels in the cell (not including the soma)
+ unsigned max_depth = 5;
+
+ // The following parameters are described as ranges.
+ // The first value is at the soma, and the last value is used on the last level.
+ // Values at levels in between are found by linear interpolation.
+ std::array<double,2> branch_probs = {1.0, 0.5}; // Probability of a branch occuring.
+ std::array<unsigned,2> compartments = {20, 2}; // Compartment count on a branch.
+ std::array<double,2> lengths = {200, 20}; // Length of branch in μm.
+
+ // The number of synapses per cell.
+ unsigned synapses = 1;
+};
+
+cell_parameters parse_cell_parameters(nlohmann::json& json) {
+ cell_parameters params;
+ sup::param_from_json(params.max_depth, "depth", json);
+ sup::param_from_json(params.branch_probs, "branch-probs", json);
+ sup::param_from_json(params.compartments, "compartments", json);
+ sup::param_from_json(params.lengths, "lengths", json);
+ sup::param_from_json(params.synapses, "synapses", json);
+
+ return params;
+}
+
+// Helper used to interpolate in branch_cell.
+template <typename T>
+double interp(const std::array<T,2>& r, unsigned i, unsigned n) {
+ double p = i * 1./(n-1);
+ double r0 = r[0];
+ double r1 = r[1];
+ return r[0] + p*(r1-r0);
+}
+
+arb::cable_cell branch_cell(arb::cell_gid_type gid, const cell_parameters& params) {
+ arb::segment_tree tree;
+
+ // Add soma.
+ double srad = 12.6157/2.0; // soma radius
+ int stag = 1; // soma tag
+ tree.append(arb::mnpos, {0, 0,-srad, srad}, {0, 0, srad, srad}, stag); // For area of 500 μm².
+
+ std::vector<std::vector<unsigned>> levels;
+ levels.push_back({0});
+
+ // Standard mersenne_twister_engine seeded with gid.
+ std::mt19937 gen(gid);
+ std::uniform_real_distribution<double> dis(0, 1);
+
+ double drad = 0.5; // Diameter of 1 μm for each dendrite cable.
+ int dtag = 3; // Dendrite tag.
+
+ double dist_from_soma = srad; // Start dendrite at the edge of the soma.
+ for (unsigned i=0; i<params.max_depth; ++i) {
+ // Branch prob at this level.
+ double bp = interp(params.branch_probs, i, params.max_depth);
+ // Length at this level.
+ double l = interp(params.lengths, i, params.max_depth);
+ // Number of compartments at this level.
+ unsigned nc = std::round(interp(params.compartments, i, params.max_depth));
+
+ std::vector<unsigned> sec_ids;
+ for (unsigned sec: levels[i]) {
+ for (unsigned j=0; j<2; ++j) {
+ if (dis(gen)<bp) {
+ auto z = dist_from_soma;
+ auto dz = l/nc;
+ auto p = sec;
+ for (unsigned k=0; k<nc; ++k) {
+ p = tree.append(p, {0,0,z+(k+1)*dz, drad}, dtag);
+ }
+ sec_ids.push_back(p);
+ }
+ }
+ }
+ if (sec_ids.empty()) {
+ break;
+ }
+ levels.push_back(sec_ids);
+
+ dist_from_soma += l;
+ }
+
+ arb::label_dict labels;
+
+ using arb::reg::tagged;
+ labels.set("soma", tagged(stag));
+ labels.set("dend", tagged(dtag));
+
+ arb::decor decor;
+
+ decor.paint("soma"_lab, arb::density("hh"));
+ decor.paint("dend"_lab, arb::density("pas"));
+
+ decor.set_default(arb::axial_resistivity{100}); // [Ω·cm]
+
+ // Add spike threshold detector at the soma.
+ decor.place(arb::mlocation{0,0}, arb::threshold_detector{10}, "detector");
+
+ // Add a synapse to the mid point of the first dendrite.
+ decor.place(arb::mlocation{0, 0.5}, arb::synapse("expsyn"), "primary_syn");
+
+ // Add additional synapses that will not be connected to anything.
+
+ if (params.synapses > 1) {
+ decor.place(arb::ls::uniform("dend"_lab, 0, params.synapses - 2, gid), arb::synapse("expsyn"), "extra_syns");
+ }
+
+ // Make a CV between every sample in the sample tree.
+ decor.set_default(arb::cv_policy_every_segment());
+
+ arb::cable_cell cell(arb::morphology(tree), labels, decor);
+
+ return cell;
+}
diff --git a/example/plasticity/plasticity.cpp b/example/plasticity/plasticity.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..4b97d458bec2ea8ea60b12a811fcf56a9a47affb
--- /dev/null
+++ b/example/plasticity/plasticity.cpp
@@ -0,0 +1,104 @@
+#include <any>
+#include <iostream>
+#include <iomanip>
+#include <unordered_map>
+
+#include <arborio/label_parse.hpp>
+
+#include <arbor/spike_source_cell.hpp>
+#include <arbor/cable_cell.hpp>
+#include <arbor/simulation.hpp>
+#include <arbor/util/any_cast.hpp>
+#include <arbor/util/any_ptr.hpp>
+
+using namespace arborio::literals;
+
+// Fan out network with n members
+// - one spike source at gid 0
+// - n-1 passive, soma-only cable cells
+// Starts out disconnected, but new connections from the source to the cable cells may be added.
+struct recipe: public arb::recipe {
+ const std::string // Textuals labels for
+ syn = "synapse", // exponential synapse on cable cells
+ det = "detector", // spike detector on cable cells
+ src = "src"; // spike source
+ const double // Parameters
+ r_soma = 12.6157/2.0, // soma radius; A=500 μm²
+ f_spike = 0.25, // spike source interval
+ weight = 5.5, // connection weight
+ delay = 0.05; // delay
+ arb::locset center = "(location 0 0.5)"_ls; // Soma center
+ arb::region all = "(all)"_reg; // Whole cell
+ arb::cell_size_type n_ = 0; // Cell count
+
+ mutable std::unordered_map<arb::cell_gid_type, std::vector<arb::cell_connection>> connected; // lookup table for connections
+ // Required but uninteresting methods
+ recipe(arb::cell_size_type n): n_{n} {}
+ arb::cell_size_type num_cells() const override { return n_; }
+ arb::cell_kind get_cell_kind(arb::cell_gid_type gid) const override { return 0 == gid ? arb::cell_kind::spike_source : arb::cell_kind::cable; }
+ std::any get_global_properties(arb::cell_kind kind) const override {
+ if (kind == arb::cell_kind::spike_source) return {};
+ arb::cable_cell_global_properties ccp;
+ ccp.default_parameters = arb::neuron_parameter_defaults;
+ return ccp;
+ }
+ // For printing Um
+ std::vector<arb::probe_info> get_probes(arb::cell_gid_type gid) const override {
+ if (gid == 0) return {};
+ return {arb::cable_probe_membrane_voltage{center}};
+ }
+ // Look up the (potential) connection to this cell
+ std::vector<arb::cell_connection> connections_on(arb::cell_gid_type gid) const override { return connected[gid]; }
+ // Connect cell `to` to the spike source
+ void add_connection(arb::cell_gid_type to) { assert(to > 0); connected[to] = {arb::cell_connection({0, src}, {syn}, weight, delay)}; }
+ // Return the cell at gid
+ arb::util::unique_any get_cell_description(arb::cell_gid_type gid) const override {
+ // source at gid 0
+ if (gid == 0) return arb::spike_source_cell{src, arb::regular_schedule(f_spike)};
+ // all others are receiving cable cells; single CV w/ HH
+ arb::segment_tree tree; tree.append(arb::mnpos, {-r_soma, 0, 0, r_soma}, {r_soma, 0, 0, r_soma}, 1);
+ auto decor = arb::decor{};
+ decor.paint(all, arb::density("hh", {{"gl", 5}}));
+ decor.place(center, arb::synapse("expsyn"), syn);
+ decor.place(center, arb::threshold_detector{-10.0}, det);
+ decor.set_default(arb::cv_policy_every_segment());
+ return arb::cable_cell({tree}, {}, decor);
+ }
+};
+
+void sampler(arb::probe_metadata pm, std::size_t n, const arb::sample_record* samples) {
+ auto* loc = arb::util::any_cast<const arb::mlocation*>(pm.meta);
+ std::cout << std::fixed << std::setprecision(4);
+ for (std::size_t i = 0; i<n; ++i) {
+ auto* value = arb::util::any_cast<const double*>(samples[i].data);
+ std::cout << "| " << samples[i].time << " | " << loc->pos << " | " << *value << " |\n";
+ }
+}
+
+void spike_cb(const std::vector<arb::spike>& spikes) {
+ for(const auto& spike: spikes) std::cout << " * " << spike.source << "@" << spike.time << '\n';
+}
+
+void print_header(double from, double to) {
+ std::cout << "\n"
+ << "Running simulation from " << from << "ms to " << to << "ms\n"
+ << "Spikes are marked: *\n"
+ << "\n"
+ << "| Time/ms | Position/um | Um/mV |\n"
+ << "|---------+-------------+----------|\n";
+}
+
+const double dt = 0.05;
+
+int main(int argc, char** argv) {
+ auto rec = recipe(3);
+ rec.add_connection(1);
+ auto sim = arb::simulation(rec);
+ sim.add_sampler(arb::all_probes, arb::regular_schedule(dt), sampler, arb::sampling_policy::exact);
+ sim.set_local_spike_callback(spike_cb);
+ print_header(0, 1);
+ sim.run(1.0, dt);
+ rec.add_connection(2);
+ print_header(1, 2);
+ sim.run(2.0, dt);
+}
diff --git a/example/plasticity/readme.md b/example/plasticity/readme.md
new file mode 100644
index 0000000000000000000000000000000000000000..e2cac0b56af041769108f5b0c7a127c945651f04
--- /dev/null
+++ b/example/plasticity/readme.md
@@ -0,0 +1,22 @@
+# Plasticity Example
+
+A miniapp that demonstrates how to use mutable connectivity. The simulation is
+run in two parts with different connectivity.
+
+All cells will print out their _generated_ spikes and (cable cells only!) their
+membrane potential at the soma. The topology is this
+```
+ +--------------+
+ | 0 Spike |
+ | Source |
+ | f=1/0.0125ms |
+ +--------------+
+ / \
+ / \ <--- This connection will not be present for the first part
+ / \ and enabled for the second part of the simulation.
+ / \
+ +---------+ +---------+
+ | 1 Cable | | 2 Cable |
+ | Cell | | Cell |
+ +---------+ +---------+
+```
diff --git a/python/context.cpp b/python/context.cpp
index 1de01e62aa803c731f4c18a644349c7b3db2f69b..2f0e6ffd6f2ce2e0caafcaa3c73ba02e2976b4fc 100644
--- a/python/context.cpp
+++ b/python/context.cpp
@@ -92,7 +92,7 @@ context_shim create_context(unsigned threads, pybind11::object gpu, pybind11::ob
return context_shim(arb::make_context(alloc, c->comm));
}
#endif
- return context_shim(arb::make_context(alloc));
+ return context_shim{arb::make_context(alloc)};
}
std::ostream& operator<<(std::ostream& o, const proc_allocation_shim& alloc) {
diff --git a/python/context.hpp b/python/context.hpp
index 8fff4bc3b2a002148f5cc65de85f00d12158982c..3ba8bfe4e576643c8621cea72cca3f0357684c03 100644
--- a/python/context.hpp
+++ b/python/context.hpp
@@ -6,7 +6,7 @@ namespace pyarb {
struct context_shim {
arb::context context;
- context_shim(arb::context&& c): context(std::move(c)) {}
+ context_shim(arb::context c): context{c} {}
};
} // namespace pyarb
diff --git a/python/example/plasticity.py b/python/example/plasticity.py
new file mode 100644
index 0000000000000000000000000000000000000000..ee20edccfccdca73dd1e8f0f5ec41e1ed9065f2b
--- /dev/null
+++ b/python/example/plasticity.py
@@ -0,0 +1,111 @@
+import arbor as A
+
+
+class recipe(A.recipe):
+ def __init__(self, n):
+ # Call our parent; needed for proper initialization.
+ A.recipe.__init__(self)
+ # Cell count; should be at least 3 for the example to work.
+ assert n >= 3
+ self.cells = n
+ # Uniform weights and delays for the connectivity.
+ self.weight = 0.75
+ self.delay = 0.1
+ # Track the connections from the source cell to a given gid.
+ self.connected = set()
+
+ def num_cells(self):
+ return self.cells
+
+ def cell_kind(self, gid):
+ # Cell 0 is the spike source, all others cable cells.
+ if gid == 0:
+ return A.cell_kind.spike_source
+ else:
+ return A.cell_kind.cable
+
+ def global_properties(self, kind):
+ # Cable cells return the NRN defaults.
+ if kind == A.cell_kind.cable:
+ return A.neuron_cable_properties()
+ # Spike source cells have nothing to report.
+ return None
+
+ def cell_description(self, gid):
+ # Cell 0 is reserved for the spike source, spiking every 0.0125ms.
+ if gid == 0:
+ return A.spike_source_cell("source", A.regular_schedule(0.0125))
+ # All cells >= 1 are cable cells w/ a simple, soma-only morphology
+ # comprising two segments of radius r=3um and length l=3um *each*.
+ #
+ # +-- --+ ^
+ # | * | 3um
+ # +-- --+ v
+ # < 6um >
+ #
+ # * Detectors and Synapses are here, at the midpoint.
+ r = 3
+ tree = A.segment_tree()
+ tree.append(A.mnpos, A.mpoint(-r, 0, 0, r), A.mpoint(r, 0, 0, r), tag=1)
+ # ... and a synapse/detector pair
+ decor = A.decor()
+ # - just have a leaky membrane here.
+ decor.paint("(all)", A.density("pas"))
+ # - synapse to receive incoming spikes from the source cell.
+ decor.place("(location 0 0.5)", A.synapse("expsyn"), "synapse")
+ # - detector for reporting spikes on the cable cells.
+ decor.place("(location 0 0.5)", A.spike_detector(-10.0), "detector")
+ # return the cable cell description
+ return A.cable_cell(tree, A.label_dict(), decor)
+
+ def connections_on(self, gid):
+ # If we have added a connection to this cell, return it, else nothing.
+ if gid in self.connected:
+ # Connect spike source to synapse(s)
+ return [A.connection((0, "source"), "synapse", self.weight, self.delay)]
+ return []
+
+ def add_connection_to_spike_source(self, to):
+ """Add a connection from the spike source at gid=0 to the cable cell
+ gid=<to>. Note that we try to minimize the information stored here,
+ which is important with large networks. Also we cannot connect the
+ source back to itself.
+ """
+ assert to != 0
+ self.connected.add(to)
+
+
+# Make an unconnected network with 2 cable cells and one spike source,
+rec = recipe(3)
+
+# but before setting up anything, connect cable cell gid=1 to spike source gid=0
+# and make the simulation of the simple network
+#
+# spike_source <gid=0> ----> cable_cell <gid=1>
+#
+# cable_cell <gid=2>
+#
+# Note that the connection is just _recorded_ in the recipe, the actual connectivity
+# is set up in the simulation construction.
+rec.add_connection_to_spike_source(1)
+sim = A.simulation(rec)
+sim.record(A.spike_recording.all)
+
+# then run the simulation for a bit
+sim.run(0.25, 0.025)
+
+# update the simulation to
+#
+# spike_source <gid=0> ----> cable_cell <gid=1>
+# \
+# ----> cable_cell <gid=2>
+rec.add_connection_to_spike_source(2)
+sim.update(rec)
+
+# and run the simulation for another bit.
+sim.run(0.5, 0.025)
+
+# When finished, print spike times and locations.
+print("spikes:")
+for sp in sim.spikes():
+ print(" ", sp)
diff --git a/python/simulation.cpp b/python/simulation.cpp
index 647522d0b37987183f98a0f862d79d80bb375ee9..03b80bd055e7cc96dffdc59e003e079aecad94ff 100644
--- a/python/simulation.cpp
+++ b/python/simulation.cpp
@@ -94,6 +94,10 @@ public:
sim_->set_binning_policy(policy, bin_interval);
}
+ void update(std::shared_ptr<py_recipe>& rec) {
+ sim_->update(py_recipe_shim(rec));
+ }
+
void record(spike_recording policy) {
auto spike_recorder = [this](const std::vector<arb::spike>& spikes) {
auto old_size = spike_record_.size();
@@ -215,6 +219,10 @@ void register_simulation(pybind11::module& m, pyarb_global_ptr global_ptr) {
"recipe"_a,
pybind11::arg_v("context", pybind11::none(), "Execution context"),
pybind11::arg_v("domains", pybind11::none(), "Domain decomposition"))
+ .def("update", &simulation_shim::update,
+ "Rebuild the connection table from recipe::connections_on and the event"
+ "generators based on recipe::event_generators.",
+ "recipe"_a)
.def("reset", &simulation_shim::reset,
pybind11::call_guard<pybind11::gil_scoped_release>(),
"Reset the state of the simulation to its initial state.")
@@ -227,7 +235,7 @@ void register_simulation(pybind11::module& m, pyarb_global_ptr global_ptr) {
"tfinal"_a, "dt"_a=0.025)
.def("set_binning_policy", &simulation_shim::set_binning_policy,
"Set the binning policy for event delivery, and the binning time interval if applicable [ms].",
- "policy"_a, "bin_interval"_a)
+ "policy"_a, "bin_interval"_a)
.def("record", &simulation_shim::record,
"Disable or enable local or global spike recording.")
.def("spikes", &simulation_shim::spikes,
diff --git a/scripts/run_cpp_examples.sh b/scripts/run_cpp_examples.sh
index 36c70669204872503052897a89ba508edd1fdb3e..efbcf8157c83146ef4adc4dfb6325009f93eb32b 100755
--- a/scripts/run_cpp_examples.sh
+++ b/scripts/run_cpp_examples.sh
@@ -26,7 +26,7 @@ check () {
fi
}
-for ex in bench brunel gap_junctions generators lfp ring single-cell "probe-demo v"
+for ex in bench brunel gap_junctions generators lfp ring single-cell "probe-demo v" plasticity
do
echo " - $ex"
dir=`echo $ex | tr ' ' '_'`
diff --git a/scripts/run_python_examples.sh b/scripts/run_python_examples.sh
index ff26908335b911feca7fa2150e16d78ab0082467..4277d7b71f8a63a91aa1be7b8d2fa40def3ad6cc 100755
--- a/scripts/run_python_examples.sh
+++ b/scripts/run_python_examples.sh
@@ -30,3 +30,4 @@ $PREFIX python python/example/network_ring.py
$PREFIX python python/example/network_ring_gpu.py # by default, gpu_id=None
$PREFIX python python/example/network_two_cells_gap_junctions.py
$PREFIX python python/example/diffusion.py
+$PREFIX python python/example/plasticity.py
diff --git a/test/unit-distributed/test_communicator.cpp b/test/unit-distributed/test_communicator.cpp
index eeefa66789ae543844d3b9392ec7d1e7f6880e65..e3ac5250e360110d795f3c4f0318f88fefbde5d7 100644
--- a/test/unit-distributed/test_communicator.cpp
+++ b/test/unit-distributed/test_communicator.cpp
@@ -644,10 +644,10 @@ TEST(communicator, all2all)
for (auto i: util::make_span(0, n_global)) {
for (unsigned j = 0; j < n_local; ++j) {
auto c = connections[i*n_local+j];
- EXPECT_EQ(i, c.source().gid);
- EXPECT_EQ(0u, c.source().index);
- EXPECT_EQ(i, c.destination());
- EXPECT_LT(c.index_on_domain(), n_local);
+ EXPECT_EQ(i, c.source.gid);
+ EXPECT_EQ(0u, c.source.index);
+ EXPECT_EQ(i, c.destination);
+ EXPECT_LT(c.index_on_domain, n_local);
}
}
@@ -689,7 +689,7 @@ TEST(communicator, mini_network)
// sort connections by source then target
auto connections = C.connections();
util::sort(connections, [](const connection& lhs, const connection& rhs) {
- return std::forward_as_tuple(lhs.source(), lhs.index_on_domain(), lhs.destination()) < std::forward_as_tuple(rhs.source(), rhs.index_on_domain(), rhs.destination());
+ return std::forward_as_tuple(lhs.source, lhs.index_on_domain, lhs.destination) < std::forward_as_tuple(rhs.source, rhs.index_on_domain, rhs.destination);
});
// Expect one set of 22 connections from every rank: these have been sorted.
@@ -701,9 +701,9 @@ TEST(communicator, mini_network)
std::vector<cell_gid_type> ex_source_gids(22u, i*3 + 1);
for (unsigned j = 0; j < 22u; ++j) {
auto c = connections[i*22 + j];
- EXPECT_EQ(ex_source_gids[j], c.source().gid);
- EXPECT_EQ(ex_source_lids[j], c.source().index);
- EXPECT_EQ(ex_target_lids[i%2][j], c.destination());
+ EXPECT_EQ(ex_source_gids[j], c.source.gid);
+ EXPECT_EQ(ex_source_lids[j], c.source.index);
+ EXPECT_EQ(ex_target_lids[i%2][j], c.destination);
}
}
}
diff --git a/test/unit/test_probe.cpp b/test/unit/test_probe.cpp
index a122552a39a5fca9e9487992081c3d9e71081e5b..61f5fef2b2e1cf84ba09870af7ecd19666f7a787 100644
--- a/test/unit/test_probe.cpp
+++ b/test/unit/test_probe.cpp
@@ -99,7 +99,7 @@ static morphology make_stick_morphology() {
}
template <typename Backend>
-void run_v_i_probe_test(const context& ctx) {
+void run_v_i_probe_test(context ctx) {
using fvm_cell = typename backend_access<Backend>::fvm_cell;
auto deref = [](const arb_value_type* p) { return backend_access<Backend>::deref(p); };
@@ -198,7 +198,7 @@ void run_v_i_probe_test(const context& ctx) {
}
template <typename Backend>
-void run_v_cell_probe_test(const context& ctx) {
+void run_v_cell_probe_test(context ctx) {
using fvm_cell = typename backend_access<Backend>::fvm_cell;
// Take the per-cable voltage over a Y-shaped cell with and without
@@ -258,7 +258,7 @@ void run_v_cell_probe_test(const context& ctx) {
}
template <typename Backend>
-void run_expsyn_g_probe_test(const context& ctx) {
+void run_expsyn_g_probe_test(context ctx) {
using fvm_cell = typename backend_access<Backend>::fvm_cell;
auto deref = [](const arb_value_type* p) { return backend_access<Backend>::deref(p); };
@@ -353,7 +353,7 @@ void run_expsyn_g_probe_test(const context& ctx) {
}
template <typename Backend>
-void run_expsyn_g_cell_probe_test(const context& ctx) {
+void run_expsyn_g_cell_probe_test(context ctx) {
using fvm_cell = typename backend_access<Backend>::fvm_cell;
auto deref = [](const auto* p) { return backend_access<Backend>::deref(p); };
@@ -483,7 +483,7 @@ void run_expsyn_g_cell_probe_test(const context& ctx) {
}
template <typename Backend>
-void run_ion_density_probe_test(const context& ctx) {
+void run_ion_density_probe_test(context ctx) {
using fvm_cell = typename backend_access<Backend>::fvm_cell;
auto deref = [](const arb_value_type* p) { return backend_access<Backend>::deref(p); };
@@ -656,7 +656,7 @@ void run_ion_density_probe_test(const context& ctx) {
}
template <typename Backend>
-void run_partial_density_probe_test(const context& ctx) {
+void run_partial_density_probe_test(context ctx) {
using fvm_cell = typename backend_access<Backend>::fvm_cell;
auto deref = [](const arb_value_type* p) { return backend_access<Backend>::deref(p); };
@@ -769,7 +769,7 @@ void run_partial_density_probe_test(const context& ctx) {
}
template <typename Backend>
-void run_axial_and_ion_current_sampled_probe_test(const context& ctx) {
+void run_axial_and_ion_current_sampled_probe_test(context ctx) {
// On a passive cable in steady-state, the capacitive membrane current will be zero,
// and the axial currents should balance the stimulus and ionic membrane currents in any CV.
//
@@ -937,7 +937,7 @@ auto run_simple_sampler(
}
template <typename Backend>
-void run_multi_probe_test(const context& ctx) {
+void run_multi_probe_test(context ctx) {
// Construct and run thorugh simple sampler a probe defined over
// cell terminal points; check metadata and values.
@@ -971,7 +971,7 @@ void run_multi_probe_test(const context& ctx) {
}
template <typename Backend>
-void run_v_sampled_probe_test(const context& ctx) {
+void run_v_sampled_probe_test(context ctx) {
soma_cell_builder builder(12.6157/2.0);
builder.add_branch(0, 200, 1.0/2, 1.0/2, 1, "dend");
builder.add_branch(0, 200, 1.0/2, 1.0/2, 1, "dend");
@@ -1015,7 +1015,7 @@ void run_v_sampled_probe_test(const context& ctx) {
template <typename Backend>
-void run_total_current_probe_test(const context& ctx) {
+void run_total_current_probe_test(context ctx) {
// Model two passive Y-shaped cells with a similar but not identical
// time constant Ï„.
//
@@ -1147,7 +1147,7 @@ void run_total_current_probe_test(const context& ctx) {
template <typename Backend>
-void run_stimulus_probe_test(const context& ctx) {
+void run_stimulus_probe_test(context ctx) {
// Model two simple stick cable cells, 3 CVs each, and stimuli on cell 0, cv 1
// and cell 1, cv 2. Run both cells in the same cell group.
@@ -1192,7 +1192,7 @@ void run_stimulus_probe_test(const context& ctx) {
}
template <typename Backend>
-void run_exact_sampling_probe_test(const context& ctx) {
+void run_exact_sampling_probe_test(context ctx) {
// As the exact sampling implementation interacts with the event delivery
// implementation within in cable cell groups, construct a somewhat
// elaborate model with 4 cells and a gap junction between cell 1 and 3.