diff --git a/arbor/fvm_layout.cpp b/arbor/fvm_layout.cpp
index 19a68d2eb26a990b8858b359216a2c75c8fd8b2f..dc6364325399b3f7a604a5d95a0317e4334bc241 100644
--- a/arbor/fvm_layout.cpp
+++ b/arbor/fvm_layout.cpp
@@ -1,3 +1,4 @@
+#include <algorithm>
#include <set>
#include <stdexcept>
#include <unordered_set>
@@ -19,6 +20,7 @@
#include "util/piecewise.hpp"
#include "util/rangeutil.hpp"
#include "util/transform.hpp"
+#include "util/unique.hpp"
namespace arb {
@@ -30,6 +32,7 @@ using util::pw_elements;
using util::pw_element;
using util::sort;
using util::sort_by;
+using util::stable_sort_by;
using util::value_by_key;
namespace {
@@ -124,43 +127,52 @@ pw_elements<U> pw_over_cable(const mcable_map<T>& mm, mcable cable, U dflt_value
// Construct cv_geometry for cell from locset describing CV boundary points.
cv_geometry cv_geometry_from_ends(const cable_cell& cell, const locset& lset) {
- struct mloc_hash {
- std::size_t operator()(const mlocation& loc) const { return loc.branch ^ std::hash<double>()(loc.pos); }
- };
- using mlocation_map = std::unordered_map<mlocation, fvm_size_type, mloc_hash>;
auto pop = [](auto& vec) { auto h = vec.back(); return vec.pop_back(), h; };
cv_geometry geom;
const auto& mp = cell.provider();
const auto& m = mp.morphology();
- if (mp.morphology().empty()) {
+ if (m.empty()) {
geom.cell_cv_divs = {0, 0};
return geom;
}
- auto canon = [&m](mlocation loc) { return canonical(m, loc); };
- auto origin = [&canon](mcable cab) { return canon(mlocation{cab.branch, cab.prox_pos}); };
- auto terminus = [&canon](mcable cab) { return canon(mlocation{cab.branch, cab.dist_pos}); };
-
mlocation_list locs = thingify(lset, mp);
- std::vector<msize_t> n_cv_cables;
- std::vector<mlocation> next_cv_head;
- next_cv_head.push_back({mnpos, 0});
+ // Filter out root, terminal locations and repeated locations so as to
+ // avoid trivial CVs outside of fork points. (This is not necessary for
+ // correctness, but is for the convenience of specification by lset.)
- mcable_list cables, all_cables;
- std::vector<msize_t> branches;
+ auto neither_root_nor_terminal = [&m](mlocation x) {
+ return !(x.pos==0 && x.branch==(m.branch_children(mnpos).size()>1u? mnpos: 0)) // root?
+ && !(x.pos==1 && m.branch_children(x.branch).empty()); // terminal?
+ };
+ locs.erase(std::partition(locs.begin(), locs.end(), neither_root_nor_terminal), locs.end());
+ util::sort(locs);
+ util::unique_in_place(locs);
- mlocation_map head_count;
- unsigned extra_cv_count = 0;
+ // Collect cables constituting each CV, maintaining a stack of CV
+ // proximal 'head' points, and recursing down branches in the morphology
+ // within each CV.
+
+ constexpr fvm_index_type no_parent = -1;
+ std::vector<std::pair<mlocation, fvm_index_type>> next_cv_head; // head loc, parent cv index
+ next_cv_head.emplace_back(mlocation{mnpos, 0}, no_parent);
+
+ mcable_list cables;
+ std::vector<msize_t> branches;
+ geom.cv_cables_divs.push_back(0);
+ fvm_index_type cv_index = 0;
while (!next_cv_head.empty()) {
- mlocation h = pop(next_cv_head);
+ auto next = pop(next_cv_head);
+ mlocation h = next.first;
cables.clear();
branches.clear();
branches.push_back(h.branch);
+ geom.cv_parent.push_back(next.second);
while (!branches.empty()) {
msize_t b = pop(branches);
@@ -178,7 +190,7 @@ cv_geometry cv_geometry_from_ends(const cable_cell& cell, const locset& lset) {
// Otherwise, recurse over child branches.
if (it!=locs.end() && it->branch==b) {
cables.push_back({b, b==h.branch? h.pos: 0, it->pos});
- next_cv_head.push_back(*it);
+ next_cv_head.emplace_back(*it, cv_index);
}
else {
if (b!=mnpos) {
@@ -190,69 +202,33 @@ cv_geometry cv_geometry_from_ends(const cable_cell& cell, const locset& lset) {
}
}
- auto empty_cable = [](mcable c) { return c.prox_pos==c.dist_pos; };
- cables.erase(std::remove_if(cables.begin(), cables.end(), empty_cable), cables.end());
-
- if (!cables.empty()) {
- if (++head_count[origin(cables.front())]==2) {
- ++extra_cv_count;
- }
-
- n_cv_cables.push_back(cables.size());
- sort(cables);
- append(all_cables, std::move(cables));
- }
- }
-
- geom.cv_cables.reserve(all_cables.size()+extra_cv_count);
- geom.cv_parent.reserve(n_cv_cables.size()+extra_cv_count);
- geom.cv_cables_divs.reserve(n_cv_cables.size()+extra_cv_count+1);
- geom.cv_cables_divs.push_back(0);
-
- mlocation_map parent_map;
-
- unsigned all_cables_index = 0;
- unsigned cv_index = 0;
-
- // Multiple CVs meeting at (0,0)?
- mlocation root{0, 0};
- unsigned n_top_children = value_by_key(head_count, root).value_or(0);
- if (n_top_children>1) {
- // Add initial trical CV.
- geom.cv_parent.push_back(mnpos);
- geom.cv_cables.push_back(mcable{0, 0, 0});
+ sort(cables);
+ append(geom.cv_cables, std::move(cables));
geom.cv_cables_divs.push_back(geom.cv_cables.size());
- parent_map[root] = cv_index++;
+ ++cv_index;
}
- for (auto n_cables: n_cv_cables) {
- mlocation head = origin(all_cables[all_cables_index]);
- msize_t parent_cv = value_by_key(parent_map, head).value_or(mnpos);
+ auto n_cv = cv_index;
+ arb_assert(n_cv>0);
+ arb_assert(geom.cv_parent.front()==-1);
+ arb_assert(util::all_of(util::subrange_view(geom.cv_parent, 1, n_cv),
+ [](auto v) { return v!=no_parent; }));
- auto cables = util::subrange_view(all_cables, all_cables_index, all_cables_index+n_cables);
- std::copy(cables.begin(), cables.end(), std::back_inserter(geom.cv_cables));
+ // Construct CV children mapping by sorting CV indices by parent.
+ assign(geom.cv_children, util::make_span(1, n_cv));
+ stable_sort_by(geom.cv_children, [&geom](auto cv) { return geom.cv_parent[cv]; });
- geom.cv_parent.push_back(parent_cv);
- geom.cv_cables_divs.push_back(geom.cv_cables.size());
+ geom.cv_children_divs.reserve(n_cv+1);
+ geom.cv_children_divs.push_back(0);
- auto this_cv = cv_index++;
- for (auto cable: cables) {
- mlocation term = terminus(cable);
-
- unsigned n_children = value_by_key(head_count, term).value_or(0);
- if (n_children>1) {
- // Add trivial CV for lindep.
- geom.cv_parent.push_back(this_cv);
- geom.cv_cables.push_back(mcable{cable.branch, 1., 1.});
- geom.cv_cables_divs.push_back((fvm_index_type)geom.cv_cables.size());
- parent_map[term] = cv_index++;
- }
- else {
- parent_map[term] = this_cv;
- }
- }
+ auto b = geom.cv_children.begin();
+ auto e = geom.cv_children.end();
+ auto from = b;
- all_cables_index += n_cables;
+ for (fvm_index_type cv = 0; cv<n_cv; ++cv) {
+ from = std::partition_point(from, e,
+ [cv, &geom](auto i) { return geom.cv_parent[i]<=cv; });
+ geom.cv_children_divs.push_back(from-b);
}
// Fill cv/cell mapping for single cell (index 0).
@@ -263,17 +239,14 @@ cv_geometry cv_geometry_from_ends(const cable_cell& cell, const locset& lset) {
geom.branch_cv_map.resize(1);
std::vector<pw_elements<fvm_size_type>>& bmap = geom.branch_cv_map.back();
- for (auto cv: util::make_span(geom.size())) {
+ for (auto cv: util::make_span(n_cv)) {
for (auto cable: geom.cables(cv)) {
if (cable.branch>=bmap.size()) {
bmap.resize(cable.branch+1);
}
// Ordering of CV ensures CV cables on any given branch are found sequentially.
- // Omit empty cables.
- if (cable.prox_pos<cable.dist_pos) {
- bmap[cable.branch].push_back(cable.prox_pos, cable.dist_pos, cv);
- }
+ bmap[cable.branch].push_back(cable.prox_pos, cable.dist_pos, cv);
}
}
@@ -327,6 +300,9 @@ cv_geometry& append(cv_geometry& geom, const cv_geometry& right) {
append_divs(geom.cv_cables_divs, right.cv_cables_divs);
append_offset(geom.cv_parent, geom_n_cv, right.cv_parent);
+ append_offset(geom.cv_children, geom_n_cv, right.cv_children);
+ append_divs(geom.cv_children_divs, right.cv_children_divs);
+
append_offset(geom.cv_to_cell, geom_n_cell, right.cv_to_cell);
append_divs(geom.cell_cv_divs, right.cell_cv_divs);
@@ -345,6 +321,8 @@ fvm_cv_discretization& append(fvm_cv_discretization& dczn, const fvm_cv_discreti
append(dczn.temperature_K, right.temperature_K);
append(dczn.diam_um, right.diam_um);
+ append(dczn.axial_resistivity, right.axial_resistivity);
+
return dczn;
}
@@ -372,13 +350,21 @@ fvm_cv_discretization fvm_cv_discretize(const cable_cell& cell, const cable_cell
double dflt_potential = *(dflt.init_membrane_potential | global_dflt.init_membrane_potential);
double dflt_temperature = *(dflt.temperature_K | global_dflt.temperature_K);
+ D.axial_resistivity.resize(1);
+ msize_t n_branch = D.geometry.n_branch(0);
+ D.axial_resistivity[0].reserve(n_branch);
+ for (msize_t i = 0; i<n_branch; ++i) {
+ D.axial_resistivity[0].push_back(pw_over_cable(cell.region_assignments().get<axial_resistivity>(),
+ mcable{i, 0., 1.}, dflt_resistivity));
+ }
+
const auto& embedding = cell.embedding();
for (auto i: count_along(D.geometry.cv_parent)) {
auto cv_cables = D.geometry.cables(i);
// Computing face_conductance:
//
- // Flux between adjacemt CVs is computed as if there were no membrane currents, and with the CV voltage
+ // Flux between adjacent CVs is computed as if there were no membrane currents, and with the CV voltage
// values taken to be exact at a reference point in each CV:
// * If the CV is unbranched, the reference point is taken to be the CV midpoint.
// * If the CV is branched, the reference point is taken to be closest branch point to
@@ -408,8 +394,7 @@ fvm_cv_discretization fvm_cv_discretize(const cable_cell& cell, const cable_cell
}
mcable span{bid, parent_refpt, cv_refpt};
- double resistance = embedding.integrate_ixa(bid,
- pw_over_cable(cell.region_assignments().get<axial_resistivity>(), span, dflt_resistivity));
+ double resistance = embedding.integrate_ixa(span, D.axial_resistivity[0].at(bid));
D.face_conductance[i] = 100/resistance; // 100 scales to µS.
}
@@ -718,7 +703,7 @@ fvm_mechanism_data fvm_build_mechanism_data(const cable_cell_global_properties&
std::sort(in.param_value.begin(), in.param_value.end());
in.target_index = pm.lid;
- in.cv = D.geometry.location_cv(cell_idx, pm.loc);
+ in.cv = D.geometry.location_cv(cell_idx, pm.loc, cv_prefer::cv_nonempty);
sl.push_back(std::move(in));
}
@@ -781,7 +766,8 @@ fvm_mechanism_data fvm_build_mechanism_data(const cable_cell_global_properties&
const auto& stimuli = cell.stimuli();
std::vector<size_type> stimuli_cv;
- assign_by(stimuli_cv, stimuli, [&D, cell_idx](auto& p) { return D.geometry.location_cv(cell_idx, p.loc); });
+ assign_by(stimuli_cv, stimuli, [&D, cell_idx](auto& p) {
+ return D.geometry.location_cv(cell_idx, p.loc, cv_prefer::cv_nonempty); });
std::vector<size_type> cv_order;
assign(cv_order, count_along(stimuli));
diff --git a/arbor/fvm_layout.hpp b/arbor/fvm_layout.hpp
index 1c5566f81482c94f8941ed2f5bf49f7d3deaa418..434ba6f5e80363941fe43f9e064e8480e5a3529d 100644
--- a/arbor/fvm_layout.hpp
+++ b/arbor/fvm_layout.hpp
@@ -4,12 +4,12 @@
#include <utility>
#include <vector>
-#include <arbor/fvm_types.hpp>
#include <arbor/cable_cell.hpp>
#include <arbor/mechanism.hpp>
#include <arbor/mechinfo.hpp>
#include <arbor/mechcat.hpp>
#include <arbor/recipe.hpp>
+#include <arbor/util/optional.hpp>
#include "execution_context.hpp"
#include "util/piecewise.hpp"
@@ -19,15 +19,52 @@
namespace arb {
// CV geometry as determined by per-cell CV boundary points.
+//
+// Details of CV cable representation:
+//
+// * The extent of the cables of a control volume corresponds to the
+// closure of the control volume on the morphology tree.
+//
+// * Every fork in the morphology tree 'belongs' to exactly
+// one CV. A fork belongs to a CV if and only if there is
+// a cable in the CV for each branch of the fork, that
+// includes the fork point.
+//
+// * If a CV has more than one cable covering a fork point, then
+// that fork point must belong to that CV.
+//
+// These requirements are a consequence of the CV geometry being determined
+// by a collection of CV boundarary locations.
+
+namespace cv_prefer {
+ // Enum for resolving which CV to return on location look-up, if the
+ // location is on a CV boundary.
+
+ enum type {
+ // Prefer more proximal CV:
+ cv_proximal,
+ // Prefer more distal CV:
+ cv_distal,
+ // Prefer distal CV unless it has zero extent on location branch.
+ // This should be used for placing point processes on CVs.
+ cv_nonempty,
+ // Prefer distal CV unless the proximal CV has zero extent on location branch.
+ // This should be used for determing to which CV a fork point belongs.
+ cv_empty
+ };
+}
struct cv_geometry {
using size_type = fvm_size_type;
using index_type = fvm_index_type;
std::vector<mcable> cv_cables; // CV unbranched sections, partitioned by CV.
- std::vector<index_type> cv_cables_divs; // Partitions cv_cables by CV index.
+ std::vector<index_type> cv_cables_divs; // Partitions cv_cables by CV index.
std::vector<index_type> cv_parent; // Index of CV parent or size_type(-1) for a cell root CV.
+ std::vector<index_type> cv_children; // CV child indices, partitioned by CV, and then in order.
+ std::vector<index_type> cv_children_divs; // Paritions cv_children by CV index.
+
std::vector<index_type> cv_to_cell; // Maps CV index to cell index.
std::vector<index_type> cell_cv_divs; // Partitions CV indices by cell.
@@ -39,6 +76,11 @@ struct cv_geometry {
return util::subrange_view(cv_cables, partn[cv_index]);
}
+ auto children(size_type cv_index) const {
+ auto partn = util::partition_view(cv_children_divs);
+ return util::subrange_view(cv_children, partn[cv_index]);
+ }
+
std::pair<index_type, index_type> cell_cv_interval(size_type cell_idx) const {
auto partn = util::partition_view(cell_cv_divs);
return partn[cell_idx];
@@ -68,8 +110,43 @@ struct cv_geometry {
return cell_cv_divs.empty()? 0: cell_cv_divs.size()-1;
}
- size_type location_cv(size_type cell_idx, mlocation loc) const {
- return cell_cv_divs.at(cell_idx)+branch_cv_map.at(cell_idx).at(loc.branch)(loc.pos).second;
+ size_type n_branch(size_type cell_idx) const {
+ return branch_cv_map.at(cell_idx).size();
+ }
+
+ 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;
+ };
+
+ 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;
+
+ // index_of() should have returned right-most matching interval.
+ arb_assert(i==i_max || loc.pos<pw_cv_offset.interval(i+1).first);
+
+ using namespace cv_prefer;
+ switch (prefer) {
+ case cv_distal:
+ break;
+ case cv_proximal:
+ if (loc.pos==cv_prox && i>0) --i;
+ break;
+ case cv_nonempty:
+ if (zero_extent(i)) {
+ if (i>0 && !zero_extent(i-1)) --i;
+ else if (i<i_max && !zero_extent(i+1)) ++i;
+ }
+ break;
+ case cv_empty:
+ if (loc.pos==cv_prox && i>0 && zero_extent(i-1)) --i;
+ break;
+ }
+
+ index_type cv_base = cell_cv_divs.at(cell_idx);
+ return cv_base+pw_cv_offset[i].second;
}
};
@@ -86,6 +163,14 @@ cv_geometry cv_geometry_from_ends(const cable_cell& cell, const locset& lset);
// diam_um is taken to be the diameter of a CV with constant diameter and same
// extent which has the same surface area (i.e. cv_area/(Ï€L) where L is the
// total length of the cables comprising the CV.)
+//
+// The bulk conductivity over the morphology is recorded here as well, for
+// calculating voltage and axial current interpolating probes.
+//
+// For the computation of inter-CV conductances and voltage interpolation, it
+// is assumed that the CV voltage is exact at every fork point that belongs
+// to the CV, or in the absence of any internal forks, is exact at the
+// midpoint of an unbranched CV.
struct fvm_cv_discretization {
using size_type = fvm_size_type;
@@ -98,12 +183,16 @@ struct fvm_cv_discretization {
size_type size() const { return geometry.size(); }
size_type n_cell() const { return geometry.n_cell(); }
+ // Following members have one element per CV.
std::vector<value_type> face_conductance; // [µS]
std::vector<value_type> cv_area; // [µm²]
std::vector<value_type> cv_capacitance; // [pF]
std::vector<value_type> init_membrane_potential; // [mV]
std::vector<value_type> temperature_K; // [K]
std::vector<value_type> diam_um; // [µm]
+
+ // For each cell, one piece-wise constant value per branch.
+ std::vector<std::vector<pw_constant_fn>> axial_resistivity; // [Ω·cm]
};
// Combine two fvm_cv_geometry groups in-place.
@@ -111,7 +200,6 @@ struct fvm_cv_discretization {
fvm_cv_discretization& append(fvm_cv_discretization&, const fvm_cv_discretization&);
// Construct fvm_cv_discretization from one or more cells.
-
fvm_cv_discretization fvm_cv_discretize(const cable_cell& cell, const cable_cell_parameter_set& global_dflt);
fvm_cv_discretization fvm_cv_discretize(const std::vector<cable_cell>& cells, const cable_cell_parameter_set& global_defaults, const arb::execution_context& ctx={});
diff --git a/arbor/fvm_lowered_cell_impl.hpp b/arbor/fvm_lowered_cell_impl.hpp
index 23ecd5035d56caeea5f8c8900959f8ef7456047e..0f54936bb7220964aa3c8d1ef0b1d577ff52ea95 100644
--- a/arbor/fvm_lowered_cell_impl.hpp
+++ b/arbor/fvm_lowered_cell_impl.hpp
@@ -520,7 +520,7 @@ void fvm_lowered_cell_impl<B>::initialize(
cell_gid_type gid = gids[cell_idx];
for (auto entry: cells[cell_idx].detectors()) {
- detector_cv.push_back(D.geometry.location_cv(cell_idx, entry.loc));
+ detector_cv.push_back(D.geometry.location_cv(cell_idx, entry.loc, cv_prefer::cv_empty));
detector_threshold.push_back(entry.item.threshold);
}
@@ -528,14 +528,16 @@ void fvm_lowered_cell_impl<B>::initialize(
probe_info pi = rec.get_probe({gid, j});
auto where = any_cast<cell_probe_address>(pi.address);
- auto cv = D.geometry.location_cv(cell_idx, where.location);
+ fvm_size_type cv;
probe_handle handle;
switch (where.kind) {
case cell_probe_address::membrane_voltage:
+ cv = D.geometry.location_cv(cell_idx, where.location, cv_prefer::cv_empty);
handle = state_->voltage.data()+cv;
break;
case cell_probe_address::membrane_current:
+ cv = D.geometry.location_cv(cell_idx, where.location, cv_prefer::cv_nonempty);
handle = state_->current_density.data()+cv;
break;
default:
@@ -568,7 +570,7 @@ std::vector<fvm_gap_junction> fvm_lowered_cell_impl<B>::fvm_gap_junctions(
const auto& cell_gj = cells[cell_idx].gap_junction_sites();
for (auto gj : cell_gj) {
- auto cv = D.geometry.location_cv(cell_idx, gj.loc);
+ auto cv = D.geometry.location_cv(cell_idx, gj.loc, cv_prefer::cv_nonempty);
gid_to_cvs[gids[cell_idx]].push_back(cv);
}
}
diff --git a/arbor/include/arbor/morph/embed_pwlin.hpp b/arbor/include/arbor/morph/embed_pwlin.hpp
index 0338186bb31fa1b45cde63100ec2f7522c019cb1..0a352a5708029cfe2d5ee7e40ef81f3caa8d3aa7 100644
--- a/arbor/include/arbor/morph/embed_pwlin.hpp
+++ b/arbor/include/arbor/morph/embed_pwlin.hpp
@@ -36,15 +36,21 @@ struct embed_pwlin {
// Computed length of mcable.
double integrate_length(mcable c) const;
double integrate_length(mlocation proxmal, mlocation distal) const;
+
+ double integrate_length(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(mcable c) const;
double integrate_area(mlocation proxmal, mlocation distal) const;
+
+ double integrate_area(mcable c, const pw_constant_fn&) const;
double integrate_area(msize_t bid, const pw_constant_fn&) const;
// Integrated inverse cross-sectional area of given mcable.
double integrate_ixa(mcable c) const;
+
+ double integrate_ixa(mcable c, const pw_constant_fn&) const;
double integrate_ixa(msize_t bid, const pw_constant_fn&) const;
// Length of whole branch.
diff --git a/arbor/morph/embed_pwlin.cpp b/arbor/morph/embed_pwlin.cpp
index eaf430a6f42bafc9fde9f7ae8d0c28eb30f059ad..ae3d0f2117a11f5d4425e4977c75361cc581e552 100644
--- a/arbor/morph/embed_pwlin.cpp
+++ b/arbor/morph/embed_pwlin.cpp
@@ -53,6 +53,36 @@ double integrate(const branch_pw_ratpoly<p, q>& f, unsigned bid, const pw_consta
return accum;
}
+// 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.
+
+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;
+
+ for (msize_t i = 0; i<g.size(); ++i) {
+ std::pair<double, double> interval = g.interval(i);
+
+ 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;
+}
+
template <typename operation>
mcable_list data_cmp(const branch_pw_ratpoly<1, 0>& f, unsigned bid, double val, operation op) {
mcable_list L;
@@ -109,19 +139,19 @@ double embed_pwlin::directed_projection(arb::mlocation loc) const {
return interpolate(data_->directed_projection, loc.branch, loc.pos);
}
-double embed_pwlin::integrate_length(msize_t bid, const pw_constant_fn& g) const {
- return integrate(data_->length, bid, g);
-}
+// Point to point integration:
-double embed_pwlin::integrate_area(msize_t bid, const pw_constant_fn& g) const {
- return integrate(data_->area, bid, g);
+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);
}
-double embed_pwlin::integrate_ixa(msize_t bid, const pw_constant_fn& g) const {
- return integrate(data_->ixa, bid, g);
+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);
}
-// Cable versions of integration methods:
+// Integrate over cable:
double embed_pwlin::integrate_length(mcable c) const {
return integrate_length(c.branch, pw_constant_fn{{c.prox_pos, c.dist_pos}, {1.}});
@@ -135,16 +165,32 @@ double embed_pwlin::integrate_ixa(mcable c) const {
return integrate_ixa(c.branch, pw_constant_fn{{c.prox_pos, c.dist_pos}, {1.}});
}
-// Point to point integration:
+// Integrate piecewise function over a branch:
-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);
+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(mlocation proximal, mlocation distal) const {
- return interpolate(data_->area, distal.branch, distal.pos) -
- interpolate(data_->area, proximal.branch, proximal.pos);
+double embed_pwlin::integrate_area(msize_t bid, const pw_constant_fn& g) const {
+ return integrate(data_->area, bid, g);
+}
+
+double embed_pwlin::integrate_ixa(msize_t bid, const pw_constant_fn& g) const {
+ return integrate(data_->ixa, bid, g);
+}
+
+// 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);
+}
+
+double embed_pwlin::integrate_area(mcable c, const pw_constant_fn& g) const {
+ return integrate(data_->area, c, g);
+}
+
+double embed_pwlin::integrate_ixa(mcable c, const pw_constant_fn& g) const {
+ return integrate(data_->ixa, c, g);
}
// Subregions defined by geometric inequalities:
diff --git a/arbor/morph/locset.cpp b/arbor/morph/locset.cpp
index 54ad51e090af6acedec4b91e6ddb9a1afd378d15..ad6a9a4e6a4c33ede2c0806e7b646acb11c706db 100644
--- a/arbor/morph/locset.cpp
+++ b/arbor/morph/locset.cpp
@@ -16,6 +16,7 @@
#include "util/transform.hpp"
#include "util/span.hpp"
#include "util/strprintf.hpp"
+#include "util/unique.hpp"
namespace arb {
namespace ls {
@@ -203,21 +204,6 @@ locset most_distal(region reg) {
return locset(most_distal_{std::move(reg)});
}
-template <typename X>
-void unique_in_place(std::vector<X>& v) {
- if (v.empty()) return;
-
- auto write = v.begin();
- auto read = write;
-
- while (++read!=v.end()) {
- if (*read==*write) continue;
- if (++write!=read) *write = std::move(*read);
- }
-
- v.erase(++write, v.end());
-}
-
mlocation_list thingify_(const most_distal_& n, const mprovider& p) {
mlocation_list L;
@@ -244,7 +230,7 @@ mlocation_list thingify_(const most_distal_& n, const mprovider& p) {
}
util::sort(L);
- unique_in_place(L);
+ util::unique_in_place(L);
return L;
}
diff --git a/arbor/util/unique.hpp b/arbor/util/unique.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..e55c122619157b8acb6aaa8efbdde3820c1d62eb
--- /dev/null
+++ b/arbor/util/unique.hpp
@@ -0,0 +1,30 @@
+#pragma once
+
+// Eliminated successive identical elements in-place in a SequenceConteiner.
+// (Can wrap or be replaced by std::unique in C++17.)
+//
+// If we ever implement a unique_view, it should go here too.
+
+#include <utility>
+
+namespace arb {
+namespace util {
+
+template <typename SeqContainer, typename EqPred = std::equal_to<>>
+void unique_in_place(SeqContainer& c, EqPred eq = EqPred{}) {
+ if (c.empty()) return;
+
+ auto end = c.end();
+ auto write = c.begin();
+ auto read = write;
+
+ while (++read!=end) {
+ if (eq(*read, *write)) continue;
+ if (++write!=read) *write = std::move(*read);
+ }
+
+ c.erase(++write, end);
+}
+
+} // namespace util
+} // namespace arb
diff --git a/test/unit/CMakeLists.txt b/test/unit/CMakeLists.txt
index f97bb2d3127688e8cefc639301a524b32552cdc5..d867d294b15d406f4f933ca16e181084b1507887 100644
--- a/test/unit/CMakeLists.txt
+++ b/test/unit/CMakeLists.txt
@@ -144,6 +144,7 @@ set(unit_sources
test_tree.cpp
test_transform.cpp
test_uninitialized.cpp
+ test_unique.cpp
test_unique_any.cpp
test_variant.cpp
test_vector.cpp
diff --git a/test/unit/test_cv_geom.cpp b/test/unit/test_cv_geom.cpp
index 0a0fdbaefa4f2df2f6ab6a94df9bb00f0ce4c192..63dca9639ff22219ac872fe958a3e787f12b19d7 100644
--- a/test/unit/test_cv_geom.cpp
+++ b/test/unit/test_cv_geom.cpp
@@ -11,16 +11,68 @@
#include "common.hpp"
#include "common_morphologies.hpp"
+#include "morph_pred.hpp"
#include "../common_cells.hpp"
using namespace arb;
using util::make_span;
+::testing::AssertionResult verify_cv_children(const cv_geometry& g) {
+ unsigned visited_children = 0;
+ for (unsigned i = 0; i<g.size(); ++i) {
+ if (!util::is_sorted(g.children(i))) {
+ return ::testing::AssertionFailure() << "CV " << i
+ << " has unsorted sequence of child CVs";
+ }
+
+ for (auto cv: g.children(i)) {
+ if ((fvm_index_type)i != g.cv_parent.at(cv)) {
+ return ::testing::AssertionFailure() << "CV " << i
+ << " has child CV " << cv << " which has parent " << g.cv_parent.at(cv);
+ }
+ ++visited_children;
+ }
+ }
+
+ if (g.cv_children.size()!=visited_children) {
+ return ::testing::AssertionFailure() << "geometry child CV count " << g.cv_children.size()
+ << " does not equal number of visited children " << visited_children;
+ }
+
+ unsigned n_nonempty_cells = 0;
+ for (auto c: util::make_span(g.n_cell())) {
+ if (!g.cell_cvs(c).empty()) {
+ ++n_nonempty_cells;
+ }
+ }
+ if (g.cv_children.size()!=g.size()-n_nonempty_cells) {
+ return ::testing::AssertionFailure() << "child CV count " << g.cv_children.size()
+ << " plus root CV count " << n_nonempty_cells
+ << " does not equal total number of CVs " << g.size();
+ }
+ return ::testing::AssertionSuccess();
+}
+
+namespace arb {
+namespace cv_prefer {
+std::ostream& operator<<(std::ostream& out, ::arb::cv_prefer::type p) {
+ switch (p) {
+ case cv_proximal: return out << "cv_proximal";
+ case cv_distal: return out << "cv_distal";
+ case cv_empty: return out << "cv_empty";
+ case cv_nonempty: return out << "cv_nonempty";
+ default: return out;
+ }
+}
+}
+}
+
TEST(cv_geom, empty) {
using namespace common_morphology;
cable_cell empty_cell{m_empty};
cv_geometry geom = cv_geometry_from_ends(empty_cell, ls::nil());
+ EXPECT_TRUE(verify_cv_children(geom));
EXPECT_TRUE(geom.cv_parent.empty());
EXPECT_TRUE(geom.cv_cables.empty());
@@ -44,11 +96,13 @@ TEST(cv_geom, trivial) {
cable_cell cell{p.second};
auto& m = cell.morphology();
-
// Equivalent ways of specifying one CV comprising whole cell:
cv_geometry geom1 = cv_geometry_from_ends(cell, ls::nil());
cv_geometry geom2 = cv_geometry_from_ends(cell, ls::terminal());
+ EXPECT_TRUE(verify_cv_children(geom1));
+ EXPECT_TRUE(verify_cv_children(geom2));
+
EXPECT_EQ(1u, geom1.size());
EXPECT_EQ(geom1.cv_cables, geom2.cv_cables);
@@ -56,6 +110,9 @@ TEST(cv_geom, trivial) {
cv_geometry geom3 = cv_geometry_from_ends(cell, ls::root());
cv_geometry geom4 = cv_geometry_from_ends(cell, join(ls::root(), ls::terminal()));
+ EXPECT_TRUE(verify_cv_children(geom3));
+ EXPECT_TRUE(verify_cv_children(geom4));
+
EXPECT_EQ(geom3.cv_cables, geom4.cv_cables);
if (m.branch_children(mnpos).size()==1) {
EXPECT_EQ(geom1.cv_cables, geom4.cv_cables);
@@ -76,7 +133,8 @@ TEST(cv_geom, one_cv_per_branch) {
cable_cell cell{p.second};
auto& m = cell.morphology();
- cv_geometry geom = cv_geometry_from_ends(cell, ls::on_branches(0));
+ cv_geometry geom = cv_geometry_from_ends(cell, sum(ls::on_branches(0), ls::on_branches(1)));
+ EXPECT_TRUE(verify_cv_children(geom));
// Expect trivial CVs at every fork point, and single-cable CVs for each branch.
std::vector<unsigned> seen_branches(m.num_branches(), 0);
@@ -84,10 +142,10 @@ TEST(cv_geom, one_cv_per_branch) {
for (auto i: make_span(geom.size())) {
auto cables = geom.cables(i);
- ASSERT_EQ(1u, cables.size());
auto c = cables.front();
if (c.prox_pos==c.dist_pos) {
+ EXPECT_LT(1u, cables.size());
if (c.branch==0 && c.prox_pos==0) {
EXPECT_TRUE(n_branch_child(mnpos)>1);
}
@@ -95,8 +153,11 @@ TEST(cv_geom, one_cv_per_branch) {
EXPECT_EQ(1., c.prox_pos);
EXPECT_TRUE(n_branch_child(c.branch)>1);
}
+ // Cables in trivial CV should be the same as those in the extent over the point.
+ EXPECT_TRUE(testing::seq_eq(mextent{m, mcable_list{c}}.cables(), cables));
}
else {
+ ASSERT_EQ(1u, cables.size());
++seen_branches[c.branch];
EXPECT_EQ(1., seen_branches[c.branch]);
EXPECT_EQ(0., c.prox_pos);
@@ -104,21 +165,9 @@ TEST(cv_geom, one_cv_per_branch) {
// Confirm parent CV is fork CV:
if (i>0) {
- mlocation pfork = canonical(m, mlocation{c.branch, 0.});
-
- auto pcables = geom.cables(geom.cv_parent[i]);
- ASSERT_EQ(1u, pcables.size());
-
- mcable p = pcables.front();
- EXPECT_EQ(pfork.branch, p.branch);
- EXPECT_EQ(p.prox_pos, p.dist_pos);
-
- if (p.branch==0) {
- EXPECT_TRUE(p.prox_pos==0 || p.prox_pos==1);
- }
- else {
- EXPECT_EQ(1., p.prox_pos);
- }
+ mextent fork_ext{m, mcable_list{{c.branch, 0}}};
+ mcable_list pcables = util::assign_from(geom.cables(geom.cv_parent[i]));
+ ASSERT_TRUE(testing::cablelist_eq(fork_ext.cables(), pcables));
}
}
}
@@ -139,6 +188,7 @@ TEST(cv_geom, midpoints) {
auto& m = cell.morphology();
cv_geometry geom = cv_geometry_from_ends(cell, ls::on_branches(0.5));
+ EXPECT_TRUE(verify_cv_children(geom));
// Expect CVs to be either: covering fork points, with one cable per branch
// at the fork (for a multiple-root-branch morphology, this would be treating
@@ -221,7 +271,8 @@ TEST(cv_geom, weird) {
// o
//
// CV 0 will comprise branches 0 and 2; CV 1 branches 1, 3, 5;
- // and CV 2 branch 4.
+ // and CV 2 branch 4. CV 0 will also cover the fork point (0,1);
+ // CV 1 will cover the fork point (1, 1).
using C = mcable;
using testing::seq_eq;
@@ -229,12 +280,13 @@ TEST(cv_geom, weird) {
cable_cell cell{common_morphology::m_reg_b6};
cv_geometry geom = cv_geometry_from_ends(cell, mlocation_list{{1, 0}, {4,0}});
+ EXPECT_TRUE(verify_cv_children(geom));
ASSERT_EQ(3u, geom.size());
- mcable_list expected0 = {C{0u, 0., 1.}, C{2u, 0., 1.}};
+ mcable_list expected0 = {C{0u, 0., 1.}, C{1u, 0., 0.}, C{2u, 0., 1.}};
EXPECT_TRUE(seq_eq(expected0, geom.cables(0)));
- mcable_list expected1 = {C{1u, 0., 1.}, C{3u, 0., 1.}, C{5u, 0., 1.}};
+ mcable_list expected1 = {C{1u, 0., 1.}, C{3u, 0., 1.}, C{4u, 0., 0.}, C{5u, 0., 1.}};
EXPECT_TRUE(seq_eq(expected1, geom.cables(1)));
mcable_list expected2 = {C{4u, 0., 1.}};
@@ -245,48 +297,143 @@ TEST(cv_geom, location_cv) {
using namespace common_morphology;
cable_cell cell{m_reg_b6};
+ auto& m = cell.morphology();
+
+ auto cv_extent = [&m](const cv_geometry& geom, auto cv) {
+ mcable_list cl;
+ util::assign(cl, geom.cables(cv));
+ return mextent{m, cl};
+ };
// Two CVs per branch, plus trivial CV at forks.
cv_geometry geom = cv_geometry_from_ends(cell,
- join(ls::on_branches(0.0), ls::on_branches(0.5), ls::on_branches(1.)));
+ join(ls::on_branches(0.), ls::on_branches(0.5), ls::on_branches(1.)));
- // Confirm CVs are all only one cable, and either trivial or half a branch.
+ // Confirm CVs are either trivial or a single cable covering half a branch.
for (auto cv: geom.cell_cvs(0)) {
auto cables = geom.cables(cv);
- ASSERT_EQ(1u, cables.size());
+ if (cables.size()==1u) {
+ // Half branch cable.
+ mcable cable = cables.front();
+ ASSERT_TRUE((cable.prox_pos==0 && cable.dist_pos==0.5 ) ||
+ (cable.prox_pos==0.5 && cable.dist_pos==1.));
+ }
+ else {
+ // Trivial CV over fork point.
+ mcable cable0 = cables.front();
+ ASSERT_TRUE(cable0.prox_pos==cable0.dist_pos);
- mcable cable = cables.front();
- ASSERT_TRUE(cable.prox_pos==cable.dist_pos ||
- (cable.prox_pos==0 && cable.dist_pos==0.5 )||
- (cable.prox_pos==0.5 && cable.dist_pos==1.));
+ mextent fork_ext{m, mcable_list{cable0}};
+ mcable_list clist = util::assign_from(cables);
+ ASSERT_TRUE(testing::cablelist_eq(fork_ext.cables(), clist));
+ }
}
- // Where ambiguous, CV found should be one with a non-trivial cable containing the location.
- for (auto bid: util::make_span(cell.morphology().num_branches())) {
- for (double pos: {0., 0.3, 0.5, 0.7, 1.}) {
- mlocation loc{bid, pos};
+ // For positions strictly within a CV extent, CV preference should make no difference.
+ for (auto prefer: {cv_prefer::cv_distal, cv_prefer::cv_proximal,
+ cv_prefer::cv_nonempty, cv_prefer::cv_empty}) {
+ SCOPED_TRACE(prefer);
+ for (auto bid: util::make_span(m.num_branches())) {
+ for (double pos: {0.3, 0.7}) {
+ mlocation loc{bid, pos};
+ SCOPED_TRACE(loc);
+
+ auto cv = geom.location_cv(0, loc, prefer);
+ ASSERT_TRUE(cv_extent(geom, cv).intersects(loc));
+
+ ASSERT_EQ(1u, geom.cables(cv).size());
+ mcable cable = geom.cables(cv).front();
+ EXPECT_TRUE(cable.branch==loc.branch);
+ EXPECT_TRUE(cable.prox_pos<cable.dist_pos);
+ }
+ }
+ }
+
+ // For positions in the middle of a branch, we should get distal CV unless
+ // CV prerence is `cv_proximal`.
+ for (auto prefer: {cv_prefer::cv_distal, cv_prefer::cv_proximal,
+ cv_prefer::cv_nonempty, cv_prefer::cv_empty}) {
+ SCOPED_TRACE(prefer);
+ for (auto bid: util::make_span(m.num_branches())) {
+ mlocation loc{bid, 0.5};
SCOPED_TRACE(loc);
- auto cv = geom.location_cv(0, loc);
+ auto cv = geom.location_cv(0, loc, prefer);
+ ASSERT_TRUE(cv_extent(geom, cv).intersects(loc));
+ ASSERT_EQ(1u, geom.cables(cv).size());
mcable cable = geom.cables(cv).front();
- EXPECT_TRUE(cable.prox_pos<cable.dist_pos);
- EXPECT_TRUE(cable.prox_pos<=pos);
- EXPECT_TRUE(cable.dist_pos>=pos);
EXPECT_TRUE(cable.branch==loc.branch);
+ EXPECT_TRUE(cable.prox_pos<cable.dist_pos);
+
+ if (prefer==cv_prefer::cv_proximal) {
+ EXPECT_EQ(0., cable.prox_pos);
+ }
+ else {
+ EXPECT_EQ(0.5, cable.prox_pos);
+ }
}
}
- // CV 0 will comprise branches 0 and 2; CV 1 branches 1, 3, 5;
- // and CV 2 branch 4 (see test cv_geom.weird above).
- cv_geometry gweird = cv_geometry_from_ends(cell, mlocation_list{{1, 0}, {4,0}});
-
- // Expect location (4, 0.) to be in CV 2, but colocated locations
- // (3, 0.), (5, 0.) and (1, 1.) to be in CV 1.
- EXPECT_EQ(2u, gweird.location_cv(0, mlocation{4, 0.}));
- EXPECT_EQ(1u, gweird.location_cv(0, mlocation{3, 0.}));
- EXPECT_EQ(1u, gweird.location_cv(0, mlocation{5, 0.}));
- EXPECT_EQ(1u, gweird.location_cv(0, mlocation{1, 1.}));
+ // For the head of a non-root branch, we should get the trivial CV over the
+ // fork for `cv_proximal` or `cv_empty`; otherwise the CV over the first
+ // half of the branch.
+ for (auto prefer: {cv_prefer::cv_distal, cv_prefer::cv_proximal,
+ cv_prefer::cv_nonempty, cv_prefer::cv_empty}) {
+ SCOPED_TRACE(prefer);
+ for (auto bid: util::make_span(m.num_branches())) {
+ if (m.branch_parent(bid)==mnpos) continue;
+
+ mlocation loc{bid, 0.};
+ SCOPED_TRACE(loc);
+
+ auto cv = geom.location_cv(0, loc, prefer);
+ ASSERT_TRUE(cv_extent(geom, cv).intersects(loc));
+
+ auto cables = geom.cables(cv);
+ switch (prefer) {
+ case cv_prefer::cv_proximal:
+ case cv_prefer::cv_empty:
+ EXPECT_NE(1u, cables.size());
+ break;
+ case cv_prefer::cv_distal:
+ case cv_prefer::cv_nonempty:
+ EXPECT_EQ(1u, cables.size());
+ EXPECT_EQ(0.5, cables.front().dist_pos);
+ break;
+ }
+ }
+ }
+
+ // For the tail of a non-terminal branch, we should get the trivial CV over the
+ // fork for `cv_distal` or `cv_empty`; otherwise the CV over the second
+ // half of the branch.
+ for (auto prefer: {cv_prefer::cv_distal, cv_prefer::cv_proximal,
+ cv_prefer::cv_nonempty, cv_prefer::cv_empty}) {
+ SCOPED_TRACE(prefer);
+ for (auto bid: util::make_span(m.num_branches())) {
+ if (m.branch_children(bid).empty()) continue;
+
+ mlocation loc{bid, 1.};
+ SCOPED_TRACE(loc);
+
+ auto cv = geom.location_cv(0, loc, prefer);
+ ASSERT_TRUE(cv_extent(geom, cv).intersects(loc));
+
+ auto cables = geom.cables(cv);
+ switch (prefer) {
+ case cv_prefer::cv_proximal:
+ case cv_prefer::cv_nonempty:
+ EXPECT_EQ(1u, cables.size());
+ EXPECT_EQ(0.5, cables.front().prox_pos);
+ break;
+ case cv_prefer::cv_distal:
+ case cv_prefer::cv_empty:
+ EXPECT_NE(1u, cables.size());
+ break;
+ }
+ }
+ }
}
TEST(cv_geom, multicell) {
@@ -301,6 +448,8 @@ TEST(cv_geom, multicell) {
cv_geometry geom2 = geom;
append(geom2, geom);
+ EXPECT_TRUE(verify_cv_children(geom));
+
ASSERT_EQ(2*n_cv, geom2.size());
for (unsigned i = 0; i<n_cv; ++i) {
EXPECT_EQ(geom.cv_parent[i], geom2.cv_parent[i]);
diff --git a/test/unit/test_cv_layout.cpp b/test/unit/test_cv_layout.cpp
index 4d5f9710c8b51f1d98a1e40f86be6a4042777783..62989fd88b003265130f4a402fc1f42baabf4cfb 100644
--- a/test/unit/test_cv_layout.cpp
+++ b/test/unit/test_cv_layout.cpp
@@ -146,13 +146,14 @@ TEST(cv_layout, zero_size_cv) {
// With one CV per branch, expect reference points for face conductance
// to be at (0, 0.5); (0, 1); (1, 0.5); (2, 1); (3, 0.5); (4, 0.5); (5, 0.5).
// The first CV should be all of branch 0; the second CV should be the
- // zero-size CV at the branch point (0, 1).
+ // zero-size CV covering the branch point (0, 1).
params.discretization = cv_policy_fixed_per_branch(1);
fvm_cv_discretization D = fvm_cv_discretize(cell, params);
unsigned cv_a = 0, cv_x = 1;
ASSERT_TRUE(util::equal(mcable_list{mcable{0, 0, 1}}, D.geometry.cables(cv_a)));
- ASSERT_TRUE(util::equal(mcable_list{mcable{0, 1, 1}}, D.geometry.cables(cv_x)));
+ ASSERT_TRUE(util::equal(mcable_list{mcable{0, 1, 1}, mcable{1, 0, 0}, mcable{2, 0, 0}},
+ D.geometry.cables(cv_x)));
// Find the two CV children of CV x.
unsigned cv_b = -1, cv_c = -1;
diff --git a/test/unit/test_fvm_lowered.cpp b/test/unit/test_fvm_lowered.cpp
index 8957d1ff9c99458282ec7d820ce4cdef49ad6377..884d2df185a5c3fc8e720fe2b26c758e74e8bf89 100644
--- a/test/unit/test_fvm_lowered.cpp
+++ b/test/unit/test_fvm_lowered.cpp
@@ -1020,14 +1020,15 @@ TEST(fvm_lowered, gj_coords_complex) {
fvcell.fvm_intdom(rec, gids, cell_to_intdom);
fvm_cv_discretization D = fvm_cv_discretize(cells, neuron_parameter_defaults, context);
+ using namespace cv_prefer;
int c0_gj_cv[2];
- for (int i = 0; i<2; ++i) c0_gj_cv[i] = D.geometry.location_cv(0, c0_gj[i]);
+ for (int i = 0; i<2; ++i) c0_gj_cv[i] = D.geometry.location_cv(0, c0_gj[i], cv_nonempty);
int c1_gj_cv[4];
- for (int i = 0; i<4; ++i) c1_gj_cv[i] = D.geometry.location_cv(1, c1_gj[i]);
+ for (int i = 0; i<4; ++i) c1_gj_cv[i] = D.geometry.location_cv(1, c1_gj[i], cv_nonempty);
int c2_gj_cv[3];
- for (int i = 0; i<3; ++i) c2_gj_cv[i] = D.geometry.location_cv(2, c2_gj[i]);
+ for (int i = 0; i<3; ++i) c2_gj_cv[i] = D.geometry.location_cv(2, c2_gj[i], cv_nonempty);
std::vector<fvm_gap_junction> GJ = fvcell.fvm_gap_junctions(cells, gids, rec, D);
EXPECT_EQ(10u, GJ.size());
diff --git a/test/unit/test_unique.cpp b/test/unit/test_unique.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..15e49ea1c690e9ecceefb823aeaf52d9e131e04b
--- /dev/null
+++ b/test/unit/test_unique.cpp
@@ -0,0 +1,39 @@
+#include "../gtest.h"
+
+#include <vector>
+#include <list>
+#include <utility>
+
+#include "util/unique.hpp"
+#include "./common.hpp"
+
+namespace {
+auto same_parity = [](auto a, auto b) { return a%2 == b%2; };
+
+template <typename C, typename Eq = std::equal_to<>>
+void run_unique_in_place_test(C data, const C& expected, Eq eq = Eq{}) {
+ arb::util::unique_in_place(data, eq);
+ EXPECT_TRUE(testing::seq_eq(data, expected));
+}
+
+template <typename C>
+void run_tests() {
+ run_unique_in_place_test(C{}, C{});
+ run_unique_in_place_test(C{1, 3, 2}, C{1, 3, 2});
+ run_unique_in_place_test(C{1, 1, 1, 3, 2}, C{1, 3, 2});
+ run_unique_in_place_test(C{1, 3, 3, 3, 2}, C{1, 3, 2});
+ run_unique_in_place_test(C{1, 3, 2, 2, 2}, C{1, 3, 2});
+ run_unique_in_place_test(C{1, 1, 3, 2, 2}, C{1, 3, 2});
+ run_unique_in_place_test(C{3, 1, 3, 1, 1, 3, 1, 1}, C{3, 1, 3, 1, 3, 1});
+ run_unique_in_place_test(C{1, 2, 4, 1, 3, 1, 2, 1}, C{1, 2, 1, 2, 1}, same_parity);
+}
+}
+
+TEST(unique_in_place, vector) {
+ run_tests<std::vector<int>>();
+}
+
+TEST(unique_in_place, list) {
+ run_tests<std::list<int>>();
+}
+