diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 90c9d8f8c4446615bfd3e045cd64a082a0d1d016..457f38aa6b2e9d6c898a759e821fdf5e06326bd9 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -3,6 +3,7 @@ set(HEADERS
)
set(BASE_SOURCES
swcio.cpp
+ cell.cpp
)
add_library(cellalgo ${BASE_SOURCES} ${HEADERS})
diff --git a/src/cell.cpp b/src/cell.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..5d9cec1b22e43d265f89ec3998f3c33e566f6ebc
--- /dev/null
+++ b/src/cell.cpp
@@ -0,0 +1,115 @@
+#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
diff --git a/src/cell.hpp b/src/cell.hpp
index df22c1aa52d8a8c7e8ff2da302c07cfb1b7d772e..bc3d49e6289099b59e24e813c1a91e6b386c16c5 100644
--- a/src/cell.hpp
+++ b/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"
- );
- }
-
- // 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)
- );
- }
- }
-
- 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"
- );
- }
-
- // 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);
- }
+ /// 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);
+ /// 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)
- {
- // 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);
- }
+ void add_cable(index_type parent, Args ...args);
- bool has_soma() const { return soma_ > -1; }
+ /// the number of segments in the cell
+ int num_segments() const;
- soma_segment* soma() {
- if(has_soma()) {
- return segments_[soma_].get()->as_soma();
- }
- return nullptr;
- }
+ bool has_soma() const;
+
+ /// access pointer to the soma
+ soma_segment* soma();
+
+ /// the volume of the cell
+ value_type volume() const;
+
+ /// the surface area of the cell
+ value_type area() const;
+
+ 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
+
diff --git a/src/cell_tree.hpp b/src/cell_tree.hpp
index c09f647bde9317dedb4ac924585c1fcf54463062..16f7e43f8141375e64855b8bf99d8dd63eea7b59 100644
--- a/src/cell_tree.hpp
+++ b/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;
};
+
diff --git a/src/math.hpp b/src/math.hpp
index e28c135d9e03d00516d36ad4f6e1254d40fa5f54..c32c029eb5f03c9f1978b83caa836a8623bb9dc6 100644
--- a/src/math.hpp
+++ b/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;
}
diff --git a/src/segment.hpp b/src/segment.hpp
index dfb4cf49195f8d1f0774c901e62cf6db1ff03ccd..7710ed9714aa112dea731cf27bfd648027bcca16 100644
--- a/src/segment.hpp
+++ b/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_;
diff --git a/src/tree.hpp b/src/tree.hpp
index 6b53aa6f329490a3c431817445f9b2be4dda25b1..eca2fca6c7e8dd198974128784a6cc2b03b5db3b 100644
--- a/src/tree.hpp
+++ b/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,19 +236,26 @@ class tree {
}
void set_ranges(int nnode) {
- auto nchild = nnode - 1;
- // data_ is partitioned as follows:
- // data_ = [children_[nchild], child_index_[nnode+1], parents_[nnode]]
- assert(data_.size() == unsigned(nchild + (nnode+1) + nnode));
- children_ = data_(0, nchild);
- child_index_ = data_(nchild, nchild+nnode+1);
- parents_ = data_(nchild+nnode+1, memory::end);
-
- // check that arrays have appropriate size
- // this should be moved into a unit test
- assert(children_.size() == unsigned(nchild));
- assert(child_index_.size() == unsigned(nnode+1));
- assert(parents_.size() == unsigned(nnode));
+ if(nnode) {
+ auto nchild = nnode - 1;
+ // data_ is partitioned as follows:
+ // data_ = [children_[nchild], child_index_[nnode+1], parents_[nnode]]
+ assert(data_.size() == unsigned(nchild + (nnode+1) + nnode));
+ children_ = data_(0, nchild);
+ child_index_ = data_(nchild, nchild+nnode+1);
+ parents_ = data_(nchild+nnode+1, memory::end);
+
+ // check that arrays have appropriate size
+ // this should be moved into a unit test
+ assert(children_.size() == unsigned(nchild));
+ 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
@@ -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);
};
+
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index d246268247d147eb7e608296af52bab191b616f7..3bfe22c59526798ea1108ad1862a7e418a734182 100644
--- a/tests/CMakeLists.txt
+++ b/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
diff --git a/tests/test_cell.cpp b/tests/test_cell.cpp
index dfc3481712bab208104c9cbedb2de90ecb02c83a..d9cb96b8e816d1942438d98ff13a5ee3e1bba940 100644
--- a/tests/test_cell.cpp
+++ b/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);
+
+ }
+}
+
diff --git a/tests/test_run.cpp b/tests/test_run.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..9ac528f7a5768c08b8043d019f9207b7339506e7
--- /dev/null
+++ b/tests/test_run.cpp
@@ -0,0 +1,17 @@
+#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);
+}