Snapshot, doesnt pass

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 <typename T>
  mlocation_map<T>& get_location_map(const T&) {
@@ -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);
  }
@@ -164,74 +164,107 @@ 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; }
+ }
  }
- 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; }
+ }
  }
 }
 
 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())) {}

arbor/fvm_layout.cpp

@@ -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();
 

arbor/morph/embed_pwlin.cpp

@@ -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));
 }
@@ -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;
@@ -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;
 
@@ -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:

arbor/util/piecewise.hpp

@@ -228,6 +228,9 @@ struct pw_elements {
  counter<pw_size_type>& inner() { return c_; }
  const counter<pw_size_type>& 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<X>* pw_;
  counter<pw_size_type> c_;
@@ -250,6 +253,11 @@ struct pw_elements {
  counter<pw_size_type>& inner() { return c_; }
  const counter<pw_size_type>& 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<X>* pw_;
  counter<pw_size_type> 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<X>{*this, i}; }
  auto operator[](size_type i) const & { return value_type{extent(i), value(i)}; }
@@ -474,6 +483,9 @@ struct pw_elements<void> {
  counter<pw_size_type>& inner() { return c_; }
  const counter<pw_size_type>& 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<void>* pw_;
  counter<pw_size_type> c_;
@@ -486,12 +498,10 @@ struct pw_elements<void> {
  pw_elements() = default;
 
  template <typename VertexSeq>
- pw_elements(const VertexSeq& vs) {
- assign(vs);
- }
-
- pw_elements(std::initializer_list<double> vs) {
- assign(vs);
+ pw_elements(const VertexSeq& vs) { assign(vs); }
+ pw_elements(std::initializer_list<double> 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<void> {
  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<double, double> 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 <typename F>
+ 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 <typename F>
+ 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};
  }