start working on compartment model

src/compartment.hpp

+#pragma once
+
+#include <cmath>
+
+#include <vector>
+
+#include "point.hpp"
+
+/*
+    We start with a high-level description of the cell
+    - list of branches of the cell
+        - soma, dendrites, axons
+        - spatial locations if provided
+            - bare minimum spatial information required is length and radius
+              at each end for each of the branches, and a soma radius
+        - model properties of each branch
+            - mechanisms
+            - clamps
+            - synapses
+        - list of compartments if they have been provided
+
+    This description is not used for solving the system
+    From the description we can then build a cell solver
+    - e.g. the FVM formulation
+    - e.g. Green's functions
+
+*/
+
+namespace nestmc {
+
+template <typename T>
+T constexpr pi() {
+    return 3.141592653589793238462643383279502;
+}
+
+template <typename T>
+class compartment {
+    using value_type = T;
+    using point_type = point<value_type>;
+
+    constexpr compartment()
+    :   p1_{point_type()},
+        p2_{point_type()},
+        radius1_{std::numeric_limits<value_type>::quiet_NaN()},
+        radius2_{std::numeric_limits<value_type>::quiet_NaN()}
+    {}
+
+    constexpr compartment(
+        point_type const& p1,
+        point_type const& p2,
+        value_type r1,
+        value_type r2
+    )
+    :   p1_{p1},
+        p2_{p2},
+        radius1_{r1},
+        radius2_{r2}
+    {}
+
+    value_type length() const {
+        return norm(p1_-p2_);
+    }
+
+    value_type area() const {
+        return volume_frustrum(length(), radius1_, radius2_);
+    }
+
+    value_type volume() const {
+        return volume_frustrum(length(), radius1_, radius2_);
+    }
+
+    constexpr point_type midpoint() const {
+        return 0.5*(p1_+p2_);
+    }
+
+    constexpr point_type midradius() const {
+        return 0.5*(radius1_+radius2_);
+    }
+
+    private :
+
+    value_type area_frustrum(
+        value_type L,
+        value_type r1,
+        value_type r2
+    ) const
+    {
+        auto dr = r1 - r2;
+        return pi<double>() * (r1+r2) * std::sqrt(L*L + dr*dr);
+    }
+
+    value_type volume_frustrum(
+        value_type L,
+        value_type r1,
+        value_type r2
+    ) const
+    {
+        auto meanr = (r1+r2) / 2.;
+        return pi<double>() * meanr * meanr * L;
+    }
+
+    point_type p1_;
+    point_type p2_;
+    value_type radius1_;
+    value_type radius2_;
+};
+
+class segment {
+
+    private :
+
+    double length_;
+    double radius_start_;
+    double radius_end_;
+    std::vector<compartment> compartments_;
+    segmentKind kind_;
+};
+
+class abstract_cell {
+    abstract_cell() = default;
+
+    private :
+};
+
+} // namespace nestmc 
+
+

src/math.hpp

+#pragma once
+
+#include <cmath>
+
+namespace nestmc {
+
+    template <typename T>
+    T constexpr pi() {
+        return 3.141592653589793238462643383279502;
+    }
+
+    template <typename T>
+    T area_frustrum(T L, T r1, T r2)
+    {
+        auto dr = r1 - r2;
+        return pi<double>() * (r1+r2) * sqrt(L*L + dr*dr);
+    }
+
+    template <typename T>
+    T volume_frustrum(T L, T r1, T r2)
+    {
+        auto meanr = (r1+r2) / 2.;
+        return pi<double>() * meanr * meanr * L;
+    }
+
+} // namespace nestmc

src/point.hpp

+#pragma once
+
+#include <limits>
+#include <ostream>
+
+template <typename T>
+struct point {
+    using value_type = T;
+    value_type x;
+    value_type y;
+    value_type z;
+
+    constexpr point(T a, T b, T c)
+    : x(a), y(b), z(c)
+    {}
+
+    // initialize to NaN by default
+    constexpr point()
+    : x(std::numeric_limits<T>::quiet_NaN()),
+      y(std::numeric_limits<T>::quiet_NaN()),
+      z(std::numeric_limits<T>::quiet_NaN())
+    {}
+};
+
+template <typename T>
+constexpr point<T>  operator+ (
+    point<T> const& lhs,
+    point<T> const& rhs)
+{
+    return point<T>(lhs.x + rhs.x, lhs.y + rhs.y, lhs.z + rhs.z);
+}
+
+template <typename T>
+constexpr point<T>  operator- (
+    point<T> const& lhs,
+    point<T> const& rhs)
+{
+    return point<T>(lhs.x - rhs.x, lhs.y - rhs.y, lhs.z - rhs.z);
+}
+
+template <typename T>
+constexpr point<T> operator* (
+    T alpha,
+    point<T> const& p)
+{
+    return point<T>(alpha*p.x, alpha*p.y, alpha*p.z);
+}
+
+template <typename T>
+T norm(point<T> const& p)
+{
+    return sqrt(p.x*p.x + p.y*p.y + p.z*p.z);
+}
+
+template <typename T>
+constexpr T dot(
+    point<T> const& lhs,
+    point<T> const& rhs)
+{
+    return lhs.x*rhs.x + lhs.y*rhs.y + lhs.z*rhs.z;
+}
+
+template <typename T>
+std::ostream& operator << (std::ostream& o, point<T> const& p) {
+    return o << "[" << p.x << ", " << p.y << ", " << p.z << "]";
+}
+
+

src/segment.hpp

+#pragma once
+
+#include <cmath>
+
+#include <vector>
+
+#include "math.hpp"
+#include "point.hpp"
+#include "util.hpp"
+
+/*
+    We start with a high-level description of the cell
+    - list of branches of the cell
+        - soma, dendrites, axons
+        - spatial locations if provided
+            - bare minimum spatial information required is length and radius
+              at each end for each of the branches, and a soma radius
+        - model properties of each branch
+            - mechanisms
+            - clamps
+            - synapses
+        - list of compartments if they have been provided
+
+    This description is not used for solving the system
+    From the description we can then build a cell solver
+    - e.g. the FVM formulation
+    - e.g. Green's functions
+
+*/
+
+namespace nestmc {
+
+enum class segmentKind {
+    soma,
+    dendrite,
+    axon
+};
+
+class segment {
+    public:
+
+    using value_type = double;
+    using point_type = point<value_type>;
+
+    segmentKind kind() const {
+        return kind_;
+    }
+
+    virtual value_type volume() const = 0;
+    virtual value_type area()  const = 0;
+
+    virtual ~segment() = default;
+
+    protected:
+
+    segmentKind kind_;
+};
+
+class spherical_segment : public segment
+{
+    public :
+
+    using base = segment;
+    using base::kind_;
+    using base::value_type;
+    using base::point_type;
+
+    spherical_segment() = delete;
+
+    spherical_segment(value_type r)
+    :   radius_{r}
+    {
+        kind_ = segmentKind::soma;
+    }
+
+    spherical_segment(value_type r, point_type const &c)
+    :   spherical_segment(r)
+    {
+        center_ = c;
+        kind_ = segmentKind::soma;
+    }
+
+    value_type volume() const override
+    {
+        return 4./3. * pi<value_type>() * radius_ * radius_ * radius_;
+    }
+
+    value_type area() const override
+    {
+        return 4. * pi<value_type>() * radius_ * radius_;
+    }
+
+    virtual ~spherical_segment() = default;
+
+    private :
+
+    // store the center and radius of the soma
+    // note that the center may be undefined
+    value_type radius_;
+    point_type center_;
+};
+
+class frustrum_segment : public segment
+{
+    public :
+
+    using segment::kind_;
+    using base = segment;
+    using base::value_type;
+    using base::point_type;
+
+    frustrum_segment() = delete;
+
+    frustrum_segment(
+        segmentKind k,
+        value_type r1,
+        value_type r2,
+        value_type len
+    ) {
+        r1_ = r1;
+        r2_ = r2;
+        length_ = len;
+        kind_ = k;
+        assert(k==segmentKind::dendrite || k ==segmentKind::axon);
+    }
+
+    frustrum_segment(
+        segmentKind k,
+        value_type r1,
+        value_type r2,
+        point_type const& p1,
+        point_type const& p2
+    )
+    : frustrum_segment(k, r1, r2, norm(p1-p2))
+    {
+        p1_ = p1;
+        p2_ = p2;
+    }
+
+    value_type volume() const override
+    {
+        return volume_frustrum(length_, r1_, r2_);
+    }
+    value_type area()   const override
+    {
+        return area_frustrum(length_, r1_, r2_);
+    }
+
+    virtual ~frustrum_segment() = default;
+
+    private :
+
+    value_type length_;
+    value_type r1_;
+    value_type r2_;
+    point_type p1_;
+    point_type p2_;
+};
+
+using segment_ptr = std::unique_ptr<segment>;
+
+template <typename T, typename... Args>
+segment_ptr make_segment(Args&&... args) {
+    return segment_ptr(new T(std::forward<Args>(args)...));
+}
+
+} // namespace nestmc
+

src/util.hpp

@@ -11,6 +11,7 @@ using memory::util::blue;
 using memory::util::cyan;
 */
 
+#include <memory>
 #include <ostream>
 #include <vector>
 
@@ -26,3 +27,12 @@ operator << (std::ostream &o, std::vector<T>const& v)
     return o;
 }
 
+namespace util {
+    // just because we aren't using C++14, doesn't mean we shouldn't go
+    // without make_unique
+    template <typename T, typename... Args>
+    std::unique_ptr<T> make_unique(Args&&... args) {
+        return std::unique_ptr<T>(new T(std::forward<Args>(args) ...));
+    }
+}
+