From c00ea33b908e759b5eee5881920b3e54fc7f5834 Mon Sep 17 00:00:00 2001 From: Thorsten Hater <24411438+thorstenhater@users.noreply.github.com> Date: Mon, 17 Oct 2022 17:01:41 +0200 Subject: [PATCH 01/24] Only accept recipes. --- arbor/communication/communicator.cpp | 2 +- arbor/communication/communicator.hpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/arbor/communication/communicator.cpp b/arbor/communication/communicator.cpp index 70fa08ba2e..5be35d8bb1 100644 --- a/arbor/communication/communicator.cpp +++ b/arbor/communication/communicator.cpp @@ -31,7 +31,7 @@ communicator::communicator(const recipe& rec, distributed_{ctx.distributed}, thread_pool_{ctx.thread_pool} {} -void communicator::update_connections(const connectivity& rec, +void communicator::update_connections(const recipe& rec, const domain_decomposition& dom_dec, const label_resolution_map& source_resolution_map, const label_resolution_map& target_resolution_map) { diff --git a/arbor/communication/communicator.hpp b/arbor/communication/communicator.hpp index c9eab08294..bf575e7061 100644 --- a/arbor/communication/communicator.hpp +++ b/arbor/communication/communicator.hpp @@ -68,7 +68,7 @@ class ARB_ARBOR_API communicator { void reset(); - void update_connections(const connectivity& rec, + void update_connections(const recipe& rec, const domain_decomposition& dom_dec, const label_resolution_map& source_resolution_map, const label_resolution_map& target_resolution_map); From cacf3d178e5bc51a214aedecdb8d768af0223174 Mon Sep 17 00:00:00 2001 From: Thorsten Hater <24411438+thorstenhater@users.noreply.github.com> Date: Wed, 19 Oct 2022 10:14:31 +0200 Subject: [PATCH 02/24] Stop communicator update-connections from concatenating a huge vector. --- arbor/CMakeLists.txt | 2 +- arbor/benchmark_cell_group.cpp | 17 ++++--- arbor/benchmark_cell_group.hpp | 2 + ...mc_cell_group.cpp => cable_cell_group.cpp} | 30 ++++++++---- ...mc_cell_group.hpp => cable_cell_group.hpp} | 9 ++-- arbor/cell_group_factory.cpp | 4 +- arbor/communication/communicator.cpp | 46 +++++++++++++++++-- arbor/communication/communicator.hpp | 1 - arbor/fvm_lowered_cell_impl.hpp | 11 ++--- arbor/include/arbor/simulation.hpp | 2 +- arbor/label_resolution.hpp | 2 +- arbor/lif_cell_group.cpp | 21 +++++---- arbor/lif_cell_group.hpp | 2 + arbor/simulation.cpp | 18 ++------ arbor/spike_source_cell_group.cpp | 18 ++++---- arbor/spike_source_cell_group.hpp | 3 ++ test/unit-distributed/test_communicator.cpp | 24 +++++----- test/unit/CMakeLists.txt | 2 +- ...ll_group.cpp => test_cable_cell_group.cpp} | 18 ++++---- 19 files changed, 142 insertions(+), 90 deletions(-) rename arbor/{mc_cell_group.cpp => cable_cell_group.cpp} (95%) rename arbor/{mc_cell_group.hpp => cable_cell_group.hpp} (91%) rename test/unit/{test_mc_cell_group.cpp => test_cable_cell_group.cpp} (86%) diff --git a/arbor/CMakeLists.txt b/arbor/CMakeLists.txt index b58ebff44c..48036c44a1 100644 --- a/arbor/CMakeLists.txt +++ b/arbor/CMakeLists.txt @@ -27,7 +27,7 @@ set(arbor_sources io/serialize_hex.cpp label_resolution.cpp lif_cell_group.cpp - mc_cell_group.cpp + cable_cell_group.cpp mechcat.cpp mechinfo.cpp memory/gpu_wrappers.cpp diff --git a/arbor/benchmark_cell_group.cpp b/arbor/benchmark_cell_group.cpp index 1819edefed..328d2957b9 100644 --- a/arbor/benchmark_cell_group.cpp +++ b/arbor/benchmark_cell_group.cpp @@ -15,6 +15,12 @@ namespace arb { +cell_size_type ARB_ARBOR_API get_sources(cell_label_range& src, const benchmark_cell& c) { + src.add_cell(); + src.add_label(c.source, {0, 1}); + return 1; +} + benchmark_cell_group::benchmark_cell_group(const std::vector& gids, const recipe& rec, cell_label_range& cg_sources, @@ -29,14 +35,11 @@ benchmark_cell_group::benchmark_cell_group(const std::vector& gid cells_.reserve(gids_.size()); for (auto gid: gids_) { - cells_.push_back(util::any_cast(rec.get_cell_description(gid))); - } - - for (const auto& c: cells_) { - cg_sources.add_cell(); + auto cell = util::any_cast(rec.get_cell_description(gid)); + cells_.push_back(cell); cg_targets.add_cell(); - cg_sources.add_label(c.source, {0, 1}); - cg_targets.add_label(c.target, {0, 1}); + cg_targets.add_label(cell.target, {0, 1}); + get_sources(cg_sources, cell); } benchmark_cell_group::reset(); diff --git a/arbor/benchmark_cell_group.hpp b/arbor/benchmark_cell_group.hpp index c106836a16..40bc3f2d4b 100644 --- a/arbor/benchmark_cell_group.hpp +++ b/arbor/benchmark_cell_group.hpp @@ -40,5 +40,7 @@ class benchmark_cell_group: public cell_group { std::vector gids_; }; +cell_size_type ARB_ARBOR_API get_sources(cell_label_range& src, const benchmark_cell& c); + } // namespace arb diff --git a/arbor/mc_cell_group.cpp b/arbor/cable_cell_group.cpp similarity index 95% rename from arbor/mc_cell_group.cpp rename to arbor/cable_cell_group.cpp index f3ac374b68..c8fbc1fbb2 100644 --- a/arbor/mc_cell_group.cpp +++ b/arbor/cable_cell_group.cpp @@ -16,7 +16,7 @@ #include "event_binner.hpp" #include "fvm_lowered_cell.hpp" #include "label_resolution.hpp" -#include "mc_cell_group.hpp" +#include "cable_cell_group.hpp" #include "profile/profiler_macro.hpp" #include "sampler_map.hpp" #include "util/filter.hpp" @@ -30,7 +30,17 @@ 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& gids, +cell_size_type ARB_ARBOR_API get_sources(cell_label_range& src, const cable_cell& c) { + src.add_cell(); + cell_size_type count = 0; + for (const auto& [label, range]: c.detector_ranges()) { + src.add_label(label, range); + count += range.end - range.begin; + } + return count; +} + +cable_cell_group::cable_cell_group(const std::vector& gids, const recipe& rec, cell_label_range& cg_sources, cell_label_range& cg_targets, @@ -38,7 +48,7 @@ mc_cell_group::mc_cell_group(const std::vector& gids, gids_(gids), lowered_(std::move(lowered)) { // Default to no binning of events - mc_cell_group::set_binning_policy(binning_kind::none, 0); + cable_cell_group::set_binning_policy(binning_kind::none, 0); // Build lookup table for gid to local index. for (auto i: util::count_along(gids_)) { @@ -70,7 +80,7 @@ mc_cell_group::mc_cell_group(const std::vector& gids, spike_sources_.shrink_to_fit(); } -void mc_cell_group::reset() { +void cable_cell_group::reset() { spikes_.clear(); sample_events_.clear(); @@ -85,7 +95,7 @@ void mc_cell_group::reset() { lowered_->reset(); } -void mc_cell_group::set_binning_policy(binning_kind policy, time_type bin_interval) { +void cable_cell_group::set_binning_policy(binning_kind policy, time_type bin_interval) { binners_.clear(); binners_.resize(gids_.size(), event_binner(policy, bin_interval)); } @@ -389,7 +399,7 @@ void run_samples( std::visit([&](auto& x) {run_samples(x, sc, raw_times, raw_samples, sample_records, scratch); }, sc.pdata_ptr->info); } -void mc_cell_group::advance(epoch ep, time_type dt, const event_lane_subrange& event_lanes) { +void cable_cell_group::advance(epoch ep, time_type dt, const event_lane_subrange& event_lanes) { time_type tstart = lowered_->time(); // Bin and collate deliverable events from event lanes. @@ -555,7 +565,7 @@ void mc_cell_group::advance(epoch ep, time_type dt, const event_lane_subrange& e } } -void mc_cell_group::add_sampler(sampler_association_handle h, cell_member_predicate probeset_ids, +void cable_cell_group::add_sampler(sampler_association_handle h, cell_member_predicate probeset_ids, schedule sched, sampler_function fn, sampling_policy policy) { std::lock_guard guard(sampler_mex_); @@ -569,17 +579,17 @@ void mc_cell_group::add_sampler(sampler_association_handle h, cell_member_predic } } -void mc_cell_group::remove_sampler(sampler_association_handle h) { +void cable_cell_group::remove_sampler(sampler_association_handle h) { std::lock_guard guard(sampler_mex_); sampler_map_.erase(h); } -void mc_cell_group::remove_all_samplers() { +void cable_cell_group::remove_all_samplers() { std::lock_guard guard(sampler_mex_); sampler_map_.clear(); } -std::vector mc_cell_group::get_probe_metadata(cell_member_type probeset_id) const { +std::vector cable_cell_group::get_probe_metadata(cell_member_type probeset_id) const { // Probe associations are fixed after construction, so we do not need to grab the mutex. std::optional maybe_tag = util::value_by_key(probe_map_.tag, probeset_id); diff --git a/arbor/mc_cell_group.hpp b/arbor/cable_cell_group.hpp similarity index 91% rename from arbor/mc_cell_group.hpp rename to arbor/cable_cell_group.hpp index 45b6c540e8..f20d30aa0d 100644 --- a/arbor/mc_cell_group.hpp +++ b/arbor/cable_cell_group.hpp @@ -12,6 +12,7 @@ #include #include #include +#include #include "backends/event.hpp" #include "cell_group.hpp" @@ -24,11 +25,11 @@ namespace arb { -class ARB_ARBOR_API mc_cell_group: public cell_group { +class ARB_ARBOR_API cable_cell_group: public cell_group { public: - mc_cell_group() = default; + cable_cell_group() = default; - mc_cell_group(const std::vector& gids, + cable_cell_group(const std::vector& gids, const recipe& rec, cell_label_range& cg_sources, cell_label_range& cg_targets, @@ -105,4 +106,6 @@ class ARB_ARBOR_API mc_cell_group: public cell_group { std::vector target_handle_divisions_; }; +cell_size_type ARB_ARBOR_API get_sources(cell_label_range& src, const cable_cell& c); + } // namespace arb diff --git a/arbor/cell_group_factory.cpp b/arbor/cell_group_factory.cpp index e7fb7006a9..d6ced713da 100644 --- a/arbor/cell_group_factory.cpp +++ b/arbor/cell_group_factory.cpp @@ -9,7 +9,7 @@ #include "execution_context.hpp" #include "fvm_lowered_cell.hpp" #include "lif_cell_group.hpp" -#include "mc_cell_group.hpp" +#include "cable_cell_group.hpp" #include "spike_source_cell_group.hpp" namespace arb { @@ -27,7 +27,7 @@ ARB_ARBOR_API cell_group_factory cell_kind_implementation( switch (ck) { case cell_kind::cable: return [bk, ctx, seed](const gid_vector& gids, const recipe& rec, cell_label_range& cg_sources, cell_label_range& cg_targets) { - return make_cell_group(gids, rec, cg_sources, cg_targets, make_fvm_lowered_cell(bk, ctx, seed)); + return make_cell_group(gids, rec, cg_sources, cg_targets, make_fvm_lowered_cell(bk, ctx, seed)); }; case cell_kind::spike_source: diff --git a/arbor/communication/communicator.cpp b/arbor/communication/communicator.cpp index 5be35d8bb1..f7446e9ca1 100644 --- a/arbor/communication/communicator.cpp +++ b/arbor/communication/communicator.cpp @@ -1,5 +1,6 @@ #include #include +#include #include #include @@ -17,6 +18,10 @@ #include "util/partition.hpp" #include "util/rangeutil.hpp" #include "util/span.hpp" +#include "spike_source_cell_group.hpp" +#include "cable_cell_group.hpp" +#include "benchmark_cell_group.hpp" +#include "lif_cell_group.hpp" #include "communication/communicator.hpp" @@ -31,9 +36,35 @@ communicator::communicator(const recipe& rec, distributed_{ctx.distributed}, thread_pool_{ctx.thread_pool} {} +// This is a bit nasty, as we basically reimplement things from all cell kinds ... +cell_label_range get_sources(cell_gid_type gid, const recipe& rec) { + auto cell = rec.get_cell_description(gid); + auto kind = rec.get_cell_kind(gid); + cell_label_range result; + if (kind == cell_kind::lif) { + get_sources(result, util::any_cast(cell)); + } + else if (kind == cell_kind::spike_source) { + get_sources(result, util::any_cast(cell)); + } + else if (kind == cell_kind::benchmark) { + get_sources(result, util::any_cast(cell)); + } + else if (kind == cell_kind::cable) { + auto c = util::any_cast(cell); + result.add_cell(); + for (const auto& [label, range]: c.detector_ranges()) { + result.add_label(label, range); + } + } + else { + throw arbor_internal_error("Unknown cell kind"); + } + return result; +} + void communicator::update_connections(const recipe& 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(); @@ -102,11 +133,20 @@ void communicator::update_connections(const recipe& rec, 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); + auto source_labels = std::unordered_map>{}; for (const auto& cell: gid_infos) { auto index = cell.index_on_domain; - auto source_resolver = resolver(&source_resolution_map); + auto source_resolvers = std::unordered_map{}; for (const auto& c: cell.conns) { - auto src_lid = source_resolver.resolve(c.source); + auto sgid = c.source.gid; + if (!source_labels.count(sgid)) { + auto sources = cell_labels_and_gids{get_sources(sgid, rec), {sgid}}; + source_labels.emplace(sgid, std::make_unique(sources)); + } + if (!source_resolvers.count(sgid)) { + source_resolvers.emplace(sgid, source_labels.at(sgid).get()); + } + auto src_lid = source_resolvers.at(sgid).resolve(c.source); auto tgt_lid = target_resolver.resolve({cell.gid, c.dest}); auto offset = offsets[*src_domain]++; ++src_domain; diff --git a/arbor/communication/communicator.hpp b/arbor/communication/communicator.hpp index bf575e7061..1088bbccf1 100644 --- a/arbor/communication/communicator.hpp +++ b/arbor/communication/communicator.hpp @@ -70,7 +70,6 @@ class ARB_ARBOR_API communicator { void update_connections(const recipe& rec, const domain_decomposition& dom_dec, - const label_resolution_map& source_resolution_map, const label_resolution_map& target_resolution_map); private: diff --git a/arbor/fvm_lowered_cell_impl.hpp b/arbor/fvm_lowered_cell_impl.hpp index 09902f5c27..ea6fd01db0 100644 --- a/arbor/fvm_lowered_cell_impl.hpp +++ b/arbor/fvm_lowered_cell_impl.hpp @@ -35,6 +35,7 @@ #include "util/rangeutil.hpp" #include "util/strprintf.hpp" #include "util/transform.hpp" +#include "cable_cell_group.hpp" namespace arb { @@ -402,18 +403,12 @@ fvm_initialization_data fvm_lowered_cell_impl::initialize( auto gid = gids[i]; const auto& c = cells[i]; - fvm_info.source_data.add_cell(); fvm_info.target_data.add_cell(); fvm_info.gap_junction_data.add_cell(); - unsigned count = 0; - for (const auto& [label, range]: c.detector_ranges()) { - fvm_info.source_data.add_label(label, range); - count+=(range.end - range.begin); - } - fvm_info.num_sources[gid] = count; + fvm_info.num_sources[gid] = get_sources(fvm_info.source_data, c); - count = 0; + unsigned count = 0; for (const auto& [label, range]: c.synapse_ranges()) { fvm_info.target_data.add_label(label, range); count+=(range.end - range.begin); diff --git a/arbor/include/arbor/simulation.hpp b/arbor/include/arbor/simulation.hpp index 4acdadafca..def7e2960d 100644 --- a/arbor/include/arbor/simulation.hpp +++ b/arbor/include/arbor/simulation.hpp @@ -45,7 +45,7 @@ class ARB_ARBOR_API simulation { static simulation_builder create(recipe const &); - void update(const connectivity& rec); + void update(const recipe& rec); void reset(); diff --git a/arbor/label_resolution.hpp b/arbor/label_resolution.hpp index dbfe2014a9..f087f4dd62 100644 --- a/arbor/label_resolution.hpp +++ b/arbor/label_resolution.hpp @@ -114,7 +114,7 @@ struct ARB_ARBOR_API resolver { resolver(const label_resolution_map* label_map): label_map_(label_map) {} cell_lid_type resolve(const cell_global_label_type& iden); -private: +// private: using state_variant = std::variant; state_variant construct_state(lid_selection_policy pol); diff --git a/arbor/lif_cell_group.cpp b/arbor/lif_cell_group.cpp index 0fa71789ec..5c112fe4ec 100644 --- a/arbor/lif_cell_group.cpp +++ b/arbor/lif_cell_group.cpp @@ -6,7 +6,13 @@ #include "util/rangeutil.hpp" #include "util/span.hpp" -using namespace arb; +namespace arb { + +cell_size_type ARB_ARBOR_API get_sources(cell_label_range& src, const lif_cell& c) { + src.add_cell(); + src.add_label(c.source, {0, 1}); + return 1; +} // Constructor containing gid of first cell in a group and a container of all cells. lif_cell_group::lif_cell_group(const std::vector& gids, const recipe& rec, cell_label_range& cg_sources, cell_label_range& cg_targets): @@ -24,14 +30,11 @@ lif_cell_group::lif_cell_group(const std::vector& gids, const rec last_time_updated_.resize(gids_.size()); for (auto lid: util::make_span(gids_.size())) { - cells_.push_back(util::any_cast(rec.get_cell_description(gids_[lid]))); - } - - for (const auto& c: cells_) { - cg_sources.add_cell(); + auto cell = util::any_cast(rec.get_cell_description(gids_[lid])); + cells_.push_back(cell); + get_sources(cg_sources, cell); cg_targets.add_cell(); - cg_sources.add_label(c.source, {0, 1}); - cg_targets.add_label(c.target, {0, 1}); + cg_targets.add_label(cell.target, {0, 1}); } } @@ -120,3 +123,5 @@ void lif_cell_group::advance_cell(time_type tfinal, time_type dt, cell_gid_type // This is the last time a cell was updated. last_time_updated_[lid] = t; } + +} // namespace arb diff --git a/arbor/lif_cell_group.hpp b/arbor/lif_cell_group.hpp index b89e442f86..6d9bed8785 100644 --- a/arbor/lif_cell_group.hpp +++ b/arbor/lif_cell_group.hpp @@ -55,4 +55,6 @@ class ARB_ARBOR_API lif_cell_group: public cell_group { std::vector last_time_updated_; }; +cell_size_type ARB_ARBOR_API get_sources(cell_label_range& src, const lif_cell& c); + } // namespace arb diff --git a/arbor/simulation.cpp b/arbor/simulation.cpp index 887462d12e..eeceded1f9 100644 --- a/arbor/simulation.cpp +++ b/arbor/simulation.cpp @@ -92,7 +92,7 @@ class simulation_state { public: simulation_state(const recipe& rec, const domain_decomposition& decomp, context ctx, arb_seed_type seed); - void update(const connectivity& rec); + void update(const recipe& rec); void reset(); @@ -118,7 +118,6 @@ class simulation_state { 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: @@ -200,7 +199,6 @@ simulation_state::simulation_state( // Generate the cell groups in parallel, with one task per cell group. auto num_groups = decomp.num_groups(); cell_groups_.resize(num_groups); - std::vector cg_sources(num_groups); std::vector cg_targets(num_groups); foreach_group_index( [&](cell_group_ptr& group, int i) { @@ -208,27 +206,21 @@ simulation_state::simulation_state( cell_label_range sources, targets; auto factory = cell_kind_implementation(group_info.kind, group_info.backend, *ctx_, seed); group = factory(group_info.gids, rec, sources, targets); - - cg_sources[i] = cell_labels_and_gids(std::move(sources), group_info.gids); cg_targets[i] = cell_labels_and_gids(std::move(targets), group_info.gids); }); - cell_labels_and_gids local_sources, local_targets; + cell_labels_and_gids local_targets; for(const auto& i: util::make_span(num_groups)) { - 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); - - source_resolution_map_ = label_resolution_map(std::move(global_sources)); target_resolution_map_ = label_resolution_map(std::move(local_targets)); communicator_ = communicator(rec, ddc_, *ctx_); update(rec); epoch_.reset(); } -void simulation_state::update(const connectivity& rec) { - communicator_.update_connections(rec, ddc_, source_resolution_map_, target_resolution_map_); +void simulation_state::update(const recipe& rec) { + communicator_.update_connections(rec, ddc_, target_resolution_map_); // Use half minimum delay of the network for max integration interval. t_interval_ = communicator_.min_delay()/2; @@ -530,7 +522,7 @@ void simulation::reset() { impl_->reset(); } -void simulation::update(const connectivity& rec) { impl_->update(rec); } +void simulation::update(const recipe& rec) { impl_->update(rec); } time_type simulation::run(time_type tfinal, time_type dt) { if (dt <= 0.0) { diff --git a/arbor/spike_source_cell_group.cpp b/arbor/spike_source_cell_group.cpp index 2f4176f3d8..b77013650d 100644 --- a/arbor/spike_source_cell_group.cpp +++ b/arbor/spike_source_cell_group.cpp @@ -13,6 +13,12 @@ namespace arb { +cell_size_type ARB_ARBOR_API get_sources(cell_label_range& src, const spike_source_cell& c) { + src.add_cell(); + src.add_label(c.source, {0, 1}); + return 1; +} + spike_source_cell_group::spike_source_cell_group( const std::vector& gids, const recipe& rec, @@ -28,16 +34,10 @@ spike_source_cell_group::spike_source_cell_group( time_sequences_.reserve(gids.size()); for (auto gid: gids_) { - cg_sources.add_cell(); + auto cell = util::any_cast(rec.get_cell_description(gid)); cg_targets.add_cell(); - try { - auto cell = util::any_cast(rec.get_cell_description(gid)); - time_sequences_.emplace_back(cell.seqs); - cg_sources.add_label(cell.source, {0, 1}); - } - catch (std::bad_any_cast& e) { - throw bad_cell_description(cell_kind::spike_source, gid); - } + get_sources(cg_sources, cell); + time_sequences_.emplace_back(cell.seqs); } } diff --git a/arbor/spike_source_cell_group.hpp b/arbor/spike_source_cell_group.hpp index 23a9c4b36a..5c59a7bd00 100644 --- a/arbor/spike_source_cell_group.hpp +++ b/arbor/spike_source_cell_group.hpp @@ -8,6 +8,7 @@ #include #include #include +#include #include "cell_group.hpp" #include "epoch.hpp" @@ -43,5 +44,7 @@ class ARB_ARBOR_API spike_source_cell_group: public cell_group { std::vector> time_sequences_; }; +cell_size_type ARB_ARBOR_API get_sources(cell_label_range& src, const spike_source_cell& c); + } // namespace arb diff --git a/test/unit-distributed/test_communicator.cpp b/test/unit-distributed/test_communicator.cpp index 89f98dc14c..d29aa42938 100644 --- a/test/unit-distributed/test_communicator.cpp +++ b/test/unit-distributed/test_communicator.cpp @@ -13,7 +13,7 @@ #include "execution_context.hpp" #include "fvm_lowered_cell.hpp" #include "lif_cell_group.hpp" -#include "mc_cell_group.hpp" +#include "cable_cell_group.hpp" #include "util/filter.hpp" #include "util/rangeutil.hpp" #include "util/span.hpp" @@ -356,7 +356,7 @@ namespace { } cell_kind get_cell_kind(cell_gid_type gid) const override { - return gid%3 != 1? cell_kind::cable: cell_kind::lif; + return cell_kind::cable; } std::vector connections_on(cell_gid_type gid) const override { @@ -509,7 +509,7 @@ TEST(communicator, ring) const auto D = partition_load_balance(R, g_context); // set up source and target label->lid resolvers - // from mc_cell_group and lif_cell_group + // from cable_cell_group and lif_cell_group std::vector mc_gids, lif_gids; for (auto g: D.groups()) { if (g.kind == cell_kind::cable) { @@ -520,7 +520,7 @@ TEST(communicator, ring) } } cell_label_range mc_srcs, mc_tgts, lif_srcs, lif_tgts; - auto mc_group = mc_cell_group(mc_gids, R, mc_srcs, mc_tgts, make_fvm_lowered_cell(backend_kind::multicore, *g_context)); + auto mc_group = cable_cell_group(mc_gids, R, mc_srcs, mc_tgts, make_fvm_lowered_cell(backend_kind::multicore, *g_context)); auto lif_group = lif_cell_group(lif_gids, R, lif_srcs, lif_tgts); auto local_sources = cell_labels_and_gids(mc_srcs, mc_gids); @@ -532,7 +532,7 @@ TEST(communicator, ring) // construct the communicator auto C = communicator(R, D, *g_context); - C.update_connections(R, D, label_resolution_map(global_sources), label_resolution_map(local_targets)); + C.update_connections(R, D, label_resolution_map(local_targets)); // every cell fires EXPECT_TRUE(test_ring(D, C, [](cell_gid_type g){return true;})); // last cell in each domain fires @@ -628,18 +628,17 @@ TEST(communicator, all2all) const auto D = partition_load_balance(R, g_context); // set up source and target label->lid resolvers - // from mc_cell_group + // from cable_cell_group std::vector mc_gids; for (auto g: D.groups()) { mc_gids.insert(mc_gids.end(), g.gids.begin(), g.gids.end()); } cell_label_range local_sources, local_targets; - auto mc_group = mc_cell_group(mc_gids, R, local_sources, local_targets, make_fvm_lowered_cell(backend_kind::multicore, *g_context)); - auto global_sources = g_context->distributed->gather_cell_labels_and_gids({local_sources, mc_gids}); + auto mc_group = cable_cell_group(mc_gids, R, local_sources, local_targets, make_fvm_lowered_cell(backend_kind::multicore, *g_context)); // construct the communicator auto C = communicator(R, D, *g_context); - C.update_connections(R, D, label_resolution_map(global_sources), label_resolution_map({local_targets, mc_gids})); + C.update_connections(R, D, label_resolution_map({local_targets, mc_gids})); auto connections = C.connections(); for (auto i: util::make_span(0, n_global)) { @@ -675,18 +674,17 @@ TEST(communicator, mini_network) const auto D = partition_load_balance(R, g_context); // set up source and target label->lid resolvers - // from mc_cell_group + // from cable_cell_group std::vector gids; for (auto g: D.groups()) { gids.insert(gids.end(), g.gids.begin(), g.gids.end()); } cell_label_range local_sources, local_targets; - auto mc_group = mc_cell_group(gids, R, local_sources, local_targets, make_fvm_lowered_cell(backend_kind::multicore, *g_context)); - auto global_sources = g_context->distributed->gather_cell_labels_and_gids({local_sources, gids}); + auto mc_group = cable_cell_group(gids, R, local_sources, local_targets, make_fvm_lowered_cell(backend_kind::multicore, *g_context)); // construct the communicator auto C = communicator(R, D, *g_context); - C.update_connections(R, D, label_resolution_map(global_sources), label_resolution_map({local_targets, gids})); + C.update_connections(R, D, label_resolution_map({local_targets, gids})); // sort connections by source then target auto connections = C.connections(); diff --git a/test/unit/CMakeLists.txt b/test/unit/CMakeLists.txt index 9270ae444b..9adc117c0c 100644 --- a/test/unit/CMakeLists.txt +++ b/test/unit/CMakeLists.txt @@ -94,7 +94,7 @@ set(unit_sources test_math.cpp test_matrix.cpp test_mcable_map.cpp - test_mc_cell_group.cpp + test_cable_cell_group.cpp test_mechanisms.cpp test_mech_temp_diam.cpp test_mechcat.cpp diff --git a/test/unit/test_mc_cell_group.cpp b/test/unit/test_cable_cell_group.cpp similarity index 86% rename from test/unit/test_mc_cell_group.cpp rename to test/unit/test_cable_cell_group.cpp index 3744968227..34245c8a60 100644 --- a/test/unit/test_mc_cell_group.cpp +++ b/test/unit/test_cable_cell_group.cpp @@ -6,7 +6,7 @@ #include "epoch.hpp" #include "fvm_lowered_cell.hpp" -#include "mc_cell_group.hpp" +#include "cable_cell_group.hpp" #include "util/rangeutil.hpp" #include "common.hpp" @@ -36,19 +36,19 @@ namespace { } ACCESS_BIND( - std::vector mc_cell_group::*, + std::vector cable_cell_group::*, private_spike_sources_ptr, - &mc_cell_group::spike_sources_) + &cable_cell_group::spike_sources_) -TEST(mc_cell_group, get_kind) { +TEST(cable_cell_group, get_kind) { cable_cell cell = make_cell(); cell_label_range srcs, tgts; - mc_cell_group group{{0}, cable1d_recipe({cell}), srcs, tgts, lowered_cell()}; + cable_cell_group group{{0}, cable1d_recipe({cell}), srcs, tgts, lowered_cell()}; EXPECT_EQ(cell_kind::cable, group.get_cell_kind()); } -TEST(mc_cell_group, test) { +TEST(cable_cell_group, test) { cable_cell cell = make_cell(); auto rec = cable1d_recipe({cell}); rec.nernst_ion("na"); @@ -56,7 +56,7 @@ TEST(mc_cell_group, test) { rec.nernst_ion("k"); cell_label_range srcs, tgts; - mc_cell_group group{{0}, rec, srcs, tgts, lowered_cell()}; + cable_cell_group group{{0}, rec, srcs, tgts, lowered_cell()}; group.advance(epoch(0, 0., 50.), 0.01, {}); // Model is expected to generate 4 spikes as a result of the @@ -64,7 +64,7 @@ TEST(mc_cell_group, test) { EXPECT_EQ(4u, group.spikes().size()); } -TEST(mc_cell_group, sources) { +TEST(cable_cell_group, sources) { // Make twenty cells, with an extra detector on gids 0, 3 and 17 // to make things more interesting. std::vector cells; @@ -86,7 +86,7 @@ TEST(mc_cell_group, sources) { rec.nernst_ion("k"); cell_label_range srcs, tgts; - mc_cell_group group{gids, rec, srcs, tgts, lowered_cell()}; + cable_cell_group group{gids, rec, srcs, tgts, lowered_cell()}; // Expect group sources to be lexicographically sorted by source id // with gids in cell group's range and indices starting from zero. From adbec3cface8f39dc03164802411eb6ec0b87e1c Mon Sep 17 00:00:00 2001 From: Thorsten Hater <24411438+thorstenhater@users.noreply.github.com> Date: Fri, 27 Jan 2023 09:56:02 +0100 Subject: [PATCH 03/24] Fix review commits. --- arbor/communication/communicator.cpp | 6 +----- arbor/label_resolution.hpp | 1 - 2 files changed, 1 insertion(+), 6 deletions(-) diff --git a/arbor/communication/communicator.cpp b/arbor/communication/communicator.cpp index f7446e9ca1..0903413473 100644 --- a/arbor/communication/communicator.cpp +++ b/arbor/communication/communicator.cpp @@ -51,11 +51,7 @@ cell_label_range get_sources(cell_gid_type gid, const recipe& rec) { get_sources(result, util::any_cast(cell)); } else if (kind == cell_kind::cable) { - auto c = util::any_cast(cell); - result.add_cell(); - for (const auto& [label, range]: c.detector_ranges()) { - result.add_label(label, range); - } + get_sources(result, util::any_cast(cell)); } else { throw arbor_internal_error("Unknown cell kind"); diff --git a/arbor/label_resolution.hpp b/arbor/label_resolution.hpp index f087f4dd62..bb69d6fc17 100644 --- a/arbor/label_resolution.hpp +++ b/arbor/label_resolution.hpp @@ -114,7 +114,6 @@ struct ARB_ARBOR_API resolver { resolver(const label_resolution_map* label_map): label_map_(label_map) {} cell_lid_type resolve(const cell_global_label_type& iden); -// private: using state_variant = std::variant; state_variant construct_state(lid_selection_policy pol); From 6273fc0573bea1247a3a6c16496c1aedd3dc28e3 Mon Sep 17 00:00:00 2001 From: Thorsten Hater <24411438+thorstenhater@users.noreply.github.com> Date: Tue, 31 Jan 2023 15:31:26 +0100 Subject: [PATCH 04/24] Remove an obsolete hierarchy. --- arbor/communication/communicator.cpp | 2 +- arbor/include/arbor/profile/profiler.hpp | 11 ++++++ arbor/include/arbor/recipe.hpp | 43 +++++------------------- arbor/simulation.cpp | 2 ++ 4 files changed, 23 insertions(+), 35 deletions(-) diff --git a/arbor/communication/communicator.cpp b/arbor/communication/communicator.cpp index 0903413473..b483c75026 100644 --- a/arbor/communication/communicator.cpp +++ b/arbor/communication/communicator.cpp @@ -130,9 +130,9 @@ void communicator::update_connections(const recipe& rec, auto src_domain = src_domains.begin(); auto target_resolver = resolver(&target_resolution_map); auto source_labels = std::unordered_map>{}; + auto source_resolvers = std::unordered_map{}; for (const auto& cell: gid_infos) { auto index = cell.index_on_domain; - auto source_resolvers = std::unordered_map{}; for (const auto& c: cell.conns) { auto sgid = c.source.gid; if (!source_labels.count(sgid)) { diff --git a/arbor/include/arbor/profile/profiler.hpp b/arbor/include/arbor/profile/profiler.hpp index 810600ebd0..871a118ed1 100644 --- a/arbor/include/arbor/profile/profiler.hpp +++ b/arbor/include/arbor/profile/profiler.hpp @@ -1,5 +1,8 @@ #pragma once +#include +#include + #include #include #include @@ -13,6 +16,14 @@ namespace arb { namespace profile { +inline +ARB_ARBOR_API void print_hwm_mem_at(const std::string& here, + const std::string& tag="") { + rusage usage; + getrusage(RUSAGE_SELF, &usage); + std::cout << "RSS:" << here << tag << ":" << usage.ru_maxrss << '\n'; +} + // type used for region identifiers using region_id_type = std::size_t; diff --git a/arbor/include/arbor/recipe.hpp b/arbor/include/arbor/recipe.hpp index 1f116b970e..07bff569bd 100644 --- a/arbor/include/arbor/recipe.hpp +++ b/arbor/include/arbor/recipe.hpp @@ -61,41 +61,8 @@ struct gap_junction_connection { peer(std::move(peer)), local(std::move(local)), weight(g) {} }; -struct ARB_ARBOR_API has_gap_junctions { - virtual std::vector gap_junctions_on(cell_gid_type) const { - return {}; - } - virtual ~has_gap_junctions() {} -}; - -struct ARB_ARBOR_API has_synapses { - virtual std::vector connections_on(cell_gid_type) const { - return {}; - } - virtual ~has_synapses() {} -}; - -struct ARB_ARBOR_API has_probes { - virtual std::vector get_probes(cell_gid_type gid) const { - return {}; - } - virtual ~has_probes() {} -}; - -struct ARB_ARBOR_API has_generators { - virtual std::vector event_generators(cell_gid_type) const { - return {}; - } - virtual ~has_generators() {} -}; - -// 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 { +struct ARB_ARBOR_API recipe { // 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. @@ -104,6 +71,14 @@ struct ARB_ARBOR_API recipe: public has_gap_junctions, has_probes, connectivity 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{}; }; + // list of gap junctions incoming to this gid + virtual std::vector gap_junctions_on(cell_gid_type) const { return {}; } + // list of synapses incoming to this gid + virtual std::vector connections_on(cell_gid_type) const { return {}; } + // list of probes on this gid + virtual std::vector get_probes(cell_gid_type gid) const { return {}; } + // list of generators on this gid + virtual std::vector event_generators(cell_gid_type) const { return {}; } virtual ~recipe() {} }; diff --git a/arbor/simulation.cpp b/arbor/simulation.cpp index eeceded1f9..152810364e 100644 --- a/arbor/simulation.cpp +++ b/arbor/simulation.cpp @@ -196,6 +196,7 @@ simulation_state::simulation_state( task_system_(ctx->thread_pool), local_spikes_({thread_private_spike_store(ctx->thread_pool), thread_private_spike_store(ctx->thread_pool)}) { + profile::print_hwm_mem_at(__FUNCTION__, ":Start"); // Generate the cell groups in parallel, with one task per cell group. auto num_groups = decomp.num_groups(); cell_groups_.resize(num_groups); @@ -217,6 +218,7 @@ simulation_state::simulation_state( communicator_ = communicator(rec, ddc_, *ctx_); update(rec); epoch_.reset(); + profile::print_hwm_mem_at(__FUNCTION__, ":Stop"); } void simulation_state::update(const recipe& rec) { From 29f79e366bf22a2e5dfc5ad735ae3fcc4d07fd6b Mon Sep 17 00:00:00 2001 From: Thorsten Hater <24411438+thorstenhater@users.noreply.github.com> Date: Wed, 1 Feb 2023 13:42:58 +0100 Subject: [PATCH 05/24] Add high-water-mark memory to profiling. --- arbor/include/arbor/profile/profiler.hpp | 13 ++------ arbor/profile/profiler.cpp | 38 +++++++++++++++++++----- arbor/simulation.cpp | 2 -- example/ring/ring.cpp | 2 +- 4 files changed, 35 insertions(+), 20 deletions(-) diff --git a/arbor/include/arbor/profile/profiler.hpp b/arbor/include/arbor/profile/profiler.hpp index 871a118ed1..a61cc52f7c 100644 --- a/arbor/include/arbor/profile/profiler.hpp +++ b/arbor/include/arbor/profile/profiler.hpp @@ -1,8 +1,5 @@ #pragma once -#include -#include - #include #include #include @@ -16,13 +13,6 @@ namespace arb { namespace profile { -inline -ARB_ARBOR_API void print_hwm_mem_at(const std::string& here, - const std::string& tag="") { - rusage usage; - getrusage(RUSAGE_SELF, &usage); - std::cout << "RSS:" << here << tag << ":" << usage.ru_maxrss << '\n'; -} // type used for region identifiers using region_id_type = std::size_t; @@ -38,6 +28,9 @@ struct profile { // the accumulated time spent in each region. std::vector times; + // Maximum Resident Set Size up to this region's end. + std::vector max_rss; + // the number of threads for which profiling information was recorded. std::size_t num_threads; diff --git a/arbor/profile/profiler.cpp b/arbor/profile/profiler.cpp index fccb33ba87..578cf80f58 100644 --- a/arbor/profile/profiler.cpp +++ b/arbor/profile/profiler.cpp @@ -10,6 +10,8 @@ #include "util/span.hpp" #include "util/rangeutil.hpp" +#include + namespace arb { namespace profile { @@ -46,6 +48,7 @@ namespace { // Holds the accumulated number of calls and time spent in a region. struct profile_accumulator { std::size_t count=0; + long max_rss = 0; double time=0.; }; @@ -127,11 +130,12 @@ struct profile_node { std::string name; double time = 0; region_id_type count = npos; + long max_rss = -1; std::vector children; profile_node() = default; - profile_node(std::string n, double t, region_id_type c): - name(std::move(n)), time(t), count(c) {} + profile_node(std::string n, double t, region_id_type c, long rss): + name(std::move(n)), time(t), count(c), max_rss{rss} {} profile_node(std::string n): name(std::move(n)), time(0), count(npos) {} }; @@ -154,14 +158,26 @@ void recorder::enter(region_id_type index) { } void recorder::leave() { + long max_rss = 0; + + rusage usage; + getrusage(RUSAGE_SELF, &usage); + max_rss = usage.ru_maxrss; +#ifdef __APPLE__ + max_rss /= 1024; // Because MacOS >= 10.11 reports in B, not kB +#endif + + // calculate the elapsed time before any other steps, to increase accuracy. auto delta = timer_type::toc(start_time_); if (index_==npos) { throw std::runtime_error("recorder::leave without matching recorder::enter"); } - accumulators_[index_].count++; - accumulators_[index_].time += delta; + auto& acc = accumulators_[index_]; + acc.count++; + acc.time += delta; + acc.max_rss = std::max(acc.max_rss, max_rss); index_ = npos; } @@ -235,13 +251,14 @@ profile profiler::results() const { profile p; p.names = region_names_; - + p.max_rss = std::vector(nregions); p.times = std::vector(nregions); p.counts = std::vector(nregions); for (auto& r: recorders_) { auto& accumulators = r.accumulators(); for (auto i: make_span(0, accumulators.size())) { p.times[i] += accumulators[i].time; + p.max_rss[i] = std::max(p.max_rss[i], accumulators[i].max_rss); p.counts[i] += accumulators[i].count; } } @@ -257,9 +274,11 @@ profile profiler::results() const { std::swap(p.counts[i], p.counts.back()); std::swap(p.times[i], p.times.back()); std::swap(p.names[i], p.names.back()); + std::swap(p.max_rss[i], p.max_rss.back()); p.counts.pop_back(); p.times.pop_back(); p.names.pop_back(); + p.max_rss.pop_back(); } return p; @@ -295,7 +314,7 @@ profile_node make_profile_tree(const profile& p) { node = &(*child); } } - node->children.emplace_back(names[idx].back(), p.times[idx], p.counts[idx]); + node->children.emplace_back(names[idx].back(), p.times[idx], p.counts[idx], p.max_rss[idx]); } sort_profile_tree(tree); @@ -312,6 +331,7 @@ struct prof_line { std::string time; std::string thread; std::string percent; + std::string max_rss; }; void print_lines(std::vector& lines, @@ -338,6 +358,7 @@ void print_lines(std::vector& lines, res.thread = buf; snprintf(buf, std::size(buf), "%.1f", float(proportion)); res.percent = buf; + res.max_rss = (n.max_rss == -1) ? "-" : std::to_string(n.max_rss); lines.push_back(res); // print each of the children in turn for (auto& c: n.children) print_lines(lines, c, wall_time, nthreads, thresh, indent + " "); @@ -348,7 +369,7 @@ void print(std::ostream& os, float wall_time, unsigned nthreads, float thresh) { - std::vector lines{{"REGION", "CALLS", "THREAD", "WALL", "\%"}}; + std::vector lines{{"REGION", "CALLS", "THREAD/s", "WALL/s", "\%", "MAX_RSS/kB"}}; print_lines(lines, n, wall_time, nthreads, thresh, ""); // fixing up lengths here std::size_t max_len_name = 0; @@ -356,12 +377,14 @@ void print(std::ostream& os, std::size_t max_len_thread = 0; std::size_t max_len_time = 0; std::size_t max_len_percent = 0; + std::size_t max_len_rss = 0; for (const auto& line: lines) { max_len_name = std::max(max_len_name, line.name.size()); max_len_count = std::max(max_len_count, line.count.size()); max_len_thread = std::max(max_len_thread, line.thread.size()); max_len_time = std::max(max_len_time, line.time.size()); max_len_percent = std::max(max_len_percent, line.percent.size()); + max_len_rss = std::max(max_len_rss, line.max_rss.size()); } auto lpad = [](const std::string& s, std::size_t n) { return std::string(n - s.size(), ' ') + s + " "; }; @@ -372,6 +395,7 @@ void print(std::ostream& os, << lpad(line.thread, max_len_thread) << lpad(line.time, max_len_time) << lpad(line.percent, max_len_percent) + << lpad(line.max_rss, max_len_rss) << '\n'; }; diff --git a/arbor/simulation.cpp b/arbor/simulation.cpp index 152810364e..eeceded1f9 100644 --- a/arbor/simulation.cpp +++ b/arbor/simulation.cpp @@ -196,7 +196,6 @@ simulation_state::simulation_state( task_system_(ctx->thread_pool), local_spikes_({thread_private_spike_store(ctx->thread_pool), thread_private_spike_store(ctx->thread_pool)}) { - profile::print_hwm_mem_at(__FUNCTION__, ":Start"); // Generate the cell groups in parallel, with one task per cell group. auto num_groups = decomp.num_groups(); cell_groups_.resize(num_groups); @@ -218,7 +217,6 @@ simulation_state::simulation_state( communicator_ = communicator(rec, ddc_, *ctx_); update(rec); epoch_.reset(); - profile::print_hwm_mem_at(__FUNCTION__, ":Stop"); } void simulation_state::update(const recipe& rec) { diff --git a/example/ring/ring.cpp b/example/ring/ring.cpp index cf8967537a..4cf951c442 100644 --- a/example/ring/ring.cpp +++ b/example/ring/ring.cpp @@ -47,7 +47,7 @@ struct ring_params { std::string name = "default"; unsigned num_cells = 100; double min_delay = 10; - double duration = 1000; + double duration = 10; cell_parameters cell; }; From a90fc042b2d2143a3d1a9b8656792a42984bce29 Mon Sep 17 00:00:00 2001 From: Thorsten Hater <24411438+thorstenhater@users.noreply.github.com> Date: Wed, 1 Feb 2023 14:27:17 +0100 Subject: [PATCH 06/24] Add first-seen memory mark and iron out printing. --- arbor/include/arbor/profile/profiler.hpp | 1 + arbor/profile/profiler.cpp | 65 ++++++++++++++++-------- example/ring/ring.cpp | 12 ++--- 3 files changed, 50 insertions(+), 28 deletions(-) diff --git a/arbor/include/arbor/profile/profiler.hpp b/arbor/include/arbor/profile/profiler.hpp index a61cc52f7c..e575c0bfde 100644 --- a/arbor/include/arbor/profile/profiler.hpp +++ b/arbor/include/arbor/profile/profiler.hpp @@ -30,6 +30,7 @@ struct profile { // Maximum Resident Set Size up to this region's end. std::vector max_rss; + std::vector fst_rss; // the number of threads for which profiling information was recorded. std::size_t num_threads; diff --git a/arbor/profile/profiler.cpp b/arbor/profile/profiler.cpp index 578cf80f58..900f6bee10 100644 --- a/arbor/profile/profiler.cpp +++ b/arbor/profile/profiler.cpp @@ -48,7 +48,8 @@ namespace { // Holds the accumulated number of calls and time spent in a region. struct profile_accumulator { std::size_t count=0; - long max_rss = 0; + long max_rss = -1; + long fst_rss = -1; double time=0.; }; @@ -131,11 +132,12 @@ struct profile_node { double time = 0; region_id_type count = npos; long max_rss = -1; + long fst_rss = -1; std::vector children; profile_node() = default; - profile_node(std::string n, double t, region_id_type c, long rss): - name(std::move(n)), time(t), count(c), max_rss{rss} {} + profile_node(std::string n, double t, region_id_type c, long mrss, long frss): + name(std::move(n)), time(t), count(c), max_rss{mrss}, fst_rss{frss} {} profile_node(std::string n): name(std::move(n)), time(0), count(npos) {} }; @@ -178,6 +180,7 @@ void recorder::leave() { acc.count++; acc.time += delta; acc.max_rss = std::max(acc.max_rss, max_rss); + if (acc.fst_rss < 0 && acc.max_rss > 0) acc.fst_rss = acc.max_rss; index_ = npos; } @@ -252,6 +255,7 @@ profile profiler::results() const { profile p; p.names = region_names_; p.max_rss = std::vector(nregions); + p.fst_rss = std::vector(nregions); p.times = std::vector(nregions); p.counts = std::vector(nregions); for (auto& r: recorders_) { @@ -259,6 +263,7 @@ profile profiler::results() const { for (auto i: make_span(0, accumulators.size())) { p.times[i] += accumulators[i].time; p.max_rss[i] = std::max(p.max_rss[i], accumulators[i].max_rss); + p.fst_rss[i] = std::max(p.fst_rss[i], accumulators[i].fst_rss); p.counts[i] += accumulators[i].count; } } @@ -275,10 +280,12 @@ profile profiler::results() const { std::swap(p.times[i], p.times.back()); std::swap(p.names[i], p.names.back()); std::swap(p.max_rss[i], p.max_rss.back()); + std::swap(p.fst_rss[i], p.fst_rss.back()); p.counts.pop_back(); p.times.pop_back(); p.names.pop_back(); p.max_rss.pop_back(); + p.fst_rss.pop_back(); } return p; @@ -314,7 +321,7 @@ profile_node make_profile_tree(const profile& p) { node = &(*child); } } - node->children.emplace_back(names[idx].back(), p.times[idx], p.counts[idx], p.max_rss[idx]); + node->children.emplace_back(names[idx].back(), p.times[idx], p.counts[idx], p.max_rss[idx], p.fst_rss[idx]); } sort_profile_tree(tree); @@ -332,6 +339,7 @@ struct prof_line { std::string thread; std::string percent; std::string max_rss; + std::string fst_rss; }; void print_lines(std::vector& lines, @@ -355,10 +363,11 @@ void print_lines(std::vector& lines, snprintf(buf, std::size(buf), "%.3f", float(n.time)); res.time = buf; snprintf(buf, std::size(buf), "%.3f", float(per_thread_time)); - res.thread = buf; + res.thread = std::string{buf}; snprintf(buf, std::size(buf), "%.1f", float(proportion)); - res.percent = buf; + res.percent = std::string{buf}; res.max_rss = (n.max_rss == -1) ? "-" : std::to_string(n.max_rss); + res.fst_rss = (n.fst_rss == -1) ? "-" : std::to_string(n.fst_rss); lines.push_back(res); // print each of the children in turn for (auto& c: n.children) print_lines(lines, c, wall_time, nthreads, thresh, indent + " "); @@ -369,7 +378,8 @@ void print(std::ostream& os, float wall_time, unsigned nthreads, float thresh) { - std::vector lines{{"REGION", "CALLS", "THREAD/s", "WALL/s", "\%", "MAX_RSS/kB"}}; + std::vector lines{{"REGION", "CALLS", "THREAD/s", "WALL/s", "\%", "MAX_RSS/kB", "1st_RSS/kB"}, + {"------", "-----", "--------", "------", "-----", "----------", "----------"}}; print_lines(lines, n, wall_time, nthreads, thresh, ""); // fixing up lengths here std::size_t max_len_name = 0; @@ -377,25 +387,36 @@ void print(std::ostream& os, std::size_t max_len_thread = 0; std::size_t max_len_time = 0; std::size_t max_len_percent = 0; - std::size_t max_len_rss = 0; + std::size_t max_len_mrss = 0; + std::size_t max_len_frss = 0; for (const auto& line: lines) { - max_len_name = std::max(max_len_name, line.name.size()); - max_len_count = std::max(max_len_count, line.count.size()); - max_len_thread = std::max(max_len_thread, line.thread.size()); - max_len_time = std::max(max_len_time, line.time.size()); + max_len_name = std::max(max_len_name, line.name.size()); + max_len_count = std::max(max_len_count, line.count.size()); + max_len_time = std::max(max_len_time, line.time.size()); + max_len_thread = std::max(max_len_thread, line.thread.size()); max_len_percent = std::max(max_len_percent, line.percent.size()); - max_len_rss = std::max(max_len_rss, line.max_rss.size()); + max_len_mrss = std::max(max_len_mrss, line.max_rss.size()); + max_len_frss = std::max(max_len_frss, line.fst_rss.size()); } - auto lpad = [](const std::string& s, std::size_t n) { return std::string(n - s.size(), ' ') + s + " "; }; - auto rpad = [](const std::string& s, std::size_t n) { return s + std::string(n - s.size(), ' ') + " "; }; - - for (const auto& line: lines) os << rpad(line.name, max_len_name) - << lpad(line.count, max_len_count) - << lpad(line.thread, max_len_thread) - << lpad(line.time, max_len_time) - << lpad(line.percent, max_len_percent) - << lpad(line.max_rss, max_len_rss) + auto lpad = [](const std::string& s, std::size_t n) { + auto pad = (n >= s.size()) ? n - s.size() : 0ul; + return std::string(pad, ' ') + s; + }; + + auto rpad = [](const std::string& s, std::size_t n) { + auto pad = (n >= s.size()) ? n - s.size() : 0ul; + return s + std::string(pad, ' '); + }; + + auto sep = " "; + for (const auto& line: lines) os << rpad(line.name, max_len_name) << sep + << lpad(line.count, max_len_count) << sep + << lpad(line.thread, max_len_thread) << sep + << lpad(line.time, max_len_time) << sep + << lpad(line.percent, max_len_percent) << sep + << lpad(line.max_rss, max_len_mrss) << sep + << lpad(line.fst_rss, max_len_frss) << sep << '\n'; }; diff --git a/example/ring/ring.cpp b/example/ring/ring.cpp index 4cf951c442..db7131e058 100644 --- a/example/ring/ring.cpp +++ b/example/ring/ring.cpp @@ -125,7 +125,7 @@ class ring_recipe: public arb::recipe { }; int main(int argc, char** argv) { - try { + // try { bool root = true; arb::proc_allocation resources; @@ -227,11 +227,11 @@ int main(int argc, char** argv) { auto report = arb::profile::make_meter_report(meters, context); std::cout << report; - } - catch (std::exception& e) { - std::cerr << "exception caught in ring miniapp: " << e.what() << "\n"; - return 1; - } + // } + // catch (std::exception& e) { + // std::cerr << "exception caught in ring miniapp: " << e.what() << "\n"; + // return 1; + // } return 0; } From bf1e65c1bc1502054f5a0762a3d601464266bac3 Mon Sep 17 00:00:00 2001 From: Thorsten Hater <24411438+thorstenhater@users.noreply.github.com> Date: Wed, 1 Feb 2023 14:29:21 +0100 Subject: [PATCH 07/24] Appease black. --- doc/scripts/make_images.py | 2 -- example/lfp/neuron_lfp_example.py | 1 - python/example/brunel.py | 2 -- python/example/calcium_stdp.py | 1 - python/example/gap_junctions.py | 1 - python/example/network_two_cells_gap_junctions.py | 3 +-- python/example/ou_lif/ou_lif.py | 3 --- python/example/probe_lfpykit.py | 13 ++++--------- python/example/single_cell_bluepyopt_l5pc.py | 1 - python/example/single_cell_cable.py | 1 - python/example/single_cell_detailed_recipe.py | 1 - python/test/unit/test_multiple_connections.py | 1 - .../unit_distributed/test_domain_decompositions.py | 1 - 13 files changed, 5 insertions(+), 26 deletions(-) diff --git a/doc/scripts/make_images.py b/doc/scripts/make_images.py index 2e6ee0c9e6..af9f4cfef6 100644 --- a/doc/scripts/make_images.py +++ b/doc/scripts/make_images.py @@ -170,7 +170,6 @@ def morph_image(morphs, methods, filename, **kwargs): # Handling this case would make rendering regions more complex, but would # not bee too hard to support. def label_image(morphology, labels, filename, **kwargs): - loc_sc = kwargs.get("loc_sc", 2) sc = kwargs.get("sc", 20) drawroot = kwargs.get("drawroot", True) @@ -294,7 +293,6 @@ def label_image(morphology, labels, filename, **kwargs): def generate(path=""): - morph_image( [inputs.branch_morph2], ["segments"], diff --git a/example/lfp/neuron_lfp_example.py b/example/lfp/neuron_lfp_example.py index 4b8b931566..66d0339b15 100755 --- a/example/lfp/neuron_lfp_example.py +++ b/example/lfp/neuron_lfp_example.py @@ -195,7 +195,6 @@ def collect_geometry(self): self.diam = diamvec def simulate(self): - neuron.h.finitialize(self.v_init) neuron.h.fcurrent() diff --git a/python/example/brunel.py b/python/example/brunel.py index ba8258e2c8..2af659ebb0 100755 --- a/python/example/brunel.py +++ b/python/example/brunel.py @@ -45,7 +45,6 @@ def __init__( poiss_lambda, seed=42, ): - arbor.recipe.__init__(self) # Make sure that in_degree_prop in the interval (0, 1] @@ -115,7 +114,6 @@ def event_generators(self, gid): if __name__ == "__main__": - parser = argparse.ArgumentParser(description="Brunel model miniapp.") parser.add_argument( "-n", diff --git a/python/example/calcium_stdp.py b/python/example/calcium_stdp.py index 6f28718416..890871f36c 100644 --- a/python/example/calcium_stdp.py +++ b/python/example/calcium_stdp.py @@ -107,7 +107,6 @@ def event_generators(self, gid): def run(time_lag): - # Time between stimuli T = 1000.0 / f diff --git a/python/example/gap_junctions.py b/python/example/gap_junctions.py index 50da0589cc..425295813e 100644 --- a/python/example/gap_junctions.py +++ b/python/example/gap_junctions.py @@ -19,7 +19,6 @@ def make_cable_cell(gid): - # Build a segment tree tree = arbor.segment_tree() diff --git a/python/example/network_two_cells_gap_junctions.py b/python/example/network_two_cells_gap_junctions.py index 181833f963..06221bb801 100755 --- a/python/example/network_two_cells_gap_junctions.py +++ b/python/example/network_two_cells_gap_junctions.py @@ -113,7 +113,6 @@ def gap_junctions_on(self, gid): if __name__ == "__main__": - parser = argparse.ArgumentParser( description="Two cells connected via a gap junction" ) @@ -197,7 +196,7 @@ def gap_junctions_on(self, gid): si_gj_R = 1 / si_gj_g # indicate the expected equilibrium potentials - for (i, j) in [[0, 1], [1, 0]]: + for i, j in [[0, 1], [1, 0]]: weighted_potential = args.Vms[i] + ( (args.Vms[j] - args.Vms[i]) * (si_gj_R + cell_R) ) / (2 * cell_R + si_gj_R) diff --git a/python/example/ou_lif/ou_lif.py b/python/example/ou_lif/ou_lif.py index c308ec18f1..dde3240480 100644 --- a/python/example/ou_lif/ou_lif.py +++ b/python/example/ou_lif/ou_lif.py @@ -20,7 +20,6 @@ def make_catalogue(): def make_cell(): - # cell morphology # =============== @@ -176,7 +175,6 @@ def event_generators(self, gid): # Here, if the value of the 'weight' parameter is 1, stimulation is switched on, # whereas if it is -1, stimulation is switched off. def get_generators(self, protocol): - prot_name = protocol["scheme"] # name of the protocol (defining its structure) start_time = protocol["time_start"] # time at which the stimulus starts in s label = protocol["label"] # target synapse (mechanism label) @@ -255,7 +253,6 @@ def get_generators(self, protocol): if __name__ == "__main__": - # set up and run simulation # ========================= diff --git a/python/example/probe_lfpykit.py b/python/example/probe_lfpykit.py index 0db4733fcd..f9dff50840 100644 --- a/python/example/probe_lfpykit.py +++ b/python/example/probe_lfpykit.py @@ -162,6 +162,7 @@ def make_cable_cell(morphology, clamp_location): # Compute extracellular potentials ############################################################################### + # ## Compute extracellular potentials # First we define a couple of classes to interface the LFPykit # library (https://LFPykit.readthedocs.io, https://github.com/LFPy/LFPykit): @@ -334,15 +335,9 @@ def get_cv_polycollection(cell_geometry, V_m, vlims=[-66, -64], cmap="viridis"): inds = cell_geometry._CV_ind == i zips.append( create_polygon( - cell_geometry.x[ - inds, - ].flatten(), - cell_geometry.y[ - inds, - ].flatten(), - cell_geometry.d[ - inds, - ].flatten(), + cell_geometry.x[inds,].flatten(), + cell_geometry.y[inds,].flatten(), + cell_geometry.d[inds,].flatten(), ) ) polycol = PolyCollection(zips, edgecolors=colors, facecolors=colors, linewidths=0.0) diff --git a/python/example/single_cell_bluepyopt_l5pc.py b/python/example/single_cell_bluepyopt_l5pc.py index 2a8b6f397f..502a608464 100755 --- a/python/example/single_cell_bluepyopt_l5pc.py +++ b/python/example/single_cell_bluepyopt_l5pc.py @@ -50,7 +50,6 @@ # (6) Create a class that inherits from arbor.recipe class single_recipe(arbor.recipe): - # (6.1) Define the class constructor def __init__(self, cell, probes): # The base C++ class constructor must be called first, to ensure that diff --git a/python/example/single_cell_cable.py b/python/example/single_cell_cable.py index 348fad9101..2a8b61c68c 100755 --- a/python/example/single_cell_cable.py +++ b/python/example/single_cell_cable.py @@ -149,7 +149,6 @@ def get_tmax(data): if __name__ == "__main__": - parser = argparse.ArgumentParser(description="Cable") parser.add_argument( diff --git a/python/example/single_cell_detailed_recipe.py b/python/example/single_cell_detailed_recipe.py index cefb4dc977..436f80f641 100644 --- a/python/example/single_cell_detailed_recipe.py +++ b/python/example/single_cell_detailed_recipe.py @@ -73,7 +73,6 @@ # (5) Create a class that inherits from arbor.recipe class single_recipe(arbor.recipe): - # (5.1) Define the class constructor def __init__(self): # The base C++ class constructor must be called first, to ensure that diff --git a/python/test/unit/test_multiple_connections.py b/python/test/unit/test_multiple_connections.py index b01c2fd3aa..965f96c35b 100644 --- a/python/test/unit/test_multiple_connections.py +++ b/python/test/unit/test_multiple_connections.py @@ -28,7 +28,6 @@ class TestMultipleConnections(unittest.TestCase): - # Constructor (overridden) def __init__(self, args): super(TestMultipleConnections, self).__init__(args) diff --git a/python/test/unit_distributed/test_domain_decompositions.py b/python/test/unit_distributed/test_domain_decompositions.py index 5535f02fc9..5f4db96739 100644 --- a/python/test/unit_distributed/test_domain_decompositions.py +++ b/python/test/unit_distributed/test_domain_decompositions.py @@ -209,7 +209,6 @@ def test_domain_decomposition_homogenous_MC(self): # 1 node with 1 cpu core, 1 gpu; assumes all cells will be placed on gpu in a single cell group @unittest.skipIf(not gpu_enabled, "GPU not enabled") def test_domain_decomposition_homogenous_GPU(self): - if mpi_enabled: comm = arb.mpi_comm() context = arb.context(threads=1, gpu_id=0, mpi=comm) From 285a734beb60a39a438ce34f1ed23c28a7d19678 Mon Sep 17 00:00:00 2001 From: Thorsten Hater <24411438+thorstenhater@users.noreply.github.com> Date: Thu, 2 Feb 2023 12:23:02 +0100 Subject: [PATCH 08/24] Fix and refactor source labels/resolver handling. --- arbor/communication/communicator.cpp | 24 ++++++++++++++---------- 1 file changed, 14 insertions(+), 10 deletions(-) diff --git a/arbor/communication/communicator.cpp b/arbor/communication/communicator.cpp index b483c75026..d4b0b05e19 100644 --- a/arbor/communication/communicator.cpp +++ b/arbor/communication/communicator.cpp @@ -121,6 +121,17 @@ void communicator::update_connections(const recipe& rec, } } + struct label_map { + std::unique_ptr map; + resolver res; + + label_map(cell_gid_type gid, const recipe& rec): + map{std::make_unique(cell_labels_and_gids{get_sources(gid, rec), {gid}})}, + res{map.get()} + {} + cell_lid_type resolve(const cell_global_label_type& lbl) { return res.resolve(lbl); } + }; + // Construct the connections. // The loop above gave the information required to construct in place // the connections as partitioned by the domain of their source gid. @@ -129,20 +140,13 @@ void communicator::update_connections(const recipe& rec, 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); - auto source_labels = std::unordered_map>{}; - auto source_resolvers = std::unordered_map{}; for (const auto& cell: gid_infos) { auto index = cell.index_on_domain; + auto sources = std::unordered_map{}; for (const auto& c: cell.conns) { auto sgid = c.source.gid; - if (!source_labels.count(sgid)) { - auto sources = cell_labels_and_gids{get_sources(sgid, rec), {sgid}}; - source_labels.emplace(sgid, std::make_unique(sources)); - } - if (!source_resolvers.count(sgid)) { - source_resolvers.emplace(sgid, source_labels.at(sgid).get()); - } - auto src_lid = source_resolvers.at(sgid).resolve(c.source); + if (!sources.count(sgid)) sources.emplace(sgid, label_map{sgid, rec}); + auto src_lid = sources.at(sgid).resolve(c.source); auto tgt_lid = target_resolver.resolve({cell.gid, c.dest}); auto offset = offsets[*src_domain]++; ++src_domain; From 41b6499b9832ba7efc89616a85d962395192d88e Mon Sep 17 00:00:00 2001 From: Thorsten Hater <24411438+thorstenhater@users.noreply.github.com> Date: Mon, 13 Feb 2023 09:19:49 +0100 Subject: [PATCH 09/24] Revert ring.cpp. --- example/ring/ring.cpp | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/example/ring/ring.cpp b/example/ring/ring.cpp index db7131e058..cf8967537a 100644 --- a/example/ring/ring.cpp +++ b/example/ring/ring.cpp @@ -47,7 +47,7 @@ struct ring_params { std::string name = "default"; unsigned num_cells = 100; double min_delay = 10; - double duration = 10; + double duration = 1000; cell_parameters cell; }; @@ -125,7 +125,7 @@ class ring_recipe: public arb::recipe { }; int main(int argc, char** argv) { - // try { + try { bool root = true; arb::proc_allocation resources; @@ -227,11 +227,11 @@ int main(int argc, char** argv) { auto report = arb::profile::make_meter_report(meters, context); std::cout << report; - // } - // catch (std::exception& e) { - // std::cerr << "exception caught in ring miniapp: " << e.what() << "\n"; - // return 1; - // } + } + catch (std::exception& e) { + std::cerr << "exception caught in ring miniapp: " << e.what() << "\n"; + return 1; + } return 0; } From 381691079132f296b7be2957b5d0ed2add89f48b Mon Sep 17 00:00:00 2001 From: Thorsten Hater <24411438+thorstenhater@users.noreply.github.com> Date: Tue, 13 Jun 2023 09:43:47 +0200 Subject: [PATCH 10/24] More mergen. --- arbor/communication/communicator.cpp | 16 ++++++++------- arbor/include/arbor/recipe.hpp | 29 +--------------------------- 2 files changed, 10 insertions(+), 35 deletions(-) diff --git a/arbor/communication/communicator.cpp b/arbor/communication/communicator.cpp index 8bc23f0aa3..3144fa9d29 100644 --- a/arbor/communication/communicator.cpp +++ b/arbor/communication/communicator.cpp @@ -174,14 +174,16 @@ void communicator::update_connections(const recipe& rec, 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; + for (const auto index: util::make_span(num_local_cells_)) { + const auto tgt_gid = gids[index]; auto sources = std::unordered_map{}; - for (const auto& c: cell.conns) { - auto sgid = c.source.gid; - if (!sources.count(sgid)) sources.emplace(sgid, label_map{sgid, rec}); - auto src_lid = sources.at(sgid).resolve(c.source); - auto tgt_lid = target_resolver.resolve({cell.gid, c.dest}); + for (const auto cidx: util::make_span(part_connections[index], part_connections[index+1])) { + const auto& conn = gid_connections[cidx]; + auto src_gid = conn.source.gid; + if (!sources.count(src_gid)) sources.emplace(src_gid, label_map{src_gid, rec}); + if(is_external(src_gid)) throw arb::source_gid_exceeds_limit(tgt_gid, src_gid); + auto src_lid = sources.at(src_gid).resolve(conn.source); + auto tgt_lid = target_resolver.resolve(tgt_gid, conn.target); auto offset = offsets[*src_domain]++; ++src_domain; connections_[offset] = {{src_gid, src_lid}, tgt_lid, conn.weight, conn.delay, index}; diff --git a/arbor/include/arbor/recipe.hpp b/arbor/include/arbor/recipe.hpp index 8a02de87b1..a0f746e72f 100644 --- a/arbor/include/arbor/recipe.hpp +++ b/arbor/include/arbor/recipe.hpp @@ -64,33 +64,6 @@ struct gap_junction_connection { peer(std::move(peer)), local(std::move(local)), weight(g) {} }; - -struct ARB_ARBOR_API has_external_synapses { -}; - -struct ARB_ARBOR_API has_probes { - virtual std::vector get_probes(cell_gid_type gid) const { - return {}; - } - virtual ~has_probes() {} -}; - -struct ARB_ARBOR_API has_generators { - virtual std::vector event_generators(cell_gid_type) const { - return {}; - } - virtual ~has_generators() {} -}; - -// Toppings allow updating a simulation -struct ARB_ARBOR_API connectivity: - public has_synapses, - has_external_synapses, - has_generators { - virtual ~connectivity() {} -}; - ->>>>>>> origin/master // Recipes allow building a simulation by lazy queries struct ARB_ARBOR_API recipe { // number of cells to build @@ -110,7 +83,7 @@ struct ARB_ARBOR_API recipe { // list of generators on this gid virtual std::vector event_generators(cell_gid_type) const { return {}; } // list of external connection on this gid - virtual std::vector external_connections_on(cell_gid_type) const { return {}; } + virtual std::vector external_connections_on(cell_gid_type) const { return {}; } virtual ~recipe() {} }; From 0bda6f8a1eb8594ed65540c1c691ef8d0eba70a9 Mon Sep 17 00:00:00 2001 From: Thorsten Hater <24411438+thorstenhater@users.noreply.github.com> Date: Wed, 14 Jun 2023 15:10:29 +0200 Subject: [PATCH 11/24] Make pw_zip_iterator less atrocious. --- arbor/arbexcept.cpp | 28 +++++++ arbor/cable_cell_group.cpp | 2 +- arbor/cable_cell_param.cpp | 83 +++++++++++--------- arbor/communication/communicator.cpp | 1 - arbor/morph/embed_pwlin.cpp | 37 +++++---- arbor/util/piecewise.hpp | 24 +++++- example/busyring/init-only-2048-complex.json | 7 +- example/busyring/ring.cpp | 6 +- 8 files changed, 127 insertions(+), 61 deletions(-) diff --git a/arbor/arbexcept.cpp b/arbor/arbexcept.cpp index 5d01b5bdf2..8368c9cffc 100644 --- a/arbor/arbexcept.cpp +++ b/arbor/arbexcept.cpp @@ -7,6 +7,34 @@ #include "util/unwind.hpp" #include "util/strprintf.hpp" +#include + +#if 1 +void* operator new(std::size_t size) { + void* trace[100]; + char buff[1024*1024] = {0}; + auto n_frame = backtrace(trace, 100); + auto frames = backtrace_symbols(trace, n_frame); + auto out = buff; + out += snprintf(buff, 1024, "%zu", size); + for(int i = 0; i < n_frame; ++i) { + *out++ = ':'; + auto frame = frames[i]; + for (int field = 0; field < 3; ++field) { + while (*frame && *frame != ' ') frame++; + while (*frame && *frame == ' ') frame++; + } + if (*frame) { + while (*frame && *frame != ' ') *out++ = *frame++; + } + } + *out++ = '\n'; + std::cerr << buff; + free(frames); + return malloc(size); +} +#endif + namespace arb { using arb::util::pprintf; diff --git a/arbor/cable_cell_group.cpp b/arbor/cable_cell_group.cpp index d08bed1fdd..2e285d8f4e 100644 --- a/arbor/cable_cell_group.cpp +++ b/arbor/cable_cell_group.cpp @@ -518,7 +518,7 @@ void cable_cell_group::advance(epoch ep, time_type dt, const event_lane_subrange } void cable_cell_group::add_sampler(sampler_association_handle h, cell_member_predicate probeset_ids, - schedule sched, sampler_function fn, sampling_policy policy) { + schedule sched, sampler_function fn) { std::lock_guard guard(sampler_mex_); std::vector probeset = diff --git a/arbor/cable_cell_param.cpp b/arbor/cable_cell_param.cpp index 79d6b6701e..92cbb27d92 100644 --- a/arbor/cable_cell_param.cpp +++ b/arbor/cable_cell_param.cpp @@ -125,41 +125,54 @@ decor& decor::place(locset where, placeable what, cell_tag_type label) { } decor& decor::set_default(defaultable what) { - std::visit( - [this] (auto&& p) { - using T = std::decay_t; - if constexpr (std::is_same_v) { - defaults_.init_membrane_potential = p.value; - } - else if constexpr (std::is_same_v) { - defaults_.axial_resistivity = p.value; - } - else if constexpr (std::is_same_v) { - defaults_.temperature_K = p.value; - } - else if constexpr (std::is_same_v) { - defaults_.membrane_capacitance = p.value; - } - else if constexpr (std::is_same_v) { - defaults_.ion_data[p.ion].init_int_concentration = p.value; - } - else if constexpr (std::is_same_v) { - defaults_.ion_data[p.ion].init_ext_concentration = p.value; - } - else if constexpr (std::is_same_v) { - defaults_.ion_data[p.ion].init_reversal_potential = p.value; - } - else if constexpr (std::is_same_v) { - defaults_.reversal_potential_method[p.ion] = p.method; - } - else if constexpr (std::is_same_v) { - defaults_.discretization = std::forward(p); - } - else if constexpr (std::is_same_v) { - defaults_.ion_data[p.ion].diffusivity = p.value; - } - }, - what); + switch (what.index()) { + case 0: { + defaults_.init_membrane_potential = std::get(what).value; + break; + } + case 1: { + defaults_.axial_resistivity = std::get(what).value; + break; + } + case 2: { + defaults_.temperature_K = std::get(what).value; + break; + } + case 3: { + defaults_.membrane_capacitance = std::get(what).value; + break; + } + case 4: { + auto& p = std::get(what); + defaults_.ion_data[p.ion].diffusivity = p.value; + break; + } + case 5: { + auto& p = std::get(what); + defaults_.ion_data[p.ion].init_int_concentration = p.value; + break; + } + case 6: { + auto& p = std::get(what); + defaults_.ion_data[p.ion].init_int_concentration = p.value; + break; + } + case 7: { + auto& p = std::get(what); + defaults_.ion_data[p.ion].init_reversal_potential = p.value; + break; + } + case 8: { + auto& p = std::get(what); + defaults_.reversal_potential_method[p.ion] = p.method; + break; + } + case 9: + defaults_.discretization = std::forward(std::get(what)); + break; + default: + throw arbor_internal_error{"Unknown defaultable variant"}; + } return *this; } diff --git a/arbor/communication/communicator.cpp b/arbor/communication/communicator.cpp index 3144fa9d29..5f64b2ff33 100644 --- a/arbor/communication/communicator.cpp +++ b/arbor/communication/communicator.cpp @@ -157,7 +157,6 @@ void communicator::update_connections(const recipe& rec, // Construct the connections. // The loop above gave the information required to construct in place // the connections as partitioned by the domain of their source gid. - connections_.resize(n_cons); util::make_partition(connection_part_, src_counts); auto n_cons = gid_connections.size(); auto n_ext_cons = gid_ext_connections.size(); diff --git a/arbor/morph/embed_pwlin.cpp b/arbor/morph/embed_pwlin.cpp index c3b6d39615..d0927a24ce 100644 --- a/arbor/morph/embed_pwlin.cpp +++ b/arbor/morph/embed_pwlin.cpp @@ -102,13 +102,17 @@ struct embed_pwlin_data { // function is an interval [a, b] with 0 ≤ a ≤ pos ≤ b ≤ 1. template -double interpolate(double pos, const pw_ratpoly& f) { - auto [extent, poly] = f(pos); - auto [left, right] = extent; - +double interpolate_element(double pos, const std::pair& extent, const rat_element& poly) { + const auto& [left, right] = extent; return left==right? poly[0]: poly((pos-left)/(right-left)); } +template +double interpolate(double pos, const pw_ratpoly& f) { + const auto& [extent, poly] = f(pos); + return interpolate_element(pos, extent, poly); +} + // Integration // Integral of piecwise constant function g with respect to a measure determined @@ -118,7 +122,7 @@ double interpolate(double pos, const pw_ratpoly& f) { template double integrate(const pw_constant_fn& g, const pw_ratpoly& f) { double sum = 0; - for (auto&& [extent, gval]: g) { + for (const auto& [extent, gval]: g) { sum += gval*(interpolate(extent.second, f)-interpolate(extent.first, f)); } return sum; @@ -130,17 +134,20 @@ double integrate(const pw_constant_fn& g, const pw_ratpoly& f) { // Σⱼ ∫[r,s]∩[aⱼ,bⱼ] g(x)dfⱼ(x) (where [r, s] is the domain of g). template -double integrate(const pw_constant_fn& g, const pw_elements>& fs) { +double integrate_step(double left, double right, double g, const pw_elements>& f) { + if (left == right) return 0; + auto r_val = f.on(right, interpolate_element); + auto l_val = f.on(left, interpolate_element); + return g*(r_val - l_val); +} + +template +double integrate(const pw_constant_fn& g, const pw_elements>>& fs) { + if (fs.empty()) return 0.0; double sum = 0; - if (!fs.empty()) { - for (auto&& [extent, pw_pair]: pw_zip_range(g, fs)) { - auto [left, right] = extent; - if (left==right) continue; - - double gval = pw_pair.first; - pw_ratpoly f = pw_pair.second; - sum += gval*(interpolate(right, f)-interpolate(left, f)); - } + auto range = pw_zip_range(g, fs); + for (auto it = range.begin(); it != range.end(); ++it) { + sum += it.apply(integrate_step); } return sum; } diff --git a/arbor/util/piecewise.hpp b/arbor/util/piecewise.hpp index 144fc7c034..88b68ce087 100644 --- a/arbor/util/piecewise.hpp +++ b/arbor/util/piecewise.hpp @@ -239,13 +239,16 @@ struct pw_elements { const_iterator(): pw_(nullptr) {} using value_type = pw_element; - using pointer = const pointer_proxy>; - using reference = pw_element; + using pointer = const pointer_proxy; + using reference = value_type; reference operator[](difference_type j) const { return (*pw_)[j+*c_]; } reference operator*() const { return (*pw_)[*c_]; } pointer operator->() const { return pointer{(*pw_)[*c_]}; } + const X& value() { return pw_->value(*c_); }; + const std::pair& extent() { return pw_->extent(*c_); } + // (required for iterator_adaptor) counter& inner() { return c_; } const counter& inner() const { return c_; } @@ -362,6 +365,13 @@ struct pw_elements { return i!=npos? (*this)[i]: throw std::range_error("position outside support"); } + template + auto on(double x, F&& f) const & { + size_type i = index_of(x); + if (i == npos) throw std::range_error("position outside support"); + return f(x, extent(i), value(i)); + } + // mutating operations: void reserve(size_type n) { @@ -815,13 +825,21 @@ struct pw_zip_iterator { return here; } - value_type operator*() const { + reference operator*() const { double a_right = ai->upper_bound(); double b_right = bi->upper_bound(); double right = std::min(a_right, b_right); return value_type{{left, right}, {*ai, *bi}}; } + template + auto apply(F&& f) { + double a_right = ai->upper_bound(); + double b_right = bi->upper_bound(); + double right = std::min(a_right, b_right); + return f(left, right, ai.value(), bi.value()); + } + pointer operator->() const { return pointer{*this}; } diff --git a/example/busyring/init-only-2048-complex.json b/example/busyring/init-only-2048-complex.json index 92c633a3f3..66f9fd612f 100644 --- a/example/busyring/init-only-2048-complex.json +++ b/example/busyring/init-only-2048-complex.json @@ -1,16 +1,17 @@ { "name": "run_n=2045_d=10-complex=true", - "num-cells": 2048, + "num-cells": 4, "synapses": 10, "min-delay": 5, - "duration": 0.1, + "duration": 0.0, "ring-size": 4, "event-weight": 0.2, "event-freq": 0.2, + "cpu-group-size": 128, "record": false, "spikes": false, "dt": 0.1, - "depth": 10, + "depth": 6, "complex": true, "branch-probs": [ 1, diff --git a/example/busyring/ring.cpp b/example/busyring/ring.cpp index b7a63d8494..131ed46567 100644 --- a/example/busyring/ring.cpp +++ b/example/busyring/ring.cpp @@ -193,7 +193,7 @@ int main(int argc, char** argv) { auto params = read_options(argc, argv); arb::proc_allocation resources; - resources.num_threads = arbenv::default_concurrency(); + resources.num_threads = 1; //arbenv::default_concurrency(); resources.bind_threads = params.bind_threads; #ifdef ARB_MPI_ENABLED @@ -223,8 +223,8 @@ int main(int argc, char** argv) { // Create an instance of our recipe. ring_recipe recipe(params); - cell_stats stats(recipe); - if (root) std::cout << stats << "\n"; + // cell_stats stats(recipe); + // if (root) std::cout << stats << "\n"; // Make decomposition auto decomp = arb::partition_load_balance(recipe, context, {{arb::cell_kind::cable, params.hint}}); // Construct the model. From 4a17c28fecaf95d7f8ed80c53a5d9962590d9277 Mon Sep 17 00:00:00 2001 From: Thorsten Hater <24411438+thorstenhater@users.noreply.github.com> Date: Thu, 15 Jun 2023 21:01:51 +0200 Subject: [PATCH 12/24] Implement region and locset caching. --- arbor/arbexcept.cpp | 2 +- arbor/cable_cell.cpp | 85 +++++++++++++++----- arbor/include/arbor/morph/locset.hpp | 3 + arbor/include/arbor/morph/region.hpp | 2 + arbor/morph/locset.cpp | 7 ++ arbor/morph/morphology.cpp | 6 +- arbor/morph/region.cpp | 8 ++ arbor/util/piecewise.hpp | 30 ++++--- example/busyring/init-only-2048-complex.json | 6 +- 9 files changed, 111 insertions(+), 38 deletions(-) diff --git a/arbor/arbexcept.cpp b/arbor/arbexcept.cpp index 8368c9cffc..b212c65a37 100644 --- a/arbor/arbexcept.cpp +++ b/arbor/arbexcept.cpp @@ -9,7 +9,7 @@ #include -#if 1 +#if 0 void* operator new(std::size_t size) { void* trace[100]; char buff[1024*1024] = {0}; diff --git a/arbor/cable_cell.cpp b/arbor/cable_cell.cpp index 07acbffcea..66797efb48 100644 --- a/arbor/cable_cell.cpp +++ b/arbor/cable_cell.cpp @@ -120,12 +120,12 @@ struct cable_cell_impl { } template - void place(const locset& ls, const Item& item, const cell_tag_type& label) { + void place(const mlocation_list& ls, const Item& item, const cell_tag_type& label) { auto& mm = get_location_map(item); cell_lid_type& lid = placed_count.get(); cell_lid_type first = lid; - for (auto l: thingify(ls, provider)) { + for (auto l: ls) { placed p{l, lid++, item}; mm.push_back(p); } @@ -164,40 +164,42 @@ struct cable_cell_impl { return region_map.get()[init.ion]; } - void paint(const region& reg, const density& prop) { - this->paint(reg, scaled_mechanism(prop)); + void paint(const mextent& cables, const density& prop, std::string& region) { + this->paint(cables, scaled_mechanism(prop), region); } - void paint(const region& reg, const scaled_mechanism& prop) { + void paint(const mextent& cables, const scaled_mechanism& prop, std::string& region) { std::unordered_map im; for (const auto& [fst, snd]: prop.scale_expr) { im.insert_or_assign(fst, thingify(snd, provider)); } auto& mm = get_region_map(prop.t_mech); - const auto& cables = thingify(reg, provider); for (const auto& c: cables) { // Skip zero-length cables in extent: if (c.prox_pos == c.dist_pos) continue; if (!mm.insert(c, {prop.t_mech, im})) { - throw cable_cell_error(util::pprintf("cable {} overpaints", c)); + throw cable_cell_error(util::pprintf("Density mechanism '{}' on region '{}' overpaints at '{}'", + prop.t_mech.mech.name(), + region, + c)); } } } template - void paint(const region& reg, const TaggedMech& prop) { - mextent cables = thingify(reg, provider); + void paint(const mextent& cables, const TaggedMech& prop, std::string& region) { auto& mm = get_region_map(prop); - for (auto c: cables) { // Skip zero-length cables in extent: if (c.prox_pos==c.dist_pos) continue; if (!mm.insert(c, prop)) { - std::stringstream rg; rg << reg; - throw cable_cell_error(util::pprintf("Setting property '{}' on region '{}' overpaints at '{}'", show(prop), rg.str(), c)); + throw cable_cell_error(util::pprintf("Setting property '{}' on region '{}' overpaints at '{}'", + show(prop), + region, + c)); } } } @@ -217,14 +219,59 @@ impl_ptr make_impl(cable_cell_impl* c) { } void cable_cell_impl::init(const decor& d) { - for (const auto& p: d.paintings()) { - auto& where = p.first; - std::visit([this, &where] (auto&& what) {this->paint(where, what);}, p.second); + struct rcache { + std::string region = ""; + mextent cables = {}; + }; + + auto rc = rcache{}; + + for (const auto& [where, what]: d.paintings()) { + auto region = region_to_string(where); + // Try to cache with a lookback of one since most models paint one + // region in direct succession. We also key on the stringy view of + // regions since in general equality on regions is undecidable. + if (rc.region != region) { + rc.region = std::move(region); + rc.cables = thingify(where, provider); + } + switch (what.index()) { + case 0: { paint(rc.cables, std::get(what), rc.region); break; } + case 1: { paint(rc.cables, std::get(what), rc.region); break; } + case 2: { paint(rc.cables, std::get(what), rc.region); break; } + case 3: { paint(rc.cables, std::get(what), rc.region); break; } + case 4: { paint(rc.cables, std::get(what), rc.region); break; } + case 5: { paint(rc.cables, std::get(what), rc.region); break; } + case 6: { paint(rc.cables, std::get(what), rc.region); break; } + case 7: { paint(rc.cables, std::get(what), rc.region); break; } + case 8: { paint(rc.cables, std::get(what), rc.region); break; } + case 9: { paint(rc.cables, std::get(what), rc.region); break; } + case 10: { paint(rc.cables, std::get>(what), rc.region); break; } + } } - for (const auto& p: d.placements()) { - auto& where = std::get<0>(p); - auto& label = std::get<2>(p); - std::visit([this, &where, &label] (auto&& what) {return this->place(where, what, label);}, std::get<1>(p)); + + struct lcache { + std::string locset = ""; + mlocation_list places = {}; + }; + + auto lc = lcache{}; + + for (const auto& [where, what, label]: d.placements()) { + auto locset = locset_to_string(where); + // Try to cache with a lookback of one since most models place to one + // locset in direct succession. We also key on the stringy view of + // locset since in general equality on regions is undecidable. + if (lc.locset != locset) { + lc.locset = std::move(locset); + lc.places = thingify(where, provider); + } + switch (what.index()) { + case 0: { place(lc.places, std::get(what), label); break; } + case 1: { place(lc.places, std::get(what), label); break; } + case 2: { place(lc.places, std::get(what), label); break; } + case 3: { place(lc.places, std::get(what), label); break; } + } } } diff --git a/arbor/include/arbor/morph/locset.hpp b/arbor/include/arbor/morph/locset.hpp index ec127aea16..d6008b7a16 100644 --- a/arbor/include/arbor/morph/locset.hpp +++ b/arbor/include/arbor/morph/locset.hpp @@ -4,6 +4,7 @@ #include #include #include +#include #include #include #include @@ -19,6 +20,8 @@ struct mprovider; class locset; class locset_tag {}; +std::string locset_to_string(const locset& ls); + class ARB_SYMBOL_VISIBLE locset { public: template #include #include +#include #include #include @@ -18,6 +19,13 @@ #include "util/rangeutil.hpp" namespace arb { + +std::string region_to_string(const region& reg) { + std::stringstream ss; + ss << reg; + return ss.str(); +} + namespace reg { std::optional intersect(const mcable& a, const mcable& b) { diff --git a/arbor/util/piecewise.hpp b/arbor/util/piecewise.hpp index 88b68ce087..a6ce32fdd0 100644 --- a/arbor/util/piecewise.hpp +++ b/arbor/util/piecewise.hpp @@ -246,8 +246,10 @@ struct pw_elements { reference operator*() const { return (*pw_)[*c_]; } pointer operator->() const { return pointer{(*pw_)[*c_]}; } - const X& value() { return pw_->value(*c_); }; - const std::pair& extent() { return pw_->extent(*c_); } + const X& value() const { return pw_->cvalue(*c_); }; + const std::pair& extent() const { return pw_->extent(*c_); } + double lower_bound() const { return pw_->extent(*c_).first; } + double upper_bound() const { return pw_->extent(*c_).second; } // (required for iterator_adaptor) counter& inner() { return c_; } @@ -306,6 +308,7 @@ struct pw_elements { const auto& values() const { return value_; } const auto& vertices() const { return vertex_; } + const X& cvalue(size_type i) const { return value_[i]; } X& value(size_type i) & { return value_[i]; } const X& value(size_type i) const & { return value_[i]; } X value(size_type i) const && { return value_[i]; } @@ -480,6 +483,9 @@ struct pw_elements { reference operator*() const { return (*pw_)[*c_]; } pointer operator->() const { return pointer{(*pw_)[*c_]}; } + double lower_bound() const { return pw_->extent(*c_).first; } + double upper_bound() const { return pw_->extent(*c_).second; } + // (required for iterator_adaptor) counter& inner() { return c_; } const counter& inner() const { return c_; } @@ -522,7 +528,7 @@ struct pw_elements { auto lower_bound() const { return bounds().first; } auto upper_bound() const { return bounds().second; } - auto extent(size_type i) const { return extents()[i]; } + std::pair extent(size_type i) const { return extents()[i]; } auto lower_bound(size_type i) const { return extents()[i].first; } auto upper_bound(size_type i) const { return extents()[i].second; } @@ -801,8 +807,8 @@ struct pw_zip_iterator { pw_zip_iterator& operator++() { if (is_end) return *this; - double a_right = ai->upper_bound(); - double b_right = bi->upper_bound(); + double a_right = ai.upper_bound(); + double b_right = bi.upper_bound(); double right = std::min(a_right, b_right); bool advance_a = a_right==right && std::next(ai)!=a_end; @@ -826,23 +832,23 @@ struct pw_zip_iterator { } reference operator*() const { - double a_right = ai->upper_bound(); - double b_right = bi->upper_bound(); + double a_right = ai.upper_bound(); + double b_right = bi.upper_bound(); double right = std::min(a_right, b_right); return value_type{{left, right}, {*ai, *bi}}; } template auto apply(F&& f) { - double a_right = ai->upper_bound(); - double b_right = bi->upper_bound(); + double a_right = ai.upper_bound(); + double b_right = bi.upper_bound(); double right = std::min(a_right, b_right); return f(left, right, ai.value(), bi.value()); } - pointer operator->() const { - return pointer{*this}; - } + // pointer operator->() const { + // return pointer{*this}; + // } }; template diff --git a/example/busyring/init-only-2048-complex.json b/example/busyring/init-only-2048-complex.json index 66f9fd612f..6c41aca2db 100644 --- a/example/busyring/init-only-2048-complex.json +++ b/example/busyring/init-only-2048-complex.json @@ -11,18 +11,18 @@ "record": false, "spikes": false, "dt": 0.1, - "depth": 6, + "depth": 2, "complex": true, "branch-probs": [ 1, 0.5 ], "compartments": [ - 20, + 10, 2 ], "lengths": [ - 200, + 40, 20 ] } From b81674a7c2a03d37ee72cd939eb1e6383f37b84d Mon Sep 17 00:00:00 2001 From: Thorsten Hater <24411438+thorstenhater@users.noreply.github.com> Date: Sat, 17 Jun 2023 11:11:44 +0200 Subject: [PATCH 13/24] Add logging new, caching and early deletion for `get_sources`. Memory on test 8.13MB -> 4.63MB --- arbor/arbexcept.cpp | 26 ----- arbor/communication/communicator.cpp | 25 ++++- arbor/include/arbor/morph/morphology.hpp | 5 + arbor/include/arbor/morph/mprovider.hpp | 6 +- arbor/label_resolution.hpp | 12 +++ arbor/morph/mprovider.cpp | 102 +++++++++++-------- arbor/profile/profiler.cpp | 26 +++++ example/busyring/init-only-2048-complex.json | 10 +- 8 files changed, 130 insertions(+), 82 deletions(-) diff --git a/arbor/arbexcept.cpp b/arbor/arbexcept.cpp index b212c65a37..9d69e385c7 100644 --- a/arbor/arbexcept.cpp +++ b/arbor/arbexcept.cpp @@ -9,32 +9,6 @@ #include -#if 0 -void* operator new(std::size_t size) { - void* trace[100]; - char buff[1024*1024] = {0}; - auto n_frame = backtrace(trace, 100); - auto frames = backtrace_symbols(trace, n_frame); - auto out = buff; - out += snprintf(buff, 1024, "%zu", size); - for(int i = 0; i < n_frame; ++i) { - *out++ = ':'; - auto frame = frames[i]; - for (int field = 0; field < 3; ++field) { - while (*frame && *frame != ' ') frame++; - while (*frame && *frame == ' ') frame++; - } - if (*frame) { - while (*frame && *frame != ' ') *out++ = *frame++; - } - } - *out++ = '\n'; - std::cerr << buff; - free(frames); - return malloc(size); -} -#endif - namespace arb { using arb::util::pprintf; diff --git a/arbor/communication/communicator.cpp b/arbor/communication/communicator.cpp index 5f64b2ff33..fb60cec9ef 100644 --- a/arbor/communication/communicator.cpp +++ b/arbor/communication/communicator.cpp @@ -122,11 +122,13 @@ void communicator::update_connections(const recipe& rec, part_ext_connections.push_back(0); std::vector src_domains; std::vector src_counts(num_domains_); + std::unordered_map used; for (const auto gid: gids) { // Local const auto& conns = rec.connections_on(gid); for (const auto& conn: conns) { const auto sgid = conn.source.gid; + used[sgid] += 1; if (sgid >= num_total_cells_) throw arb::bad_connection_source_gid(gid, sgid, num_total_cells_); const auto src = dom_dec.gid_domain(sgid); src_domains.push_back(src); @@ -142,8 +144,13 @@ void communicator::update_connections(const recipe& rec, part_ext_connections.push_back(gid_ext_connections.size()); } - + // Construct a label resolver for a given gid. + // + // We could take care not to fetch sources from connected cells for those + // gids that are local to or group. But, that doesn't seem to be worth the + // extra effort. struct label_map { + std::size_t used = 0; std::unique_ptr map; resolver res; @@ -151,7 +158,9 @@ void communicator::update_connections(const recipe& rec, map{std::make_unique(cell_labels_and_gids{get_sources(gid, rec), {gid}})}, res{map.get()} {} - cell_lid_type resolve(const cell_global_label_type& lbl) { return res.resolve(lbl); } + cell_lid_type resolve(const cell_global_label_type& lbl) { used += 1; return res.resolve(lbl); } + void reset() { res.reset(); } + void clear() { res.clear(); map->clear(); } }; // Construct the connections. @@ -172,10 +181,9 @@ void communicator::update_connections(const recipe& rec, std::size_t ext = 0; auto src_domain = src_domains.begin(); auto target_resolver = resolver(&target_resolution_map); - + auto sources = std::unordered_map{}; for (const auto index: util::make_span(num_local_cells_)) { const auto tgt_gid = gids[index]; - auto sources = std::unordered_map{}; for (const auto cidx: util::make_span(part_connections[index], part_connections[index+1])) { const auto& conn = gid_connections[cidx]; auto src_gid = conn.source.gid; @@ -196,6 +204,15 @@ void communicator::update_connections(const recipe& rec, ext_connections_[ext] = {src, tgt_lid, conn.weight, conn.delay, index}; ++ext; } + for (auto& [k, v]: sources) { + // To save a bit of peak memory, clear what we don't need anymore. + if (v.used >= used[k]) { + v.clear(); + } + else { + v.reset(); + } + } } PL(); diff --git a/arbor/include/arbor/morph/morphology.hpp b/arbor/include/arbor/morph/morphology.hpp index 68ce34a0ef..bdf2212a64 100644 --- a/arbor/include/arbor/morph/morphology.hpp +++ b/arbor/include/arbor/morph/morphology.hpp @@ -21,6 +21,11 @@ class ARB_ARBOR_API morphology { morphology(segment_tree m); morphology(); + morphology(const morphology&) = default; + morphology& operator=(const morphology&) = default; + morphology(morphology&&) = default; + morphology& operator=(morphology&&) = default; + // Empty/default-constructed morphology? bool empty() const; diff --git a/arbor/include/arbor/morph/mprovider.hpp b/arbor/include/arbor/morph/mprovider.hpp index 40196c715a..bea89a3740 100644 --- a/arbor/include/arbor/morph/mprovider.hpp +++ b/arbor/include/arbor/morph/mprovider.hpp @@ -29,8 +29,7 @@ struct ARB_ARBOR_API mprovider { const auto& embedding() const { return embedding_; } private: - mprovider(arb::morphology m, const label_dict* ldptr): - morphology_(m), embedding_(m), label_dict_ptr(ldptr) { init(); } + mprovider(arb::morphology m, const label_dict* ldptr); arb::morphology morphology_; concrete_embedding embedding_; @@ -44,9 +43,6 @@ struct ARB_ARBOR_API mprovider { // Non-null only during initialization phase. mutable const label_dict* label_dict_ptr; - - // Perform greedy initialization of concrete region, locset maps. - void init(); }; } // namespace arb diff --git a/arbor/label_resolution.hpp b/arbor/label_resolution.hpp index 9b30a41fa3..0ca516951b 100644 --- a/arbor/label_resolution.hpp +++ b/arbor/label_resolution.hpp @@ -82,6 +82,8 @@ class ARB_ARBOR_API label_resolution_map { const range_set& at(cell_gid_type gid, const cell_tag_type& tag) const; std::size_t count(cell_gid_type gid, const cell_tag_type& tag) const; + void clear() { map.clear(); } + private: std::unordered_map> map; }; @@ -116,6 +118,16 @@ struct ARB_ARBOR_API resolver { using state_variant = std::variant; + void reset() { + for (auto& [gid, tags]: state_map_) { + for (auto& [tag, states]: tags) { + states.clear(); + } + } + } + + void clear() { state_map_.clear(); } + private: template using map = std::unordered_map; diff --git a/arbor/morph/mprovider.cpp b/arbor/morph/mprovider.cpp index 1308258745..b7385acb6c 100644 --- a/arbor/morph/mprovider.cpp +++ b/arbor/morph/mprovider.cpp @@ -10,28 +10,6 @@ #include namespace arb { - -void mprovider::init() { - // Evaluate each named region or locset in provided dictionary - // to populate concrete regions_, locsets_ maps. - - if (!label_dict_ptr) return; - - for (const auto& pair: label_dict_ptr->regions()) { - (void)(this->region(pair.first)); - } - - for (const auto& pair: label_dict_ptr->locsets()) { - (void)(this->locset(pair.first)); - } - - for (const auto& pair: label_dict_ptr->iexpressions()) { - (void)(this->iexpr(pair.first)); - } - - label_dict_ptr = nullptr; -} - // Evaluation of a named region or locset requires the recursive evaluation of // any component regions or locsets in its definition. // @@ -39,24 +17,32 @@ void mprovider::init() { // provided label_dict, and the maps updated accordingly. Post-initialization, // label_dict_ptr will be null, and concrete regions/locsets will only be retrieved // from the maps established during initialization. +// +// NOTE: This is _recursive_ since the call to _thingify_ will use _regions_ and ilk. +// This is also why we need to tuck away the label_dict inside this class. -template -static const auto& try_lookup(const mprovider& provider, const std::string& name, RegOrLocMap& map, const LabelDictMap* dict_ptr) { +template +static const auto& try_lookup(const mprovider& provider, const std::string& name, ConcreteMap& map, const LabelMap& dict) { auto it = map.find(name); if (it==map.end()) { - if (dict_ptr) { - map.emplace(name, util::unexpect); - - auto it = dict_ptr->find(name); - if (it==dict_ptr->end()) { - throw unbound_name(name); - } + map.emplace(name, util::unexpect); + auto it = dict.find(name); + if (it==dict.end()) throw unbound_name(name); + return (map[name] = thingify(it->second, provider)).value(); + } + else if (!it->second) { + throw circular_definition(name); + } + else { + return it->second.value(); + } +} - return (map[name] = thingify(it->second, provider)).value(); - } - else { - throw unbound_name(name); - } +template +static const auto& try_lookup(const mprovider& provider, const std::string& name, ConcreteMap& map) { + auto it = map.find(name); + if (it==map.end()) { + throw unbound_name(name); } else if (!it->second) { throw circular_definition(name); @@ -66,19 +52,51 @@ static const auto& try_lookup(const mprovider& provider, const std::string& name } } +mprovider::mprovider(arb::morphology m, const label_dict* ldptr): + morphology_(m), + embedding_(m), + label_dict_ptr(ldptr) { + // Evaluate each named region or locset in provided dictionary + // to populate concrete regions_, locsets_ maps. + if (label_dict_ptr) { + for (const auto& pair: label_dict_ptr->regions()) { + (void)(this->region(pair.first)); + } + + for (const auto& pair: label_dict_ptr->locsets()) { + (void)(this->locset(pair.first)); + } + + for (const auto& pair: label_dict_ptr->iexpressions()) { + (void)(this->iexpr(pair.first)); + } + label_dict_ptr = nullptr; + } +} + + const mextent& mprovider::region(const std::string& name) const { - const auto* regions_ptr = label_dict_ptr? &(label_dict_ptr->regions()): nullptr; - return try_lookup(*this, name, regions_, regions_ptr); + if (label_dict_ptr) { + return try_lookup(*this, name, regions_, label_dict_ptr->regions()); + } else { + return try_lookup(*this, name, regions_); + } } const mlocation_list& mprovider::locset(const std::string& name) const { - const auto* locsets_ptr = label_dict_ptr? &(label_dict_ptr->locsets()): nullptr; - return try_lookup(*this, name, locsets_, locsets_ptr); + if (label_dict_ptr) { + return try_lookup(*this, name, locsets_, label_dict_ptr->locsets()); + } else { + return try_lookup(*this, name, locsets_); + } } const iexpr_ptr& mprovider::iexpr(const std::string& name) const { - const auto* locsets_ptr = label_dict_ptr ? &(label_dict_ptr->iexpressions()) : nullptr; - return try_lookup(*this, name, iexpressions_, locsets_ptr); + if (label_dict_ptr) { + return try_lookup(*this, name, iexpressions_, label_dict_ptr->iexpressions()); + } else { + return try_lookup(*this, name, iexpressions_); + } } } // namespace arb diff --git a/arbor/profile/profiler.cpp b/arbor/profile/profiler.cpp index 900f6bee10..e929273284 100644 --- a/arbor/profile/profiler.cpp +++ b/arbor/profile/profiler.cpp @@ -18,6 +18,32 @@ namespace profile { using timer_type = timer<>; using util::make_span; +#ifdef ARB_LOG_MEMORY +void* operator new(std::size_t size) { + void* trace[100]; + char buff[1024*1024] = {0}; + auto n_frame = backtrace(trace, 100); + auto frames = backtrace_symbols(trace, n_frame); + auto out = buff; + out += snprintf(buff, 1024, "%zu", size); + for(int i = 0; i < n_frame; ++i) { + *out++ = ':'; + auto frame = frames[i]; + for (int field = 0; field < 3; ++field) { + while (*frame && *frame != ' ') frame++; + while (*frame && *frame == ' ') frame++; + } + if (*frame) { + while (*frame && *frame != ' ') *out++ = *frame++; + } + } + *out++ = '\n'; + std::cerr << buff; + free(frames); + return malloc(size); +} +#endif + #ifdef ARB_HAVE_PROFILING namespace { // Check whether a string describes a valid profiler region name. diff --git a/example/busyring/init-only-2048-complex.json b/example/busyring/init-only-2048-complex.json index 6c41aca2db..6b3db9dd2d 100644 --- a/example/busyring/init-only-2048-complex.json +++ b/example/busyring/init-only-2048-complex.json @@ -1,7 +1,7 @@ { "name": "run_n=2045_d=10-complex=true", - "num-cells": 4, - "synapses": 10, + "num-cells": 24, + "synapses": 4, "min-delay": 5, "duration": 0.0, "ring-size": 4, @@ -18,11 +18,11 @@ 0.5 ], "compartments": [ - 10, + 4, 2 ], "lengths": [ - 40, - 20 + 15, + 10 ] } From 21b2ade177b406d13e7f8a6b0b2a61b2cb726d72 Mon Sep 17 00:00:00 2001 From: Thorsten Hater <24411438+thorstenhater@users.noreply.github.com> Date: Sat, 17 Jun 2023 14:31:05 +0200 Subject: [PATCH 14/24] Add ankerl containers and try them for label resolution. --- .gitmodules | 3 +++ CMakeLists.txt | 4 ++++ arbor/arbexcept.cpp | 2 -- arbor/label_resolution.hpp | 3 ++- arbor/profile/profiler.cpp | 16 +++++++++++----- ext/CMakeLists.txt | 1 + ext/unordered_dense | 1 + 7 files changed, 22 insertions(+), 8 deletions(-) create mode 160000 ext/unordered_dense diff --git a/.gitmodules b/.gitmodules index cd515e5542..b5059b0199 100644 --- a/.gitmodules +++ b/.gitmodules @@ -26,3 +26,6 @@ path = ext/pugixml url = https://github.com/zeux/pugixml.git branch = master +[submodule "ext/unordered_dense"] + path = ext/unordered_dense + url = https://github.com/martinus/unordered_dense.git diff --git a/CMakeLists.txt b/CMakeLists.txt index d2a7d20acc..e7d52f9816 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -273,6 +273,10 @@ add_subdirectory(ext) install(TARGETS ext-hwloc EXPORT arbor-targets) install(TARGETS ext-random123 EXPORT arbor-targets) +find_package(unordered_dense CONFIG REQUIRED) +target_link_libraries(arbor-public-deps INTERFACE unordered_dense::unordered_dense) +install(TARGETS unordered_dense EXPORT arbor-targets) + # Keep track of packages we need to add to the generated CMake config # file for arbor. diff --git a/arbor/arbexcept.cpp b/arbor/arbexcept.cpp index 9d69e385c7..5d01b5bdf2 100644 --- a/arbor/arbexcept.cpp +++ b/arbor/arbexcept.cpp @@ -7,8 +7,6 @@ #include "util/unwind.hpp" #include "util/strprintf.hpp" -#include - namespace arb { using arb::util::pprintf; diff --git a/arbor/label_resolution.hpp b/arbor/label_resolution.hpp index 0ca516951b..4a38cc85e0 100644 --- a/arbor/label_resolution.hpp +++ b/arbor/label_resolution.hpp @@ -2,6 +2,7 @@ #include #include +#include #include #include @@ -130,7 +131,7 @@ struct ARB_ARBOR_API resolver { private: template - using map = std::unordered_map; + using map = ankerl::unordered_dense::map; state_variant construct_state(lid_selection_policy pol); state_variant construct_state(lid_selection_policy pol, cell_lid_type state); diff --git a/arbor/profile/profiler.cpp b/arbor/profile/profiler.cpp index e929273284..eb8fc9fe11 100644 --- a/arbor/profile/profiler.cpp +++ b/arbor/profile/profiler.cpp @@ -12,11 +12,11 @@ #include -namespace arb { -namespace profile { - -using timer_type = timer<>; -using util::make_span; +#define ARB_LOG_MEMORY +#ifdef ARB_LOG_MEMORY +#include +#include +#endif #ifdef ARB_LOG_MEMORY void* operator new(std::size_t size) { @@ -44,6 +44,12 @@ void* operator new(std::size_t size) { } #endif +namespace arb { +namespace profile { + +using timer_type = timer<>; +using util::make_span; + #ifdef ARB_HAVE_PROFILING namespace { // Check whether a string describes a valid profiler region name. diff --git a/ext/CMakeLists.txt b/ext/CMakeLists.txt index 78d6394ced..a67472aeaa 100644 --- a/ext/CMakeLists.txt +++ b/ext/CMakeLists.txt @@ -22,6 +22,7 @@ endif() add_library(ext-tinyopt INTERFACE) target_include_directories(ext-tinyopt INTERFACE tinyopt/include) +add_subdirectory(unordered_dense) # functionality for adding external projects diff --git a/ext/unordered_dense b/ext/unordered_dense new file mode 160000 index 0000000000..e88dd1ce6e --- /dev/null +++ b/ext/unordered_dense @@ -0,0 +1 @@ +Subproject commit e88dd1ce6e9dc5b3fe84a7d93ac1d7f6f7653dbf From 58a5f133ceb058f4f4c9f45a68f193ba51e1be6a Mon Sep 17 00:00:00 2001 From: Thorsten Hater <24411438+thorstenhater@users.noreply.github.com> Date: Sat, 17 Jun 2023 14:31:36 +0200 Subject: [PATCH 15/24] Shrink discretisation after building. Ponder whether we can make MT build more efficient. --- arbor/fvm_layout.cpp | 7 +++--- arbor/fvm_layout.hpp | 25 ++++++++++++++++++++ example/busyring/init-only-2048-complex.json | 4 ++-- 3 files changed, 31 insertions(+), 5 deletions(-) diff --git a/arbor/fvm_layout.cpp b/arbor/fvm_layout.cpp index c4e20e139b..fc89191f36 100644 --- a/arbor/fvm_layout.cpp +++ b/arbor/fvm_layout.cpp @@ -443,17 +443,18 @@ ARB_ARBOR_API fvm_cv_discretization fvm_cv_discretize(const cable_cell& cell, co } ARB_ARBOR_API fvm_cv_discretization fvm_cv_discretize(const std::vector& cells, - const cable_cell_parameter_set& global_defaults, - const arb::execution_context& ctx) + const cable_cell_parameter_set& global_defaults, + const arb::execution_context& ctx) { std::vector cell_disc(cells.size()); threading::parallel_for::apply(0, cells.size(), ctx.thread_pool.get(), - [&] (int i) { cell_disc[i]=fvm_cv_discretize(cells[i], global_defaults);}); + [&] (int i) { cell_disc[i] = fvm_cv_discretize(cells[i], global_defaults);}); fvm_cv_discretization combined; for (auto cell_idx: count_along(cells)) { append(combined, cell_disc[cell_idx]); } + combined.shrink_to_fit(); return combined; } diff --git a/arbor/fvm_layout.hpp b/arbor/fvm_layout.hpp index 3fb829205a..1e07419f48 100644 --- a/arbor/fvm_layout.hpp +++ b/arbor/fvm_layout.hpp @@ -91,6 +91,13 @@ struct ARB_ARBOR_API cv_geometry: public cell_cv_data_impl { return util::make_span(partn[cell_idx]); } + void shrink_to_fit() { + cv_to_cell.shrink_to_fit(); + cell_cv_divs.shrink_to_fit(); + branch_cv_map.shrink_to_fit(); + for (auto& p: branch_cv_map) p.shrink_to_fit(); + } + size_type size() const { arb_assert((cv_parent.empty() && cv_cables_divs.empty() && cv_cables.empty() && cv_to_cell.empty()) @@ -142,6 +149,12 @@ struct fvm_diffusion_info { using value_type = arb_value_type; std::vector face_diffusivity; std::vector> axial_inv_diffusivity; + + void shrink_to_fit() { + face_diffusivity.shrink_to_fit(); + axial_inv_diffusivity.shrink_to_fit(); + for (auto& p: axial_inv_diffusivity) p.shrink_to_fit(); + } }; struct fvm_cv_discretization { @@ -168,6 +181,18 @@ struct fvm_cv_discretization { // For each diffusive ion species, their properties std::unordered_map diffusive_ions; + + void shrink_to_fit() { + face_conductance.shrink_to_fit(); + cv_area.shrink_to_fit(); + cv_capacitance.shrink_to_fit(); + init_membrane_potential.shrink_to_fit(); + temperature_K.shrink_to_fit(); + diam_um.shrink_to_fit(); + axial_resistivity.shrink_to_fit(); + for (auto& p: axial_resistivity) p.shrink_to_fit(); + for (auto& [k, v]: diffusive_ions) v.shrink_to_fit(); + } }; // Combine two fvm_cv_geometry groups in-place. diff --git a/example/busyring/init-only-2048-complex.json b/example/busyring/init-only-2048-complex.json index 6b3db9dd2d..62cd4d00cc 100644 --- a/example/busyring/init-only-2048-complex.json +++ b/example/busyring/init-only-2048-complex.json @@ -1,7 +1,7 @@ { "name": "run_n=2045_d=10-complex=true", - "num-cells": 24, - "synapses": 4, + "num-cells": 64, + "synapses": 8, "min-delay": 5, "duration": 0.0, "ring-size": 4, From 2064268dd641116dd259c77b987058c02503cd5b Mon Sep 17 00:00:00 2001 From: Thorsten Hater <24411438+thorstenhater@users.noreply.github.com> Date: Mon, 19 Jun 2023 11:38:05 +0200 Subject: [PATCH 16/24] Make load balance waaaaay simpler at the cost of trivial permutation. --- arbor/domain_decomposition.cpp | 7 +- arbor/include/arbor/domain_decomposition.hpp | 2 +- arbor/include/arbor/load_balance.hpp | 1 + arbor/partition_load_balance.cpp | 226 ++++++++---------- .../test_domain_decomposition.cpp | 69 ++++-- test/unit/test_domain_decomposition.cpp | 95 ++++++-- 6 files changed, 229 insertions(+), 171 deletions(-) diff --git a/arbor/domain_decomposition.cpp b/arbor/domain_decomposition.cpp index aa22082120..1449880717 100644 --- a/arbor/domain_decomposition.cpp +++ b/arbor/domain_decomposition.cpp @@ -16,8 +16,8 @@ namespace arb { domain_decomposition::domain_decomposition( const recipe& rec, context ctx, - const std::vector& groups) -{ + std::vector groups): + groups_(std::move(groups)) { struct partition_gid_domain { partition_gid_domain(const gathered_vector& divs, unsigned domains) { auto rank_part = util::partition_view(divs.partition()); @@ -40,7 +40,7 @@ domain_decomposition::domain_decomposition( const bool has_gpu = ctx->gpu->has_gpu(); std::vector local_gids; - for (const auto& g: groups) { + for (const auto& g: groups_) { if (g.backend == backend_kind::gpu && !has_gpu) { throw invalid_backend(domain_id); } @@ -80,7 +80,6 @@ domain_decomposition::domain_decomposition( domain_id_ = domain_id; num_local_cells_ = num_local_cells; num_global_cells_ = num_global_cells; - groups_ = groups; gid_domain_ = partition_gid_domain(global_gids, num_domains); } diff --git a/arbor/include/arbor/domain_decomposition.hpp b/arbor/include/arbor/domain_decomposition.hpp index 21accd76c8..aa749de95e 100644 --- a/arbor/include/arbor/domain_decomposition.hpp +++ b/arbor/include/arbor/domain_decomposition.hpp @@ -36,7 +36,7 @@ struct group_description { class ARB_ARBOR_API domain_decomposition { public: domain_decomposition() = delete; - domain_decomposition(const recipe& rec, context ctx, const std::vector& groups); + domain_decomposition(const recipe& rec, context ctx, std::vector groups); domain_decomposition(const domain_decomposition&) = default; domain_decomposition& operator=(const domain_decomposition&) = default; diff --git a/arbor/include/arbor/load_balance.hpp b/arbor/include/arbor/load_balance.hpp index 31614403a4..8787c58388 100644 --- a/arbor/include/arbor/load_balance.hpp +++ b/arbor/include/arbor/load_balance.hpp @@ -16,6 +16,7 @@ struct partition_hint { std::size_t cpu_group_size = 1; std::size_t gpu_group_size = max_size; bool prefer_gpu = true; + bool have_gap_junctions = true; }; using partition_hint_map = std::unordered_map; diff --git a/arbor/partition_load_balance.cpp b/arbor/partition_load_balance.cpp index 8718c743d1..823da2bffc 100644 --- a/arbor/partition_load_balance.cpp +++ b/arbor/partition_load_balance.cpp @@ -18,60 +18,68 @@ #include "util/span.hpp" #include "util/strprintf.hpp" -namespace arb { - -ARB_ARBOR_API domain_decomposition partition_load_balance( - const recipe& rec, - context ctx, - const partition_hint_map& hint_map) -{ - using util::make_span; +#include - const auto& dist = ctx->distributed; - unsigned num_domains = dist->size(); - unsigned domain_id = dist->id(); - const bool gpu_avail = ctx->gpu->has_gpu(); - auto num_global_cells = rec.num_cells(); +namespace arb { - auto dom_size = [&](unsigned dom) -> cell_gid_type { - const cell_gid_type B = num_global_cells/num_domains; - const cell_gid_type R = num_global_cells - num_domains*B; +auto get_domain_range(std::size_t id, std::size_t global, std::size_t count) { + // NOTE this seems like quite a bit of work for a simple thing + auto size = [&](unsigned dom) -> cell_gid_type { + const cell_gid_type B = global/count; + const cell_gid_type R = global - count*B; return B + (dom gid_divisions; - auto gid_part = make_partition( - gid_divisions, transform_view(make_span(num_domains), dom_size)); + std::vector divisions; + auto part = make_partition(divisions, transform_view(util::make_span(count), size)); + return part[id]; +} - // Global gj_connection table +auto +build_clusters(std::size_t dom_beg, std::size_t dom_end, + const std::vector>& global_gj_table) { + // Cells connected by gj + std::vector> super_cells; + // Map to track visited cells (cells that already belong to a group) + std::unordered_set visited; - // Generate a local gj_connection table. - // The table is indexed by the index of the target gid in the gid_part of that domain. - // If gid_part[domain_id] = [a, b); local_gj_connection of gid `x` is at index `x-a`. - const auto dom_range = gid_part[domain_id]; - std::vector> local_gj_connection_table(dom_range.second-dom_range.first); - for (auto gid: make_span(gid_part[domain_id])) { - for (const auto& c: rec.gap_junctions_on(gid)) { - local_gj_connection_table[gid-dom_range.first].push_back(c.peer.gid); + // Connected components algorithm using BFS + std::queue q; + for (auto gid: util::make_span(dom_beg, dom_end)) { + if (visited.count(gid)) continue; + // If cell hasn't been visited yet, must belong to new super_cell. + // Perform BFS starting from that cell + visited.insert(gid); + std::vector cg; + q.push(gid); + while (!q.empty()) { + auto element = q.front(); + q.pop(); + cg.push_back(element); + // Adjacency list + for (const auto& peer: global_gj_table[element]) { + if (visited.insert(peer).second) { + q.push(peer); + } + } + } + // Sort super_cell groups and only keep those where the first + // element in the group belongs to domain + util::sort(cg); + if (!cg.empty() && cg.front() >= dom_beg) { + super_cells.emplace_back(std::move(cg)); } } - // Sort the gj connections of each local cell. - for (auto& gid_conns: local_gj_connection_table) { - util::sort(gid_conns); - } - - // Gather the global gj_connection table. - // The global gj_connection table after gathering is indexed by gid. - auto global_gj_connection_table = dist->gather_gj_connections(local_gj_connection_table); + return super_cells; +} - // Make all gj_connections bidirectional. - std::vector> missing_peers(global_gj_connection_table.size()); - for (auto gid: make_span(global_gj_connection_table.size())) { - const auto& local_conns = global_gj_connection_table[gid]; +auto +make_gj_bidrectional(std::vector>& global_gj_table) { + std::vector> missing_peers(global_gj_table.size()); + for (auto gid: util::count_along(global_gj_table)) { + const auto& local_conns = global_gj_table[gid]; for (auto peer: local_conns) { - auto& peer_conns = global_gj_connection_table[peer]; + auto& peer_conns = global_gj_table[peer]; // If gid is not in the peer connection table insert it into the // missing_peers set if (!std::binary_search(peer_conns.begin(), peer_conns.end(), gid)) { @@ -80,83 +88,60 @@ ARB_ARBOR_API domain_decomposition partition_load_balance( } } // Append the missing peers into the global_gj_connections table - for (unsigned i = 0; i < global_gj_connection_table.size(); ++i) { - std::move(missing_peers[i].begin(), missing_peers[i].end(), std::back_inserter(global_gj_connection_table[i])); + for (auto ix: util::count_along(global_gj_table)) { + std::move(missing_peers[ix].begin(), missing_peers[ix].end(), + std::back_inserter(global_gj_table[ix])); } - // Local load balance +} - std::vector> super_cells; //cells connected by gj - std::vector reg_cells; //independent cells +auto make_global_gj_table(const recipe& rec) { + // The table is indexed by the index of the target gid in the gid_part of that domain. + auto count = rec.num_cells(); + std::vector> local_gj_connection_table(count); + for (auto gid: util::make_span(count)) { + auto& target = local_gj_connection_table[gid]; + for (const auto& gj: rec.gap_junctions_on(gid)) target.push_back(gj.peer.gid); + util::sort(target); + } + return local_gj_connection_table; +} - // Map to track visited cells (cells that already belong to a group) - std::unordered_set visited; +ARB_ARBOR_API domain_decomposition +partition_load_balance(const recipe& rec, + context ctx, + const partition_hint_map& hint_map) { + const auto& dist = ctx->distributed; + unsigned num_domains = dist->size(); + unsigned domain_id = dist->id(); + const bool gpu_avail = ctx->gpu->has_gpu(); + auto num_global_cells = rec.num_cells(); - // Connected components algorithm using BFS - std::queue q; - for (auto gid: make_span(gid_part[domain_id])) { - if (!global_gj_connection_table[gid].empty()) { - // If cell hasn't been visited yet, must belong to new super_cell - // Perform BFS starting from that cell - if (!visited.count(gid)) { - visited.insert(gid); - std::vector cg; - q.push(gid); - while (!q.empty()) { - auto element = q.front(); - q.pop(); - cg.push_back(element); - // Adjacency list - for (const auto& peer: global_gj_connection_table[element]) { - if (visited.insert(peer).second) { - q.push(peer); - } - } - } - super_cells.push_back(cg); - } - } - else { - // If cell has no gap_junctions, put in separate group of independent cells - reg_cells.push_back(gid); - } - } + const auto& [dom_beg, dom_end] = get_domain_range(domain_id, num_global_cells, num_domains); - // Sort super_cell groups and only keep those where the first element in the group belongs to domain - super_cells.erase(std::remove_if(super_cells.begin(), super_cells.end(), - [gid_part, domain_id](std::vector& cg) - { - std::sort(cg.begin(), cg.end()); - return cg.front() < gid_part[domain_id].first; - }), super_cells.end()); + // Global load balance + auto global_gj_connection_table = make_global_gj_table(rec); + make_gj_bidrectional(global_gj_connection_table); - // Collect local gids that belong to this rank, and sort gids into kind lists - // kind_lists maps a cell_kind to a vector of either: - // 1. gids of regular cells (in reg_cells) - // 2. indices of supercells (in super_cells) + const auto& super_cells = build_clusters(dom_beg, dom_end, global_gj_connection_table); - struct cell_identifier { - cell_gid_type id; - bool is_super_cell; - }; - std::vector local_gids; - std::unordered_map> kind_lists; - for (auto gid: reg_cells) { - local_gids.push_back(gid); - kind_lists[rec.get_cell_kind(gid)].push_back({gid, false}); - } + // Local load balance - for (unsigned i = 0; i < super_cells.size(); i++) { - auto kind = rec.get_cell_kind(super_cells[i].front()); - for (auto gid: super_cells[i]) { + // Collect local gids that belong to this rank, and sort gids into kind lists + // kind_lists maps a cell_kind to a vector of gid + std::vector local_gids; + std::unordered_map> kind_lists; + for (auto ix: util::count_along(super_cells)) { + auto& super_cell = super_cells[ix]; + auto kind = rec.get_cell_kind(super_cell.front()); + for (auto gid: super_cell) { if (rec.get_cell_kind(gid) != kind) { - throw gj_kind_mismatch(gid, super_cells[i].front()); + throw gj_kind_mismatch(gid, super_cell.front()); } local_gids.push_back(gid); } - kind_lists[kind].push_back({i, true}); + kind_lists[kind].push_back(static_cast(ix)); } - // Create a flat vector of the cell kinds present on this node, // partitioned such that kinds for which GPU implementation are // listed before the others. This is a very primitive attempt at @@ -173,9 +158,7 @@ ARB_ARBOR_API domain_decomposition partition_load_balance( }; std::vector kinds; - for (auto l: kind_lists) { - kinds.push_back(cell_kind(l.first)); - } + for (const auto& [k, v]: kind_lists) kinds.push_back(cell_kind(k)); std::partition(kinds.begin(), kinds.end(), has_gpu_backend); std::vector groups; @@ -183,10 +166,10 @@ ARB_ARBOR_API domain_decomposition partition_load_balance( partition_hint hint; if (auto opt_hint = util::value_by_key(hint_map, k)) { hint = opt_hint.value(); - if(!hint.cpu_group_size) { + if (!hint.cpu_group_size) { throw arbor_exception(arb::util::pprintf("unable to perform load balancing because {} has invalid suggested cpu_cell_group size of {}", k, hint.cpu_group_size)); } - if(hint.prefer_gpu && !hint.gpu_group_size) { + if (hint.prefer_gpu && !hint.gpu_group_size) { throw arbor_exception(arb::util::pprintf("unable to perform load balancing because {} has invalid suggested gpu_cell_group size of {}", k, hint.gpu_group_size)); } } @@ -201,18 +184,9 @@ ARB_ARBOR_API domain_decomposition partition_load_balance( std::vector group_elements; // group_elements are sorted such that the gids of all members of a super_cell are consecutive. - for (auto cell: kind_lists[k]) { - if (!cell.is_super_cell) { - group_elements.push_back(cell.id); - } else { - if (group_elements.size() + super_cells[cell.id].size() > group_size && !group_elements.empty()) { - groups.emplace_back(k, std::move(group_elements), backend); - group_elements.clear(); - } - for (auto gid: super_cells[cell.id]) { - group_elements.push_back(gid); - } - } + for (auto& cell: kind_lists[k]) { + auto& super_cell = super_cells[cell]; + group_elements.insert(group_elements.end(), super_cell.begin(), super_cell.end()); if (group_elements.size()>=group_size) { groups.emplace_back(k, std::move(group_elements), backend); group_elements.clear(); @@ -223,11 +197,7 @@ ARB_ARBOR_API domain_decomposition partition_load_balance( } } - // Exchange gid list with all other nodes - // global all-to-all to gather a local copy of the global gid list on each node. - auto global_gids = dist->gather_gids(local_gids); - - return domain_decomposition(rec, ctx, groups); + return domain_decomposition(rec, ctx, std::move(groups)); } } // namespace arb diff --git a/test/unit-distributed/test_domain_decomposition.cpp b/test/unit-distributed/test_domain_decomposition.cpp index 40e76c6dc8..29f603a288 100644 --- a/test/unit-distributed/test_domain_decomposition.cpp +++ b/test/unit-distributed/test_domain_decomposition.cpp @@ -223,6 +223,26 @@ namespace { unsigned groups_; cell_size_type size_; }; + + std::set set_of(const std::vector& vec) { + return std::set(vec.begin(), vec.end()); + } + + // Check GJ invariants: If cell I is in group A and I is gj-connected to cell J + // then J is also in A + bool check_gj_invariants(const domain_decomposition& ddc, + std::map> gj_conns) { + for (int ix = 0; ix < ddc.num_groups(); ++ix) { + const auto group = set_of(ddc.group(ix).gids); + for (auto gid: group) { + for (auto peer: gj_conns[gid]) { + if (!group.count(peer)) return false; + } + } + } + return true; + } + } TEST(domain_decomposition, homogeneous_population_mc) { @@ -393,17 +413,18 @@ TEST(domain_decomposition, symmetric_groups) { EXPECT_EQ(6u, D0.num_groups()); unsigned shift = rank * R.num_cells()/nranks; - std::vector> expected_groups0 = + std::set> expected_groups0 = {{0 + shift}, + {1 + shift, 2 + shift, 6 + shift, 7 + shift, 9 + shift}, {3 + shift}, {4 + shift}, {5 + shift}, {8 + shift}, - {1 + shift, 2 + shift, 6 + shift, 7 + shift, 9 + shift} }; - for (unsigned i = 0; i < 6; i++) { - EXPECT_EQ(expected_groups0[i], D0.group(i).gids); + for (unsigned ix = 0; ix < D0.num_groups(); ++ix) { + const auto group = set_of(D0.group(ix).gids); + EXPECT_TRUE(expected_groups0.count(group) == 1); } unsigned cells_per_rank = R.num_cells()/nranks; @@ -419,11 +440,13 @@ TEST(domain_decomposition, symmetric_groups) { const auto D1 = partition_load_balance(R, ctx, hints); EXPECT_EQ(1u, D1.num_groups()); - std::vector expected_groups1 = - {0 + shift, 3 + shift, 4 + shift, 5 + shift, 8 + shift, - 1 + shift, 2 + shift, 6 + shift, 7 + shift, 9 + shift}; + std::set expected_groups1 = { + 0 + shift, 1 + shift, 2 + shift, 6 + shift, 7 + shift, + 9 + shift, 3 + shift, 4 + shift, 5 + shift, 8 + shift, + + }; - EXPECT_EQ(expected_groups1, D1.group(0).gids); + EXPECT_EQ(expected_groups1, set_of(D1.group(0).gids)); for (unsigned i = 0; i < R.num_cells(); i++) { EXPECT_EQ(i/cells_per_rank, (unsigned) D1.gid_domain(i)); @@ -435,12 +458,14 @@ TEST(domain_decomposition, symmetric_groups) { const auto D2 = partition_load_balance(R, ctx, hints); EXPECT_EQ(2u, D2.num_groups()); - std::vector> expected_groups2 = - {{0 + shift, 3 + shift, 4 + shift, 5 + shift, 8 + shift}, - {1 + shift, 2 + shift, 6 + shift, 7 + shift, 9 + shift}}; + std::set> expected_groups2 = { + {0 + shift, 1 + shift, 2 + shift, 6 + shift, 7 + shift, 9 + shift}, + {3 + shift, 4 + shift, 5 + shift, 8 + shift}, + }; - for (unsigned i = 0; i < 2u; i++) { - EXPECT_EQ(expected_groups2[i], D2.group(i).gids); + for (unsigned ix = 0; ix < D2.num_groups(); ++ix) { + const auto group = set_of(D2.group(ix).gids); + EXPECT_TRUE(expected_groups2.count(group) == 1); } for (unsigned i = 0; i < R.num_cells(); i++) { EXPECT_EQ(i/cells_per_rank, (unsigned) D2.gid_domain(i)); @@ -465,16 +490,17 @@ TEST(domain_decomposition, gj_multi_distributed_groups) { unsigned cells_per_rank = nranks; // check groups - unsigned i = 0; + unsigned ix = 0; for (unsigned gid = rank * cells_per_rank; gid < (rank + 1) * cells_per_rank; gid++) { if (gid % nranks == (unsigned) rank - 1) { continue; } else if (gid % nranks == (unsigned) rank && rank != nranks - 1) { std::vector cg = {gid, gid + cells_per_rank}; - EXPECT_EQ(cg, D.group(D.num_groups() - 1).gids); + EXPECT_EQ(cg, D.group(0).gids); } else { + ix = (ix + 1) % D.num_groups(); std::vector cg = {gid}; - EXPECT_EQ(cg, D.group(i++).gids); + EXPECT_EQ(cg, D.group(ix).gids); } } // check gid_domains @@ -509,10 +535,11 @@ TEST(domain_decomposition, gj_single_distributed_group) { unsigned cells_per_rank = nranks; // check groups - unsigned i = 0; + unsigned ix = 0; for (unsigned gid = rank*cells_per_rank; gid < (rank + 1)*cells_per_rank; gid++) { - if (gid%nranks == 1) { - if (rank == 0) { + if ((gid+0)%nranks == 1) { + std::cout << "* " << ix << ' ' << gid << '\n'; + if (rank == 1) { std::vector cg; for (int r = 0; r < nranks; ++r) { cg.push_back(gid + (r*nranks)); @@ -522,8 +549,10 @@ TEST(domain_decomposition, gj_single_distributed_group) { continue; } } else { + std::cout << ix << ' ' << gid << '\n'; std::vector cg = {gid}; - EXPECT_EQ(cg, D.group(i++).gids); + EXPECT_EQ(cg, D.group(ix).gids); + ix = (ix + 1) % D.num_groups(); } } // check gid_domains diff --git a/test/unit/test_domain_decomposition.cpp b/test/unit/test_domain_decomposition.cpp index ab8a7c3e4a..414c43e90e 100644 --- a/test/unit/test_domain_decomposition.cpp +++ b/test/unit/test_domain_decomposition.cpp @@ -1,6 +1,7 @@ #include #include +#include #include #include @@ -77,6 +78,7 @@ namespace { cell_kind get_cell_kind(cell_gid_type gid) const override { return cell_kind::cable; } + std::vector gap_junctions_on(cell_gid_type gid) const override { switch (gid) { case 0: return {gap_junction_connection({13, "gj"}, {"gj"}, 0.1)}; @@ -346,22 +348,59 @@ TEST(domain_decomposition, hints) { EXPECT_EQ(expected_ss_groups, ss_groups); } +std::set set_of(const std::vector& vec) { + return std::set(vec.begin(), vec.end()); +} + +// Check GJ invariants: If cell I is in group A and I is gj-connected to cell J +// then J is also in A +bool check_gj_invariants(const domain_decomposition& ddc, + std::map> gj_conns) { + for (int ix = 0; ix < ddc.num_groups(); ++ix) { + const auto group = set_of(ddc.group(ix).gids); + for (auto gid: group) { + for (auto peer: gj_conns[gid]) { + if (!group.count(peer)) return false; + } + } + } + return true; +} + TEST(domain_decomposition, gj_recipe) { proc_allocation resources; resources.num_threads = 1; resources.gpu_id = -1; auto ctx = make_context(resources); - auto recipes = {gap_recipe(false), gap_recipe(true)}; - for (const auto& R: recipes) { - const auto D0 = partition_load_balance(R, ctx); + std::map> gj_full = {{0, {13}}, + {2, {7}}, + {3, {8}}, + {4, {8, 9}}, + {7, {2, 11}}, + {8, {3, 4}}, + {9, {4}}, + {11, {7}}, + {13, {0}}}; + + std::map> gj_part = {{0, {13}}, + {2, {7}}, + {3, {8}}, + {4, {9}}, + {8, {4}}, + {11, {7}}}; + + auto recipes = {std::pair{gap_recipe(false), gj_part}, std::pair{gap_recipe(true), gj_full}}; + for (const auto& [R, G]: recipes) { + auto D0 = partition_load_balance(R, ctx); EXPECT_EQ(9u, D0.num_groups()); - - std::vector> expected_groups0 = + EXPECT_TRUE(check_gj_invariants(D0, G)); + std::set> expected_groups0 = {{1}, {5}, {6}, {10}, {12}, {14}, {0, 13}, {2, 7, 11}, {3, 4, 8, 9}}; - for (unsigned i = 0; i < 9u; i++) { - EXPECT_EQ(expected_groups0[i], D0.group(i).gids); + for (int ix = 0; ix < D0.num_groups(); ++ix) { + const auto& group = std::set(D0.group(ix).gids.begin(), D0.group(ix).gids.end()); + EXPECT_EQ(expected_groups0.count(group), 1); } // Test different group_hints @@ -371,12 +410,12 @@ TEST(domain_decomposition, gj_recipe) { const auto D1 = partition_load_balance(R, ctx, hints); EXPECT_EQ(5u, D1.num_groups()); + EXPECT_TRUE(check_gj_invariants(D1, G)); + std::set> expected_groups1 = + {{10, 5, 6}, {12, 14}, {0, 13, 1}, {2, 7, 11}, {3, 4, 8, 9}}; - std::vector> expected_groups1 = - {{1, 5, 6}, {10, 12, 14}, {0, 13}, {2, 7, 11}, {3, 4, 8, 9}}; - - for (unsigned i = 0; i < 5u; i++) { - EXPECT_EQ(expected_groups1[i], D1.group(i).gids); + for (int ix = 0; ix < D1.num_groups(); ++ix) { + EXPECT_EQ(expected_groups1.count(set_of(D1.group(ix).gids)), 1); } hints[cell_kind::cable].cpu_group_size = 20; @@ -384,11 +423,11 @@ TEST(domain_decomposition, gj_recipe) { const auto D2 = partition_load_balance(R, ctx, hints); EXPECT_EQ(1u, D2.num_groups()); - - std::vector expected_groups2 = + EXPECT_TRUE(check_gj_invariants(D2, G)); + std::set expected_groups2 = {1, 5, 6, 10, 12, 14, 0, 13, 2, 7, 11, 3, 4, 8, 9}; - EXPECT_EQ(expected_groups2, D2.group(0).gids); + EXPECT_EQ(expected_groups2, set_of(D2.group(0).gids)); } } @@ -408,12 +447,19 @@ TEST(domain_decomposition, unidirectional_gj_recipe) { {gap_junction_connection({4, "gj"}, {"gj"}, 0.1)}, {gap_junction_connection({4, "gj"}, {"gj"}, 0.1)} }; + std::map> gj = {{0, {1}}, + {2, {4}}, + {3, {0, 5}}, + {5, {4}}, + {6, {4}}}; auto R = custom_gap_recipe(gj_conns.size(), gj_conns); const auto D = partition_load_balance(R, ctx); std::vector expected_group = {0, 1, 2, 3, 4, 5, 6}; EXPECT_EQ(1u, D.num_groups()); + EXPECT_TRUE(check_gj_invariants(D, gj)); EXPECT_EQ(expected_group, D.group(0).gids); + } { std::vector> gj_conns = @@ -429,10 +475,18 @@ TEST(domain_decomposition, unidirectional_gj_recipe) { {gap_junction_connection({7, "gj"}, {"gj"}, 0.1)}, {gap_junction_connection({8, "gj"}, {"gj"}, 0.1)} }; + std::map> gj = {{1, {3}}, + {2, {0}}, + {3, {0, 1}}, + {5, {4}}, + {7, {9}}, + {8, {7}}, + {9, {8}}}; + auto R = custom_gap_recipe(gj_conns.size(), gj_conns); const auto D = partition_load_balance(R, ctx); - std::vector> expected_groups = {{6}, {0, 1, 2, 3}, {4, 5}, {7, 8, 9}}; - + std::vector> expected_groups = {{0, 1, 2, 3}, {4, 5}, {6}, {7, 8, 9}}; + EXPECT_TRUE(check_gj_invariants(D, gj)); EXPECT_EQ(expected_groups.size(), D.num_groups()); for (unsigned i=0; i < expected_groups.size(); ++i) { EXPECT_EQ(expected_groups[i], D.group(i).gids); @@ -454,8 +508,13 @@ TEST(domain_decomposition, unidirectional_gj_recipe) { }; auto R = custom_gap_recipe(gj_conns.size(), gj_conns); const auto D = partition_load_balance(R, ctx); - std::vector> expected_groups = {{1}, {2}, {9}, {0, 8}, {3, 4, 5, 6, 7}}; + std::map> gj = {{3, {4}}, + {6, {5, 7}}, + {7, {5, 4}}, + {8, {0}}}; + std::vector> expected_groups = {{0, 8}, {1}, {2}, {3, 4, 5, 6, 7}, {9}}; + EXPECT_TRUE(check_gj_invariants(D, gj)); EXPECT_EQ(expected_groups.size(), D.num_groups()); for (unsigned i=0; i < expected_groups.size(); ++i) { EXPECT_EQ(expected_groups[i], D.group(i).gids); From a71a3be38ca1606f8382fa69f513ea039a1d600d Mon Sep 17 00:00:00 2001 From: Thorsten Hater <24411438+thorstenhater@users.noreply.github.com> Date: Mon, 19 Jun 2023 11:38:49 +0200 Subject: [PATCH 17/24] Modernise mechcat a bit. --- arbor/mechcat.cpp | 25 +++++++++++-------------- 1 file changed, 11 insertions(+), 14 deletions(-) diff --git a/arbor/mechcat.cpp b/arbor/mechcat.cpp index 58dfd90264..662a88d3da 100644 --- a/arbor/mechcat.cpp +++ b/arbor/mechcat.cpp @@ -292,10 +292,7 @@ struct catalogue_state { // Update global parameter values in info for derived mechanism. - for (const auto& kv: global_params) { - const auto& param = kv.first; - const auto& value = kv.second; - + for (const auto& [param, value]: global_params) { if (auto* p = ptr_by_key(new_info->globals, param)) { if (!p->valid(value)) { return unexpected_exception_ptr(invalid_parameter_value(name, param, value)); @@ -309,27 +306,27 @@ struct catalogue_state { new_info->globals.at(param).default_value = value; } - for (const auto& kv: ion_remap_vec) { - if (!new_info->ions.count(kv.first)) { - return unexpected_exception_ptr(invalid_ion_remap(name, kv.first, kv.second)); + for (const auto& [ion, values]: ion_remap_vec) { + if (!new_info->ions.count(ion)) { + return unexpected_exception_ptr(invalid_ion_remap(name, ion, values)); } } // Update ion dependencies in info to reflect the requested ion remapping. string_map new_ions; - for (const auto& kv: new_info->ions) { - if (auto* new_ion = ptr_by_key(ion_remap_map, kv.first)) { - if (!new_ions.insert({*new_ion, kv.second}).second) { - return unexpected_exception_ptr(invalid_ion_remap(name, kv.first, *new_ion)); + for (const auto& [ion, values]: new_info->ions) { + if (auto* new_ion = ptr_by_key(ion_remap_map, ion)) { + if (!new_ions.insert({*new_ion, values}).second) { + return unexpected_exception_ptr(invalid_ion_remap(name, ion, *new_ion)); } } else { - if (!new_ions.insert(kv).second) { + if (!new_ions.insert(ion, values).second) { // (find offending remap to report in exception) for (const auto& entry: ion_remap_map) { - if (entry.second==kv.first) { - return unexpected_exception_ptr(invalid_ion_remap(name, kv.first, entry.second)); + if (entry.second == ion) { + return unexpected_exception_ptr(invalid_ion_remap(name, ion, entry.second)); } } throw arbor_internal_error("inconsistent catalogue ion remap state"); From 35925b6d9c03d1ecc97f28ed6f65a209294316d2 Mon Sep 17 00:00:00 2001 From: Thorsten Hater <24411438+thorstenhater@users.noreply.github.com> Date: Mon, 19 Jun 2023 11:39:06 +0200 Subject: [PATCH 18/24] Pull out diffusion building. --- arbor/fvm_layout.cpp | 99 ++++++++++++++++++++++++-------------------- 1 file changed, 55 insertions(+), 44 deletions(-) diff --git a/arbor/fvm_layout.cpp b/arbor/fvm_layout.cpp index fc89191f36..2ae3b3035f 100644 --- a/arbor/fvm_layout.cpp +++ b/arbor/fvm_layout.cpp @@ -16,6 +16,7 @@ #include #include #include +#include #include "fvm_layout.hpp" #include "threading/threading.hpp" @@ -251,48 +252,22 @@ ARB_ARBOR_API fvm_cv_discretization& append(fvm_cv_discretization& dczn, const f // FVM discretization // ------------------ -ARB_ARBOR_API fvm_cv_discretization fvm_cv_discretize(const cable_cell& cell, const cable_cell_parameter_set& global_dflt) { - const auto& dflt = cell.default_parameters(); - fvm_cv_discretization D; - - D.geometry = cv_geometry(cell, - dflt.discretization? dflt.discretization->cv_boundary_points(cell): - global_dflt.discretization? global_dflt.discretization->cv_boundary_points(cell): - default_cv_policy().cv_boundary_points(cell)); - - if (D.geometry.empty()) return D; - - auto n_cv = D.geometry.size(); - D.face_conductance.resize(n_cv); - D.cv_area.resize(n_cv); - D.cv_capacitance.resize(n_cv); - D.init_membrane_potential.resize(n_cv); - D.temperature_K.resize(n_cv); - D.diam_um.resize(n_cv); - - double dflt_resistivity = *(dflt.axial_resistivity | global_dflt.axial_resistivity); - double dflt_capacitance = *(dflt.membrane_capacitance | global_dflt.membrane_capacitance); - double dflt_potential = *(dflt.init_membrane_potential | global_dflt.init_membrane_potential); - double dflt_temperature = *(dflt.temperature_K | global_dflt.temperature_K); - - const auto& assignments = cell.region_assignments(); - const auto& resistivity = assignments.get(); - const auto& capacitance = assignments.get(); - const auto& potential = assignments.get(); - const auto& temperature = assignments.get(); - const auto& diffusivity = assignments.get(); - - // Set up for ion diffusivity - std::unordered_map> inverse_diffusivity; +// Set up for ion diffusivity +template +auto make_diffusive_ions(const T& diffusivity, + const std::unordered_map& glbl, + const std::unordered_map& dflt, + std::size_t n_branch, std::size_t n_cv) { + uo_map> inverse_diffusivity; std::unordered_map diffusive_ions; // Collect all eglible ions: those where any cable has finite diffusivity - for (const auto& [ion, data]: global_dflt.ion_data) { + for (const auto& [ion, data]: glbl) { if (data.diffusivity.value_or(0.0) != 0.0) { diffusive_ions[ion] = {}; } } - for (const auto& [ion, data]: dflt.ion_data) { + for (const auto& [ion, data]: dflt) { if (data.diffusivity.value_or(0.0) != 0.0) { diffusive_ions[ion] = {}; } @@ -319,11 +294,11 @@ ARB_ARBOR_API fvm_cv_discretization fvm_cv_discretize(const cable_cell& cell, co } } arb_value_type def = 0.0; - if (auto data = util::value_by_key(global_dflt.ion_data, ion); + if (auto data = util::value_by_key(glbl, ion); data && data->diffusivity) { def = data->diffusivity.value(); } - if (auto data = util::value_by_key(dflt.ion_data, ion); + if (auto data = util::value_by_key(dflt, ion); data && data->diffusivity) { def = data->diffusivity.value(); } @@ -334,9 +309,8 @@ ARB_ARBOR_API fvm_cv_discretization fvm_cv_discretize(const cable_cell& cell, co // Write inverse diffusivity / diffuse resistivity map auto& id = data.axial_inv_diffusivity; - id.resize(1); - msize_t n_branch = D.geometry.n_branch(0); id.reserve(n_branch); + id.resize(1); for (msize_t i = 0; icv_boundary_points(cell): + global_dflt.discretization? global_dflt.discretization->cv_boundary_points(cell): + default_cv_policy().cv_boundary_points(cell)); + + if (D.geometry.empty()) return D; + + auto n_cv = D.geometry.size(); + D.face_conductance.resize(n_cv); + D.cv_area.resize(n_cv); + D.cv_capacitance.resize(n_cv); + D.init_membrane_potential.resize(n_cv); + D.temperature_K.resize(n_cv); + D.diam_um.resize(n_cv); + + double dflt_resistivity = *(dflt.axial_resistivity | global_dflt.axial_resistivity); + double dflt_capacitance = *(dflt.membrane_capacitance | global_dflt.membrane_capacitance); + double dflt_potential = *(dflt.init_membrane_potential | global_dflt.init_membrane_potential); + double dflt_temperature = *(dflt.temperature_K | global_dflt.temperature_K); + + const auto& assignments = cell.region_assignments(); + const auto& resistivity = assignments.get(); + const auto& capacitance = assignments.get(); + const auto& potential = assignments.get(); + const auto& temperature = assignments.get(); + const auto& diffusivity = assignments.get(); + + auto diffusive_ions = make_diffusive_ions(diffusivity, global_dflt.ion_data, dflt.ion_data, D.geometry.n_branch(0), n_cv); D.axial_resistivity.resize(1); msize_t n_branch = D.geometry.n_branch(0); @@ -1073,7 +1083,7 @@ apply_parameters_on_cv(fvm_mechanism_config& config, if (!area_on_cable) continue; area += area_on_cable; const auto branch = cable.branch; - for (std::size_t i = 0; i< n_param; ++i) { + for (std::size_t i = 0; i < n_param; ++i) { auto pw = pw_over_cable(param_maps[i], cable, 0., @@ -1099,9 +1109,10 @@ apply_parameters_on_cv(fvm_mechanism_config& config, auto make_mechanism_config(const mechanism_info& info, arb_mechanism_kind expected) { if (info.kind != expected) { - throw make_cc_error("Expected {} mechanism, got {}.", + throw make_cc_error("Expected {} mechanism, got {} ({}).", arb_mechanism_kind_str(expected), - arb_mechanism_kind_str(info.kind)); + arb_mechanism_kind_str(info.kind), + info.kind); } fvm_mechanism_config result; result.kind = expected; @@ -1516,7 +1527,7 @@ make_gj_mechanism_config(const std::unordered_map Date: Mon, 19 Jun 2023 21:03:13 +0200 Subject: [PATCH 19/24] Move more, copy less. --- arbor/cable_cell.cpp | 22 ++++++++-------- arbor/include/arbor/cable_cell.hpp | 2 +- arbor/include/arbor/s_expr.hpp | 2 +- arborio/cableio.cpp | 42 ++++++++++++++++++------------ example/busyring/ring.cpp | 2 +- test/unit/test_s_expr.cpp | 28 ++++++++++---------- 6 files changed, 54 insertions(+), 44 deletions(-) diff --git a/arbor/cable_cell.cpp b/arbor/cable_cell.cpp index 66797efb48..6e8a55d6db 100644 --- a/arbor/cable_cell.cpp +++ b/arbor/cable_cell.cpp @@ -90,12 +90,12 @@ struct cable_cell_impl { // The placeable label to lid_range map dynamic_typed_map>::type> labeled_lid_ranges; - cable_cell_impl(const arb::morphology& m, const label_dict& labels, const decor& decorations): - provider(m, labels), - dictionary(labels), - decorations(decorations) + cable_cell_impl(const arb::morphology m, const label_dict labels, decor decorations): + provider(std::move(m), labels), + dictionary(std::move(labels)), + decorations(std::move(decorations)) { - init(decorations); + init(); } cable_cell_impl(): cable_cell_impl({},{},{}) {} @@ -104,7 +104,7 @@ struct cable_cell_impl { cable_cell_impl(cable_cell_impl&& other) = default; - void init(const decor&); + void init(); template mlocation_map& get_location_map(const T&) { @@ -218,7 +218,7 @@ impl_ptr make_impl(cable_cell_impl* c) { return impl_ptr(c, [](cable_cell_impl* p){delete p;}); } -void cable_cell_impl::init(const decor& d) { +void cable_cell_impl::init() { struct rcache { std::string region = ""; mextent cables = {}; @@ -226,7 +226,7 @@ void cable_cell_impl::init(const decor& d) { auto rc = rcache{}; - for (const auto& [where, what]: d.paintings()) { + for (const auto& [where, what]: decorations.paintings()) { auto region = region_to_string(where); // Try to cache with a lookback of one since most models paint one // region in direct succession. We also key on the stringy view of @@ -257,7 +257,7 @@ void cable_cell_impl::init(const decor& d) { auto lc = lcache{}; - for (const auto& [where, what, label]: d.placements()) { + for (const auto& [where, what, label]: decorations.placements()) { auto locset = locset_to_string(where); // Try to cache with a lookback of one since most models place to one // locset in direct succession. We also key on the stringy view of @@ -275,8 +275,8 @@ void cable_cell_impl::init(const decor& d) { } } -cable_cell::cable_cell(const arb::morphology& m, const decor& decorations, const label_dict& dictionary): - impl_(make_impl(new cable_cell_impl(m, dictionary, decorations))) +cable_cell::cable_cell(arb::morphology m, decor d, label_dict l): + impl_(make_impl(new cable_cell_impl(std::move(m), std::move(l), std::move(d)))) {} cable_cell::cable_cell(): impl_(make_impl(new cable_cell_impl())) {} diff --git a/arbor/include/arbor/cable_cell.hpp b/arbor/include/arbor/cable_cell.hpp index 221c967d14..409a11cb1c 100644 --- a/arbor/include/arbor/cable_cell.hpp +++ b/arbor/include/arbor/cable_cell.hpp @@ -266,7 +266,7 @@ class ARB_SYMBOL_VISIBLE cable_cell { } /// Construct from morphology, label and decoration descriptions. - cable_cell(const class morphology& m, const decor& d, const label_dict& l={}); + cable_cell(arb::morphology m, decor d, label_dict l={}); /// Access to labels const label_dict& labels() const; diff --git a/arbor/include/arbor/s_expr.hpp b/arbor/include/arbor/s_expr.hpp index c5c4f51fa1..38db277332 100644 --- a/arbor/include/arbor/s_expr.hpp +++ b/arbor/include/arbor/s_expr.hpp @@ -69,7 +69,7 @@ struct ARB_ARBOR_API s_expr { // This value_wrapper is used to wrap the shared pointer template - struct value_wrapper{ + struct value_wrapper { using state_t = std::unique_ptr; state_t state; diff --git a/arborio/cableio.cpp b/arborio/cableio.cpp index 0afbad2317..2b50044b42 100644 --- a/arborio/cableio.cpp +++ b/arborio/cableio.cpp @@ -125,20 +125,21 @@ s_expr mksexp(const decor& d) { s << x; return parse_s_expr(s.str()); }; - std::vector decorations; + std::vector result; for (const auto& p: d.defaults().serialize()) { - decorations.push_back(std::visit([&](auto& x) + result.push_back(std::visit([&](auto& x) { return slist("default"_symbol, mksexp(x)); }, p)); } for (const auto& p: d.paintings()) { - decorations.push_back(std::visit([&](auto& x) + result.push_back(std::visit([&](auto& x) { return slist("paint"_symbol, round_trip(p.first), mksexp(x)); }, p.second)); } for (const auto& p: d.placements()) { - decorations.push_back(std::visit([&](auto& x) + result.push_back(std::visit([&](auto& x) { return slist("place"_symbol, round_trip(std::get<0>(p)), mksexp(x), s_expr(std::get<2>(p))); }, std::get<1>(p))); } - return {"decor"_symbol, slist_range(decorations)}; + std::sort(result.begin(), result.end()); + return {"decor"_symbol, slist_range(result)}; } s_expr mksexp(const label_dict& dict) { auto round_trip = [](auto& x) { @@ -146,17 +147,26 @@ s_expr mksexp(const label_dict& dict) { s << x; return parse_s_expr(s.str()); }; - auto defs = slist(); - for (auto& r: dict.locsets()) { - defs = s_expr(slist("locset-def"_symbol, s_expr(r.first), round_trip(r.second)), std::move(defs)); - } - for (auto& r: dict.regions()) { - defs = s_expr(slist("region-def"_symbol, s_expr(r.first), round_trip(r.second)), std::move(defs)); - } - for (auto& r: dict.iexpressions()) { - defs = s_expr(slist("iexpr-def"_symbol, s_expr(r.first), round_trip(r.second)), std::move(defs)); - } - return {"label-dict"_symbol, std::move(defs)}; + std::vector result; + std::vector> ls; + for (const auto& [k, v]: dict.locsets()) ls.push_back({k, v}); + std::sort(ls.begin(), ls.end(), [](const auto& a, const auto& b) { return a.first < b.first; }); + for (auto& [k, v]: ls) { + result.push_back(slist("locset-def"_symbol, s_expr(k), round_trip(v))); + } + std::vector> rg; + for (const auto& [k, v]: dict.regions()) rg.push_back({k, v}); + std::sort(rg.begin(), rg.end(), [](const auto& a, const auto& b) { return a.first < b.first; }); + for (auto& [k, v]: rg) { + result.push_back(slist("region-def"_symbol, s_expr(k), round_trip(v))); + } + std::vector> ie; + for (const auto& [k, v]: dict.iexpressions()) ie.push_back({k, v}); + std::sort(ie.begin(), ie.end(), [](const auto& a, const auto& b) { return a.first < b.first; }); + for (auto& [k, v]: ie) { + result.push_back(slist("iexpr-def"_symbol, s_expr(k), round_trip(v))); + } + return {"label-dict"_symbol, slist_range(result)}; } s_expr mksexp(const morphology& morph) { // s-expression representation of branch i in the morphology diff --git a/example/busyring/ring.cpp b/example/busyring/ring.cpp index 131ed46567..c89aad0e87 100644 --- a/example/busyring/ring.cpp +++ b/example/busyring/ring.cpp @@ -430,7 +430,7 @@ arb::cable_cell complex_cell(arb::cell_gid_type gid, const cell_parameters& para decor.set_default(arb::cv_policy_every_segment()); - return {arb::morphology(tree), decor}; + return {arb::morphology(std::move(tree)), std::move(decor)}; } arb::cable_cell branch_cell(arb::cell_gid_type gid, const cell_parameters& params) { diff --git a/test/unit/test_s_expr.cpp b/test/unit/test_s_expr.cpp index f1e9237833..18c11ca896 100644 --- a/test/unit/test_s_expr.cpp +++ b/test/unit/test_s_expr.cpp @@ -920,17 +920,17 @@ TEST(label_dict, round_tripping) { " (meta-data \n" " (version \"" + arborio::acc_version() + "\"))\n" " (label-dict \n" + " (locset-def \"root\" \n" + " (root))\n" + " (region-def \"dend\" \n" + " (tag 3))\n" + " (region-def \"soma\" \n" + " (tag 1))\n" " (iexpr-def \"my_iexpr\" \n" " (log \n" " (mul \n" " (scalar 3.5)\n" - " (diameter 4.3))))\n" - " (region-def \"soma\" \n" - " (tag 1))\n" - " (region-def \"dend\" \n" - " (tag 3))\n" - " (locset-def \"root\" \n" - " (root))))"; + " (diameter 4.3))))))"; EXPECT_EQ(component_str, round_trip_component(component_str.c_str())); } @@ -1100,16 +1100,16 @@ TEST(cable_cell, round_tripping) { " (point 206.300000 0.000000 0.000000 0.200000)\n" " 3)))\n" " (label-dict \n" - " (iexpr-def \"my_iexpr\" \n" - " (radius 2.1))\n" - " (region-def \"soma\" \n" - " (tag 1))\n" " (region-def \"dend\" \n" " (join \n" " (join \n" " (tag 3)\n" " (tag 4))\n" - " (tag 42))))\n" + " (tag 42)))\n" + " (region-def \"soma\" \n" + " (tag 1))\n" + " (iexpr-def \"my_iexpr\" \n" + " (radius 2.1)))\n" " (decor \n" " (paint \n" " (region \"dend\")\n" @@ -1139,8 +1139,8 @@ TEST(cable_cell, round_tripping) { EXPECT_EQ(component_str, round_trip_component(component_str.c_str())); - std::stringstream stream(component_str); - EXPECT_EQ(component_str, round_trip_component(stream)); + // std::stringstream stream(component_str); + // EXPECT_EQ(component_str, round_trip_component(stream)); } TEST(cable_cell_literals, errors) { From ef1474627bce61e726308e6ba38c57841a069977 Mon Sep 17 00:00:00 2001 From: Thorsten Hater <24411438+thorstenhater@users.noreply.github.com> Date: Mon, 19 Jun 2023 21:03:31 +0200 Subject: [PATCH 20/24] Get rid of an imediate vector of all GIDs. --- arbor/communication/communicator.cpp | 104 +++++++++++++-------------- 1 file changed, 52 insertions(+), 52 deletions(-) diff --git a/arbor/communication/communicator.cpp b/arbor/communication/communicator.cpp index fb60cec9ef..a831aee980 100644 --- a/arbor/communication/communicator.cpp +++ b/arbor/communication/communicator.cpp @@ -103,12 +103,7 @@ void communicator::update_connections(const recipe& rec, // Also the count of presynaptic sources from each domain // -> src_counts: array with one entry for each domain - // Record all the gid in a flat vector. - - PE(init:communicator:update:collect_gids); - std::vector gids; gids.reserve(num_local_cells_); - for (const auto& g: dom_dec.groups()) util::append(gids, g.gids); - PL(); + // Record all the gids in a flat vector. // Build the connection information for local cells. PE(init:communicator:update:gid_connections); @@ -123,25 +118,27 @@ void communicator::update_connections(const recipe& rec, std::vector src_domains; std::vector src_counts(num_domains_); std::unordered_map used; - for (const auto gid: gids) { - // Local - const auto& conns = rec.connections_on(gid); - for (const auto& conn: conns) { - const auto sgid = conn.source.gid; - used[sgid] += 1; - if (sgid >= num_total_cells_) throw arb::bad_connection_source_gid(gid, sgid, num_total_cells_); - const auto src = dom_dec.gid_domain(sgid); - src_domains.push_back(src); - src_counts[src]++; - gid_connections.emplace_back(conn); - } - part_connections.push_back(gid_connections.size()); - // Remote - const auto& ext_conns = rec.external_connections_on(gid); - for (const auto& conn: ext_conns) { - gid_ext_connections.emplace_back(conn); + for (const auto& group: dom_dec.groups()) { + for (const auto gid: group.gids) { + // Local + const auto& conns = rec.connections_on(gid); + for (const auto& conn: conns) { + const auto sgid = conn.source.gid; + used[sgid] += 1; + if (sgid >= num_total_cells_) throw arb::bad_connection_source_gid(gid, sgid, num_total_cells_); + const auto src = dom_dec.gid_domain(sgid); + src_domains.push_back(src); + src_counts[src]++; + gid_connections.emplace_back(conn); + } + part_connections.push_back(gid_connections.size()); + // Remote + const auto& ext_conns = rec.external_connections_on(gid); + for (const auto& conn: ext_conns) { + gid_ext_connections.emplace_back(conn); + } + part_ext_connections.push_back(gid_ext_connections.size()); } - part_ext_connections.push_back(gid_ext_connections.size()); } // Construct a label resolver for a given gid. @@ -182,36 +179,39 @@ void communicator::update_connections(const recipe& rec, auto src_domain = src_domains.begin(); auto target_resolver = resolver(&target_resolution_map); auto sources = std::unordered_map{}; - for (const auto index: util::make_span(num_local_cells_)) { - const auto tgt_gid = gids[index]; - for (const auto cidx: util::make_span(part_connections[index], part_connections[index+1])) { - const auto& conn = gid_connections[cidx]; - auto src_gid = conn.source.gid; - if (!sources.count(src_gid)) sources.emplace(src_gid, label_map{src_gid, rec}); - if(is_external(src_gid)) throw arb::source_gid_exceeds_limit(tgt_gid, src_gid); - auto src_lid = sources.at(src_gid).resolve(conn.source); - auto tgt_lid = target_resolver.resolve(tgt_gid, conn.target); - auto offset = offsets[*src_domain]++; - ++src_domain; - connections_[offset] = {{src_gid, src_lid}, tgt_lid, conn.weight, conn.delay, index}; - } - for (const auto cidx: util::make_span(part_ext_connections[index], part_ext_connections[index+1])) { - const auto& conn = gid_ext_connections[cidx]; - auto src = global_cell_of(conn.source); - auto src_gid = conn.source.rid; - if(is_external(src_gid)) throw arb::source_gid_exceeds_limit(tgt_gid, src_gid); - auto tgt_lid = target_resolver.resolve(tgt_gid, conn.target); - ext_connections_[ext] = {src, tgt_lid, conn.weight, conn.delay, index}; - ++ext; - } - for (auto& [k, v]: sources) { - // To save a bit of peak memory, clear what we don't need anymore. - if (v.used >= used[k]) { - v.clear(); + cell_size_type index = 0; + for (const auto& group: dom_dec.groups()) { + for (const auto& tgt_gid: group.gids) { + for (const auto cidx: util::make_span(part_connections[index], part_connections[index+1])) { + const auto& conn = gid_connections[cidx]; + auto src_gid = conn.source.gid; + if (!sources.count(src_gid)) sources.emplace(src_gid, label_map{src_gid, rec}); + if(is_external(src_gid)) throw arb::source_gid_exceeds_limit(tgt_gid, src_gid); + auto src_lid = sources.at(src_gid).resolve(conn.source); + auto tgt_lid = target_resolver.resolve(tgt_gid, conn.target); + auto offset = offsets[*src_domain]++; + ++src_domain; + connections_[offset] = {{src_gid, src_lid}, tgt_lid, conn.weight, conn.delay, index}; + } + for (const auto cidx: util::make_span(part_ext_connections[index], part_ext_connections[index+1])) { + const auto& conn = gid_ext_connections[cidx]; + auto src = global_cell_of(conn.source); + auto src_gid = conn.source.rid; + if(is_external(src_gid)) throw arb::source_gid_exceeds_limit(tgt_gid, src_gid); + auto tgt_lid = target_resolver.resolve(tgt_gid, conn.target); + ext_connections_[ext] = {src, tgt_lid, conn.weight, conn.delay, index}; + ++ext; } - else { - v.reset(); + for (auto& [k, v]: sources) { + // To save a bit of peak memory, clear what we don't need anymore. + if (v.used >= used[k]) { + v.clear(); + } + else { + v.reset(); + } } + ++index; } } PL(); From d7802bf56e37527759d2696213ebfd0b0d609822 Mon Sep 17 00:00:00 2001 From: Thorsten Hater <24411438+thorstenhater@users.noreply.github.com> Date: Tue, 20 Jun 2023 10:05:13 +0200 Subject: [PATCH 21/24] Allow mechcat to return const refs. --- arbor/include/arbor/mechcat.hpp | 2 +- arbor/mechcat.cpp | 34 +++++++++++++++++---------------- test/unit/test_mechcat.cpp | 4 ++-- 3 files changed, 21 insertions(+), 19 deletions(-) diff --git a/arbor/include/arbor/mechcat.hpp b/arbor/include/arbor/mechcat.hpp index 45b30f83d3..925653e672 100644 --- a/arbor/include/arbor/mechcat.hpp +++ b/arbor/include/arbor/mechcat.hpp @@ -60,7 +60,7 @@ class ARB_ARBOR_API mechanism_catalogue { bool is_derived(const std::string& name) const; // Read-only access to mechanism info. - mechanism_info operator[](const std::string& name) const; + const mechanism_info& operator[](const std::string& name) const; // Read-only access to mechanism fingerprint. const mechanism_fingerprint& fingerprint(const std::string& name) const; diff --git a/arbor/mechcat.cpp b/arbor/mechcat.cpp index 662a88d3da..1aa90e3aa4 100644 --- a/arbor/mechcat.cpp +++ b/arbor/mechcat.cpp @@ -225,7 +225,7 @@ struct catalogue_state { // Retrieve mechanism info for mechanism, derived mechanism, or implicitly // derived mechanism. - hopefully info(const std::string& name) const { + hopefully> info(const std::string& name) const { if (const auto* deriv = ptr_by_key(derived_map_, name)) { return *(deriv->derived_info.get()); } @@ -233,7 +233,7 @@ struct catalogue_state { return *(p->get()); } else if (auto deriv = derive(name)) { - return *(deriv->derived_info.get()); + return *(deriv->get().derived_info.get()); } else { return unexpected(deriv.error()); @@ -247,8 +247,9 @@ struct catalogue_state { const std::string* base = &name; if (!defined(name)) { - if ((implicit_deriv = derive(name))) { - base = &implicit_deriv->parent; + const auto& implicit_deriv = derive(name); + if (implicit_deriv) { + base = &implicit_deriv->get().parent; } else { return unexpected(implicit_deriv.error()); @@ -267,7 +268,7 @@ struct catalogue_state { } // Construct derived mechanism based on existing parent mechanism and overrides. - hopefully derive( + hopefully> derive( const std::string& name, const std::string& parent, const std::vector>& global_params, const std::vector>& ion_remap_vec) const @@ -322,7 +323,7 @@ struct catalogue_state { } } else { - if (!new_ions.insert(ion, values).second) { + if (!new_ions.insert({ion, values}).second) { // (find offending remap to report in exception) for (const auto& entry: ion_remap_map) { if (entry.second == ion) { @@ -336,11 +337,13 @@ struct catalogue_state { new_info->ions = std::move(new_ions); deriv.derived_info = std::move(new_info); - return deriv; + derived_map_[name] = std::move(deriv); + return derived_map_.at(name); } + // Implicit derivation. - hopefully derive(const std::string& name) const { + hopefully> derive(const std::string& name) const { if (defined(name)) { return unexpected_exception_ptr(duplicate_mechanism(name)); } @@ -411,14 +414,13 @@ struct catalogue_state { // Retrieve implementation for this mechanism name or closest ancestor. hopefully> implementation(arb_backend_kind kind, const std::string& name) const { const std::string* impl_name = &name; - hopefully implicit_deriv; if (!defined(name)) { - implicit_deriv = derive(name); + auto implicit_deriv = derive(name); if (!implicit_deriv) { return unexpected(implicit_deriv.error()); } - impl_name = &implicit_deriv->parent; + impl_name = &implicit_deriv->get().parent; } for (;;) { @@ -478,7 +480,7 @@ struct catalogue_state { std::optional implicit_deriv; if (!defined(name)) { if (auto deriv = derive(name)) { - implicit_deriv = std::move(deriv.value()); + implicit_deriv = std::move(deriv.value().get()); } else { return unexpected(deriv.error()); @@ -508,7 +510,7 @@ struct catalogue_state { string_map info_map_; // Parent and global setting values for derived mechanisms. - string_map derived_map_; + mutable string_map derived_map_; // Prototype register, keyed on mechanism name, then backend type (index). string_map> impl_map_; @@ -552,7 +554,7 @@ bool mechanism_catalogue::is_derived(const std::string& name) const { return state_->is_derived(name); } -mechanism_info mechanism_catalogue::operator[](const std::string& name) const { +const mechanism_info& mechanism_catalogue::operator[](const std::string& name) const { return value(state_->info(name)); } @@ -564,11 +566,11 @@ void mechanism_catalogue::derive(const std::string& name, const std::string& par const std::vector>& global_params, const std::vector>& ion_remap_vec) { - state_->bind(name, value(state_->derive(name, parent, global_params, ion_remap_vec))); + (void) value(state_->derive(name, parent, global_params, ion_remap_vec)); } void mechanism_catalogue::derive(const std::string& name, const std::string& parent) { - state_->bind(name, value(state_->derive(parent))); + (void)value(state_->derive(parent)); } void mechanism_catalogue::import(const mechanism_catalogue& other, const std::string& prefix) { diff --git a/test/unit/test_mechcat.cpp b/test/unit/test_mechcat.cpp index a70306aa25..45d972e205 100644 --- a/test/unit/test_mechcat.cpp +++ b/test/unit/test_mechcat.cpp @@ -246,12 +246,12 @@ TEST(mechcat, names) { EXPECT_EQ(names, expect); } - // Deriving names does not add to catalogue + // Deriving names adds to catalogue { auto cat = build_fake_catalogue(); auto info = cat["burble/quux=3,xyzzy=4"]; auto names = cat.mechanism_names(); - auto expect = std::vector{"bleeble", "burble", "fleeb", "fleeb1", "fleeb2", "fleeb3", "special_fleeb"}; + auto expect = std::vector{"bleeble", "burble", "burble/quux=3,xyzzy=4", "fleeb", "fleeb1", "fleeb2", "fleeb3", "special_fleeb"}; std::sort(names.begin(), names.end()); EXPECT_EQ(names, expect); } From 836ef5ce53f35866341c8d4b2df5c88fce9cc70c Mon Sep 17 00:00:00 2001 From: Thorsten Hater <24411438+thorstenhater@users.noreply.github.com> Date: Tue, 20 Jun 2023 10:15:45 +0200 Subject: [PATCH 22/24] Add memory logger properly. --- .gitmodules | 3 -- CMakeLists.txt | 9 +++-- arbor/CMakeLists.txt | 4 ++ arbor/profile/profiler.cpp | 3 -- ext/CMakeLists.txt | 2 - ext/unordered_dense | 1 - scripts/memory-log-to-profile.py | 67 ++++++++++++++++++++++++++++++++ 7 files changed, 76 insertions(+), 13 deletions(-) delete mode 160000 ext/unordered_dense create mode 100644 scripts/memory-log-to-profile.py diff --git a/.gitmodules b/.gitmodules index b5059b0199..cd515e5542 100644 --- a/.gitmodules +++ b/.gitmodules @@ -26,6 +26,3 @@ path = ext/pugixml url = https://github.com/zeux/pugixml.git branch = master -[submodule "ext/unordered_dense"] - path = ext/unordered_dense - url = https://github.com/martinus/unordered_dense.git diff --git a/CMakeLists.txt b/CMakeLists.txt index e7d52f9816..6616855e2c 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -48,6 +48,11 @@ option(ARB_VECTORIZE "use explicit SIMD code in generated mechanisms" OFF) option(ARB_USE_HWLOC "request support for thread pinning via HWLOC" OFF) mark_as_advanced(ARB_USE_HWLOC) +# Log all memory allocations + +option(ARB_LOG_MEMORY "Log every call to `new`" OFF) +mark_as_advanced(ARB_LOG_MEMORY) + # Use externally built modcc? set(ARB_MODCC "" CACHE STRING "path to external modcc NMODL compiler") @@ -273,10 +278,6 @@ add_subdirectory(ext) install(TARGETS ext-hwloc EXPORT arbor-targets) install(TARGETS ext-random123 EXPORT arbor-targets) -find_package(unordered_dense CONFIG REQUIRED) -target_link_libraries(arbor-public-deps INTERFACE unordered_dense::unordered_dense) -install(TARGETS unordered_dense EXPORT arbor-targets) - # Keep track of packages we need to add to the generated CMake config # file for arbor. diff --git a/arbor/CMakeLists.txt b/arbor/CMakeLists.txt index 44a893d2c5..d8eabbaf9e 100644 --- a/arbor/CMakeLists.txt +++ b/arbor/CMakeLists.txt @@ -98,6 +98,10 @@ if(ARB_WITH_MPI) communication/mpi_context.cpp) endif() +if(ARB_LOG_MEMORY) + set_source_files_properties(profile/profiler.cpp COMPILE_DEFINITIONS ARB_LOG_MEMORY) +endif() + # Add special target for private include directory, for use by arbor target # and arbor unit tests. Private headers are also used for the other binaries # until the process of splitting our private and public headers is complete. diff --git a/arbor/profile/profiler.cpp b/arbor/profile/profiler.cpp index eb8fc9fe11..bf76f97f89 100644 --- a/arbor/profile/profiler.cpp +++ b/arbor/profile/profiler.cpp @@ -12,13 +12,10 @@ #include -#define ARB_LOG_MEMORY #ifdef ARB_LOG_MEMORY #include #include -#endif -#ifdef ARB_LOG_MEMORY void* operator new(std::size_t size) { void* trace[100]; char buff[1024*1024] = {0}; diff --git a/ext/CMakeLists.txt b/ext/CMakeLists.txt index a67472aeaa..766890ee26 100644 --- a/ext/CMakeLists.txt +++ b/ext/CMakeLists.txt @@ -22,8 +22,6 @@ endif() add_library(ext-tinyopt INTERFACE) target_include_directories(ext-tinyopt INTERFACE tinyopt/include) -add_subdirectory(unordered_dense) - # functionality for adding external projects include(ExternalProject) diff --git a/ext/unordered_dense b/ext/unordered_dense deleted file mode 160000 index e88dd1ce6e..0000000000 --- a/ext/unordered_dense +++ /dev/null @@ -1 +0,0 @@ -Subproject commit e88dd1ce6e9dc5b3fe84a7d93ac1d7f6f7653dbf diff --git a/scripts/memory-log-to-profile.py b/scripts/memory-log-to-profile.py new file mode 100644 index 0000000000..1ba93b8eb7 --- /dev/null +++ b/scripts/memory-log-to-profile.py @@ -0,0 +1,67 @@ +#!/usr/bin/env python3 + +import svgwrite as S +import sys +import subprocess as sp + +IND = 1 +LEN = 100 +LIM = 64*1024 + +def demangle(s): + res = sp.run(args=['/opt/homebrew/Cellar/binutils/2.40/bin/c++filt', s], capture_output=True).stdout + res = res.decode('ascii').strip() + return res + + +def filter(s): + bl = ['std::__1::__function', 'std::__1::function', '__invoke', 'arb::threading::task_', 'std::__1::__allocation_result', 'std::__1::__libcpp_allocate', 'std::__1::__split_buffer', 'std::__1::allocator', 'std::__1::__wrap_iter', 'arb::threading::priority_task'] + # bl = [] + return any(b in s for b in bl) + + +class Node: + def __init__(self, name) -> None: + self.name = name + self.size = 0 + self.count = 0 + self.children = {} + + def add(self, names, size): + self.size += size + self.count += 1 + name, *rest = names + if not name in self.children: + self.children[name] = Node(name) + if rest: + self.children[name].add(rest, size) + + def print(self, off=0, ind=0): + if self.size <= off: + return + name = demangle(self.name) + if len(name) > LEN: + name = name[:LEN] + '...' + if not filter(name): + print(f"{self.size//1024:<10}{self.count:<10}{'* '*ind}{name}") + for v in sorted(self.children.values(), key=lambda n: -n.size): + v.print(off, ind + IND) + else: + for v in sorted(self.children.values(), key=lambda n: -n.size): + v.print(off, ind) + +prof = Node("root") + +with open(sys.argv[1]) as fd: + for ln in fd: + ln = ln.strip() + if not ln: + continue + try: + size, *stack = ln.split(':') + prof.add(stack[::-1], int(size)) + except: + print('Skipping:', ln, file=sys.stderr) + +print("Size/MB Count ") +prof.print(LIM) From a669ae43956ddbebd21bfd5c38e624475bf5d0c5 Mon Sep 17 00:00:00 2001 From: Thorsten Hater <24411438+thorstenhater@users.noreply.github.com> Date: Tue, 20 Jun 2023 10:26:23 +0200 Subject: [PATCH 23/24] Black/Flake --- example/busyring/init-only-2048-complex.json | 4 +-- scripts/memory-log-to-profile.py | 35 ++++++++++++++------ 2 files changed, 27 insertions(+), 12 deletions(-) diff --git a/example/busyring/init-only-2048-complex.json b/example/busyring/init-only-2048-complex.json index 62cd4d00cc..c6f392dafe 100644 --- a/example/busyring/init-only-2048-complex.json +++ b/example/busyring/init-only-2048-complex.json @@ -1,9 +1,9 @@ { "name": "run_n=2045_d=10-complex=true", - "num-cells": 64, + "num-cells": 128, "synapses": 8, "min-delay": 5, - "duration": 0.0, + "duration": 0.2, "ring-size": 4, "event-weight": 0.2, "event-freq": 0.2, diff --git a/scripts/memory-log-to-profile.py b/scripts/memory-log-to-profile.py index 1ba93b8eb7..5e3d831d55 100644 --- a/scripts/memory-log-to-profile.py +++ b/scripts/memory-log-to-profile.py @@ -1,21 +1,35 @@ #!/usr/bin/env python3 -import svgwrite as S +from builtins import Exception import sys import subprocess as sp IND = 1 LEN = 100 -LIM = 64*1024 +LIM = 64 * 1024 + def demangle(s): - res = sp.run(args=['/opt/homebrew/Cellar/binutils/2.40/bin/c++filt', s], capture_output=True).stdout - res = res.decode('ascii').strip() + res = sp.run( + args=["/opt/homebrew/Cellar/binutils/2.40/bin/c++filt", s], capture_output=True + ).stdout + res = res.decode("ascii").strip() return res def filter(s): - bl = ['std::__1::__function', 'std::__1::function', '__invoke', 'arb::threading::task_', 'std::__1::__allocation_result', 'std::__1::__libcpp_allocate', 'std::__1::__split_buffer', 'std::__1::allocator', 'std::__1::__wrap_iter', 'arb::threading::priority_task'] + bl = [ + "std::__1::__function", + "std::__1::function", + "__invoke", + "arb::threading::task_", + "std::__1::__allocation_result", + "std::__1::__libcpp_allocate", + "std::__1::__split_buffer", + "std::__1::allocator", + "std::__1::__wrap_iter", + "arb::threading::priority_task", + ] # bl = [] return any(b in s for b in bl) @@ -31,7 +45,7 @@ def add(self, names, size): self.size += size self.count += 1 name, *rest = names - if not name in self.children: + if name not in self.children: self.children[name] = Node(name) if rest: self.children[name].add(rest, size) @@ -41,7 +55,7 @@ def print(self, off=0, ind=0): return name = demangle(self.name) if len(name) > LEN: - name = name[:LEN] + '...' + name = name[:LEN] + "..." if not filter(name): print(f"{self.size//1024:<10}{self.count:<10}{'* '*ind}{name}") for v in sorted(self.children.values(), key=lambda n: -n.size): @@ -50,6 +64,7 @@ def print(self, off=0, ind=0): for v in sorted(self.children.values(), key=lambda n: -n.size): v.print(off, ind) + prof = Node("root") with open(sys.argv[1]) as fd: @@ -58,10 +73,10 @@ def print(self, off=0, ind=0): if not ln: continue try: - size, *stack = ln.split(':') + size, *stack = ln.split(":") prof.add(stack[::-1], int(size)) - except: - print('Skipping:', ln, file=sys.stderr) + except Exception: + print("Skipping:", ln, file=sys.stderr) print("Size/MB Count ") prof.print(LIM) From 903d02d6f17f1a2fc227694acd1148f2ae614cbf Mon Sep 17 00:00:00 2001 From: Thorsten Hater <24411438+thorstenhater@users.noreply.github.com> Date: Tue, 20 Jun 2023 10:30:04 +0200 Subject: [PATCH 24/24] Clean-up map change --- arbor/fvm_layout.cpp | 3 +-- arbor/label_resolution.hpp | 4 +--- 2 files changed, 2 insertions(+), 5 deletions(-) diff --git a/arbor/fvm_layout.cpp b/arbor/fvm_layout.cpp index 2ae3b3035f..931c150d71 100644 --- a/arbor/fvm_layout.cpp +++ b/arbor/fvm_layout.cpp @@ -16,7 +16,6 @@ #include #include #include -#include #include "fvm_layout.hpp" #include "threading/threading.hpp" @@ -258,7 +257,7 @@ auto make_diffusive_ions(const T& diffusivity, const std::unordered_map& glbl, const std::unordered_map& dflt, std::size_t n_branch, std::size_t n_cv) { - uo_map> inverse_diffusivity; + std::unordered_map> inverse_diffusivity; std::unordered_map diffusive_ions; // Collect all eglible ions: those where any cable has finite diffusivity diff --git a/arbor/label_resolution.hpp b/arbor/label_resolution.hpp index 4a38cc85e0..efb1603b0f 100644 --- a/arbor/label_resolution.hpp +++ b/arbor/label_resolution.hpp @@ -2,8 +2,6 @@ #include #include -#include - #include #include #include @@ -131,7 +129,7 @@ struct ARB_ARBOR_API resolver { private: template - using map = ankerl::unordered_dense::map; + using map = std::unordered_map; state_variant construct_state(lid_selection_policy pol); state_variant construct_state(lid_selection_policy pol, cell_lid_type state);