From 9b35ba5266975cf951890492e086fbf53da84dfe Mon Sep 17 00:00:00 2001
From: Sam Yates <halfflat@gmail.com>
Date: Tue, 9 Nov 2021 21:55:34 +0800
Subject: [PATCH] Deal with zero radius points in a morphology (#1719)
* Add test that is sensitive to small radius loss of precision in ixa computation.
* Add test that checks for consistent behaviour when there is an (isolated) zero radius point in the morphology.
* Special case p==1 q==1 case in ratelem so that it can interpolate ixa in the presence of non-finite interpolants.
* Change naming in pw_element/pw_elements: 'element' refers to the extent+value pairs comprising a 'pw_elements' object; 'value' refers to the value associated with an element; 'extent' refers to the closed interval which is the support of the element.
* Allow values (but not extents) in a pw_elements object to be mutable, using proxies for iterator access.
* Write piecewise zip (pw_zip) in terms of a lazy pw_zip_view. Rename zip functions.
* Document pw_elements more thoroughly.
* Simplify and document embed_pwlin.cpp routines, expressing everything in terms of simple interpolate/integrate operations.
* Represent integrated inverse cross-sectional area in the embedding by multiple piecewise-rational functions over a branch, each of which contribute separately to an ixa-based integration, so that precision loss associated with small (or zero) radii can be avoided.
* Comment embed_pwlin.cpp more thoroughly.
Fixes #1526
---
arbor/fvm_layout.cpp | 6 +-
arbor/fvm_layout.hpp | 8 +-
arbor/morph/embed_pwlin.cpp | 342 +++++++-----
arbor/morph/place_pwlin.cpp | 7 +-
arbor/util/iterutil.hpp | 7 +-
arbor/util/piecewise.hpp | 819 +++++++++++++++++++++--------
arbor/util/ratelem.hpp | 50 +-
test/unit/test_morph_embedding.cpp | 70 +++
test/unit/test_piecewise.cpp | 220 +++++---
test/unit/test_ratelem.cpp | 32 +-
10 files changed, 1089 insertions(+), 472 deletions(-)
diff --git a/arbor/fvm_layout.cpp b/arbor/fvm_layout.cpp
index 438f3313..d926192c 100644
--- a/arbor/fvm_layout.cpp
+++ b/arbor/fvm_layout.cpp
@@ -118,7 +118,7 @@ pw_elements<U> pw_over_cable(const mcable_map<T>& mm, mcable cable, U dflt_value
}
if (cable.prox_pos!=0 || cable.dist_pos!=1) {
- pw = zip(pw, pw_elements<>({cable.prox_pos, cable.dist_pos}));
+ pw = pw_zip_with(pw, pw_elements<>({cable.prox_pos, cable.dist_pos}));
}
return pw;
}
@@ -1315,8 +1315,8 @@ fvm_mechanism_data fvm_build_mechanism_data(
const mcable_map<init_ext_concentration>& econc_on_cable = initial_econc_map[ion];
const mcable_map<init_reversal_potential>& rvpot_on_cable = initial_rvpot_map[ion];
- auto pw_times = [](const pw_elements<double>& a, const pw_elements<double>& b) {
- return zip(a, b, [](double left, double right, pw_element<double> a, pw_element<double> b) { return a.element*b.element; });
+ auto pw_times = [](const pw_elements<double>& pwa, const pw_elements<double>& pwb) {
+ return pw_zip_with(pwa, pwb, [](std::pair<double, double>, double a, double b) { return a*b; });
};
for (auto i: count_along(config.cv)) {
diff --git a/arbor/fvm_layout.hpp b/arbor/fvm_layout.hpp
index dc9ffa82..8d1356fb 100644
--- a/arbor/fvm_layout.hpp
+++ b/arbor/fvm_layout.hpp
@@ -116,15 +116,15 @@ struct cv_geometry {
size_type location_cv(size_type cell_idx, mlocation loc, cv_prefer::type prefer) const {
auto& pw_cv_offset = branch_cv_map.at(cell_idx).at(loc.branch);
auto zero_extent = [&pw_cv_offset](auto j) {
- return pw_cv_offset.interval(j).first==pw_cv_offset.interval(j).second;
+ return pw_cv_offset.extent(j).first==pw_cv_offset.extent(j).second;
};
auto i = pw_cv_offset.index_of(loc.pos);
auto i_max = pw_cv_offset.size()-1;
- auto cv_prox = pw_cv_offset.interval(i).first;
+ auto cv_prox = pw_cv_offset.extent(i).first;
// index_of() should have returned right-most matching interval.
- arb_assert(i==i_max || loc.pos<pw_cv_offset.interval(i+1).first);
+ arb_assert(i==i_max || loc.pos<pw_cv_offset.extent(i+1).first);
using namespace cv_prefer;
switch (prefer) {
@@ -145,7 +145,7 @@ struct cv_geometry {
}
index_type cv_base = cell_cv_divs.at(cell_idx);
- return cv_base+pw_cv_offset[i].element;
+ return cv_base+pw_cv_offset.value(i);
}
};
diff --git a/arbor/morph/embed_pwlin.cpp b/arbor/morph/embed_pwlin.cpp
index 03393acc..dfe84c73 100644
--- a/arbor/morph/embed_pwlin.cpp
+++ b/arbor/morph/embed_pwlin.cpp
@@ -14,118 +14,72 @@
namespace arb {
+// Each integrable or interpolated quantity is represented by
+// a sequence of spans covering each branch, where a span is itself
+// a piecewise polynomial or rational function.
+//
+// The piecewise functions are represented by util::pw_elements<util::rat_element<p, q>>
+// objects. util::rat_element describes an order (p, q) rational function in terms of
+// the values of that function at p+q+1 equally spaced points along the element, including
+// the two endpoints.
+
using util::rat_element;
+using util::pw_elements;
-template <unsigned p, unsigned q>
-using pw_ratpoly = util::pw_elements<rat_element<p, q>>;
+// Represents piecewise rational function over a subset of a branch.
+// When the function represents an integral, the integral between two points
+// in the domain of a pw_ratpoly can be computed by taking the difference
+// in the interpolated value at the two points.
template <unsigned p, unsigned q>
-using branch_pw_ratpoly = std::vector<pw_ratpoly<p, q>>;
+using pw_ratpoly = util::pw_elements<rat_element<p, q>>;
-// Special case handling required for degenerate branches of length zero:
-template <typename Elem>
-static bool is_degenerate(const util::pw_elements<Elem>& pw) {
- return pw.bounds().second==0;
-}
+// One branch can be covered by more than one pw_ratpoly (of the same degree).
+// When the function represents an integral, the integral between two points
+// that do not lie in the same domain of a pw_ratpoly must be summed by
+// considering the contribution of each pw_ratpoly element separately. Multiple
+// pw_ratpoly elements over the same branch are required to avoid cases of loss
+// of precision and singularities.
template <unsigned p, unsigned q>
-double interpolate(const branch_pw_ratpoly<p, q>& f, unsigned bid, double pos) {
- const auto& pw = f.at(bid);
- if (is_degenerate(pw)) pos = 0;
+using branch_pw_spans = util::pw_elements<pw_ratpoly<p, q>>;
- auto [bounds, element] = pw(pos);
+struct embed_pwlin_data {
+ // Vectors below are all indexed by branch index.
- if (bounds.first==bounds.second) return element[0];
- else {
- double x = (pos-bounds.first)/(bounds.second-bounds.first);
- return element(x);
- }
-}
+ // Length and directed projection are piecewise linear, continuous, and
+ // monotonically increasing. They are represented by a single piecewise
+ // linear function per branch.
-// Length, area, and ixa are polynomial or rational polynomial functions of branch position,
-// continuos and monotonically increasing with respect to distance from root.
-//
-// Integration wrt a piecewise constant function is performed by taking the difference between
-// interpolated values at the end points of each constant interval.
+ std::vector<pw_ratpoly<1, 0>> length; // [µm]
+ std::vector<pw_ratpoly<1, 0>> directed_projection; // [µm]
-template <unsigned p, unsigned q>
-double integrate(const branch_pw_ratpoly<p, q>& f, unsigned bid, const pw_constant_fn& g) {
- double accum = 0;
- for (msize_t i = 0; i<g.size(); ++i) {
- std::pair<double, double> interval = g.interval(i);
- accum += g.element(i)*(interpolate(f, bid, interval.second)-interpolate(f, bid, interval.first));
- }
- return accum;
-}
+ // Radius is piecewise linear, but not necessarily continuous. Where the
+ // morphology describes different radii at the same point, the interpolated
+ // radius is regarded as being right-continuous, and corresponds to the
+ // value described by the last segment added that covers the point.
-// Performance note: when integrating over a cable within a branch, the code effectively
-// performs a linear search for the starting interval. This can be replaced with a binary
-// search for a small increase in code complexity.
+ std::vector<pw_ratpoly<1, 0>> radius; // [µm]
-template <unsigned p, unsigned q>
-double integrate(const branch_pw_ratpoly<p, q>& f, mcable c, const pw_constant_fn& g) {
- msize_t bid = c.branch;
- double accum = 0;
+ // Morphological surface area between points on a branch is given as
+ // the difference between the interpolated value of a piecewise quadratic
+ // function over the branch. The function is monotonically increasing,
+ // but not necessarily continuous. Where the value jumps, the function
+ // is interpreted as being right-continuous.
- for (msize_t i = 0; i<g.size(); ++i) {
- std::pair<double, double> interval = g.interval(i);
+ std::vector<pw_ratpoly<2, 0>> area; // [µm²]
- if (interval.second<c.prox_pos) {
- continue;
- }
- else if (interval.first>=c.dist_pos) {
- break;
- }
- else {
- interval.first = std::max(interval.first, c.prox_pos);
- interval.second = std::min(interval.second, c.dist_pos);
-
- if (interval.first<interval.second) {
- accum += g.element(i)*(interpolate(f, bid, interval.second)-interpolate(f, bid, interval.first));
- }
- }
- }
- return accum;
-}
+ // Integrated inverse cross-sectional area (ixa) between points on a branch
+ // is used to compute the conductance between points on the morphology.
+ // It is represented by an order (1, 1) rational function (i.e. a linear
+ // function divided by a linear function).
+ //
+ // Because a small radius can lead to very large (or infinite, if the
+ // radius is zero) values, one branch may be covered by multiple
+ // piecewise functions, and computing ixa may require summing the
+ // contributions from more than one.
-template <typename operation>
-mcable_list data_cmp(const branch_pw_ratpoly<1, 0>& f, unsigned bid, double val, operation op) {
- mcable_list L;
- const auto& pw = f.at(bid);
- for (const auto& piece: pw) {
- auto extents = piece.interval;
- auto left_val = piece.element(0);
- auto right_val = piece.element(1);
-
- if (!op(left_val, val) && !op(right_val, val)) {
- continue;
- }
- if (op(left_val, val) && op(right_val, val)) {
- L.push_back({bid, extents.first, extents.second});
- continue;
- }
-
- auto cable_loc = (val - left_val)/(right_val - left_val);
- auto edge = math::lerp(extents.first, extents.second, cable_loc);
-
- if (op(left_val, val)) {
- L.push_back({bid, extents.first, edge});
- continue;
- }
- if (!op(left_val, val)) {
- L.push_back({bid, edge, extents.second});
- continue;
- }
- }
- return L;
-}
-
-struct embed_pwlin_data {
- branch_pw_ratpoly<1, 0> length; // [µm]
- branch_pw_ratpoly<1, 0> directed_projection; // [µm]
- branch_pw_ratpoly<1, 0> radius; // [µm]
- branch_pw_ratpoly<2, 0> area; // [µm²]
- branch_pw_ratpoly<1, 1> ixa; // [1/µm]
+ std::vector<branch_pw_spans<1, 1>> ixa; // [1/µm]
explicit embed_pwlin_data(msize_t n_branch):
length(n_branch),
@@ -136,24 +90,88 @@ struct embed_pwlin_data {
{}
};
+// Interpolation
+
+// Value at pos on a branch interpolated from a piecewise rational function
+// which covers pos. pos is in [0, 1], and the bounds of the piecewise rational
+// function is an interval [a, b] with 0 ≤ a ≤ pos ≤ b ≤ 1.
+
+template <unsigned p, unsigned q>
+double interpolate(double pos, const pw_ratpoly<p, q>& f) {
+ auto [extent, poly] = f(pos);
+ auto [left, right] = extent;
+
+ return left==right? poly[0]: poly((pos-left)/(right-left));
+}
+
+// Integration
+
+// Integral of piecwise constant function g with respect to a measure determined
+// by a single piecewise monotonic right continuous rational function f, over
+// an interval [r, s]: ∫[r,s] g(x) df(x) (where [r, s] is the domain of g).
+
+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) {
+ sum += gval*(interpolate(extent.second, f)-interpolate(extent.first, f));
+ }
+ return sum;
+}
+
+// Integral of piecewise constant function g with respect to a measure determined
+// by a contiguous sequence of piecewise monotonic right continuous rational
+// functions fi with support [aⱼ, bⱼ], over an interval [r, s] with a₀ ≤ r ≤ s ≤ bₙ:
+// Σⱼ ∫[r,s]∩[aⱼ,bⱼ] g(x)dfⱼ(x) (where [r, s] is the domain of g).
+
+template <unsigned p, unsigned q>
+double integrate(const pw_constant_fn& g, const pw_elements<pw_ratpoly<p, q>>& fs) {
+ double sum = 0;
+ 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<p, q> f = pw_pair.second;
+ sum += gval*(interpolate(right, f)-interpolate(left, f));
+ }
+ return sum;
+}
+
+// Implementation of public embed_pwlin methods:
+
double embed_pwlin::radius(mlocation loc) const {
- return interpolate(data_->radius, loc.branch, loc.pos);
+ return interpolate(loc.pos, data_->radius.at(loc.branch));
}
double embed_pwlin::directed_projection(arb::mlocation loc) const {
- return interpolate(data_->directed_projection, loc.branch, loc.pos);
+ return interpolate(loc.pos, data_->directed_projection.at(loc.branch));
}
// Point to point integration:
double embed_pwlin::integrate_length(mlocation proximal, mlocation distal) const {
- return interpolate(data_->length, distal.branch, distal.pos) -
- interpolate(data_->length, proximal.branch, proximal.pos);
+ return interpolate(distal.pos, data_->length.at(distal.branch)) -
+ interpolate(proximal.pos, data_->length.at(proximal.branch));
}
double embed_pwlin::integrate_area(mlocation proximal, mlocation distal) const {
- return interpolate(data_->area, distal.branch, distal.pos) -
- interpolate(data_->area, proximal.branch, proximal.pos);
+ return interpolate(distal.pos, data_->area.at(distal.branch)) -
+ interpolate(proximal.pos, data_->area.at(proximal.branch));
+}
+
+// Integrate piecewise function over a branch:
+
+double embed_pwlin::integrate_length(msize_t bid, const pw_constant_fn& g) const {
+ return integrate(g, data_->length.at(bid));
+}
+
+double embed_pwlin::integrate_area(msize_t bid, const pw_constant_fn& g) const {
+ return integrate(g, data_->area.at(bid));
+}
+
+double embed_pwlin::integrate_ixa(msize_t bid, const pw_constant_fn& g) const {
+ return integrate(g, data_->ixa.at(bid));
}
// Integrate over cable:
@@ -170,36 +188,65 @@ double embed_pwlin::integrate_ixa(mcable c) const {
return integrate_ixa(c.branch, pw_constant_fn{{c.prox_pos, c.dist_pos}, {1.}});
}
-// Integrate piecewise function over a branch:
-
-double embed_pwlin::integrate_length(msize_t bid, const pw_constant_fn& g) const {
- return integrate(data_->length, bid, g);
-}
-
-double embed_pwlin::integrate_area(msize_t bid, const pw_constant_fn& g) const {
- return integrate(data_->area, bid, g);
-}
+// Integrate piecewise function over a cable:
-double embed_pwlin::integrate_ixa(msize_t bid, const pw_constant_fn& g) const {
- return integrate(data_->ixa, bid, g);
+static pw_constant_fn restrict(const pw_constant_fn& g, double left, double right) {
+ return pw_zip_with(g, pw_elements<void>{{left, right}});
}
-// Integrate piecewise function over a cable:
-
double embed_pwlin::integrate_length(mcable c, const pw_constant_fn& g) const {
- return integrate(data_->length, c, g);
+ return integrate_length(c.branch, restrict(g, c.prox_pos, c.dist_pos));
}
double embed_pwlin::integrate_area(mcable c, const pw_constant_fn& g) const {
- return integrate(data_->area, c, g);
+ return integrate_area(c.branch, restrict(g, c.prox_pos, c.dist_pos));
}
double embed_pwlin::integrate_ixa(mcable c, const pw_constant_fn& g) const {
- return integrate(data_->ixa, c, g);
+ return integrate_ixa(c.branch, restrict(g, c.prox_pos, c.dist_pos));
}
// Subregions defined by geometric inequalities:
+// Given a piecewise linear function f over a branch, a comparison operation op and a threshold v,
+// determine the subset of the branch where the op(f(x), v) is true.
+//
+// Functions are supplied for each branch via a branch-index vector of pw_ratpoly<1, 0> functions;
+// supplied branch id is the index into this vector, and is used to construct the cables
+// corresponding to the matching subset.
+
+template <typename operation>
+mcable_list data_cmp(const std::vector<pw_ratpoly<1, 0>>& f_on_branch, msize_t bid, double val, operation op) {
+ mcable_list L;
+ for (auto&& piece: f_on_branch.at(bid)) {
+ auto [left, right] = piece.extent;
+ auto left_val = piece.value(0);
+ auto right_val = piece.value(1);
+
+ if (!op(left_val, val) && !op(right_val, val)) {
+ continue;
+ }
+ if (op(left_val, val) && op(right_val, val)) {
+ L.push_back({bid, left, right});
+ continue;
+ }
+
+ auto cable_loc = (val - left_val)/(right_val - left_val);
+ auto edge = math::lerp(left, right, cable_loc);
+
+ if (op(left_val, val)) {
+ L.push_back({bid, left, edge});
+ continue;
+ }
+ if (!op(left_val, val)) {
+ L.push_back({bid, edge, right});
+ continue;
+ }
+ }
+ return L;
+}
+
+
mcable_list embed_pwlin::radius_cmp(msize_t bid, double val, comp_op op) const {
switch (op) {
case comp_op::lt: return data_cmp(data_->radius, bid, val, [](auto l, auto r){return l < r;});
@@ -275,11 +322,13 @@ embed_pwlin::embed_pwlin(const arb::morphology& m) {
segment_cables_[seg.id] = mcable{bid, pos0, pos1};
}
- double length_0 = parent==mnpos? 0: data_->length[parent].back().element[1];
+ double length_0 = parent==mnpos? 0: data_->length[parent].back().value[1];
data_->length[bid].push_back(0., 1, rat_element<1, 0>(length_0, length_0+branch_length));
- double area_0 = parent==mnpos? 0: data_->area[parent].back().element[2];
- double ixa_0 = parent==mnpos? 0: data_->ixa[parent].back().element[2];
+ double area_0 = parent==mnpos? 0: data_->area[parent].back().value[2];
+
+ double ixa_last = 0;
+ pw_ratpoly<1, 1> ixa_pw;
for (auto i: util::count_along(segments)) {
auto prox = segments[i].prox;
@@ -304,11 +353,52 @@ embed_pwlin::embed_pwlin(const arb::morphology& m) {
data_->area[bid].push_back(p0, p1, rat_element<2, 0>(area_0, area_half, area_1));
area_0 = area_1;
- // (Test for positive dx explicitly in case r0 is zero.)
- double ixa_half = ixa_0 + (dx>0? dx/(pi*r0*(r0+r1)): 0);
- double ixa_1 = ixa_0 + (dx>0? dx/(pi*r0*r1): 0);
- data_->ixa[bid].push_back(p0, p1, rat_element<1, 1>(ixa_0, ixa_half, ixa_1));
- ixa_0 = ixa_1;
+ if (dx>0) {
+ double ixa_0 = 0;
+ double ixa_half = dx/(pi*r0*(r0+r1));
+ double ixa_1 = dx/(pi*r0*r1);
+
+ if (r0==0 && r1==0) {
+ // ixa is just not defined on this segment; represent with all
+ // infinite node values.
+ ixa_0 = INFINITY;
+ ixa_half = INFINITY;
+ ixa_1 = INFINITY;
+ }
+ else if (r0==0) {
+ ixa_0 = -INFINITY;
+ ixa_half = -dx/(pi*r1*r1);
+ ixa_1 = 0;
+ }
+
+ // Start a new piecewise function if last ixa value is
+ // large compared to ixa_1 - ixa_0 (leading to loss
+ // of precision), or if ixa_0 is already non-zero (something
+ // above has declared we should start anew).
+
+ constexpr double max_ixa_ratio = 1e5;
+ if (!ixa_pw.empty() && ixa_0 == 0 && ixa_last<max_ixa_ratio*ixa_1) {
+ // Extend last pw representation on the branch.
+ ixa_0 += ixa_last;
+ ixa_half += ixa_last;
+ ixa_1 += ixa_last;
+ }
+ else {
+ // Start a new pw representation on the branch:
+ if (!ixa_pw.empty()) {
+ auto [left, right] = ixa_pw.bounds();
+ data_->ixa[bid].push_back(left, right, ixa_pw);
+ }
+ ixa_pw.clear();
+ }
+ ixa_pw.push_back(p0, p1, rat_element<1, 1>(ixa_0, ixa_half, ixa_1));
+ ixa_last = ixa_1;
+ }
+ }
+ // push last ixa pw function on the branch, if nonempty.
+ if (!ixa_pw.empty()) {
+ auto [left, right] = ixa_pw.bounds();
+ data_->ixa[bid].push_back(left, right, ixa_pw);
}
arb_assert((data_->radius[bid].size()>0));
diff --git a/arbor/morph/place_pwlin.cpp b/arbor/morph/place_pwlin.cpp
index 8bec27ca..9814aeb2 100644
--- a/arbor/morph/place_pwlin.cpp
+++ b/arbor/morph/place_pwlin.cpp
@@ -52,8 +52,8 @@ mpoint place_pwlin::at(mlocation loc) const {
const auto& pw_index = data_->segment_index.at(loc.branch);
double pos = is_degenerate(pw_index)? 0: loc.pos;
- auto [bounds, index] = pw_index(pos);
- return interpolate_segment(bounds, data_->segments.at(index), pos);
+ auto index = pw_index(pos);
+ return interpolate_segment(index.extent, data_->segments.at(index), pos);
}
std::vector<mpoint> place_pwlin::all_at(mlocation loc) const {
@@ -61,7 +61,8 @@ std::vector<mpoint> place_pwlin::all_at(mlocation loc) const {
const auto& pw_index = data_->segment_index.at(loc.branch);
double pos = is_degenerate(pw_index)? 0: loc.pos;
- for (auto [bounds, index]: util::make_range(pw_index.equal_range(pos))) {
+ for (auto index: util::make_range(pw_index.equal_range(pos))) {
+ auto bounds = index.extent;
auto seg = data_->segments.at(index);
// Add both ends of zero length segment, if they differ.
diff --git a/arbor/util/iterutil.hpp b/arbor/util/iterutil.hpp
index 1556a50f..14b6cde0 100644
--- a/arbor/util/iterutil.hpp
+++ b/arbor/util/iterutil.hpp
@@ -91,8 +91,11 @@ decltype(auto) back(Seq& seq) {
template <typename V>
struct pointer_proxy {
V v;
- pointer_proxy(const V& v): v(v) {}
- pointer_proxy(V&& v): v(std::move(v)) {}
+
+ template <typename... Args>
+ pointer_proxy(Args&&... args): v(std::forward<Args>(args)...) {}
+
+ V* operator->() { return &v; }
const V* operator->() const { return &v; }
};
diff --git a/arbor/util/piecewise.hpp b/arbor/util/piecewise.hpp
index f250c1b7..8f029b4f 100644
--- a/arbor/util/piecewise.hpp
+++ b/arbor/util/piecewise.hpp
@@ -1,15 +1,105 @@
#pragma once
-// Create/manipulate 1-d piece-wise defined objects.
+// Create/manipulate 1-d piecewise defined objects.
//
-// Using vectors everywhere here for ease; consider making
-// something more container/sequence-generic later.
+// A `pw_element<A>` describes a _value_ of type `A` and an _extent_ of
+// type std::pair<double, double>. An extent with value `{l, r}` has
+// l ≤ r and represents a closed interval [l, r] of the real line.
+//
+// For void `A`, `pw_element<void>` holds only the extent of an element.
+//
+// Once constructed, the value of a `pw_element<A>` (with `A` not void) is
+// mutable, and can be assigned directly with `operator=`, but the extent is
+// constant. A `pw_element<A>` can be implicitly converted to a value of type
+// `A`.
+//
+// A `pw_element<A>` can be bound via a structured binding to an extent
+// and value, or just an extent if A is void.
+//
+// A `pw_elements<A>` object comprises a contiguous sequence Eáµ¢ of
+// `pw_element<A>` elements. If Eᵢ has extent [lᵢ, rᵢ] and Eᵢ₊₠has extent
+// [lᵢ₊â‚, rᵢ₊â‚] then ráµ¢ must equal lᵢ₊â‚.
+//
+// The `value_type` of `pw_elements<A>` is the type of the elements, i.e.
+// `pw_element<A>`, as `pw_elements<A>` presents as a container.
+// To avoid ambiguity the type `A` is termed the _codomain_ of a `pw_elements<A>`
+// object.
+//
+// When the codomain is `void`, which is also the default for the type
+// parameter `A`, `pw_element<>` contains only the extent, and `pw_elements<>`
+// does not hold any values for the elements.
+//
+//
+// Construction:
+//
+// A `pw_elements<A>` object can be constructed from a sequence of _vertices_
+// and a sequence of _values_ (if `A` is not `void`). A vertex sequence of
+// doubles xâ‚€, xâ‚, …, xâ‚™ and value sequence vâ‚€, vâ‚, …, vₙ₋₠defines n elements
+// with extents [xâ‚€, xâ‚], [xâ‚, xâ‚‚], etc.
+//
+// A default constructed `pw_elements<A>` has no elements.
+//
+// A `pw_element<A>` can be appended to a `pw_elements<A>` object with:
+//
+// pw_elements<A>::push_back(pw_element<A>)
+//
+// The lower bound of the element must equal the upper bound of the existing
+// sequence, if it is not empty.
+//
+// The extent and value can be given explicitly with:
+//
+// template <typename U&&>
+// pw_elements<A>::push_back(double left, double right, U&& value)
+//
+// template <typename U&&>
+// pw_elements<A>::push_back(double right, U&& value)
+//
+// where the lower bound of the extent of the new element can be omitted
+// if the sequence is non-empty. When `A` is void, these methods are
+// instead:
+//
+// pw_elements<>::push_back(double left, double right)
+// pw_elements<>::push_back(double right)
+//
+//
+// Conversion:
+//
+// A `pw_elements<A>` object can be explicitly constructed from a `pw_elements<B>`
+// object if the values of type B can be converted to values of type A.
+//
+// A `pw_elements<void>` object can be explicitly constucted from any
+// `pw_elements<B>` object, where it takes just the element extents.
+//
+//
+// Element access:
+//
+// Elements are retrievable by index by `operator[]`. The value of the ith
+// element is given by `value(i)`, and its extent as a pair of doubles by
+// `extent(i)`.
+//
+// The method `vertices()` returns the vector of nodes defining the element
+// extents, while for non-void codomains, `values()` returns the vector of
+// element values.
+//
+// Elements can be looked up by their position with `operator()(double x)` and
+// `equal_range(double x)`. The element (or element proxy, see below) retrieved
+// for `x` is the right-most element whose extent includes `x`; this makes
+// `pw_elements<A>` act like a right-continuous function of position.
+// `equal_range(x)` returns a pair of iterators that correspond to the sequence
+// of elements whose extents all include `x`.
+//
+// Elements that are returned by a non-const lvalue `pw_elements<A>` via
+// `operator[]`, `operator()`, or via an iterator are represented by a proxy
+// object that allows in-place mutation of the element's value.
+#include <cmath>
#include <initializer_list>
#include <iterator>
#include <type_traits>
#include <vector>
+#include "util/iterutil.hpp"
+#include "util/transform.hpp"
#include "util/meta.hpp"
#include "util/partition.hpp"
@@ -19,83 +109,161 @@ namespace util {
using pw_size_type = unsigned;
constexpr pw_size_type pw_npos = -1;
-// Generic random access const iterator for a collection
-// providing operator[].
+template <typename X = void>
+struct pw_element {
+ pw_element():
+ extent(NAN, NAN),
+ value()
+ {}
-template <typename T>
-struct indexed_const_iterator {
- using size_type = decltype(std::size(std::declval<T>()));
- using difference_type = std::make_signed_t<size_type>;
+ pw_element(std::pair<double, double> extent, X value):
+ extent(std::move(extent)),
+ value(std::move(value))
+ {}
- using value_type = decltype(std::declval<T>()[0]);
- struct pointer {
- value_type v;
- const value_type* operator->() const { return &v; }
- };
+ pw_element(const pw_element&) = default;
+ pw_element(pw_element&&) = default;
- using reference = value_type;
- using iterator_category = std::random_access_iterator_tag;
+ operator X() const { return value; }
+ pw_element& operator=(X x) { value = std::move(x); return *this; };
- const T* ptr_ = nullptr;
- size_type i_ = 0;
+ double lower_bound() const { return extent.first; }
+ double upper_bound() const { return extent.second; }
- bool operator==(indexed_const_iterator x) const { return ptr_ == x.ptr_ && i_ == x.i_; }
- bool operator!=(indexed_const_iterator x) const { return !(*this==x); }
- bool operator<=(indexed_const_iterator x) const { return i_<=x.i_; }
- bool operator<(indexed_const_iterator x) const { return i_<x.i_; }
- bool operator>=(indexed_const_iterator x) const { return i_>=x.i_; }
- bool operator>(indexed_const_iterator x) const { return i_>x.i_; }
+ const std::pair<double, double> extent;
+ X value;
+};
- difference_type operator-(indexed_const_iterator x) const { return i_-x.i_; }
- indexed_const_iterator& operator++() { return ++i_, *this; }
- indexed_const_iterator& operator--() { return --i_, *this; }
- indexed_const_iterator operator++(int) { auto x = *this; return ++i_, x; }
- indexed_const_iterator operator--(int) { auto x = *this; return --i_, x; }
+template <>
+struct pw_element<void> {
+ pw_element():
+ extent(NAN, NAN)
+ {}
- indexed_const_iterator operator+(difference_type n) { return indexed_const_iterator{ptr_, i_+n}; }
- indexed_const_iterator operator-(difference_type n) { return indexed_const_iterator{ptr_, i_-n}; }
- indexed_const_iterator& operator+=(difference_type n) { return i_+=n, *this; }
- indexed_const_iterator& operator-=(difference_type n) { return i_-=n, *this; }
+ pw_element(std::pair<double, double> extent):
+ extent(std::move(extent))
+ {}
- friend indexed_const_iterator operator+(difference_type n, indexed_const_iterator x) {
- indexed_const_iterator r(std::move(x));
- return r+=n;
- }
+ pw_element(const pw_element&) = default;
+ pw_element(pw_element&&) = default;
- reference operator*() const { return (*ptr_)[i_]; }
- pointer operator->() const { return pointer{(*ptr_)[i_]}; }
-};
+ double lower_bound() const { return extent.first; }
+ double upper_bound() const { return extent.second; }
+ const std::pair<double, double> extent;
+};
template <typename X = void>
+struct pw_elements;
+
+// Proxy object for mutable iterators into a pw_elements<X> object.
+template <typename X>
+struct pw_element_proxy {
+ pw_element_proxy(pw_elements<X>& pw, pw_size_type i):
+ extent(pw.extent(i)), value(pw.value(i)) {}
+
+ operator pw_element<X>() const { return pw_element<X>{extent, value}; }
+ 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<double, double> extent;
+ X& value;
+};
+
+// Compute indices into vertex set corresponding to elements that cover a point x:
+
+namespace {
+std::pair<std::ptrdiff_t, std::ptrdiff_t> equal_range_indices(const std::vector<double>& vertices, double x) {
+ if (vertices.empty()) return {0, 0};
+
+ auto eq = std::equal_range(vertices.begin(), vertices.end(), x);
+
+ // Let n be the number of elements, indexed from 0 to n-1, with
+ // vertices indexed from 0 to n. Observe:
+ // * eq.first points to least vertex v_i ≥ x.
+ // or else to vertices.end() if v < x for all vertices v.
+ // * eq.second points to vertices.end() if the last vertex v_n ≤ x,
+ // or else to the least vertex v_k > x.
+ //
+ // Elements then correspond to the index range [b, e), where:
+ // * b=0 if i=0, else b=i-1, as v_i will be the upper vertex for
+ // the first element whose (closed) support contains x.
+ // * e=k if k<n, since v_k will be the upper vertex for the
+ // the last element (index k-1) whose support contains x.
+ // Otherwise, if k==n or eq.second is vertices.end(), the
+ // last element (index n-1) contains x, and so e=n.
+
+ if (eq.first==vertices.end()) return {0, 0};
+ if (eq.first>vertices.begin()) --eq.first;
+ if (eq.second==vertices.end()) --eq.second;
+
+ return std::make_pair(eq.first-vertices.begin(), eq.second-vertices.begin());
+}
+} // anonymous namespace
+
+template <typename X>
struct pw_elements {
using size_type = pw_size_type;
+ using difference_type = std::make_signed_t<pw_size_type>;
static constexpr size_type npos = pw_npos;
- struct value_type {
- std::pair<double, double> interval;
- X element;
+ using value_type = pw_element<X>;
+ using codomain = X;
- bool operator==(const value_type& other) const { return interval==other.interval && element==other.element; }
- bool operator!=(const value_type& other) const { return interval!=other.interval || element!=other.element; }
+ struct iterator: iterator_adaptor<iterator, counter<pw_size_type>> {
+ using typename iterator_adaptor<iterator, counter<pw_size_type>>::difference_type;
+ iterator(pw_elements<X>& pw, pw_size_type index): pw_(&pw), c_(index) {}
+ iterator(): pw_(nullptr) {}
+
+ using value_type = pw_element<X>;
+ using pointer = pointer_proxy<pw_element_proxy<X>>;
+ using reference = pw_element_proxy<X>;
+
+ reference operator[](difference_type j) { return reference{*pw_, j+*c_}; }
+ reference operator*() { return reference{*pw_, *c_}; }
+ value_type operator*() const { return (*pw_)[*c_]; }
+ pointer operator->() { return pointer{*pw_, *c_}; }
+
+ // (required for iterator_adaptor)
+ counter<pw_size_type>& inner() { return c_; }
+ const counter<pw_size_type>& inner() const { return c_; }
+
+ protected:
+ pw_elements<X>* pw_;
+ counter<pw_size_type> c_;
};
- using const_iterator = indexed_const_iterator<pw_elements<X>>;
- using iterator = const_iterator;
+ struct const_iterator: iterator_adaptor<const_iterator, counter<pw_size_type>> {
+ using typename iterator_adaptor<const_iterator, counter<pw_size_type>>::difference_type;
+ const_iterator(const pw_elements<X>& pw, pw_size_type index): pw_(&pw), c_(index) {}
+ const_iterator(): pw_(nullptr) {}
- // Consistency requirements:
- // 1. empty() || element_.size()+1 = vertex_.size()
- // 2. vertex_[i]<=vertex_[j] for i<=j.
+ using value_type = pw_element<X>;
+ using pointer = const pointer_proxy<pw_element<X>>;
+ using reference = pw_element<X>;
- std::vector<double> vertex_;
- std::vector<X> element_;
+ reference operator[](difference_type j) const { return (*pw_)[j+*c_]; }
+ reference operator*() const { return (*pw_)[*c_]; }
+ pointer operator->() const { return pointer{(*pw_)[*c_]}; }
+
+ // (required for iterator_adaptor)
+ counter<pw_size_type>& inner() { return c_; }
+ const counter<pw_size_type>& inner() const { return c_; }
+
+ protected:
+ const pw_elements<X>* pw_;
+ counter<pw_size_type> c_;
+ };
// Ctors and assignment:
pw_elements() = default;
- template <typename VSeq, typename ESeq>
- pw_elements(const VSeq& vs, const ESeq& es) {
+ template <typename VertexSeq, typename ValueSeq>
+ pw_elements(const VertexSeq& vs, const ValueSeq& es) {
assign(vs, es);
}
@@ -108,7 +276,7 @@ struct pw_elements {
template <typename Y>
explicit pw_elements(const pw_elements<Y>& from):
- vertex_(from.vertex_), element_(from.element_.begin(), from.element_.end())
+ vertex_(from.vertex_), value_(from.value_.begin(), from.value_.end())
{}
pw_elements& operator=(pw_elements&&) = default;
@@ -116,90 +284,107 @@ struct pw_elements {
// Access:
- auto intervals() const { return util::partition_view(vertex_); }
- auto interval(size_type i) const { return intervals()[i]; }
+ auto extents() const { return util::partition_view(vertex_); }
+
+ auto bounds() const { return extents().bounds(); }
+ auto lower_bound() const { return bounds().first; }
+ auto upper_bound() const { return bounds().second; }
- auto bounds() const { return intervals().bounds(); }
+ auto 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; }
- size_type size() const { return element_.size(); }
+ size_type size() const { return value_.size(); }
bool empty() const { return size()==0; }
bool operator==(const pw_elements& x) const {
- return vertex_==x.vertex_ && element_==x.element_;
+ return vertex_==x.vertex_ && value_==x.value_;
}
bool operator!=(const pw_elements& x) const { return !(*this==x); }
- const auto& elements() const { return element_; }
+ const auto& values() const { return value_; }
const auto& vertices() const { return vertex_; }
- X& element(size_type i) & { return element_[i]; }
- const X& element(size_type i) const & { return element_[i]; }
- value_type operator[](size_type i) const { return value_type{interval(i), element(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]; }
- const_iterator cbegin() const { return const_iterator{this, 0}; }
+ 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)}; }
+ auto operator[](size_type i) const && { return value_type{extent(i), value(i)}; }
+
+ auto front() & { return pw_element_proxy<X>{*this, 0}; }
+ auto front() const & { return value_type{extent(0), value(0)}; }
+ auto front() const && { return value_type{extent(0), value(0)}; }
+
+ auto back() & { return pw_element_proxy<X>{*this, size()-1}; }
+ auto back() const & { return value_type{extent(size()-1), value(size()-1)}; }
+ auto back() const && { return value_type{extent(size()-1), value(size()-1)}; }
+
+ iterator begin() { return iterator{*this, 0}; }
+ iterator end() { return iterator{*this, size()}; }
+
+ const_iterator cbegin() const { return const_iterator{*this, 0}; }
const_iterator begin() const { return cbegin(); }
- const_iterator cend() const { return const_iterator{this, size()}; }
+ const_iterator cend() const { return const_iterator{*this, size()}; }
const_iterator end() const { return cend(); }
- value_type front() const { return (*this)[0]; }
- value_type back() const { return (*this)[size()-1]; }
// Return index of right-most element whose corresponding closed interval contains x.
size_type index_of(double x) const {
if (empty()) return npos;
- auto partn = intervals();
+ auto partn = extents();
if (x == partn.bounds().second) return size()-1;
else return partn.index(x);
}
- // Return iterator pair spanning elements whose corresponding closed intervals contain x.
- std::pair<iterator, iterator> equal_range(double x) const {
- auto eq = std::equal_range(vertex_.begin(), vertex_.end(), x);
+ // Return iterator pair spanning elements whose corresponding closed extents contain x.
+ std::pair<iterator, iterator> equal_range(double x) {
+ auto qi = equal_range_indices(vertex_, x);
+ return {begin()+qi.first, begin()+qi.second};
+ }
- // Let n be the number of elements, indexed from 0 to n-1, with
- // vertices indexed from 0 to n. Observe:
- // * eq.first points to least vertex v_i ≥ x.
- // * eq.second points to vertex_.end() if the last vertex v_n ≤ x,
- // or else to the least vertex v_k > x.
- //
- // Elements then correspond to the index range [b, e), where:
- // * b=0 if i=0, else b=i-1, as v_i will be the upper vertex for
- // the first element whose (closed) support contains x.
- // * e=k if k<n, since v_k will be the upper vertex for the
- // the last element (index k-1) whose support contains x.
- // Otherwise, if k==n or eq.second is vertex_.end(), the
- // last element (index n-1) contains x, and so e=n.
+ std::pair<const_iterator, const_iterator> equal_range(double x) const {
+ auto qi = equal_range_indices(vertex_, x);
+ return {cbegin()+qi.first, cbegin()+qi.second};
+ }
- if (eq.first==vertex_.end()) return {end(), end()};
- if (eq.first>vertex_.begin()) --eq.first;
- if (eq.second==vertex_.end()) --eq.second;
+ auto operator()(double x) const && {
+ size_type i = index_of(x);
+ return i!=npos? (*this)[i]: throw std::range_error("position outside support");
+ }
- return {begin()+(eq.first-vertex_.begin()), begin()+(eq.second-vertex_.begin())};
+ auto operator()(double x) & {
+ size_type i = index_of(x);
+ return i!=npos? (*this)[i]: throw std::range_error("position outside support");
}
- value_type operator()(double x) const {
+ auto operator()(double x) const & {
size_type i = index_of(x);
- if (i==npos) {
- throw std::range_error("position outside support");
- }
- return (*this)[i];
+ return i!=npos? (*this)[i]: throw std::range_error("position outside support");
}
// mutating operations:
void reserve(size_type n) {
vertex_.reserve(n+1);
- element_.reserve(n);
+ value_.reserve(n);
}
void clear() {
vertex_.clear();
- element_.clear();
+ value_.clear();
+ }
+
+ void push_back(pw_element<X> elem) {
+ double left = elem.lower_bound();
+ double right = elem.upper_bound();
+ push_back(left, right, std::move(elem.value));
}
template <typename U>
- void push_back(double left, double right, U&& elem) {
+ void push_back(double left, double right, U&& v) {
if (!empty() && left!=vertex_.back()) {
throw std::runtime_error("noncontiguous element");
}
@@ -208,20 +393,19 @@ struct pw_elements {
throw std::runtime_error("inverted element");
}
- // Extend element_ first in case a conversion/copy/move throws.
-
- element_.push_back(std::forward<U>(elem));
+ // Extend value_ first in case a conversion/copy/move throws.
+ value_.push_back(std::forward<U>(v));
if (vertex_.empty()) vertex_.push_back(left);
vertex_.push_back(right);
}
template <typename U>
- void push_back(double right, U&& elem) {
+ void push_back(double right, U&& v) {
if (empty()) {
throw std::runtime_error("require initial left vertex for element");
}
- push_back(vertex_.back(), right, elem);
+ push_back(vertex_.back(), right, std::forward<U>(v));
}
void assign(std::initializer_list<double> vs, std::initializer_list<X> es) {
@@ -230,15 +414,15 @@ struct pw_elements {
}
template <typename Seq1, typename Seq2>
- void assign(const Seq1& vertices, const Seq2& elements) {
+ void assign(const Seq1& vertices, const Seq2& values) {
using std::begin;
using std::end;
auto vi = begin(vertices);
auto ve = end(vertices);
- auto ei = begin(elements);
- auto ee = end(elements);
+ auto ei = begin(values);
+ auto ee = end(values);
if (ei==ee) { // empty case
if (vi!=ve) {
@@ -248,13 +432,13 @@ struct pw_elements {
return;
}
- double left = *vi++;
if (vi==ve) {
throw std::runtime_error("vertex list too short");
}
clear();
+ double left = *vi++;
double right = *vi++;
push_back(left, right, *ei++);
@@ -270,91 +454,144 @@ struct pw_elements {
throw std::runtime_error("vertex list too long");
}
}
+
+private:
+ // Consistency requirements:
+ // 1. empty() || value_.size()+1 = vertex_.size()
+ // 2. vertex_[i]<=vertex_[j] for i<=j.
+
+ std::vector<double> vertex_;
+ std::vector<X> value_;
};
-// With X = void, present the element intervals only,
-// keeping othewise the same interface.
+// With X = void, present the element intervals only, keeping otherwise the
+// same interface.
-template <> struct pw_elements<void> {
+template <>
+struct pw_elements<void> {
using size_type = pw_size_type;
static constexpr size_type npos = pw_npos;
+ using difference_type = std::make_signed_t<pw_size_type>;
- std::vector<double> vertex_;
+ using value_type = pw_element<void>;
+ using codomain = void;
+
+ struct const_iterator: iterator_adaptor<const_iterator, counter<pw_size_type>> {
+ using typename iterator_adaptor<const_iterator, counter<pw_size_type>>::difference_type;
+ const_iterator(const pw_elements<void>& pw, pw_size_type index): pw_(&pw), c_(index) {}
+ const_iterator(): pw_(nullptr) {}
+
+ using value_type = pw_element<void>;
+ using pointer = const pointer_proxy<pw_element<void>>;
+ using reference = pw_element<void>;
+
+ reference operator[](difference_type j) const { return (*pw_)[j+*c_]; }
+ reference operator*() const { return (*pw_)[*c_]; }
+ pointer operator->() const { return pointer{(*pw_)[*c_]}; }
- struct value_type {
- std::pair<double, double> interval;
+ // (required for iterator_adaptor)
+ counter<pw_size_type>& inner() { return c_; }
+ const counter<pw_size_type>& inner() const { return c_; }
- bool operator==(const value_type& other) const { return interval==other.interval; }
- bool operator!=(const value_type& other) const { return interval!=other.interval; }
+ protected:
+ const pw_elements<void>* pw_;
+ counter<pw_size_type> c_;
};
- using const_iterator = indexed_const_iterator<pw_elements<void>>;
using iterator = const_iterator;
- // ctors and assignment:
+ // Ctors and assignment:
+
+ pw_elements() = default;
- template <typename VSeq>
- explicit pw_elements(const VSeq& vs) { assign(vs); }
+ template <typename VertexSeq>
+ pw_elements(const VertexSeq& vs) {
+ assign(vs);
+ }
- pw_elements(std::initializer_list<double> vs) { assign(vs); }
+ pw_elements(std::initializer_list<double> vs) {
+ assign(vs);
+ }
- pw_elements() = default;
pw_elements(pw_elements&&) = default;
pw_elements(const pw_elements&) = default;
template <typename Y>
explicit pw_elements(const pw_elements<Y>& from):
- vertex_(from.vertex_) {}
+ vertex_(from.vertices()) {}
pw_elements& operator=(pw_elements&&) = default;
pw_elements& operator=(const pw_elements&) = default;
- // access:
+ // Access:
+
+ auto extents() const { return util::partition_view(vertex_); }
- auto intervals() const { return util::partition_view(vertex_); }
- auto interval(size_type i) const { return intervals()[i]; }
- value_type operator[](size_type i) const { return value_type{interval(i)}; }
+ auto bounds() const { return extents().bounds(); }
+ auto lower_bound() const { return bounds().first; }
+ auto upper_bound() const { return bounds().second; }
- auto bounds() const { return intervals().bounds(); }
+ auto 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; }
- size_type size() const { return vertex_.empty()? 0: vertex_.size()-1; }
+ size_type size() const { return empty()? 0: vertex_.size()-1; }
bool empty() const { return vertex_.empty(); }
- bool operator==(const pw_elements& x) const { return vertex_==x.vertex_; }
+ bool operator==(const pw_elements& x) const {
+ return vertex_==x.vertex_;
+ }
+
bool operator!=(const pw_elements& x) const { return !(*this==x); }
const auto& vertices() const { return vertex_; }
+ void value(size_type i) const {}
+ pw_element<void> operator[](size_type i) const { return value_type{extent(i)}; }
+
+ pw_element<void> front() const { return value_type{extent(0)}; }
+ pw_element<void> back() const { return value_type{extent(size()-1)}; }
+
+ const_iterator cbegin() const { return const_iterator{*this, 0}; }
+ const_iterator begin() const { return cbegin(); }
+ const_iterator cend() const { return const_iterator{*this, size()}; }
+ const_iterator end() const { return cend(); }
+
+ // Return index of right-most element whose corresponding closed interval contains x.
size_type index_of(double x) const {
if (empty()) return npos;
- auto partn = intervals();
+ auto partn = extents();
if (x == partn.bounds().second) return size()-1;
else return partn.index(x);
}
- // Return iterator pair spanning elements whose corresponding closed intervals contain x.
- std::pair<iterator, iterator> equal_range(double x) const {
- auto eq = std::equal_range(vertex_.begin(), vertex_.end(), x);
-
- if (eq.first==vertex_.end()) return {end(), end()};
- if (eq.first>vertex_.begin()) --eq.first;
- if (eq.second==vertex_.end()) --eq.second;
-
- return {begin()+(eq.first-vertex_.begin()), begin()+(eq.second-vertex_.begin())};
+ // Return iterator pair spanning elements whose corresponding closed extents contain x.
+ std::pair<const_iterator,const_iterator> equal_range(double x) const {
+ auto qi = equal_range_indices(vertex_, x);
+ return {cbegin()+qi.first, cbegin()+qi.second};
}
- const_iterator cbegin() const { return const_iterator{this, 0}; }
- const_iterator begin() const { return cbegin(); }
- const_iterator cend() const { return const_iterator{this, size()}; }
- const_iterator end() const { return cend(); }
- value_type front() const { return (*this)[0]; }
- value_type back() const { return (*this)[size()-1]; }
+ auto operator()(double x) const {
+ size_type i = index_of(x);
+ return i!=npos? (*this)[i]: throw std::range_error("position outside support");
+ }
// mutating operations:
- void reserve(size_type n) { vertex_.reserve(n+1); }
- void clear() { vertex_.clear(); }
+ void reserve(size_type n) {
+ vertex_.reserve(n+1);
+ }
+
+ void clear() {
+ vertex_.clear();
+ }
+
+ void push_back(const pw_element<void>& elem) {
+ double left = elem.lower_bound();
+ double right = elem.upper_bound();
+ push_back(left, right);
+ }
void push_back(double left, double right) {
if (!empty() && left!=vertex_.back()) {
@@ -373,11 +610,12 @@ template <> struct pw_elements<void> {
if (empty()) {
throw std::runtime_error("require initial left vertex for element");
}
- vertex_.push_back(right);
+
+ push_back(vertex_.back(), right);
}
void assign(std::initializer_list<double> vs) {
- assign(util::make_range(vs.begin(), vs.end()));
+ assign(make_range(vs.begin(), vs.end()));
}
template <typename Seq1>
@@ -408,66 +646,79 @@ template <> struct pw_elements<void> {
push_back(right);
}
}
-};
-template <typename X>
-using pw_element = typename pw_elements<X>::value_type;
+private:
+ // Consistency requirements:
+ // 1. empty() || value_.size()+1 = vertex_.size()
+ // 2. vertex_[i]<=vertex_[j] for i<=j.
-namespace impl {
- template <typename A, typename B>
- struct piecewise_pairify {
- std::pair<A, B> operator()(
- double, double,
- const pw_element<A>& a_elem,
- const pw_element<B>& b_elem) const
- {
- return {a_elem.element, b_elem.element};
- }
- };
+ std::vector<double> vertex_;
+};
- template <typename X>
- struct piecewise_pairify<X, void> {
- X operator()(
- double, double,
- const pw_element<X>& a_elem,
- const pw_element<void>& b_elem) const
- {
- return a_elem.element;
+// The piecewise map applies a transform to the values in a piecewise
+// object, preserving the extents. If the original piecewise object
+// is `pw_elements<void>`, then the function is called with zero arguments.
+//
+// If the mapping function returns void, return a `pw_elements<void>` object
+// with the same extents as the original piecewise object.
+
+template <typename X, typename Fn>
+auto pw_map(const pw_elements<X>& pw, Fn&& fn) {
+ if constexpr (std::is_void_v<X>) {
+ using Out = std::invoke_result_t<Fn&&>;
+ if constexpr (std::is_void_v<Out>) {
+ // Evalate fn for side effects.
+ for (const auto& elem: pw) fn();
+ return pw;
}
- };
-
- template <typename X>
- struct piecewise_pairify<void, X> {
- X operator()(
- double, double,
- const pw_element<void>& a_elem,
- const pw_element<X>& b_elem) const
- {
- return b_elem.element;
+ else {
+ auto mapped = util::transform_view(pw, [&fn](auto&&) { return fn(); });
+ return pw_elements<Out>(pw.vertices(), std::vector<Out>(mapped.begin(), mapped.end()));
}
- };
-
- template <>
- struct piecewise_pairify<void, void> {
- void operator()(
- double, double,
- const pw_element<void>&,
- const pw_element<void>&) const {}
- };
+ }
+ else {
+ using Out = std::invoke_result_t<Fn&&, X>;
+ if constexpr (std::is_void_v<Out>) {
+ // Evalate fn for side effects.
+ for (const auto& v: pw.values()) fn(v);
+ return pw_elements<void>(pw);
+ }
+ else {
+ auto mapped = util::transform_view(pw, [&fn](auto&& elem) { return fn(elem.value); });
+ return pw_elements<Out>(pw.vertices(), std::vector<Out>(mapped.begin(), mapped.end()));
+ }
+ }
}
-// Zip combines successive elements from two pw_elements sequences where the
-// elements overlap. Let A_i and B_i denote the ordered elements from each of
-// the sequences A and B, and Z_i denote the elements from the resulting
-// sequence Z.
+// The piecewise zip combines successive elements from two `pw_elements`
+// sequences where the elements overlap.
+//
+// * `pw_zip_view` performs a lazy zip, represented by a range of
+// `pw_zip_iterator` objects. The iterators dereference to objects of
+// type `pw_element<std::pair<pw_element<A>, pw_element<B>>`.
+//
+// * `pw_zip` performs a strict zip, returning a piecewise object of type
+// `pw_elements<std::pair<pw_element<A>, pw_element<B>>`.
+//
+// * `pw_zip_with` performs a map composed with a zip, where the map is
+// given by a function which takes three arguments: the extent of
+// the zipped element as `std::pair<double, double>`; the element
+// from the first sequence of type `<pw_element<A>`; and the element
+// from the second sequence of type `pw_element<B>`. It is equivalent
+// in action to performing a `pw_zip` followed by a `pw_map`.
+//
+// By default, `pw_zip_with` uses `pw_pairify` to combine elements,
+// which returns a tuple of the non-void values from each pair of elements
+// in the zip.
//
-// * The support (`bounds()`) of the zipped result is the intersections of the
-// support of the two sequences.
+// The elements of the zip are determined as follows:
//
-// * Each element Z_k in the zip corresponds to the intersection of some
+// Let A_i and B_i denote the ordered elements from each of the sequences A and
+// B, and Z_i denote the elements forming the resulting sequence Z.
+//
+// * Each element Z_k in the zip corresponds to the intersection of an
// element A_i and B_j. The extent of Z_k is the intersection of the
-// extents of A_i and B_j, and its value is determined by the supplied
-// `combine` function.
+// extents of A_i and B_j, and its value is the pair {A_i, B_j}.
//
// * For every element in A_i in A, if A_i intersects with an element of
// B, then there will be an Z_k corresponding to the intersection of A_i
@@ -477,53 +728,157 @@ namespace impl {
// not repeat. If Z_k is derived from A_i and B_j, then Z_(k+1) is derived
// from A_(i+1) and B_(j+1) if possible, or else from A_i and B_(j+1) or
// A_(i+1) and B_j.
-//
-// The Combine functional takes four arguments: double left, double right,
-// pw_elements<A>::value_type, pw_elements<B>::value_type b. The default
-// combine functional returns std::pair<A, B>, unless one of A and B is void.
-//
-// TODO: Consider making a lazy `zip_view` version of zip.
-template <typename A, typename B, typename Combine = impl::piecewise_pairify<A, B>>
-auto zip(const pw_elements<A>& a, const pw_elements<B>& b, Combine combine = {})
-{
- using Out = decltype(combine(0., 0., a.front(), b.front()));
- constexpr bool is_void = std::is_void_v<Out>;
+template <typename A, typename B>
+struct pw_zip_iterator {
+ using value_type = pw_element<std::pair<pw_element<A>, pw_element<B>>>;
+ using pointer = pointer_proxy<value_type>;
+ using reference = value_type;
+ using difference_type = std::ptrdiff_t;
+
+ // Proxy iterator so not _really_ a forward iterator.
+ using iterator_category = std::forward_iterator_tag;
+
+ bool is_end = true;
+ typename pw_elements<A>::const_iterator ai, a_end;
+ typename pw_elements<B>::const_iterator bi, b_end;
+ double left;
+
+ pw_zip_iterator() = default;
+ pw_zip_iterator(const pw_elements<A>& a, const pw_elements<B>& b) {
+ double lmax = std::max(a.lower_bound(), b.lower_bound());
+ double rmin = std::min(a.upper_bound(), b.upper_bound());
+
+ is_end = rmin<lmax;
+ if (!is_end) {
+ ai = a.equal_range(lmax).first;
+ a_end = a.equal_range(rmin).second;
+ bi = b.equal_range(lmax).first;
+ b_end = b.equal_range(rmin).second;
+ left = lmax;
+ }
+ }
- pw_elements<Out> z;
- if (a.empty() || b.empty()) return z;
+ pw_zip_iterator(const pw_zip_iterator& other):
+ is_end(other.is_end),
+ ai(other.ai),
+ a_end(other.a_end),
+ bi(other.bi),
+ b_end(other.b_end),
+ left(other.left)
+ {}
- double lmax = std::max(a.bounds().first, b.bounds().first);
- double rmin = std::min(a.bounds().second, b.bounds().second);
- if (rmin<lmax) return z;
+ bool operator==(pw_zip_iterator other) const {
+ if (is_end && other.is_end) return true;
+ return is_end==other.is_end && ai==other.ai && a_end==other.a_end && bi==other.bi && b_end==other.b_end && left==other.left;
+ }
- auto ai = a.equal_range(lmax).first;
- auto bi = b.equal_range(lmax).first;
+ bool operator!=(pw_zip_iterator other) const {
+ return !(*this==other);
+ }
- auto a_end = a.equal_range(rmin).second;
- auto b_end = b.equal_range(rmin).second;
+ pw_zip_iterator& operator++() {
+ if (is_end) return *this;
- double left = lmax;
- double a_right = ai->interval.second;
- double b_right = bi->interval.second;
- for (;;) {
+ double a_right = ai->upper_bound();
+ double b_right = bi->upper_bound();
double right = std::min(a_right, b_right);
- if constexpr (is_void) {
- z.push_back(left, right);
+
+ bool advance_a = a_right==right && std::next(ai)!=a_end;
+ bool advance_b = b_right==right && std::next(bi)!=b_end;
+ if (!advance_a && !advance_b) {
+ is_end = true;
}
else {
- z.push_back(left, right, combine(left, right, *ai, *bi));
+ if (advance_a) ++ai;
+ if (advance_b) ++bi;
+ left = right;
}
- bool advance_a = a_right==right && std::next(ai)!=a_end;
- bool advance_b = b_right==right && std::next(bi)!=b_end;
- if (!advance_a && !advance_b) break;
+ return *this;
+ }
+
+ pw_zip_iterator operator++(int) {
+ pw_zip_iterator here = *this;
+ ++*this;
+ return here;
+ }
+
+ value_type 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}};
+ }
+
+ pointer operator->() const {
+ return pointer{*this};
+ }
+};
- if (advance_a) a_right = (++ai)->interval.second;
- if (advance_b) b_right = (++bi)->interval.second;
- left = right;
+template <typename A, typename B>
+range<pw_zip_iterator<A, B>> pw_zip_range(const pw_elements<A>& a, const pw_elements<B>& b) {
+ return {pw_zip_iterator<A, B>{a, b}, pw_zip_iterator<A, B>{}};
+}
+
+template <typename A, typename B>
+auto pw_zip(const pw_elements<A>& a, const pw_elements<B>& b) {
+ pw_elements<std::pair<pw_element<A>, pw_element<B>>> out;
+ for (auto elem: pw_zip_range(a, b)) {
+ out.push_back(elem);
+ }
+ return out;
+}
+
+// `pw_pairify` is a functional intended for use with `pw_zip_with` that takes
+// an extent, which is ignored, and two elements of types `pw_element<A>` and
+// `pw_element<B>` and returns void if both A and B are void, or the pair of
+// their values if neither A nor B is void, or else the single non-void value
+// of type A or B.
+
+namespace impl {
+ template <typename A, typename B>
+ std::pair<A, B> pw_pairify_(const pw_element<A>& a_elem, const pw_element<B>& b_elem) {
+ return {a_elem.value, b_elem.value};
+ };
+
+ template <typename A>
+ A pw_pairify_(const pw_element<A>& a_elem, const pw_element<void>&) {
+ return a_elem.value;
+ };
+
+ template <typename B>
+ B pw_pairify_(const pw_element<void>&, const pw_element<B>& b_elem) {
+ return b_elem.value;
+ };
+
+ inline void pw_pairify_(const pw_element<void>&, const pw_element<void>&) {}
+}
+
+struct pw_pairify {
+ template <typename A, typename B>
+ auto operator()(std::pair<double, double>, const pw_element<A>& a_elem, const pw_element<B>& b_elem) const {
+ return impl::pw_pairify_(a_elem, b_elem);
+ }
+};
+
+template <typename A, typename B, typename Fn = pw_pairify>
+auto pw_zip_with(const pw_elements<A>& a, const pw_elements<B>& b, Fn&& fn = Fn{}) {
+ using Out = decltype(fn(std::pair<double, double>{}, a.front(), b.front()));
+ pw_elements<Out> out;
+
+ for (auto elem: pw_zip_range(a, b)) {
+ if constexpr (std::is_void_v<Out>) {
+ fn(elem.extent, elem.value.first, elem.value.second);
+ out.push_back(elem.extent.first, elem.extent.second);
+ }
+ else {
+ out.push_back(elem.extent.first, elem.extent.second,
+ fn(elem.extent, elem.value.first, elem.value.second));
+ }
}
- return z;
+ return out;
}
} // namespace util
diff --git a/arbor/util/ratelem.hpp b/arbor/util/ratelem.hpp
index cbb2bfac..396de714 100644
--- a/arbor/util/ratelem.hpp
+++ b/arbor/util/ratelem.hpp
@@ -53,7 +53,6 @@ void rat_interpolate(std::array<double, m>& h, const std::array<double, m+1>& g,
h[i] = ook*((x - i)*g[i+1] + (i+k - x)*g[i]);
}
else {
- // Using h[i] = k/(g[i+1]/(x - i) + g[i]/(i+k - x)) is more robust to
// singularities, but the expense should not be necessary if we stick
// to strictly monotonic elements for rational polynomials.
h[i] = k*g[i]*g[i+1]/(g[i]*(x - i) + g[i+1]*(i+k - x));
@@ -62,7 +61,7 @@ void rat_interpolate(std::array<double, m>& h, const std::array<double, m+1>& g,
}
template <unsigned a, unsigned c, unsigned k, bool upper>
-double rat_eval(const std::array<double, 1+a+c>& g, const std::array<double, 2+a+c>& p, double x) {
+double recursive_rat_eval(const std::array<double, 1+a+c>& g, const std::array<double, 2+a+c>& p, double x) {
if constexpr (a==0 && c==0) {
return g[0];
}
@@ -73,29 +72,49 @@ double rat_eval(const std::array<double, 1+a+c>& g, const std::array<double, 2+a
h[i] = p[i+1] + k/((x - i)/(g[i+1]-p[i+1]) + (i+k - x)/(g[i]-p[i+1]));
}
- return rat_eval<0, c-1, k+1, upper>(h, g, x);
+ return recursive_rat_eval<0, c-1, k+1, upper>(h, g, x);
}
else {
std::array<double, a+c> h;
rat_interpolate<k, upper>(h, g, x);
- return rat_eval<a-1, c, k+1, upper>(h, g, x);
+ return recursive_rat_eval<a-1, c, k+1, upper>(h, g, x);
}
}
+template <unsigned p, unsigned q>
+double rat_eval(const std::array<double, 1+p+q>& data, double x) {
+ // upper => interpolate polynomials first;
+ // !upper => interpolate reciprocal polynomials first;
+ // a is the number of (reciprocal) polynomial interpolation steps;
+ // c is the number of 'rhombus' interpolation steps.
+
+ constexpr bool upper = p>=q;
+ constexpr unsigned a = upper? p-q+(q>0): q-p+(p>0);
+ constexpr unsigned c = p+q-a;
-template <unsigned a, unsigned c, unsigned k, bool upper>
-double rat_eval(const std::array<double, 1+a+c>& g, double x) {
if constexpr (a==0 && c==0) {
- return g[0];
+ return data[0];
+ }
+ else if constexpr (a==1 && c==1) {
+ // Special case when p==1 and q==1 so that we can have
+ // well-defined behaviour in the case where the interpolants
+ // are strictly monotonic, but might take non-finite values.
+ // (Used for ixa interpolation in embed_pwlin.cpp).
+
+ double y0 = data[0];
+ double y1 = data[1];
+ double y2 = data[2];
+
+ return y1+(2*x-1)/(x/(y2-y1)+(1-x)/(y1-y0));
}
else {
+ x *= p+q;
std::array<double, a+c> h;
- rat_interpolate<k, upper>(h, g, x);
- return rat_eval<a-1, c, k+1, upper>(h, g, x);
+ rat_interpolate<1, upper>(h, data, x);
+ return recursive_rat_eval<a-1, c, 2, upper>(h, data, x);
}
}
-
} // namespace impl
template <unsigned p, unsigned q>
@@ -132,16 +151,7 @@ struct rat_element {
// Rational interpolation at x.
double operator()(double x) const {
- // upper j => interpolate polynomials first;
- // !upper => interpolate reciprocal polynomials first;
- // a is the number of (reciprocal) polynomial interpolation steps;
- // c is the number of 'rhombus' interpolation steps.
-
- constexpr bool upper = p>=q;
- constexpr unsigned a = upper? p-q+(q>0): q-p+(p>0);
- constexpr unsigned c = p+q-a;
-
- return impl::rat_eval<a, c, 1, upper>(data_, x*(p+q));
+ return impl::rat_eval<p, q>(data_, x);
}
// Node values.
diff --git a/test/unit/test_morph_embedding.cpp b/test/unit/test_morph_embedding.cpp
index 00345b08..ee8d0f73 100644
--- a/test/unit/test_morph_embedding.cpp
+++ b/test/unit/test_morph_embedding.cpp
@@ -286,3 +286,73 @@ TEST(embedding, area_0_length_segment) {
double expected_a2 = expected_a1 + pi*(20*20-10*10);
EXPECT_TRUE(near_relative(a2, expected_a2, reltol));
}
+
+TEST(embedding, small_radius) {
+ using testing::near_relative;
+ constexpr double pi = math::pi<double>;
+ constexpr double reltol = 1e-10;
+
+ segment_tree t;
+ t.append(mnpos, { 0, 0, 0, 10}, {10, 0, 0, 0.00001}, 0);
+ t.append(0, {10, 0, 0, 40}, {20, 0, 0, 40}, 0);
+
+ embedding em{morphology(t)};
+
+ // Integrated inverse cross-sectional area in segment 1
+ // corresponding to cable(0, 0.5, 1):
+ // radius râ‚€ = 40, râ‚ = 40, length = 10.
+
+ double expected_ixa = 10/(40*40*pi);
+ double computed_ixa = em.integrate_ixa(mcable{0, 0.5, 1.0});
+ ASSERT_FALSE(std::isnan(computed_ixa));
+ EXPECT_TRUE(near_relative(expected_ixa, computed_ixa, reltol));
+}
+
+TEST(embedding, zero_radius) {
+ using testing::near_relative;
+ constexpr double pi = math::pi<double>;
+ constexpr double reltol = 1e-10;
+
+ segment_tree t;
+ t.append(mnpos, { 0, 0, 0, 10}, {10, 0, 0, 0}, 0);
+ t.append(0, {10, 0, 0, 40}, {20, 0, 0, 40}, 0);
+
+ embedding em{morphology(t)};
+
+ // Integrated inverse cross-sectional area in segment 1
+ // corresponding to cable(0, 0.5, 1):
+ // radius râ‚€ = 40, râ‚ = 40, length = 10.
+
+ double expected_ixa = 10/(40*40*pi);
+ double computed_ixa = em.integrate_ixa(mcable{0, 0.5, 1.0});
+ ASSERT_FALSE(std::isnan(computed_ixa));
+ EXPECT_TRUE(near_relative(expected_ixa, computed_ixa, reltol));
+
+ // Integrating over the zero radius point though should give us
+ // INFINITY.
+
+ double infinite_ixa = em.integrate_ixa(mcable{0, 0.25, 0.75});
+ ASSERT_FALSE(std::isnan(infinite_ixa));
+ EXPECT_TRUE(std::isinf(infinite_ixa));
+
+ // Integrating to the zero radius point should also give us
+ // INFINITY, because we integrate over closed intervals.
+
+ double also_infinite_ixa = em.integrate_ixa(mcable{0, 0.25, 0.5});
+ ASSERT_FALSE(std::isnan(also_infinite_ixa));
+ EXPECT_TRUE(std::isinf(also_infinite_ixa));
+
+ // Should be able to integrate ixa over a tree that starts
+ // with a zero radius.
+
+ segment_tree t2;
+ t2.append(mnpos, { 0, 0, 0, 0}, {10, 0, 0, 10}, 0);
+ embedding em2{morphology(t2)};
+
+ expected_ixa = 5/(10*5*pi);
+ computed_ixa = em2.integrate_ixa(mcable{0, 0.5, 1.0});
+ ASSERT_FALSE(std::isnan(computed_ixa));
+ EXPECT_TRUE(near_relative(expected_ixa, computed_ixa, reltol));
+
+ EXPECT_TRUE(std::isinf(em2.integrate_ixa(mcable{0, 0, 0.5})));
+}
diff --git a/test/unit/test_piecewise.cpp b/test/unit/test_piecewise.cpp
index f05b540a..0f65cf89 100644
--- a/test/unit/test_piecewise.cpp
+++ b/test/unit/test_piecewise.cpp
@@ -82,33 +82,33 @@ TEST(piecewise, assign) {
ASSERT_EQ(4u, p.size());
- EXPECT_EQ(10, p[0].element);
- EXPECT_EQ( 8, p[1].element);
- EXPECT_EQ( 9, p[2].element);
- EXPECT_EQ( 4, p[3].element);
+ EXPECT_EQ(10, p[0].value);
+ EXPECT_EQ( 8, p[1].value);
+ EXPECT_EQ( 9, p[2].value);
+ EXPECT_EQ( 4, p[3].value);
using dp = std::pair<double, double>;
- EXPECT_EQ(dp(1.0, 1.5), p.interval(0));
- EXPECT_EQ(dp(1.5, 2.0), p.interval(1));
- EXPECT_EQ(dp(2.0, 2.5), p.interval(2));
- EXPECT_EQ(dp(2.5, 3.0), p.interval(3));
+ EXPECT_EQ(dp(1.0, 1.5), p.extent(0));
+ EXPECT_EQ(dp(1.5, 2.0), p.extent(1));
+ EXPECT_EQ(dp(2.0, 2.5), p.extent(2));
+ EXPECT_EQ(dp(2.5, 3.0), p.extent(3));
pw_elements<int> q1(p);
pw_elements<int> q2;
q2 = p;
pw_elements<int> q3(p);
- q3.assign(p.vertices(), p.elements());
+ q3.assign(p.vertices(), p.values());
EXPECT_EQ((std::vector<double>{1.0, 1.5, 2.0, 2.5, 3.0}), p.vertices());
- EXPECT_EQ((std::vector<int>{10, 8, 9, 4}), p.elements());
+ EXPECT_EQ((std::vector<int>{10, 8, 9, 4}), p.values());
EXPECT_EQ(q1.vertices(), p.vertices());
EXPECT_EQ(q2.vertices(), p.vertices());
EXPECT_EQ(q3.vertices(), p.vertices());
- EXPECT_EQ(q1.elements(), p.elements());
- EXPECT_EQ(q2.elements(), p.elements());
- EXPECT_EQ(q3.elements(), p.elements());
+ EXPECT_EQ(q1.values(), p.values());
+ EXPECT_EQ(q2.values(), p.values());
+ EXPECT_EQ(q3.values(), p.values());
q3.assign({}, {});
EXPECT_TRUE(q3.empty());
@@ -131,10 +131,10 @@ TEST(piecewise, assign_void) {
ASSERT_EQ(4u, p.size());
using dp = std::pair<double, double>;
- EXPECT_EQ(dp(1.0, 1.5), p.interval(0));
- EXPECT_EQ(dp(1.5, 2.0), p.interval(1));
- EXPECT_EQ(dp(2.0, 2.5), p.interval(2));
- EXPECT_EQ(dp(2.5, 3.0), p.interval(3));
+ EXPECT_EQ(dp(1.0, 1.5), p.extent(0));
+ EXPECT_EQ(dp(1.5, 2.0), p.extent(1));
+ EXPECT_EQ(dp(2.0, 2.5), p.extent(2));
+ EXPECT_EQ(dp(2.5, 3.0), p.extent(3));
pw_elements<> q1(p);
pw_elements<> q2;
@@ -168,14 +168,14 @@ TEST(piecewise, access) {
p.assign(v, x);
for (unsigned i = 0; i<4; ++i) {
- EXPECT_EQ(v[i], p[i].interval.first);
- EXPECT_EQ(v[i+1], p[i].interval.second);
+ EXPECT_EQ(v[i], p[i].extent.first);
+ EXPECT_EQ(v[i+1], p[i].extent.second);
- EXPECT_EQ(v[i], p.interval(i).first);
- EXPECT_EQ(v[i+1], p.interval(i).second);
+ EXPECT_EQ(v[i], p.extent(i).first);
+ EXPECT_EQ(v[i+1], p.extent(i).second);
- EXPECT_EQ(x[i], p[i].element);
- EXPECT_EQ(x[i], p.element(i));
+ EXPECT_EQ(x[i], p[i].value);
+ EXPECT_EQ(x[i], p.value(i));
}
EXPECT_EQ(p[0], p.front());
@@ -191,12 +191,16 @@ TEST(piecewise, bounds) {
pw_elements<int> p{{1., 1.5, 2., 2.5, 3.}, {10, 8, 9, 4}};
EXPECT_EQ(1., p.bounds().first);
+ EXPECT_EQ(1., p.lower_bound());
EXPECT_EQ(3., p.bounds().second);
+ EXPECT_EQ(3., p.upper_bound());
pw_elements<> v{{1., 1.5, 2., 2.5, 3.}};
EXPECT_EQ(1., v.bounds().first);
+ EXPECT_EQ(1., p.lower_bound());
EXPECT_EQ(3., v.bounds().second);
+ EXPECT_EQ(3., p.upper_bound());
}
TEST(piecewise, index_of) {
@@ -240,21 +244,21 @@ TEST(piecewise, equal_range) {
auto er1 = p.equal_range(1.0);
ASSERT_EQ(1, er1.second-er1.first);
- EXPECT_EQ(10, er1.first->element);
+ EXPECT_EQ(10, er1.first->value);
auto er2 = p.equal_range(2.0);
ASSERT_EQ(2, er2.second-er2.first);
auto iter = er2.first;
- EXPECT_EQ(10, iter++->element);
- EXPECT_EQ(9, iter->element);
+ EXPECT_EQ(10, iter++->value);
+ EXPECT_EQ(9, iter->value);
auto er3_5 = p.equal_range(3.5);
ASSERT_EQ(1, er3_5.second-er3_5.first);
- EXPECT_EQ(8, er3_5.first->element);
+ EXPECT_EQ(8, er3_5.first->value);
auto er4 = p.equal_range(4.0);
ASSERT_EQ(1, er4.second-er4.first);
- EXPECT_EQ(8, er4.first->element);
+ EXPECT_EQ(8, er4.first->value);
auto er5 = p.equal_range(5.0);
ASSERT_EQ(er5.first, er5.second);
@@ -269,27 +273,26 @@ TEST(piecewise, equal_range) {
auto er1 = p.equal_range(1.0);
ASSERT_EQ(2, er1.second-er1.first);
auto iter = er1.first;
- EXPECT_EQ(10, iter++->element);
- EXPECT_EQ(11, iter++->element);
+ EXPECT_EQ(10, iter++->value);
+ EXPECT_EQ(11, iter++->value);
auto er2 = p.equal_range(2.0);
ASSERT_EQ(4, er2.second-er2.first);
iter = er2.first;
- EXPECT_EQ(11, iter++->element);
- EXPECT_EQ(12, iter++->element);
- EXPECT_EQ(13, iter++->element);
- EXPECT_EQ(14, iter++->element);
+ EXPECT_EQ(11, iter++->value);
+ EXPECT_EQ(12, iter++->value);
+ EXPECT_EQ(13, iter++->value);
+ EXPECT_EQ(14, iter++->value);
auto er3 = p.equal_range(3.0);
ASSERT_EQ(2, er3.second-er3.first);
iter = er3.first;
- EXPECT_EQ(14, iter++->element);
- EXPECT_EQ(15, iter++->element);
+ EXPECT_EQ(14, iter++->value);
+ EXPECT_EQ(15, iter++->value);
auto er5 = p.equal_range(5.0);
ASSERT_EQ(er5.first, er5.second);
}
-
}
TEST(piecewise, push) {
@@ -302,14 +305,14 @@ TEST(piecewise, push) {
q.clear();
q.push_back(1.1, 3.1, 4);
q.push_back(3.1, 4.3, 5);
- EXPECT_EQ(dp(1.1, 3.1), q.interval(0));
- EXPECT_EQ(dp(3.1, 4.3), q.interval(1));
- EXPECT_EQ(4, q[0].element);
- EXPECT_EQ(5, q[1].element);
+ EXPECT_EQ(dp(1.1, 3.1), q.extent(0));
+ EXPECT_EQ(dp(3.1, 4.3), q.extent(1));
+ EXPECT_EQ(4, q[0].value);
+ EXPECT_EQ(5, q[1].value);
q.push_back(7.2, 6);
- EXPECT_EQ(dp(4.3, 7.2), q.interval(2));
- EXPECT_EQ(6, q[2].element);
+ EXPECT_EQ(dp(4.3, 7.2), q.extent(2));
+ EXPECT_EQ(6, q[2].value);
// Supplied left side doesn't match current right.
EXPECT_THROW(q.push_back(7.4, 9.1, 7), std::runtime_error);
@@ -329,12 +332,53 @@ TEST(piecewise, push_void) {
EXPECT_EQ(3u, p.size());
EXPECT_EQ((std::vector<double>{0.1, 0.2, 0.4, 0.5}), p.vertices());
- EXPECT_EQ(dp(0.2,0.4), p.interval(1));
+ EXPECT_EQ(dp(0.2,0.4), p.extent(1));
// Supplied left side doesn't match current right.
EXPECT_THROW(p.push_back(0.7, 0.9), std::runtime_error);
}
+TEST(piecewise, mutate) {
+ pw_elements<int> p({1., 2., 3., 3., 4., 5.}, {10, 11, 12, 13, 14});
+ ASSERT_EQ(10, p.value(0));
+ ASSERT_EQ(11, p.value(1));
+
+ p.value(0) = 20;
+ EXPECT_EQ((std::vector<int>{20, 11, 12, 13, 14}), p.values());
+
+ p.front() = 30;
+ EXPECT_EQ((std::vector<int>{30, 11, 12, 13, 14}), p.values());
+
+ p[0] = 40;
+ EXPECT_EQ((std::vector<int>{40, 11, 12, 13, 14}), p.values());
+
+ for (auto&& elem: util::make_range(p.equal_range(3.))) {
+ elem = 7;
+ }
+ EXPECT_EQ((std::vector<int>{40, 7, 7, 7, 14}), p.values());
+
+ p(3.) = 50; // right most element intersecting 3.
+ EXPECT_EQ((std::vector<int>{40, 7, 7, 50, 14}), p.values());
+
+ p(3.).value = 60;
+ EXPECT_EQ((std::vector<int>{40, 7, 7, 60, 14}), p.values());
+
+ ASSERT_TRUE(std::is_const_v<decltype(p[0].extent)>);
+}
+
+TEST(piecewise, map) {
+ double xx[5] = {1, 2.25, 3.25, 3.5, 4.};
+ pw_elements<int> p(xx, (int [4]){3, 4, 5, 6});
+
+ auto void_fn = [](auto) {};
+ pw_elements<> void_expected(xx);
+ EXPECT_EQ(void_expected, pw_map(p, void_fn));
+
+ auto str_fn = [](int j) { return std::string(j, '.'); };
+ pw_elements<std::string> str_expected(xx, (const char* [4]){"...", "....", ".....", "......"});
+ EXPECT_EQ(str_expected, pw_map(p, str_fn));
+}
+
TEST(piecewise, zip) {
pw_elements<int> p03;
p03.assign((double [3]){0., 1.5, 3.}, (int [2]){10, 11});
@@ -343,42 +387,57 @@ TEST(piecewise, zip) {
p14.assign((double [5]){1, 2.25, 3.25, 3.5, 4.}, (int [4]){3, 4, 5, 6});
using ii = std::pair<int, int>;
- pw_elements<ii> p03_14 = zip(p03, p14);
+ using pipi = std::pair<pw_element<int>, pw_element<int>>;
+
+ // Zipped elements are pairs of pw_element.
+ pw_elements<pipi> p03_14_pw = pw_zip(p03, p14);
+ EXPECT_EQ(1., p03_14_pw.bounds().first);
+ EXPECT_EQ(3., p03_14_pw.bounds().second);
+
+ EXPECT_EQ((std::vector<double>{1., 1.5, 2.25, 3.}), p03_14_pw.vertices());
+ EXPECT_EQ((std::vector<pipi>{{p03[0], p14[0]}, {p03[1], p14[0]}, {p03[1], p14[1]}}),
+ p03_14_pw.values());
+
+ // To get pairs of just the element values, use pw_zip_with (with the default
+ // pw_pairify map).
+
+ pw_elements<ii> p03_14 = pw_zip_with(p03, p14);
EXPECT_EQ(1., p03_14.bounds().first);
EXPECT_EQ(3., p03_14.bounds().second);
EXPECT_EQ((std::vector<double>{1., 1.5, 2.25, 3.}), p03_14.vertices());
- EXPECT_EQ((std::vector<ii>{ii(10, 3), ii(11, 3), ii(11, 4)}), p03_14.elements());
+ EXPECT_EQ((std::vector<ii>{{10, 3}, {11, 3}, {11, 4}}), p03_14.values());
- pw_elements<ii> p14_03 = zip(p14, p03);
+ pw_elements<ii> p14_03 = pw_zip_with(p14, p03);
EXPECT_EQ(p03_14.vertices(), p14_03.vertices());
- std::vector<ii> flipped = util::assign_from(util::transform_view(p14_03.elements(),
+ std::vector<ii> flipped = util::assign_from(util::transform_view(p14_03.values(),
[](ii p) { return ii{p.second, p.first}; }));
- EXPECT_EQ(p03_14.elements(), flipped);
+ EXPECT_EQ(p03_14.values(), flipped);
pw_elements<> v03;
v03.assign((double [3]){0., 1.5, 3.});
- EXPECT_EQ((std::vector<int>{3, 3, 4}), zip(v03, p14).elements());
- EXPECT_EQ((std::vector<int>{3, 3, 4}), zip(p14, v03).elements());
+ EXPECT_EQ((std::vector<int>{3, 3, 4}), pw_zip_with(v03, p14).values());
+ EXPECT_EQ((std::vector<int>{3, 3, 4}), pw_zip_with(p14, v03).values());
- EXPECT_EQ((std::vector<double>{1., 1.5, 2.25, 3.}), zip(v03, p14).vertices());
- EXPECT_EQ((std::vector<double>{1., 1.5, 2.25, 3.}), zip(p14, v03).vertices());
+ EXPECT_EQ((std::vector<double>{1., 1.5, 2.25, 3.}), pw_zip_with(v03, p14).vertices());
+ EXPECT_EQ((std::vector<double>{1., 1.5, 2.25, 3.}), pw_zip_with(p14, v03).vertices());
- auto project = [](double l, double r, pw_element<void>, const pw_element<int>& b) -> double {
- double b_width = b.interval.second-b.interval.first;
- return b.element*(r-l)/b_width;
+ auto project = [](std::pair<double, double> extent, pw_element<void>, const pw_element<int>& b) -> double {
+ auto [l, r] = extent;
+ double b_width = b.extent.second-b.extent.first;
+ return b.value*(r-l)/b_width;
};
pw_elements<void> vxx; // elements cover bounds of p14
vxx.assign((double [6]){0.2, 1.7, 1.95, 2.325, 2.45, 4.9});
- pw_elements<double> pxx = zip(vxx, p14, project);
- double p14_sum = util::sum(util::transform_view(p14, [](auto v) { return v.element; }));
- double pxx_sum = util::sum(util::transform_view(pxx, [](auto v) { return v.element; }));
- EXPECT_DOUBLE_EQ(p14_sum, pxx_sum);
+ pw_elements<double> pxx = pw_zip_with(vxx, p14, project);
+ double p14_sum = util::sum(util::transform_view(p14, [](auto&& v) { return v.value; }));
+ double pxx_sum = util::sum(util::transform_view(pxx, [](auto&& v) { return v.value; }));
+ EXPECT_DOUBLE_EQ(p14_sum, pxx_sum);
}
TEST(piecewise, zip_zero_length_elements) {
@@ -398,71 +457,70 @@ TEST(piecewise, zip_zero_length_elements) {
using ii = std::pair<int, int>;
{
- pw_elements<ii> zz = zip(p03a, p03b);
+ pw_elements<ii> zz = pw_zip_with(p03a, p03b);
EXPECT_EQ(0., zz.bounds().first);
EXPECT_EQ(3., zz.bounds().second);
std::vector<double> expected_vertices = {0, 0, 0, 1, 1.5, 3, 3, 3};
- std::vector<ii> expected_elements = {ii(10, 20), ii(11, 21), ii(11,22), ii(11,23), ii(12, 23), ii(13,24), ii(13,25)};
+ std::vector<ii> expected_values = {ii(10, 20), ii(11, 21), ii(11,22), ii(11,23), ii(12, 23), ii(13,24), ii(13,25)};
EXPECT_EQ(expected_vertices, zz.vertices());
- EXPECT_EQ(expected_elements, zz.elements());
+ EXPECT_EQ(expected_values, zz.values());
- pw_elements<ii> yy = zip(p03b, p03a);
- flip(expected_elements);
+ pw_elements<ii> yy = pw_zip_with(p03b, p03a);
+ flip(expected_values);
EXPECT_EQ(expected_vertices, yy.vertices());
- EXPECT_EQ(expected_elements, yy.elements());
+ EXPECT_EQ(expected_values, yy.values());
}
{
- pw_elements<ii> zz = zip(p03a, p33);
+ pw_elements<ii> zz = pw_zip_with(p03a, p33);
EXPECT_EQ(3., zz.bounds().first);
EXPECT_EQ(3., zz.bounds().second);
std::vector<double> expected_vertices = {3, 3, 3};
- std::vector<ii> expected_elements = {ii(12, 30), ii(13, 31)};
+ std::vector<ii> expected_values = {ii(12, 30), ii(13, 31)};
EXPECT_EQ(expected_vertices, zz.vertices());
- EXPECT_EQ(expected_elements, zz.elements());
+ EXPECT_EQ(expected_values, zz.values());
- pw_elements<ii> yy = zip(p33, p03a);
- flip(expected_elements);
+ pw_elements<ii> yy = pw_zip_with(p33, p03a);
+ flip(expected_values);
EXPECT_EQ(expected_vertices, yy.vertices());
- EXPECT_EQ(expected_elements, yy.elements());
-
+ EXPECT_EQ(expected_values, yy.values());
}
{
- pw_elements<ii> zz = zip(p03a, p14);
+ pw_elements<ii> zz = pw_zip_with(p03a, p14);
EXPECT_EQ(1., zz.bounds().first);
EXPECT_EQ(3., zz.bounds().second);
std::vector<double> expected_vertices = {1, 1.5, 2, 3, 3};
- std::vector<ii> expected_elements = {ii(11, 40), ii(12, 40), ii(12, 41), ii(13, 41)};
+ std::vector<ii> expected_values = {ii(11, 40), ii(12, 40), ii(12, 41), ii(13, 41)};
EXPECT_EQ(expected_vertices, zz.vertices());
- EXPECT_EQ(expected_elements, zz.elements());
+ EXPECT_EQ(expected_values, zz.values());
- pw_elements<ii> yy = zip(p14, p03a);
- flip(expected_elements);
+ pw_elements<ii> yy = pw_zip_with(p14, p03a);
+ flip(expected_values);
EXPECT_EQ(expected_vertices, yy.vertices());
- EXPECT_EQ(expected_elements, yy.elements());
+ EXPECT_EQ(expected_values, yy.values());
}
{
// Check void version too!
pw_elements<> v03a(p03a), v03b(p03b);
- pw_elements<> zz = zip(v03a, v03b);
+ pw_elements<> zz = pw_zip_with(v03a, v03b);
EXPECT_EQ(0., zz.bounds().first);
EXPECT_EQ(3., zz.bounds().second);
std::vector<double> expected_vertices = {0, 0, 0, 1, 1.5, 3, 3, 3};
EXPECT_EQ(expected_vertices, zz.vertices());
- pw_elements<> yy = zip(v03b, v03a);
+ pw_elements<> yy = pw_zip_with(v03b, v03a);
EXPECT_EQ(expected_vertices, yy.vertices());
}
}
diff --git a/test/unit/test_ratelem.cpp b/test/unit/test_ratelem.cpp
index 2f21777f..d16e1fb3 100644
--- a/test/unit/test_ratelem.cpp
+++ b/test/unit/test_ratelem.cpp
@@ -1,6 +1,6 @@
#include "../gtest.h"
-#include <list>
+#include <cmath>
#include "util/ratelem.hpp"
#include "util/rangeutil.hpp"
@@ -115,3 +115,33 @@ TYPED_TEST(ratelem_pq, interpolate_monotonic) {
}
}
}
+
+TEST(ratelem, p1q1singular) {
+ // Check special case behaviour for p==1, q==1 when interpolants
+ // are strictly monotonic but possibly infinite.
+
+ auto f1 = [](double x) { return (1-x)/x; };
+ rat_element<1, 1> r1(f1);
+
+ for (unsigned i = 0; i<=4; ++i) {
+ double x = (double)i/4.0;
+ EXPECT_DOUBLE_EQ(f1(x), r1(x));
+ }
+
+ auto f2 = [](double x) { return x/(1-x); };
+ rat_element<1, 1> r2(f2);
+
+ for (unsigned i = 0; i<=4; ++i) {
+ double x = (double)i/4.0;
+ EXPECT_DOUBLE_EQ(f2(x), r2(x));
+ }
+
+ // With p==1, q==1, all infinite node values should
+ // give us NaN when we try to evaluate.
+
+ rat_element<1, 1> nope(INFINITY, INFINITY, INFINITY);
+ for (unsigned i = 0; i<=4; ++i) {
+ double x = (double)i/4.0;
+ EXPECT_TRUE(std::isnan(nope(x)));
+ }
+}
--
GitLab