From 55f93dfe0e37ec3e7deb1cf072bba8cb7d8a18a6 Mon Sep 17 00:00:00 2001 From: Thorsten Hater <24411438+thorstenhater@users.noreply.github.com> Date: Mon, 14 Aug 2023 11:06:28 +0200 Subject: [PATCH 1/6] Use index/get instead of std::visit. --- arbor/cable_cell_param.cpp | 104 ++++++++++++++++++++++--------------- 1 file changed, 61 insertions(+), 43 deletions(-) diff --git a/arbor/cable_cell_param.cpp b/arbor/cable_cell_param.cpp index 909ec7faba..b3bede9cc3 100644 --- a/arbor/cable_cell_param.cpp +++ b/arbor/cable_cell_param.cpp @@ -125,49 +125,67 @@ 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) { - if (p.value.type() != iexpr_type::scalar) throw cable_cell_error{"Default values cannot have a scale."}; - defaults_.init_membrane_potential = *p.value.get_scalar(); - } - else if constexpr (std::is_same_v) { - if (p.value.type() != iexpr_type::scalar) throw cable_cell_error{"Default values cannot have a scale."}; - defaults_.axial_resistivity = *p.value.get_scalar(); - } - else if constexpr (std::is_same_v) { - if (p.value.type() != iexpr_type::scalar) throw cable_cell_error{"Default values cannot have a scale."}; - defaults_.temperature_K = *p.value.get_scalar(); - } - else if constexpr (std::is_same_v) { - if (p.value.type() != iexpr_type::scalar) throw cable_cell_error{"Default values cannot have a scale."}; - defaults_.membrane_capacitance = *p.value.get_scalar(); - } - else if constexpr (std::is_same_v) { - if (p.value.type() != iexpr_type::scalar) throw cable_cell_error{"Default values cannot have a scale."}; - defaults_.ion_data[p.ion].init_int_concentration = *p.value.get_scalar(); - } - else if constexpr (std::is_same_v) { - if (p.value.type() != iexpr_type::scalar) throw cable_cell_error{"Default values cannot have a scale."}; - defaults_.ion_data[p.ion].init_ext_concentration = *p.value.get_scalar(); - } - else if constexpr (std::is_same_v) { - if (p.value.type() != iexpr_type::scalar) throw cable_cell_error{"Default values cannot have a scale."}; - defaults_.ion_data[p.ion].init_reversal_potential = *p.value.get_scalar(); - } - 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) { - if (p.value.type() != iexpr_type::scalar) throw cable_cell_error{"Default values cannot have a scale."}; - defaults_.ion_data[p.ion].diffusivity = p.value.get_scalar(); - } - }, - what); + // NOTE: the index/get approach is considerably faster than std::visit. + switch (what.index()) { + case 0: { + auto& p = std::get(what); + if (p.value.type() != iexpr_type::scalar) throw cable_cell_error{"Default values cannot have a scale."}; + defaults_.init_membrane_potential = *p.value.get_scalar(); + break; + } + case 1: { + auto& p = std::get(what); + if (p.value.type() != iexpr_type::scalar) throw cable_cell_error{"Default values cannot have a scale."}; + defaults_.axial_resistivity = *p.value.get_scalar(); + break; + } + case 2: { + auto& p = std::get(what); + if (p.value.type() != iexpr_type::scalar) throw cable_cell_error{"Default values cannot have a scale."}; + defaults_.temperature_K = *p.value.get_scalar(); + break; + } + case 3: { + auto& p = std::get(what); + if (p.value.type() != iexpr_type::scalar) throw cable_cell_error{"Default values cannot have a scale."}; + defaults_.membrane_capacitance = *p.value.get_scalar(); + break; + } + case 4: { + auto& p = std::get(what); + if (p.value.type() != iexpr_type::scalar) throw cable_cell_error{"Default values cannot have a scale."}; + defaults_.ion_data[p.ion].diffusivity = p.value.get_scalar(); + break; + } + case 5: { + auto& p = std::get(what); + if (p.value.type() != iexpr_type::scalar) throw cable_cell_error{"Default values cannot have a scale."}; + defaults_.ion_data[p.ion].init_ext_concentration = *p.value.get_scalar(); + break; + } + case 6: { + auto& p = std::get(what); + if (p.value.type() != iexpr_type::scalar) throw cable_cell_error{"Default values cannot have a scale."}; + defaults_.ion_data[p.ion].init_int_concentration = p.value.get_scalar(); + break; + } + case 7: { + auto& p = std::get(what); + if (p.value.type() != iexpr_type::scalar) throw cable_cell_error{"Default values cannot have a scale."}; + defaults_.ion_data[p.ion].init_reversal_potential = *p.value.get_scalar(); + 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; } From adae932a331788bee6a60a301e52a1b2eab847ce Mon Sep 17 00:00:00 2001 From: Thorsten Hater <24411438+thorstenhater@users.noreply.github.com> Date: Mon, 14 Aug 2023 17:04:14 +0200 Subject: [PATCH 2/6] Snapshot, doesnt pass --- arbor/cable_cell.cpp | 101 +++++++++++++++++++++++------------ arbor/fvm_layout.cpp | 1 - arbor/morph/embed_pwlin.cpp | 52 +++++++++++++++--- arbor/util/piecewise.hpp | 41 +++++++++++--- arbor/util/pw_over_cable.hpp | 21 +++----- example/busyring/ring.cpp | 5 +- 6 files changed, 154 insertions(+), 67 deletions(-) diff --git a/arbor/cable_cell.cpp b/arbor/cable_cell.cpp index 3afbc82400..d6927768c9 100644 --- a/arbor/cable_cell.cpp +++ b/arbor/cable_cell.cpp @@ -95,7 +95,7 @@ struct cable_cell_impl { dictionary(labels), decorations(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&) { @@ -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,53 +164,41 @@ 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, const std::string& reg) { + this->paint(cables, scaled_mechanism(prop), reg); } - void paint(const region& reg, const scaled_mechanism& prop) { + void paint(const mextent& cables, const scaled_mechanism& prop, const std::string& reg) { 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})) { - std::stringstream rg; rg << reg; throw cable_cell_error(util::pprintf("Setting mechanism '{}' on region '{}' overpaints at cable {}", - prop.t_mech.mech.name(), rg.str(), c)); + prop.t_mech.mech.name(), reg, c)); } } } template - void paint(const region& reg, const TaggedMech& prop) { - mextent cables = thingify(reg, provider); + void paint(const mextent& cables, const TaggedMech& prop, const std::string& reg) { 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)); - } + if (!mm.insert(c, prop)) throw cable_cell_error(util::pprintf("Setting property '{}' on region '{}' overpaints at '{}'", + show(prop), reg, c)); } } - mlocation_list concrete_locset(const locset& l) const { - return thingify(l, provider); - } - - mextent concrete_region(const region& r) const { - return thingify(r, provider); - } + mlocation_list concrete_locset(const locset& l) const { return thingify(l, provider); } + mextent concrete_region(const region& r) const { return thingify(r, provider); } }; using impl_ptr = std::unique_ptr; @@ -218,20 +206,65 @@ 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) { - for (const auto& p: d.paintings()) { - auto& where = p.first; - std::visit([this, &where] (auto&& what) {this->paint(where, what);}, p.second); +void cable_cell_impl::init() { + struct rcache { + std::string region = ""; + mextent cables = {}; + }; + + auto rc = rcache{}; + + for (const auto& [where, what]: decorations.paintings()) { + auto region = util::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]: decorations.placements()) { + auto locset = util::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; } + } } } 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))) + impl_(make_impl(new cable_cell_impl(std::move(m), std::move(dictionary), std::move(decorations)))) {} cable_cell::cable_cell(): impl_(make_impl(new cable_cell_impl())) {} diff --git a/arbor/fvm_layout.cpp b/arbor/fvm_layout.cpp index 5baa95c036..27737750de 100644 --- a/arbor/fvm_layout.cpp +++ b/arbor/fvm_layout.cpp @@ -1204,7 +1204,6 @@ make_density_mechanism_config(const region_assignment& assignments, const auto& info = data.catalogue[name]; auto config = make_mechanism_config(info, arb_mechanism_kind_density); - auto parameters = ordered_parameters(info); auto n_param = parameters.size(); diff --git a/arbor/morph/embed_pwlin.cpp b/arbor/morph/embed_pwlin.cpp index 945dffdf6b..ad063c2cc0 100644 --- a/arbor/morph/embed_pwlin.cpp +++ b/arbor/morph/embed_pwlin.cpp @@ -103,8 +103,8 @@ struct embed_pwlin_data { template double interpolate(double pos, const pw_ratpoly& f) { - auto [extent, poly] = f(pos); - auto [left, right] = extent; + const auto& [extent, poly] = f(pos); + const auto& [left, right] = extent; return left==right? poly[0]: poly((pos-left)/(right-left)); } @@ -118,7 +118,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; @@ -133,7 +133,7 @@ template double integrate(const pw_constant_fn& g, const pw_elements>& fs) { double sum = 0; if (!fs.empty()) { - for (auto&& [extent, pw_pair]: pw_zip_range(g, fs)) { + for (const auto& [extent, pw_pair]: pw_zip_range(g, fs)) { auto [left, right] = extent; if (left==right) continue; @@ -198,19 +198,55 @@ double embed_pwlin::integrate_ixa(const mcable& c) const { // Integrate piecewise function over a cable: static pw_constant_fn restrict_to(const pw_constant_fn& g, double left, double right) { - return pw_zip_with(g, pw_elements{{left, right}}); + return pw_zip_with(g, pw_elements{left, right}); } double embed_pwlin::integrate_length(const mcable& c, const pw_constant_fn& g) const { - return integrate_length(c.branch, restrict_to(g, c.prox_pos, c.dist_pos)); + if (g.empty() || c.dist_pos == c.prox_pos) return 0; + auto rn = pw_elements{c.prox_pos, c.dist_pos}; + const auto& ar = data_->length.at(c.branch); + double sum = 0; + auto pwzi = util::pw_zip_iterator(g, rn); + while (!pwzi.is_end) { + sum += pwzi.apply_left([&ar](auto l, auto r, const auto& vl) { + return vl*(interpolate(r, ar) - interpolate(l, ar)); + }); + ++pwzi; + } + return sum; } double embed_pwlin::integrate_area(const mcable& c, const pw_constant_fn& g) const { - return integrate_area(c.branch, restrict_to(g, c.prox_pos, c.dist_pos)); + if (g.empty() || c.dist_pos == c.prox_pos) return 0; + auto rn = pw_elements{c.prox_pos, c.dist_pos}; + const auto& ar = data_->area.at(c.branch); + double sum = 0; + auto pwzi = util::pw_zip_iterator(g, rn); + while (!pwzi.is_end) { + sum += pwzi.apply_left([&ar](auto l, auto r, const auto& vl) { + return vl*(interpolate(r, ar) - interpolate(l, ar)); + }); + ++pwzi; + } + return sum; } double embed_pwlin::integrate_ixa(const mcable& c, const pw_constant_fn& g) const { - return integrate_ixa(c.branch, restrict_to(g, c.prox_pos, c.dist_pos)); + if (g.empty() || c.dist_pos == c.prox_pos) return 0; + auto rn = pw_elements{c.prox_pos, c.dist_pos}; + auto pw = util::pw_zip_with(g, rn); + const auto& ix = data_->ixa.at(c.branch); + if (ix.empty()) return 0.0; + double sum = 0; + auto pwzi = util::pw_zip_iterator(pw, ix); + while (!pwzi.is_end) { + sum += pwzi.apply([](auto l, auto r, const auto& gval, const pw_ratpoly<1, 1>& f) { + return gval*(interpolate(r, f)-interpolate(l, f)); + }); + ++pwzi; + } + return sum; + } // Subregions defined by geometric inequalities: diff --git a/arbor/util/piecewise.hpp b/arbor/util/piecewise.hpp index ff83d6ea01..f9899b44d7 100644 --- a/arbor/util/piecewise.hpp +++ b/arbor/util/piecewise.hpp @@ -228,6 +228,9 @@ struct pw_elements { counter& inner() { return c_; } const counter& inner() const { return c_; } + double lower_bound() const { return pw_->extent(*c_).first; } + double upper_bound() const { return pw_->extent(*c_).second; } + protected: pw_elements* pw_; counter c_; @@ -250,6 +253,11 @@ struct pw_elements { counter& inner() { return c_; } const counter& inner() const { return c_; } + double lower_bound() const { return pw_->extent(*c_).first; } + double upper_bound() const { return pw_->extent(*c_).second; } + + const X& value() const { return pw_->cvalue(*c_); }; + protected: const pw_elements* pw_; counter c_; @@ -306,6 +314,7 @@ struct pw_elements { 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]; } + const X& cvalue(size_type i) const { return value_[i]; } auto operator[](size_type i) & { return pw_element_proxy{*this, i}; } auto operator[](size_type i) const & { return value_type{extent(i), value(i)}; } @@ -474,6 +483,9 @@ struct pw_elements { counter& inner() { return c_; } const counter& inner() const { return c_; } + double lower_bound() const { return pw_->extent(*c_).first; } + double upper_bound() const { return pw_->extent(*c_).second; } + protected: const pw_elements* pw_; counter c_; @@ -486,12 +498,10 @@ struct pw_elements { pw_elements() = default; template - pw_elements(const VertexSeq& vs) { - assign(vs); - } - - pw_elements(std::initializer_list vs) { - assign(vs); + pw_elements(const VertexSeq& vs) { assign(vs); } + pw_elements(std::initializer_list vs) { assign(vs); } + pw_elements(double lo, double hi) { + push_back(lo, hi); } pw_elements(pw_elements&&) = default; @@ -512,7 +522,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; } @@ -822,6 +832,23 @@ struct pw_zip_iterator { return value_type{{left, right}, {*ai, *bi}}; } + template + auto apply_left(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()); + } + + + 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/arbor/util/pw_over_cable.hpp b/arbor/util/pw_over_cable.hpp index ee699a9fcb..bf701b17e7 100644 --- a/arbor/util/pw_over_cable.hpp +++ b/arbor/util/pw_over_cable.hpp @@ -16,7 +16,10 @@ struct get_value { // Convert mcable_map values to a piecewise function over an mcable. // The projection gives the map from the values in the mcable_map to the values in the piecewise function. template -util::pw_elements pw_over_cable(const mcable_map& mm, mcable cable, U dflt_value, Proj projection = Proj{}) { +util::pw_elements pw_over_cable(const mcable_map& mm, + const mcable& cable, + U&& dflt_value, + Proj projection = Proj{}) { using value_type = typename mcable_map::value_type; msize_t bid = cable.branch; @@ -29,28 +32,20 @@ util::pw_elements pw_over_cable(const mcable_map& mm, mcable cable, U dflt auto map_on_branch = util::make_range( std::equal_range(mm.begin(), mm.end(), bid, [](as_branch a, as_branch b) { return a.value({cable.prox_pos, cable.dist_pos}, {dflt_value}); - } + if (map_on_branch.empty()) return {{cable.prox_pos, cable.dist_pos}, {dflt_value}}; util::pw_elements pw; pw.reserve(2*map_on_branch.size()); double pw_right = 0.0; for (const auto& [elf, els]: map_on_branch) { - if (elf.prox_pos > pw_right) { - pw.push_back(pw_right, elf.prox_pos, dflt_value); - } + if (elf.prox_pos > pw_right) pw.push_back(pw_right, elf.prox_pos, dflt_value); pw.push_back(elf.prox_pos, elf.dist_pos, projection(elf, els)); pw_right = pw.upper_bound(); } - if (pw_right < 1.0) { - pw.push_back(pw_right, 1.0, dflt_value); - } + if (pw_right < 1.0) pw.push_back(pw_right, 1.0, dflt_value); - if (cable.prox_pos!=0 || cable.dist_pos!=1) { - pw = util::pw_zip_with(pw, util::pw_elements<>({cable.prox_pos, cable.dist_pos})); - } + if (cable.prox_pos!=0 || cable.dist_pos!=1) pw = util::pw_zip_with(pw, util::pw_elements<>(cable.prox_pos, cable.dist_pos)); return pw; } } // namespace util diff --git a/example/busyring/ring.cpp b/example/busyring/ring.cpp index ad5b055655..74fea03de6 100644 --- a/example/busyring/ring.cpp +++ b/example/busyring/ring.cpp @@ -285,10 +285,7 @@ int main(int argc, char** argv) { } auto report = arb::profile::make_meter_report(meters, context); - if (root) { - std::cout << report << '\n' - << arb::profile::profiler_summary() << "\n"; - } + if (root) arb::profile::print_profiler_summary(std::cout, 5); } catch (std::exception& e) { std::cerr << "exception caught in ring miniapp: " << e.what() << "\n"; From 4f9d5ddbb32456f72d805a0f094a2800405a8be5 Mon Sep 17 00:00:00 2001 From: Thorsten Hater <24411438+thorstenhater@users.noreply.github.com> Date: Tue, 15 Aug 2023 09:47:51 +0200 Subject: [PATCH 3/6] Fix cases. --- arbor/cable_cell.cpp | 2 +- arbor/cable_cell_param.cpp | 4 ++-- arbor/fvm_layout.cpp | 2 +- arbor/include/arbor/morph/embed_pwlin.hpp | 4 ++-- arbor/util/pw_over_cable.hpp | 2 +- 5 files changed, 7 insertions(+), 7 deletions(-) diff --git a/arbor/cable_cell.cpp b/arbor/cable_cell.cpp index d6927768c9..3334d95f82 100644 --- a/arbor/cable_cell.cpp +++ b/arbor/cable_cell.cpp @@ -264,7 +264,7 @@ void cable_cell_impl::init() { } cable_cell::cable_cell(const arb::morphology& m, const decor& decorations, const label_dict& dictionary): - impl_(make_impl(new cable_cell_impl(std::move(m), std::move(dictionary), std::move(decorations)))) + impl_(make_impl(new cable_cell_impl(m, dictionary, decorations))) {} cable_cell::cable_cell(): impl_(make_impl(new cable_cell_impl())) {} diff --git a/arbor/cable_cell_param.cpp b/arbor/cable_cell_param.cpp index b3bede9cc3..d4062e6dcc 100644 --- a/arbor/cable_cell_param.cpp +++ b/arbor/cable_cell_param.cpp @@ -160,13 +160,13 @@ decor& decor::set_default(defaultable what) { case 5: { auto& p = std::get(what); if (p.value.type() != iexpr_type::scalar) throw cable_cell_error{"Default values cannot have a scale."}; - defaults_.ion_data[p.ion].init_ext_concentration = *p.value.get_scalar(); + defaults_.ion_data[p.ion].init_int_concentration = *p.value.get_scalar(); break; } case 6: { auto& p = std::get(what); if (p.value.type() != iexpr_type::scalar) throw cable_cell_error{"Default values cannot have a scale."}; - defaults_.ion_data[p.ion].init_int_concentration = p.value.get_scalar(); + defaults_.ion_data[p.ion].init_ext_concentration = p.value.get_scalar(); break; } case 7: { diff --git a/arbor/fvm_layout.cpp b/arbor/fvm_layout.cpp index 27737750de..618c913327 100644 --- a/arbor/fvm_layout.cpp +++ b/arbor/fvm_layout.cpp @@ -311,7 +311,7 @@ fvm_cv_discretize(const cable_cell& cell, const cable_cell_parameter_set& global data.end(), [](const auto& kv) { const auto& v = kv.second.value.get_scalar(); - return !v || *v != 0.0 || *v == *v; + return !v || (*v != 0.0 && *v == *v); }); if (diffusive) { // Provide a (non-sensical) default. diff --git a/arbor/include/arbor/morph/embed_pwlin.hpp b/arbor/include/arbor/morph/embed_pwlin.hpp index c282cb1f97..4fc36270fa 100644 --- a/arbor/include/arbor/morph/embed_pwlin.hpp +++ b/arbor/include/arbor/morph/embed_pwlin.hpp @@ -45,14 +45,14 @@ struct ARB_ARBOR_API embed_pwlin { // Computed length of mcable. double integrate_length(const mcable& c) const; - double integrate_length(mlocation proxmal, mlocation distal) const; + double integrate_length(mlocation prox, mlocation dist) const; double integrate_length(const mcable& c, const pw_constant_fn&) const; double integrate_length(msize_t bid, const pw_constant_fn&) const; // Membrane surface area of given mcable. double integrate_area(const mcable& c) const; - double integrate_area(mlocation proxmal, mlocation distal) const; + double integrate_area(mlocation prox, mlocation dist) const; double integrate_area(const mcable& c, const pw_constant_fn&) const; double integrate_area(msize_t bid, const pw_constant_fn&) const; diff --git a/arbor/util/pw_over_cable.hpp b/arbor/util/pw_over_cable.hpp index bf701b17e7..22d0f2d1ed 100644 --- a/arbor/util/pw_over_cable.hpp +++ b/arbor/util/pw_over_cable.hpp @@ -18,7 +18,7 @@ struct get_value { template util::pw_elements pw_over_cable(const mcable_map& mm, const mcable& cable, - U&& dflt_value, + U dflt_value, Proj projection = Proj{}) { using value_type = typename mcable_map::value_type; msize_t bid = cable.branch; From 78f821b7eb190a0e052687f73c62b305cef9b546 Mon Sep 17 00:00:00 2001 From: Thorsten Hater <24411438+thorstenhater@users.noreply.github.com> Date: Tue, 17 Oct 2023 21:24:32 +0200 Subject: [PATCH 4/6] Complain -- ie terminate -- on unknown cases. --- arbor/cable_cell.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/arbor/cable_cell.cpp b/arbor/cable_cell.cpp index 3334d95f82..1d71050402 100644 --- a/arbor/cable_cell.cpp +++ b/arbor/cable_cell.cpp @@ -235,6 +235,7 @@ void cable_cell_impl::init() { 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; } + default: throw arbor_internal_error{"Unknown paintable variant"}; } } @@ -259,6 +260,7 @@ void cable_cell_impl::init() { 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; } + default: throw arbor_internal_error{"Unknown placeable variant"}; } } } From 7ef78b6b8db61e2180604ce37e9e50515f75ed48 Mon Sep 17 00:00:00 2001 From: Thorsten Hater <24411438+thorstenhater@users.noreply.github.com> Date: Tue, 5 Mar 2024 13:44:40 +0100 Subject: [PATCH 5/6] Also add proxy for constant iterator, hoping we'll save the pw_els> case. --- arbor/util/piecewise.hpp | 28 ++++++++++++++++++++++------ 1 file changed, 22 insertions(+), 6 deletions(-) diff --git a/arbor/util/piecewise.hpp b/arbor/util/piecewise.hpp index f9899b44d7..857aabcd7a 100644 --- a/arbor/util/piecewise.hpp +++ b/arbor/util/piecewise.hpp @@ -160,16 +160,32 @@ struct pw_element_proxy { extent(pw.extent(i)), value(pw.value(i)) {} operator pw_element() const { return pw_element{extent, value}; } - operator X() const { return value; } - pw_element_proxy& operator=(X x) {value = std::move(x); return *this; }; + operator X&() const { return value; } + pw_element_proxy& operator=(X x) { value = std::move(x); return *this; }; double lower_bound() const { return extent.first; } double upper_bound() const { return extent.second; } - const std::pair extent; + const std::pair& extent; X& value; }; +template +struct pw_element_proxy { + pw_element_proxy(const pw_elements& pw, pw_size_type i): + extent(pw.extent(i)), value(pw.value(i)) {} + + operator pw_element() const { return pw_element{extent, value}; } + operator const X&() const { return value; } + + double lower_bound() const { return extent.first; } + double upper_bound() const { return extent.second; } + + const std::pair& extent; + const X& value; +}; + + // Compute indices into vertex set corresponding to elements that cover a point x: namespace { @@ -243,10 +259,10 @@ struct pw_elements { using value_type = pw_element; using pointer = const pointer_proxy>; - using reference = pw_element; + using reference = pw_element_proxy; - reference operator[](difference_type j) const { return (*pw_)[j+*c_]; } - reference operator*() const { return (*pw_)[*c_]; } + reference operator[](difference_type j) const { return reference{*pw_, j + *c_}; } + reference operator*() const { return reference{*pw_, *c_}; } pointer operator->() const { return pointer{(*pw_)[*c_]}; } // (required for iterator_adaptor) From 37b4d8164c52cd43d36327ab0afa68e144f01626 Mon Sep 17 00:00:00 2001 From: Thorsten Hater <24411438+thorstenhater@users.noreply.github.com> Date: Tue, 5 Mar 2024 14:20:13 +0100 Subject: [PATCH 6/6] Only print profile if we can. --- example/busyring/ring.cpp | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/example/busyring/ring.cpp b/example/busyring/ring.cpp index c5e7213ee8..796d458221 100644 --- a/example/busyring/ring.cpp +++ b/example/busyring/ring.cpp @@ -282,7 +282,12 @@ int main(int argc, char** argv) { } auto report = arb::profile::make_meter_report(meters, context); + std::cout << report; +#ifdef ARB_PROFILE_ENABLED if (root) arb::profile::print_profiler_summary(std::cout, 5); +#endif + + } catch (std::exception& e) { std::cerr << "exception caught in ring miniapp: " << e.what() << "\n";