Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Reduce transient memory allocations in integrate/pw_lin. #2208

Closed
wants to merge 8 commits into from
+248 −119
Closed

Reduce transient memory allocations in integrate/pw_lin. #2208

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
55f93df
Use index/get instead of std::visit.
thorstenhater Aug 14, 2023
adae932
Snapshot, doesnt pass
thorstenhater Aug 14, 2023
4f9d5dd
Fix cases.
thorstenhater Aug 15, 2023
ec90987
Merge remote-tracking branch 'origin/master' into perf/reduce-pw-lin-…
thorstenhater Oct 17, 2023
78f821b
Complain -- ie terminate -- on unknown cases.
thorstenhater Oct 17, 2023
7ef78b6
Also add proxy for constant iterator, hoping we'll save the
thorstenhater Mar 5, 2024
1fecaa1
Merge remote-tracking branch 'origin/master' into perf/reduce-pw-lin-…
thorstenhater Mar 5, 2024
37b4d81
Only print profile if we can.
thorstenhater Mar 5, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
103 changes: 69 additions & 34 deletions arbor/cable_cell.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -91,7 +91,7 @@ struct cable_cell_impl {
dictionary(labels),
decorations(decorations)
{
init(decorations);
init();
}

cable_cell_impl(): cable_cell_impl({},{},{}) {}
Expand All @@ -100,7 +100,7 @@ struct cable_cell_impl {

cable_cell_impl(cable_cell_impl&& other) = default;

void init(const decor&);
void init();

template <typename T>
mlocation_map<T>& get_location_map(const T&) {
Expand All @@ -116,12 +116,13 @@ struct cable_cell_impl {
}

template <typename Item>
void place(const locset& ls, const Item& item, const hash_type& label) {

void place(const mlocation_list& ls, const Item& item, const hash_type& label) {
auto& mm = get_location_map(item);
cell_lid_type& lid = placed_count.get<Item>();
cell_lid_type first = lid;

for (auto l: thingify(ls, provider)) {
for (auto l: ls) {
placed<Item> p{l, lid++, item};
mm.push_back(p);
}
Expand Down Expand Up @@ -160,70 +161,104 @@ struct cable_cell_impl {
return region_map.get<init_reversal_potential>()[init.ion];
}

void paint(const region& reg, const density& prop) {
this->paint(reg, scaled_mechanism<density>(prop));
void paint(const mextent& cables, const density& prop, const std::string& reg) {
this->paint(cables, scaled_mechanism<density>(prop), reg);
}

void paint(const region& reg, const scaled_mechanism<density>& prop) {
void paint(const mextent& cables, const scaled_mechanism<density>& prop, const std::string& reg) {
std::unordered_map<std::string, iexpr_ptr> 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 <typename TaggedMech>
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<cable_cell_impl, void (*)(cable_cell_impl*)>;
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<init_membrane_potential>(what), rc.region); break; }
case 1: { paint(rc.cables, std::get<axial_resistivity>(what), rc.region); break; }
case 2: { paint(rc.cables, std::get<temperature>(what), rc.region); break; }
case 3: { paint(rc.cables, std::get<membrane_capacitance>(what), rc.region); break; }
case 4: { paint(rc.cables, std::get<ion_diffusivity>(what), rc.region); break; }
case 5: { paint(rc.cables, std::get<init_int_concentration>(what), rc.region); break; }
case 6: { paint(rc.cables, std::get<init_ext_concentration>(what), rc.region); break; }
case 7: { paint(rc.cables, std::get<init_reversal_potential>(what), rc.region); break; }
case 8: { paint(rc.cables, std::get<density>(what), rc.region); break; }
case 9: { paint(rc.cables, std::get<voltage_process>(what), rc.region); break; }
case 10: { paint(rc.cables, std::get<scaled_mechanism<density>>(what), rc.region); break; }
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

For new paintable types, adding it here this might get overlooked. std::visit would offer a way to throw a compile time error for such cases, but from the comments elsewhere, I'm assuming this was changed from std::visit for performance reasons here as well. Maybe at least add a default case that throws an error at runtime for unmatched types.

default: throw arbor_internal_error{"Unknown paintable variant"};
}
}
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<i_clamp>(what), label); break; }
case 1: { place(lc.places, std::get<threshold_detector>(what), label); break; }
case 2: { place(lc.places, std::get<synapse>(what), label); break; }
case 3: { place(lc.places, std::get<junction>(what), label); break; }
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same here as with the paintings case matching.

default: throw arbor_internal_error{"Unknown placeable variant"};
}
}
}

Expand Down
107 changes: 63 additions & 44 deletions arbor/cable_cell_param.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -129,50 +129,69 @@ 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<decltype(p)>;
if constexpr (std::is_same_v<init_membrane_potential, T>) {
if (p.scale.type() != iexpr_type::scalar) throw cable_cell_error{"Default values cannot have a scale."};
defaults_.init_membrane_potential = *p.scale.get_scalar()*p.value;
}
else if constexpr (std::is_same_v<axial_resistivity, T>) {
if (p.scale.type() != iexpr_type::scalar) throw cable_cell_error{"Default values cannot have a scale."};
defaults_.axial_resistivity = *p.scale.get_scalar()*p.value;
}
else if constexpr (std::is_same_v<temperature, T>) {
if (p.scale.type() != iexpr_type::scalar) throw cable_cell_error{"Default values cannot have a scale."};
defaults_.temperature_K = *p.scale.get_scalar()*p.value;
}
else if constexpr (std::is_same_v<membrane_capacitance, T>) {
if (p.scale.type() != iexpr_type::scalar) throw cable_cell_error{"Default values cannot have a scale."};
defaults_.membrane_capacitance = *p.scale.get_scalar()*p.value;
}
else if constexpr (std::is_same_v<init_int_concentration, T>) {
if (p.scale.type() != iexpr_type::scalar) throw cable_cell_error{"Default values cannot have a scale."};
defaults_.ion_data[p.ion].init_int_concentration = *p.scale.get_scalar()*p.value;
}
else if constexpr (std::is_same_v<init_ext_concentration, T>) {
if (p.scale.type() != iexpr_type::scalar) throw cable_cell_error{"Default values cannot have a scale."};
defaults_.ion_data[p.ion].init_ext_concentration = *p.scale.get_scalar()*p.value;
}
else if constexpr (std::is_same_v<init_reversal_potential, T>) {
if (p.scale.type() != iexpr_type::scalar) throw cable_cell_error{"Default values cannot have a scale."};
defaults_.ion_data[p.ion].init_reversal_potential = *p.scale.get_scalar()*p.value;
}
else if constexpr (std::is_same_v<ion_reversal_potential_method, T>) {
defaults_.reversal_potential_method[p.ion] = p.method;
}
else if constexpr (std::is_same_v<cv_policy, T>) {
defaults_.discretization = std::forward<cv_policy>(p);
}
else if constexpr (std::is_same_v<ion_diffusivity, T>) {
if (p.scale.type() != iexpr_type::scalar) throw cable_cell_error{"Default values cannot have a scale."};
auto s = p.scale.get_scalar();
defaults_.ion_data[p.ion].diffusivity = s ? std::optional{*s*p.value} : s;
}
},
what);
// NOTE: the index/get approach is considerably faster than std::visit.
switch (what.index()) {
case 0: {
auto& p = std::get<init_membrane_potential>(what);
if (p.scale.type() != iexpr_type::scalar) throw cable_cell_error{"Default values cannot have a scale."};
defaults_.init_membrane_potential = *p.scale.get_scalar()*p.value;
break;
}
case 1: {
auto& p = std::get<axial_resistivity>(what);
if (p.scale.type() != iexpr_type::scalar) throw cable_cell_error{"Default values cannot have a scale."};
defaults_.axial_resistivity = *p.scale.get_scalar()*p.value;
break;
}
case 2: {
auto& p = std::get<temperature>(what);
if (p.scale.type() != iexpr_type::scalar) throw cable_cell_error{"Default values cannot have a scale."};
defaults_.temperature_K = *p.scale.get_scalar()*p.value;
break;
}
case 3: {
auto& p = std::get<membrane_capacitance>(what);
if (p.scale.type() != iexpr_type::scalar) throw cable_cell_error{"Default values cannot have a scale."};
defaults_.membrane_capacitance = *p.scale.get_scalar()*p.value;
break;
}
case 4: {
auto& p = std::get<ion_diffusivity>(what);
if (p.scale.type() != iexpr_type::scalar) throw cable_cell_error{"Default values cannot have a scale."};
auto s = p.scale.get_scalar();
defaults_.ion_data[p.ion].diffusivity = s ? std::optional{*s*p.value} : s;
break;
}
case 5: {
auto& p = std::get<init_int_concentration>(what);
if (p.scale.type() != iexpr_type::scalar) throw cable_cell_error{"Default values cannot have a scale."};
defaults_.ion_data[p.ion].init_int_concentration = *p.scale.get_scalar()*p.value;
break;
}
case 6: {
auto& p = std::get<init_ext_concentration>(what);
if (p.scale.type() != iexpr_type::scalar) throw cable_cell_error{"Default values cannot have a scale."};
defaults_.ion_data[p.ion].init_ext_concentration = *p.scale.get_scalar()*p.value;
break;
}
case 7: {
auto& p = std::get<init_reversal_potential>(what);
if (p.scale.type() != iexpr_type::scalar) throw cable_cell_error{"Default values cannot have a scale."};
defaults_.ion_data[p.ion].init_reversal_potential = *p.scale.get_scalar()*p.value;
break;
}
case 8: {
auto& p = std::get<ion_reversal_potential_method>(what);
defaults_.reversal_potential_method[p.ion] = p.method;
break;
}
case 9: {
defaults_.discretization = std::forward<cv_policy>(std::get<cv_policy>(what));
break;
}
default:
throw arbor_internal_error{"Unknown defaultable variant"};
}
return *this;
}

Expand Down
3 changes: 1 addition & 2 deletions arbor/fvm_layout.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -305,7 +305,7 @@ fvm_cv_discretize(const cable_cell& cell, const cable_cell_parameter_set& global
[](const auto& kv) {
const auto& [k, v] = kv;
auto s = v.scale.get_scalar();
return !s || *s*v.value != 0.0;
return !s || (*s != 0.0 && std::isfinite(*s));
});
if (diffusive) {
// Provide a (non-sensical) default.
Expand Down Expand Up @@ -1198,7 +1198,6 @@ make_density_mechanism_config(const region_assignment<density>& 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();

Expand Down
4 changes: 2 additions & 2 deletions arbor/include/arbor/morph/embed_pwlin.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down
Loading