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
Closed

Reduce transient memory allocations in integrate/pw_lin. #2208

Show file tree
Hide file tree
Changes from 3 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
99 changes: 66 additions & 33 deletions arbor/cable_cell.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -95,7 +95,7 @@ struct cable_cell_impl {
dictionary(labels),
decorations(decorations)
{
init(decorations);
init();
}

cable_cell_impl(): cable_cell_impl({},{},{}) {}
Expand All @@ -104,7 +104,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 @@ -120,12 +120,12 @@ struct cable_cell_impl {
}

template <typename Item>
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<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 @@ -164,69 +164,102 @@ 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_K>(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.

}
}
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.

}
}
}

Expand Down
104 changes: 61 additions & 43 deletions arbor/cable_cell_param.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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<decltype(p)>;
if constexpr (std::is_same_v<init_membrane_potential, T>) {
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<axial_resistivity, T>) {
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<temperature_K, T>) {
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<membrane_capacitance, T>) {
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<init_int_concentration, T>) {
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<init_ext_concentration, T>) {
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<init_reversal_potential, T>) {
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<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.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<init_membrane_potential>(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<axial_resistivity>(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<temperature_K>(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<membrane_capacitance>(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<ion_diffusivity>(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<init_int_concentration>(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 6: {
auto& p = std::get<init_ext_concentration>(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 7: {
auto& p = std::get<init_reversal_potential>(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<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 @@ -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.
Expand Down Expand Up @@ -1204,7 +1204,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
52 changes: 44 additions & 8 deletions arbor/morph/embed_pwlin.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -103,8 +103,8 @@ struct embed_pwlin_data {

template <unsigned p, unsigned q>
double interpolate(double pos, const pw_ratpoly<p, q>& 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));
}
Expand All @@ -118,7 +118,7 @@ double interpolate(double pos, const pw_ratpoly<p, q>& f) {
template <unsigned p, unsigned q>
double integrate(const pw_constant_fn& g, const pw_ratpoly<p, q>& 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;
Expand All @@ -133,7 +133,7 @@ template <unsigned p, unsigned q>
double integrate(const pw_constant_fn& g, const pw_elements<pw_ratpoly<p, q>>& 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;

Expand Down Expand Up @@ -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<void>{{left, right}});
return pw_zip_with(g, pw_elements<void>{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<void>{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<void>{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<void>{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:
Expand Down
Loading