diff --git a/src/algorithms.hpp b/src/algorithms.hpp
index 7be1406ac7581df6924c146e6b86aaa8201755bb..94e198cb4da61d9e86e1b9424645b1ad739b113b 100644
--- a/src/algorithms.hpp
+++ b/src/algorithms.hpp
@@ -6,8 +6,9 @@
#include <type_traits>
#include <vector>
-#include "util.hpp"
-#include "util/debug.hpp"
+#include <util.hpp>
+#include <util/debug.hpp>
+#include <util/meta.hpp>
/*
* Some simple wrappers around stl algorithms to improve readability of code
@@ -23,18 +24,18 @@ namespace mc {
namespace algorithms {
template <typename C>
-typename C::value_type
+typename util::sequence_traits<C>::value_type
sum(C const& c)
{
- using value_type = typename C::value_type;
- return std::accumulate(c.begin(), c.end(), value_type{0});
+ using value_type = typename util::sequence_traits<C>::value_type;
+ return std::accumulate(std::begin(c), std::end(c), value_type{0});
}
template <typename C>
-typename C::value_type
+typename util::sequence_traits<C>::value_type
mean(C const& c)
{
- return sum(c)/c.size();
+ return sum(c)/util::size(c);
}
template <typename C>
diff --git a/src/compartment.hpp b/src/compartment.hpp
index 083155ffb7c5a3ea0d9d089740b6d7a3345ae234..db548fa39b7f4958b6575449a3734004d10ae687 100644
--- a/src/compartment.hpp
+++ b/src/compartment.hpp
@@ -4,8 +4,12 @@
#include <utility>
#include <common_types.hpp>
+#include <math.hpp>
#include <util/counter.hpp>
+#include <util/iterutil.hpp>
+#include <util/partition.hpp>
#include <util/span.hpp>
+#include "util/rangeutil.hpp"
#include <util/transform.hpp>
namespace nest {
@@ -79,6 +83,206 @@ compartment_range<size_type, real_type> make_compartment_range(
compartment_maker<size_type, real_type>(num_compartments, length, radius_L, radius_R));
}
+/// Divided compartments for use with (e.g.) fvm control volume setup
+
+struct semi_compartment {
+ using value_type = double;
+ using value_pair = std::pair<value_type, value_type>;
+
+ value_type length;
+ value_type area;
+ value_type volume;
+ value_pair radii;
+
+ semi_compartment& operator+=(const semi_compartment& x) {
+ length += x.length;
+ area += x.area;
+ volume += x.volume;
+ radii.second = x.radii.second;
+ return *this;
+ }
+
+ static semi_compartment frustrum(value_type l, value_type r1, value_type r2) {
+ using namespace math;
+ return semi_compartment{
+ l,
+ area_frustrum(l, r1, r2),
+ volume_frustrum(l, r1, r2),
+ {r1, r2}
+ };
+ }
+};
+
+struct div_compartment {
+ using value_type = typename semi_compartment::value_type;
+ using value_pair = typename semi_compartment::value_pair;
+ using size_type = cell_local_size_type;
+
+ size_type index = 0;
+ semi_compartment left;
+ semi_compartment right;
+
+ div_compartment() = default;
+ div_compartment(size_type i, semi_compartment l, semi_compartment r):
+ index(i), left(std::move(l)), right(std::move(r))
+ {}
+
+ value_type length() const { return left.length+right.length; }
+ value_type area() const { return left.area+right.area; }
+ value_type volume() const { return left.volume+right.volume; }
+ value_pair radii() const { return {left.radii.first, right.radii.second}; }
+};
+
+/// Divided compartments can be made from cables with sub-segments by
+/// sampling or integrating or by approximating cable as a single frustrum.
+
+class div_compartment_by_ends {
+public:
+ using value_type = div_compartment::value_type;
+ using size_type = div_compartment::size_type;
+
+ // `lengths` must be a sequence of length at least one.
+ // `radii` must be a sequence of length `size(lengths)+1`.
+ template <typename Radii, typename Lengths>
+ div_compartment_by_ends(size_type n, const Radii& radii, const Lengths& lengths):
+ oon_(1/value_type(n)),
+ length_(algorithms::sum(lengths)),
+ ra_(util::front(radii)),
+ rb_(util::back(radii))
+ {}
+
+ div_compartment operator()(size_type i) const {
+ value_type r1 = math::lerp(ra_, rb_, i*oon_);
+ value_type rc = math::lerp(ra_, rb_, (i+0.5)*oon_);
+ value_type r2 = math::lerp(ra_, rb_, (i+1)*oon_);
+ value_type semilength = length_*oon_*0.5;
+
+ return div_compartment(
+ i,
+ semi_compartment::frustrum(semilength, r1, rc),
+ semi_compartment::frustrum(semilength, rc, r2)
+ );
+ }
+
+private:
+ value_type oon_;
+ value_type length_;
+ value_type ra_;
+ value_type rb_;
+};
+
+class div_compartment_sampler {
+public:
+ using value_type = div_compartment::value_type;
+ using size_type = div_compartment::size_type;
+
+ // `lengths` must be a sequence of length at least one.
+ // `radii` must be a sequence of length `size(lengths)+1`.
+ template <typename Radii, typename Lengths>
+ div_compartment_sampler(size_type n, const Radii& radii, const Lengths& lengths) {
+ // set up offset-to-subsegment lookup and interpolation
+ using namespace util;
+
+ segs_ = make_partition(offsets_, lengths);
+ nseg_ = size(segs_);
+ scale_ = segs_.bounds().second/n;
+ assign(radii_, radii);
+ EXPECTS(size(radii_)==size(offsets_));
+ }
+
+ div_compartment operator()(size_type i) const {
+ using namespace math;
+
+ auto r1 = radius_at(locate(scale_*i));
+ auto rc = radius_at(locate(scale_*(i+0.5)));
+ auto r2 = radius_at(locate(scale_*(i+1)));
+
+ value_type semilength = 0.5*scale_;
+ return div_compartment(
+ i,
+ semi_compartment::frustrum(semilength, r1, rc),
+ semi_compartment::frustrum(semilength, rc, r2)
+ );
+ }
+
+protected:
+ struct sub_segment_index {
+ size_type i; // index
+ value_type p; // proportion [0,1] along sub-segment
+
+ sub_segment_index(size_type i_, value_type p_): i(i_), p(p_) {}
+ bool operator<(sub_segment_index x) const {
+ return i<x.i || (i==x.i && p<x.p);
+ }
+ };
+
+ sub_segment_index locate(value_type x) const {
+ EXPECTS(x>=0);
+
+ auto i = segs_.index(x);
+ if (i==segs_.npos) {
+ i = nseg_-1;
+ }
+
+ auto seg = segs_[i];
+ if (x>=seg.second) {
+ return sub_segment_index(i, 1);
+ }
+ return sub_segment_index(i, (x-seg.first)/(seg.second-seg.first));
+ }
+
+ value_type radius_at(sub_segment_index s) const {
+ return math::lerp(radii_[s.i], radii_[s.i+1], s.p);
+ }
+
+ size_type nseg_ = 0;
+ std::vector<value_type> offsets_;
+ std::vector<value_type> radii_;
+ value_type scale_ = 0;
+ decltype(util::partition_view(offsets_)) segs_;
+};
+
+/// Overrides operator() with a more accurate method
+class div_compartment_integrator: public div_compartment_sampler {
+public:
+ template <typename Radii, typename Lengths>
+ div_compartment_integrator(size_type n, const Radii& radii, const Lengths& lengths):
+ div_compartment_sampler(n, radii, lengths) {}
+
+ div_compartment operator()(size_type i) const {
+ using namespace math;
+
+ auto sleft = locate(scale_*i);
+ auto scentre = locate(scale_*(i+0.5));
+ auto sright = locate(scale_*(i+1));
+
+ return div_compartment(i, integrate(sleft, scentre), integrate(scentre, sright));
+ }
+
+protected:
+ semi_compartment sub_segment_frustrum(sub_segment_index a, sub_segment_index b) const {
+ EXPECTS(a.i==b.i && a.p<=b.p);
+
+ auto seg = segs_[a.i];
+ auto l = (b.p-a.p)*(seg.second-seg.first);
+ return semi_compartment::frustrum(l, radius_at(a), radius_at(b));
+ }
+
+ semi_compartment integrate(sub_segment_index a, sub_segment_index b) const {
+ sub_segment_index x = std::min(b, sub_segment_index(a.i, 1));
+
+ auto s = sub_segment_frustrum(a, x);
+ while (a.i<b.i) {
+ ++a.i;
+ a.p = 0;
+ x = std::min(b, sub_segment_index(a.i, 1));
+ s += sub_segment_frustrum(a, x);
+ }
+
+ return s;
+ }
+};
+
} // namespace mc
} // namespace nest
diff --git a/src/segment.hpp b/src/segment.hpp
index eedb5095f89166deb5f71735da28d0a74a3ab0cb..dd6e8fee88dc1cf9a2ab12896a1760d6a933a95b 100644
--- a/src/segment.hpp
+++ b/src/segment.hpp
@@ -3,13 +3,13 @@
#include <cmath>
#include <vector>
-#include "algorithms.hpp"
-#include "common_types.hpp"
-#include "compartment.hpp"
-#include "math.hpp"
-#include "parameter_list.hpp"
-#include "point.hpp"
-#include "util.hpp"
+#include <algorithms.hpp>
+#include <common_types.hpp>
+#include <compartment.hpp>
+#include <math.hpp>
+#include <parameter_list.hpp>
+#include <point.hpp>
+#include <util.hpp>
namespace nest {
namespace mc {
@@ -379,6 +379,11 @@ public:
return lengths_;
}
+ std::vector<value_type> const& radii() const
+ {
+ return radii_;
+ }
+
cable_segment* as_cable() override
{
return this;
@@ -465,6 +470,18 @@ segment_ptr make_segment(Args&&... args)
return segment_ptr(new SegmentType(std::forward<Args>(args)...));
}
+/// Divided compartment adaptors for cable segments
+
+template <typename DivCompClass>
+DivCompClass div_compartments(const cable_segment* cable, unsigned ncomp) {
+ return DivCompClass(ncomp, cable->radii(), cable->lengths());
+}
+
+template <typename DivCompClass>
+DivCompClass div_compartments(const cable_segment* cable) {
+ return DivCompClass(cable->num_compartments(), cable->radii(), cable->lengths());
+}
+
} // namespace mc
} // namespace nest
diff --git a/tests/unit/test_compartments.cpp b/tests/unit/test_compartments.cpp
index 8bc034ade553675ac8f85ce19bc8bfccf5d14d5d..8701c2716f6fe40b44e0c2832a198bf1c8ca4f3b 100644
--- a/tests/unit/test_compartments.cpp
+++ b/tests/unit/test_compartments.cpp
@@ -1,12 +1,19 @@
#include <cmath>
+#include <string>
#include "gtest.h"
+#include <algorithms.hpp>
#include <compartment.hpp>
+#include <math.hpp>
#include <util.hpp>
+#include <util/span.hpp>
+#include <util/transform.hpp>
-using nest::mc::util::left;
-using nest::mc::util::right;
+using namespace nest::mc;
+using namespace nest::mc::algorithms;
+using namespace nest::mc::math;
+using namespace nest::mc::util;
// not much to test here: just test that values passed into the constructor
// are correctly stored in members
@@ -59,3 +66,250 @@ TEST(compartments, make_compartment_range)
auto rng_empty = make_compartment_range(0, 1.0, 1.0, 0.);
EXPECT_EQ(rng_empty.begin(), rng_empty.end());
}
+
+// Divided compartments
+// (FVM-friendly compartment data)
+
+template <std::size_t N>
+struct pw_cable_data {
+ double radii[N+1];
+ double lengths[N];
+
+ static constexpr std::size_t nseg() { return N; }
+ double r1() const { return radii[0]; }
+ double r2() const { return radii[N]; }
+ double length() const { return sum(lengths); }
+
+ double area() const {
+ return sum(transform_view(make_span(0, N),
+ [&](unsigned i) { return area_frustrum(lengths[i], radii[i], radii[i+1]); }));
+ }
+
+ double volume() const {
+ return sum(transform_view(make_span(0, N),
+ [&](unsigned i) { return volume_frustrum(lengths[i], radii[i], radii[i+1]); }));
+ }
+};
+
+pw_cable_data<1> cable_one = {
+ {2.0, 5.0},
+ {10.0}
+};
+
+pw_cable_data<4> cable_linear = {
+ {2.0, 3.5, 6.0, 6.5, 6.75},
+ {3.0, 5.0, 1.0, 0.5}
+};
+
+pw_cable_data<4> cable_jumble = {
+ {2.0, 6.0, 3.5, 6.75, 6.5},
+ {3.0, 5.0, 1.0, 0.5}
+};
+
+void expect_equal_divs(const div_compartment& da, const div_compartment& db) {
+ EXPECT_EQ(da.index, db.index);
+
+ double eps = std::numeric_limits<double>::epsilon();
+ double e1 = std::min(da.length(), db.length())*8*eps;
+ double e2 = std::min(da.area(), db.area())*8*eps;
+ double e3 = std::min(da.volume(), db.volume())*8*eps;
+
+ EXPECT_NEAR(da.left.length, db.left.length, e1);
+ EXPECT_NEAR(da.left.area, db.left.area, e2);
+ EXPECT_NEAR(da.left.volume, db.left.volume, e3);
+ EXPECT_NEAR(da.left.radii.first, db.left.radii.first, e1);
+ EXPECT_NEAR(da.left.radii.second, db.left.radii.second, e1);
+
+ EXPECT_NEAR(da.right.length, db.right.length, e1);
+ EXPECT_NEAR(da.right.area, db.right.area, e2);
+ EXPECT_NEAR(da.right.volume, db.right.volume, e3);
+ EXPECT_NEAR(da.right.radii.first, db.right.radii.first, e1);
+ EXPECT_NEAR(da.right.radii.second, db.right.radii.second, e1);
+}
+
+TEST(compartments, div_ends) {
+ using namespace math;
+
+ {
+ div_compartment_by_ends divcomps{1, cable_one.radii, cable_one.lengths};
+
+ auto d = divcomps(0);
+ auto r1 = cable_one.radii[0];
+ auto r2 = cable_one.radii[1];
+ auto l = cable_one.lengths[0];
+
+ EXPECT_DOUBLE_EQ(r1, d.radii().first);
+ EXPECT_DOUBLE_EQ(r2, d.radii().second);
+ EXPECT_DOUBLE_EQ(l, d.length());
+ EXPECT_DOUBLE_EQ(area_frustrum(l, r2, r1), d.area());
+ EXPECT_DOUBLE_EQ(volume_frustrum(l, r2, r1), d.volume());
+
+ auto sl = l/2.0;
+ auto rc = mean(r1, r2);
+
+ div_compartment expected{
+ 0,
+ semi_compartment{sl, area_frustrum(sl, r1, rc), volume_frustrum(sl, r1, rc), {r1, rc}},
+ semi_compartment{sl, area_frustrum(sl, rc, r2), volume_frustrum(sl, rc, r2), {rc, r2}}
+ };
+
+ SCOPED_TRACE("cable_one");
+ expect_equal_divs(expected, d);
+ }
+
+ {
+ // for a linear cable, expect this compartment maker to
+ // create consistent compartments
+
+ constexpr unsigned ncomp = 7;
+ div_compartment_by_ends divlin{ncomp, cable_linear.radii, cable_linear.lengths};
+
+ auto r1 = cable_linear.r1();
+ auto r2 = cable_linear.r2();
+ auto l = cable_linear.length();
+
+ pw_cable_data<1> one = { {r1, r2}, {l} };
+ div_compartment_by_ends divone{ncomp, one.radii, one.lengths};
+
+ for (unsigned i=0; i<ncomp; ++i) {
+ SCOPED_TRACE("cable_linear compartment "+std::to_string(i));
+ auto da = divlin(i);
+ auto db = divone(i);
+
+ EXPECT_DOUBLE_EQ(l/ncomp, da.length());
+ expect_equal_divs(da, db);
+ }
+ }
+}
+
+TEST(compartments, div_sample) {
+ using namespace math;
+
+ // expect by_ends and sampler to give same results on linear cable
+ {
+ constexpr unsigned ncomp = 7;
+ div_compartment_sampler divsampler{ncomp, cable_linear.radii, cable_linear.lengths};
+ div_compartment_by_ends divends{ncomp, cable_linear.radii, cable_linear.lengths};
+
+ auto l = cable_linear.length();
+ for (unsigned i=0; i<ncomp; ++i) {
+ SCOPED_TRACE("cable_linear compartment "+std::to_string(i));
+ auto da = divsampler(i);
+ auto db = divends(i);
+
+ EXPECT_DOUBLE_EQ(l/ncomp, da.length());
+ expect_equal_divs(da, db);
+ }
+ }
+
+ // expect (up to rounding) correct total area and volume if compartments
+ // align with sub-segments; when they don't align, expect error to decrease
+ // with ncomp
+ {
+ double area_expected = cable_jumble.area();
+ double volume_expected = cable_jumble.volume();
+ double eps = std::numeric_limits<double>::epsilon();
+
+ double common_dx = 0.5; // depends on cable_jumble.lengths;
+ // check our common_dx actually is ok
+ for (double l: cable_jumble.lengths) {
+ ASSERT_DOUBLE_EQ(l/common_dx, std::round(l/common_dx));
+ }
+
+ double length = cable_jumble.length();
+ unsigned nbase = std::round(length/common_dx);
+
+ for (unsigned m: {1u, 3u, 7u}) {
+ unsigned ncomp = m*nbase;
+ div_compartment_sampler divs{ncomp, cable_jumble.radii, cable_jumble.lengths};
+
+ double area = sum(transform_view(make_span(0, ncomp),
+ [&](unsigned i) { return divs(i).area(); }));
+
+ double volume = sum(transform_view(make_span(0, ncomp),
+ [&](unsigned i) { return divs(i).volume(); }));
+
+ double e2 = std::min(area, area_expected)*ncomp*eps;
+ double e3 = std::min(volume, volume_expected)*ncomp*eps;
+
+ SCOPED_TRACE("cable_jumble ncomp "+std::to_string(ncomp));
+
+ EXPECT_NEAR(area_expected, area, e2);
+ EXPECT_NEAR(volume_expected, volume, e3);
+ }
+
+ double coarse_area = area_frustrum(length, cable_jumble.r1(), cable_jumble.r2());
+ double coarse_volume = volume_frustrum(length, cable_jumble.r1(), cable_jumble.r2());
+ double area_error = std::abs(area_expected-coarse_area);
+ double volume_error = std::abs(volume_expected-coarse_volume);
+
+ for (unsigned m: {1u, 10u, 100u}) {
+ unsigned ncomp = m*nbase+1u;
+ div_compartment_sampler divs{ncomp, cable_jumble.radii, cable_jumble.lengths};
+
+ double area = sum(transform_view(make_span(0, ncomp),
+ [&](unsigned i) { return divs(i).area(); }));
+
+ double volume = sum(transform_view(make_span(0, ncomp),
+ [&](unsigned i) { return divs(i).volume(); }));
+
+ SCOPED_TRACE("cable_jumble ncomp "+std::to_string(ncomp));
+
+ double err = std::abs(area_expected-area);
+ EXPECT_LT(err, area_error);
+ area_error = err;
+
+ err = std::abs(volume_expected-volume);
+ EXPECT_LT(err, volume_error);
+ volume_error = err;
+ }
+ }
+}
+
+TEST(compartments, div_integrator) {
+ using namespace math;
+
+ // expect integrator and sampler to give same results on linear cable
+ {
+ constexpr unsigned ncomp = 7;
+ div_compartment_sampler divintegrator{ncomp, cable_linear.radii, cable_linear.lengths};
+ div_compartment_by_ends divends{ncomp, cable_linear.radii, cable_linear.lengths};
+
+ auto l = cable_linear.length();
+ for (unsigned i=0; i<ncomp; ++i) {
+ SCOPED_TRACE("cable_linear compartment "+std::to_string(i));
+ auto da = divintegrator(i);
+ auto db = divends(i);
+
+ EXPECT_DOUBLE_EQ(l/ncomp, da.length());
+ expect_equal_divs(da, db);
+ }
+ }
+
+ // expect integrator to give same (up to rounding) total areas and volumes
+ // as the cable
+ {
+ double area_expected = cable_jumble.area();
+ double volume_expected = cable_jumble.volume();
+ double eps = std::numeric_limits<double>::epsilon();
+
+ for (unsigned ncomp: make_span(1u, 23u)) {
+ div_compartment_integrator divs{ncomp, cable_jumble.radii, cable_jumble.lengths};
+
+ double area = sum(transform_view(make_span(0, ncomp),
+ [&](unsigned i) { return divs(i).area(); }));
+
+ double volume = sum(transform_view(make_span(0, ncomp),
+ [&](unsigned i) { return divs(i).volume(); }));
+
+ double e2 = std::min(area, area_expected)*ncomp*eps;
+ double e3 = std::min(volume, volume_expected)*ncomp*eps;
+
+ SCOPED_TRACE("cable_jumble ncomp "+std::to_string(ncomp));
+
+ EXPECT_NEAR(area_expected, area, e2);
+ EXPECT_NEAR(volume_expected, volume, e3);
+ }
+ }
+}
+