basic abstract cell

src/CMakeLists.txt

@@ -3,6 +3,7 @@ set(HEADERS
     )
 set(BASE_SOURCES
     swcio.cpp
+    cell.cpp
 )
 
 add_library(cellalgo ${BASE_SOURCES} ${HEADERS})

src/cell.cpp

+#include "cell.hpp"
+
+namespace nestmc {
+
+int cell::num_segments() const
+{
+    return segments_.size();
+}
+
+void cell::add_soma(value_type radius, point_type center)
+{
+    if(has_soma()) {
+        throw std::domain_error(
+            "attempt to add a soma to a cell that already has one"
+        );
+    }
+
+    // soma has intself as its own parent
+    soma_ = num_segments();
+    parents_.push_back(num_segments());
+
+    // add segment for the soma
+    if(center.is_set()) {
+        segments_.push_back(
+            make_segment<soma_segment>(radius, center)
+        );
+    }
+    else {
+        segments_.push_back(
+            make_segment<soma_segment>(radius)
+        );
+    }
+}
+
+// add a cable that is provided by the user as a segment_ptr
+void cell::add_cable(cell::index_type parent, segment_ptr&& cable)
+{
+    // check for a valid parent id
+    if(cable->is_soma()) {
+        throw std::domain_error(
+            "attempt to add a soma as a segment"
+        );
+    }
+
+    // check for a valid parent id
+    if(parent>num_segments()) {
+        throw std::out_of_range(
+            "parent index of cell segment is out of range"
+        );
+    }
+    segments_.push_back(std::move(cable));
+    parents_.push_back(parent);
+}
+
+
+bool cell::has_soma() const
+{
+    return soma_ > -1;
+}
+
+soma_segment* cell::soma()
+{
+    if(has_soma()) {
+        return segments_[soma_].get()->as_soma();
+    }
+    return nullptr;
+}
+
+cell::value_type cell::volume() const
+{
+    return
+       std::accumulate(
+            segments_.begin(), segments_.end(),
+            0.,
+            [](double value, segment_ptr const& seg) {
+                return seg->volume() + value;
+            }
+       );
+}
+
+cell::value_type cell::area() const
+{
+    return
+       std::accumulate(
+            segments_.begin(), segments_.end(),
+            0.,
+            [](double value, segment_ptr const& seg) {
+                return seg->area() + value;
+            }
+       );
+}
+
+std::vector<segment_ptr> const& cell::segments() const
+{
+    return segments_;
+}
+
+void cell::construct()
+{
+    if(num_segments()) {
+        tree_ = cell_tree(parents_);
+    }
+}
+
+cell_tree const& cell::graph() const
+{
+    return tree_;
+}
+
+std::vector<int> const& cell::segment_parents() const
+{
+    return parents_;
+}
+
+} // namespace nestmc

src/cell.hpp

@@ -8,85 +8,56 @@
 
 namespace nestmc {
 
-    // we probably need two cell types
-    //  1. the abstract cell type (which would be this one)
-    //  2.
+    // high-level abstract representation of a cell and its segments
     class cell {
         public:
-        using index_type = int16_t;
+
+        // types
+        using index_type = int;
         using value_type = double;
         using point_type = point<value_type>;
 
-        int num_segments() const
-        {
-            return segments_.size();
-        }
+        /// add a soma to the cell
+        /// radius must be specified
+        void add_soma(value_type radius, point_type center=point_type());
 
-        void add_soma(value_type radius, point_type center=point_type())
-        {
-            if(has_soma()) {
-                throw std::domain_error(
-                    "attempt to add a soma to a cell that already has one"
-                );
-            }
+        /// add a cable
+        /// parent is the index of the parent segment for the cable section
+        /// cable is the segment that will be moved into the cell
+        void add_cable(index_type parent, segment_ptr&& cable);
 
-            // soma has intself as its own parent
-            soma_ = num_segments();
-            parents_.push_back(num_segments());
+        /// add a cable by constructing it in place
+        /// parent is the index of the parent segment for the cable section
+        /// args are the arguments to be used to consruct the new cable
+        template <typename... Args>
+        void add_cable(index_type parent, Args ...args);
 
-            // add segment for the soma
-            if(center.is_set()) {
-                segments_.push_back(
-                    make_segment<soma_segment>(radius, center)
-                );
-            }
-            else {
-                segments_.push_back(
-                    make_segment<soma_segment>(radius)
-                );
-            }
-        }
+        /// the number of segments in the cell
+        int num_segments() const;
 
-        void add_cable(segment_ptr&& cable, index_type parent)
-        {
-            // check for a valid parent id
-            if(cable->is_soma()) {
-                throw std::domain_error(
-                    "attempt to add a soma as a segment"
-                );
-            }
+        bool has_soma() const;
 
-            // check for a valid parent id
-            if(parent>num_segments()) {
-                throw std::out_of_range(
-                    "parent index of cell segment is out of range"
-                );
-            }
-            segments_.push_back(std::move(cable));
-            parents_.push_back(parent);
-        }
+        /// access pointer to the soma
+        soma_segment* soma();
 
-        template <typename... Args>
-        void add_cable(index_type parent, Args ...args)
-        {
-            // check for a valid parent id
-            if(parent>num_segments()) {
-                throw std::out_of_range(
-                    "parent index of cell segment is out of range"
-                );
-            }
-            segments_.push_back(make_segment<cable_segment>(std::forward<Args>(args)...));
-            parents_.push_back(parent);
-        }
+        /// the volume of the cell
+        value_type volume() const;
 
-        bool has_soma() const { return soma_ > -1; }
+        /// the surface area of the cell
+        value_type area() const;
 
-        soma_segment* soma() {
-            if(has_soma()) {
-                return segments_[soma_].get()->as_soma();
-            }
-            return nullptr;
-        }
+        std::vector<segment_ptr> const& segments() const;
+
+        /// generate the internal representation of the connectivity
+        /// graph for the cell segments
+        void construct();
+
+        /// the connectivity graph for the cell segments
+        cell_tree const& graph() const;
+
+        /// return reference to array that enumerates the index of the parent of
+        /// each segment
+        std::vector<int> const& segment_parents() const;
 
         private:
 
@@ -103,12 +74,27 @@ namespace nestmc {
         int soma_ = -1;
 
         //
-        // fixed cell description, which is computed from the
-        // rough layout description above
+        // fixed cell description, which is computed from the layout description
+        // above
         //
 
-        // cell_tree that describes the connection layout
-        //cell_tree tree_;
+        cell_tree tree_;
     };
 
+    // create a cable by forwarding cable construction parameters provided by the user
+    template <typename... Args>
+    void cell::add_cable(cell::index_type parent, Args ...args)
+    {
+        // check for a valid parent id
+        if(parent>num_segments()) {
+            throw std::out_of_range(
+                "parent index of cell segment is out of range"
+            );
+        }
+        segments_.push_back(make_segment<cable_segment>(std::forward<Args>(args)...));
+        parents_.push_back(parent);
+    }
+
+
 } // namespace nestmc
+

src/cell_tree.hpp

@@ -33,6 +33,9 @@ class cell_tree {
     using index_type = memory::HostVector<int_type>;
     using index_view = index_type::view_type;
 
+    /// default empty constructor
+    cell_tree() = default;
+
     /// construct from a parent index
     cell_tree(std::vector<int> const& parent_index)
     {
@@ -68,11 +71,27 @@ class cell_tree {
       soma_(other.soma())
     { }
 
+    // assignment from rvalue
+    cell_tree& operator=(cell_tree&& other)
+    {
+        std::swap(other.tree_, tree_);
+        std::swap(other.soma_, soma_);
+        return *this;
+    }
+
+    // assignment
+    cell_tree& operator=(cell_tree const& other)
+    {
+        tree_ = other.tree_;
+        soma_ = other.soma_;
+        return *this;
+    }
+
     // move constructor
     cell_tree(cell_tree&& other)
-    : tree_(std::move(other.tree_)),
-      soma_(other.soma())
-    { }
+    {
+        *this = std::move(other);
+    }
 
     int_type soma() const {
         return soma_;
@@ -262,8 +281,14 @@ class cell_tree {
         }
     }
 
-    // storage for the tree structure of cell segments
+    //////////////////////////////////////////////////
+    // state
+    //////////////////////////////////////////////////
+
+    /// storage for the tree structure of cell segments
     tree tree_;
 
+    /// index of the soma
     int_type soma_ = 0;
 };
+

src/math.hpp

@@ -5,22 +5,26 @@
 namespace nestmc {
 namespace math {
     template <typename T>
-    T constexpr pi() {
+    T constexpr pi()
+    {
         return T(3.1415926535897932384626433832795);
     }
 
     template <typename T>
-    T constexpr mean(T a, T b) {
+    T constexpr mean(T a, T b)
+    {
         return (a+b) / T(2);
     }
 
     template <typename T>
-    T constexpr square(T a) {
+    T constexpr square(T a)
+    {
         return a*a;
     }
 
     template <typename T>
-    T constexpr cube(T a) {
+    T constexpr cube(T a)
+    {
         return a*a*a;
     }
 

src/segment.hpp

@@ -8,28 +8,15 @@
 #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 {
 
+template <typename T,
+          typename valid = typename std::is_floating_point<T>::type>
+struct segment_properties {
+    T rL = 180.0;   // resistivity [Ohm.cm]
+    T cm = 0.01;    // capacitance [F/m^2] : 10 nF/mm^2 = 0.01 F/m^2
+};
+
 enum class segmentKind {
     soma,
     dendrite,
@@ -81,6 +68,8 @@ class segment {
         return nullptr;
     }
 
+    segment_properties<value_type> properties;
+
     protected:
 
     segmentKind kind_;

src/tree.hpp

@@ -19,20 +19,34 @@ class tree {
     using index_type = memory::HostVector<int_type>;
     using index_view = index_type::view_type;
 
+    tree() = default;
+
+    tree& operator=(tree&& other) {
+        std::swap(data_, other.data_);
+        std::swap(child_index_, other.child_index_);
+        std::swap(children_, other.children_);
+        std::swap(parents_, other.parents_);
+        return *this;
+    }
+
+    tree& operator=(tree const& other) {
+        data_ = other.data_;
+        set_ranges(other.num_nodes());
+        return *this;
+    }
+
+    // copy constructors take advantage of the assignment operators
+    // defined above
     tree(tree const& other)
-    :   data_(other.data_)
     {
-        set_ranges(other.num_nodes());
+        *this = other;
     }
 
     tree(tree&& other)
-    : data_(std::move(other.data_))
     {
-        set_ranges(other.num_nodes());
+        *this = std::move(other);
     }
 
-    tree() = default;
-
     /// create the tree from a parent_index
     template <typename I>
     std::vector<I>
@@ -135,7 +149,10 @@ class tree {
         return child_index_[b+1] - child_index_[b];
     }
     size_t num_nodes() const {
-        return child_index_.size() - 1;
+        // the number of nodes is the size of the child index minus 1
+        // ... except for the case of an empty tree
+        auto sz = child_index_.size();
+        return sz ? sz - 1 : 0;
     }
 
     /// return the child index
@@ -219,6 +236,7 @@ class tree {
     }
 
     void set_ranges(int nnode) {
+        if(nnode) {
             auto nchild = nnode - 1;
             // data_ is partitioned as follows:
             // data_ = [children_[nchild], child_index_[nnode+1], parents_[nnode]]
@@ -233,6 +251,12 @@ class tree {
             assert(child_index_.size() == unsigned(nnode+1));
             assert(parents_.size()     == unsigned(nnode));
         }
+        else {
+            children_    = data_(0, 0);
+            child_index_ = data_(0, 0);
+            parents_     = data_(0, 0);
+        }
+    }
 
     /// Renumber the sub-tree with old_node as its root with new_node as
     /// the new index of old_node. This is a helper function for the
@@ -304,15 +328,24 @@ class tree {
             }
         }
         if(add_parent_as_child) {
-            new_node = add_children(new_node, old_tree.parent(old_node), old_node, p, old_tree);
+            new_node =
+                add_children(
+                    new_node, old_tree.parent(old_node), old_node, p, old_tree
+                );
         }
 
         return new_node;
     }
 
+    //////////////////////////////////////////////////
+    // state
+    //////////////////////////////////////////////////
     index_type data_;
 
-    index_view children_;
-    index_view child_index_;
-    index_view parents_;
+    // provide default parameters so that tree type can
+    // be default constructed
+    index_view children_   = data_(0,0);
+    index_view child_index_= data_(0,0);
+    index_view parents_    = data_(0,0);
 };
+

tests/CMakeLists.txt

@@ -18,6 +18,7 @@ set(TEST_SOURCES
     test_segment.cpp
     test_swcio.cpp
     test_tree.cpp
+    test_run.cpp
 
     # unit test driver
     main.cpp

tests/test_cell.cpp

@@ -57,7 +57,7 @@ TEST(cell_type, add_segment)
                 segmentKind::dendrite,
                 cable_radius, cable_radius, cable_length
             );
-        c.add_cable(std::move(seg), 0);
+        c.add_cable(0, std::move(seg));
 
         EXPECT_EQ(c.num_segments(), 2);
     }
@@ -100,3 +100,66 @@ TEST(cell_type, add_segment)
         EXPECT_EQ(c.num_segments(), 2);
     }
 }
+
+TEST(cell_type, multiple_cables)
+{
+    using namespace nestmc;
+
+    // generate a cylindrical cable segment of length 1/pi and radius 1
+    //      volume = 1
+    //      area   = 2
+    auto seg = [](segmentKind k) {
+        return make_segment<cable_segment>( k, 1.0, 1.0, 1./math::pi<double>() );
+    };
+
+    //  add a pre-defined segment
+    {
+        cell c;
+
+        auto soma_radius = std::pow(3./(4.*math::pi<double>()), 1./3.);
+
+        // cell strucure as follows
+        // left   :  segment numbering
+        // right  :  segment type (soma, axon, dendrite)
+        //
+        //          0           s
+        //         / \         / \.
+        //        1   2       d   a
+        //       / \         / \.
+        //      3   4       d   d
+
+        // add a soma
+        c.add_soma(soma_radius, {0,0,1});
+
+        // hook the dendrite and axons
+        c.add_cable(0, seg(segmentKind::dendrite));
+        c.add_cable(0, seg(segmentKind::axon));
+        c.add_cable(1, seg(segmentKind::dendrite));
+        c.add_cable(1, seg(segmentKind::dendrite));
+
+        EXPECT_EQ(c.num_segments(), 5);
+        // each of the 5 segments has volume 1 by design
+        EXPECT_EQ(c.volume(), 5.);
+        // each of the 4 cables has volume 2., and the soma has an awkward area
+        // that isn't a round number
+        EXPECT_EQ(c.area(), 8. + math::area_sphere(soma_radius));
+
+        // construct the graph
+        c.construct();
+        auto const& con = c.graph();
+
+        EXPECT_EQ(con.num_segments(), 5u);
+        EXPECT_EQ(con.parent(0), -1);
+        EXPECT_EQ(con.parent(1), 0);
+        EXPECT_EQ(con.parent(2), 0);
+        EXPECT_EQ(con.parent(3), 1);
+        EXPECT_EQ(con.parent(4), 1);
+        EXPECT_EQ(con.num_children(0), 2u);
+        EXPECT_EQ(con.num_children(1), 2u);
+        EXPECT_EQ(con.num_children(2), 0u);
+        EXPECT_EQ(con.num_children(3), 0u);
+        EXPECT_EQ(con.num_children(4), 0u);
+
+    }
+}
+

tests/test_run.cpp

+#include "gtest.h"
+
+#include "../src/cell.hpp"
+
+TEST(run, init)
+{
+    using namespace nestmc;
+
+    nestmc::cell cell;
+
+    cell.add_soma(18.8);
+    auto& props = cell.soma()->properties;
+
+    cell.construct();
+
+    EXPECT_EQ(cell.graph().num_segments(), 1u);
+}