New compartment info structure for FVM.

* Make `algorithm::sum`, `algorithm::mean` more generic,
  allowing use with array types.
* Add `div_compartment` compartment representation, that
  holds geometric information for each half of a compartment
  that will then be used in calculating control volumes.
* Add three compartmentalisation schemes/policies that
  discretize a segment into `div_compartment` objects:
    * `div_compartment_by_ends` divides based only on the
      segment end points and radii.
    * `div_compartment_sampler` forms frusta by sampling
      the segment radius at each compartment boundary
    * `div_compartment_integrator` computes the compartment
      areas and volumes exactly by summing all frustra
      in the intersection of the segment and the compartmnet
      span.

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>

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
 

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
 

tests/unit/test_compartments.cpp

 #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);
+        }
+    }
+}
+