diff --git a/src/cell.hpp b/src/cell.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..df22c1aa52d8a8c7e8ff2da302c07cfb1b7d772e
--- /dev/null
+++ b/src/cell.hpp
@@ -0,0 +1,114 @@
+#pragma once
+
+#include <vector>
+#include <stdexcept>
+
+#include "segment.hpp"
+#include "cell_tree.hpp"
+
+namespace nestmc {
+
+ // we probably need two cell types
+ // 1. the abstract cell type (which would be this one)
+ // 2.
+ class cell {
+ public:
+ using index_type = int16_t;
+ using value_type = double;
+ using point_type = point<value_type>;
+
+ int num_segments() const
+ {
+ return segments_.size();
+ }
+
+ 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);
+ }
+
+ 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);
+ }
+
+ bool has_soma() const { return soma_ > -1; }
+
+ soma_segment* soma() {
+ if(has_soma()) {
+ return segments_[soma_].get()->as_soma();
+ }
+ return nullptr;
+ }
+
+ private:
+
+ //
+ // the local description of the cell which can be modified by the user
+ // in a ad-hoc manner (adding segments, modifying segments, etc)
+ //
+
+ // storage for connections
+ std::vector<index_type> parents_;
+ // the segments
+ std::vector<segment_ptr> segments_;
+ // index of the soma
+ int soma_ = -1;
+
+ //
+ // fixed cell description, which is computed from the
+ // rough layout description above
+ //
+
+ // cell_tree that describes the connection layout
+ //cell_tree tree_;
+ };
+
+} // namespace nestmc
diff --git a/src/cell_tree.hpp b/src/cell_tree.hpp
index 0d4c496f459821c59caa80319d1536bba71b3671..c09f647bde9317dedb4ac924585c1fcf54463062 100644
--- a/src/cell_tree.hpp
+++ b/src/cell_tree.hpp
@@ -50,13 +50,15 @@ class cell_tree {
/// construct from a tree
// copy constructor
- cell_tree(tree const& t)
- : tree_(t)
+ cell_tree(tree const& t, int s)
+ : tree_(t),
+ soma_(s)
{ }
// move constructor
- cell_tree(tree&& t)
- : tree_(std::move(t))
+ cell_tree(tree&& t, int s)
+ : tree_(std::move(t)),
+ soma_(s)
{ }
/// construct from a cell tree
@@ -166,12 +168,11 @@ class cell_tree {
}
std::string to_string() const {
- std::string s;
-
- s += "[" + std::to_string(root) + ","
- + std::to_string(diameter) + "," + std::to_string(depth) + "]";
-
- return s;
+ return
+ "[" + std::to_string(root) + ","
+ + std::to_string(diameter) + ","
+ + std::to_string(depth) +
+ "]";
}
int root;
diff --git a/src/math.hpp b/src/math.hpp
index d936824966c59d18758f232bd4b351ce77cdeca1..e28c135d9e03d00516d36ad4f6e1254d40fa5f54 100644
--- a/src/math.hpp
+++ b/src/math.hpp
@@ -3,24 +3,60 @@
#include <cmath>
namespace nestmc {
-
+namespace math {
template <typename T>
T constexpr pi() {
- return 3.141592653589793238462643383279502;
+ return T(3.1415926535897932384626433832795);
+ }
+
+ template <typename T>
+ T constexpr mean(T a, T b) {
+ return (a+b) / T(2);
+ }
+
+ template <typename T>
+ T constexpr square(T a) {
+ return a*a;
+ }
+
+ template <typename T>
+ T constexpr cube(T a) {
+ return a*a*a;
}
+ // volume of a conic frustrum
template <typename T>
- T area_frustrum(T L, T r1, T r2)
+ T constexpr area_frustrum(T L, T r1, T r2)
{
- auto dr = r1 - r2;
- return pi<double>() * (r1+r2) * sqrt(L*L + dr*dr);
+ return pi<T>() * (r1+r2) * sqrt(square(L) + square(r1-r2));
}
+ // surface area of conic frustrum, not including the area of the
+ // circles at either end (i.e. just the area of the surface of rotation)
template <typename T>
- T volume_frustrum(T L, T r1, T r2)
+ T constexpr volume_frustrum(T L, T r1, T r2)
{
- auto meanr = (r1+r2) / 2.;
- return pi<double>() * meanr * meanr * L;
+ // this formulation uses one less multiplication
+ return pi<T>()/T(3) * (square(r1+r2) - r1*r2) * L;
+ //return pi<T>()/T(3) * (square(r1) + r1*r2 + square(r2)) * L;
}
+ // volume of a sphere
+ // four-thirds pi r-cubed
+ template <typename T>
+ T constexpr volume_sphere(T r)
+ {
+ return T(4)/T(3) * pi<T>() * cube(r);
+ }
+
+ // surface area of a sphere
+ // four pi r-squared
+ template <typename T>
+ T constexpr area_sphere(T r)
+ {
+ return T(4) * pi<T>() * square(r);
+ }
+
+} // namespace nestmc::math
} // namespace nestmc
+
diff --git a/src/point.hpp b/src/point.hpp
index fc3f0daea1f6d44779ecf9a36784ab3e48580eae..00ded16ca2b70e3179cc18bb26a004344972c02b 100644
--- a/src/point.hpp
+++ b/src/point.hpp
@@ -11,7 +11,7 @@ struct point {
value_type z;
constexpr point(T a, T b, T c)
- : x(a), y(b), z(c)
+ : x{a}, y{b}, z{c}
{}
// initialize to NaN by default
@@ -20,6 +20,10 @@ struct point {
y(std::numeric_limits<T>::quiet_NaN()),
z(std::numeric_limits<T>::quiet_NaN())
{}
+
+ constexpr bool is_set() const {
+ return (x==x && y==y && z==z);
+ }
};
template <typename T>
diff --git a/src/segment.hpp b/src/segment.hpp
index 6cefbd78cd45b8013b69a2f07a5ce18211fd8ed7..dfb4cf49195f8d1f0774c901e62cf6db1ff03ccd 100644
--- a/src/segment.hpp
+++ b/src/segment.hpp
@@ -36,6 +36,11 @@ enum class segmentKind {
axon
};
+// forward declarations of segment specializations
+class soma_segment;
+class cable_segment;
+
+// abstract base class for a cell segment
class segment {
public:
@@ -46,17 +51,42 @@ class segment {
return kind_;
}
+ bool is_soma() const
+ {
+ return kind_==segmentKind::soma;
+ }
+
+ bool is_dendrite() const
+ {
+ return kind_==segmentKind::dendrite;
+ }
+
+ bool is_axon() const
+ {
+ return kind_==segmentKind::axon;
+ }
+
virtual value_type volume() const = 0;
- virtual value_type area() const = 0;
+ virtual value_type area() const = 0;
virtual ~segment() = default;
+ virtual cable_segment* as_cable()
+ {
+ return nullptr;
+ }
+
+ virtual soma_segment* as_soma()
+ {
+ return nullptr;
+ }
+
protected:
segmentKind kind_;
};
-class spherical_segment : public segment
+class soma_segment : public segment
{
public :
@@ -65,16 +95,16 @@ class spherical_segment : public segment
using base::value_type;
using base::point_type;
- spherical_segment() = delete;
+ soma_segment() = delete;
- spherical_segment(value_type r)
+ soma_segment(value_type r)
: radius_{r}
{
kind_ = segmentKind::soma;
}
- spherical_segment(value_type r, point_type const &c)
- : spherical_segment(r)
+ soma_segment(value_type r, point_type const &c)
+ : soma_segment(r)
{
center_ = c;
kind_ = segmentKind::soma;
@@ -82,15 +112,28 @@ class spherical_segment : public segment
value_type volume() const override
{
- return 4./3. * pi<value_type>() * radius_ * radius_ * radius_;
+ return math::volume_sphere(radius_);
}
value_type area() const override
{
- return 4. * pi<value_type>() * radius_ * radius_;
+ return math::area_sphere(radius_);
+ }
+
+ value_type radius() const
+ {
+ return radius_;
}
- virtual ~spherical_segment() = default;
+ point_type const& center() const
+ {
+ return center_;
+ }
+
+ soma_segment* as_soma() override
+ {
+ return this;
+ }
private :
@@ -100,61 +143,121 @@ class spherical_segment : public segment
point_type center_;
};
-class frustrum_segment : public segment
+class cable_segment : public segment
{
public :
- using segment::kind_;
using base = segment;
+ using base::kind_;
using base::value_type;
using base::point_type;
- frustrum_segment() = delete;
+ // delete the default constructor
+ cable_segment() = delete;
- frustrum_segment(
+ // constructors for a cable with no location information
+ cable_segment(
+ segmentKind k,
+ std::vector<value_type> r,
+ std::vector<value_type> lens
+ ) {
+ kind_ = k;
+ assert(k==segmentKind::dendrite || k==segmentKind::axon);
+
+ radii_ = std::move(r);
+ lengths_ = std::move(lens);
+ }
+
+ cable_segment(
segmentKind k,
value_type r1,
value_type r2,
value_type len
+ )
+ : cable_segment{k, {r1, r2}, std::vector<value_type>{len}}
+ { }
+
+ // constructor that lets the user describe the cable as a
+ // seriew of radii and locations
+ cable_segment(
+ segmentKind k,
+ std::vector<value_type> r,
+ std::vector<point_type> p
) {
- r1_ = r1;
- r2_ = r2;
- length_ = len;
kind_ = k;
- assert(k==segmentKind::dendrite || k ==segmentKind::axon);
+ assert(k==segmentKind::dendrite || k==segmentKind::axon);
+
+ radii_ = std::move(r);
+ locations_ = std::move(p);
+ update_lengths();
}
- frustrum_segment(
+ // helper that lets user specify with the radius and location
+ // of just the end points of the cable
+ // i.e. describing the cable as a single frustrum
+ cable_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))
+ : cable_segment(k, {r1, r2}, {p1, p2})
+ { }
+
+ value_type volume() const override
{
- p1_ = p1;
- p2_ = p2;
+ auto sum = value_type{0};
+ for(auto i=0; i<num_sub_segments(); ++i) {
+ sum += math::volume_frustrum(lengths_[i], radii_[i], radii_[i+1]);
+ }
+ return sum;
+
}
- value_type volume() const override
+ value_type area() const override
{
- return volume_frustrum(length_, r1_, r2_);
+ auto sum = value_type{0};
+ for(auto i=0; i<num_sub_segments(); ++i) {
+ sum += math::area_frustrum(lengths_[i], radii_[i], radii_[i+1]);
+ }
+ return sum;
+ }
+
+ bool has_locations() const {
+ return locations_.size() > 0;
+ }
+
+ // the number sub-segments that define the cable segment
+ int num_sub_segments() const {
+ return radii_.size()-1;
}
- value_type area() const override
+
+ std::vector<value_type> const& lengths() const
{
- return area_frustrum(length_, r1_, r2_);
+ return lengths_;
}
- virtual ~frustrum_segment() = default;
+ cable_segment* as_cable() override
+ {
+ return this;
+ }
private :
- value_type length_;
- value_type r1_;
- value_type r2_;
- point_type p1_;
- point_type p2_;
+ void update_lengths()
+ {
+ if(locations_.size()) {
+ lengths_.resize(num_sub_segments());
+ for(auto i=0; i<num_sub_segments(); ++i) {
+ lengths_[i] = norm(locations_[i] - locations_[i+1]);
+ }
+ }
+ }
+
+ std::vector<value_type> lengths_;
+ std::vector<value_type> radii_;
+ std::vector<point_type> locations_;
};
using segment_ptr = std::unique_ptr<segment>;
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 4c921c290c5459bce2f48306843d307f2a6470c7..d246268247d147eb7e608296af52bab191b616f7 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -1,10 +1,11 @@
set(HEADERS
- ../src/swcio.hpp
- ../src/tree.hpp
+ ../src/cell.hpp
../src/cell_tree.hpp
- ../src/point.hpp
../src/math.hpp
+ ../src/point.hpp
../src/segment.hpp
+ ../src/swcio.hpp
+ ../src/tree.hpp
)
set(TEST_SOURCES
@@ -12,10 +13,11 @@ set(TEST_SOURCES
gtest-all.cpp
# unit tests
- test_swcio.cpp
- test_tree.cpp
+ test_cell.cpp
test_point.cpp
test_segment.cpp
+ test_swcio.cpp
+ test_tree.cpp
# unit test driver
main.cpp
diff --git a/tests/test_cell.cpp b/tests/test_cell.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..dfc3481712bab208104c9cbedb2de90ecb02c83a
--- /dev/null
+++ b/tests/test_cell.cpp
@@ -0,0 +1,102 @@
+#include "gtest.h"
+
+#include "../src/cell.hpp"
+
+TEST(cell_type, soma)
+{
+ // test that insertion of a soma works
+ // define with no centre point
+ {
+ nestmc::cell c;
+ auto soma_radius = 2.1;
+
+ EXPECT_EQ(c.has_soma(), false);
+ c.add_soma(soma_radius);
+ EXPECT_EQ(c.has_soma(), true);
+
+ auto s = c.soma();
+ EXPECT_EQ(s->radius(), soma_radius);
+ EXPECT_EQ(s->center().is_set(), false);
+ }
+
+ // test that insertion of a soma works
+ // define with centre point @ (0,0,1)
+ {
+ nestmc::cell c;
+ auto soma_radius = 3.2;
+
+ EXPECT_EQ(c.has_soma(), false);
+ c.add_soma(soma_radius, {0,0,1});
+ EXPECT_EQ(c.has_soma(), true);
+
+ auto s = c.soma();
+ EXPECT_EQ(s->radius(), soma_radius);
+ EXPECT_EQ(s->center().is_set(), true);
+
+ // add expression template library for points
+ //EXPECT_EQ(s->center(), point<double>(0,0,1));
+ }
+}
+
+TEST(cell_type, add_segment)
+{
+ using namespace nestmc;
+ // add a pre-defined segment
+ {
+ cell c;
+
+ auto soma_radius = 2.1;
+ auto cable_radius = 0.1;
+ auto cable_length = 8.3;
+
+ // add a soma because we need something to attach the first dendrite to
+ c.add_soma(soma_radius, {0,0,1});
+
+ auto seg =
+ make_segment<cable_segment>(
+ segmentKind::dendrite,
+ cable_radius, cable_radius, cable_length
+ );
+ c.add_cable(std::move(seg), 0);
+
+ EXPECT_EQ(c.num_segments(), 2);
+ }
+
+ // add segment on the fly
+ {
+ cell c;
+
+ auto soma_radius = 2.1;
+ auto cable_radius = 0.1;
+ auto cable_length = 8.3;
+
+ // add a soma because we need something to attach the first dendrite to
+ c.add_soma(soma_radius, {0,0,1});
+
+ c.add_cable(
+ 0,
+ segmentKind::dendrite, cable_radius, cable_radius, cable_length
+ );
+
+ EXPECT_EQ(c.num_segments(), 2);
+ }
+ {
+ cell c;
+
+ auto soma_radius = 2.1;
+ auto cable_radius = 0.1;
+ auto cable_length = 8.3;
+
+ // add a soma because we need something to attach the first dendrite to
+ c.add_soma(soma_radius, {0,0,1});
+
+ c.add_cable(
+ 0,
+ segmentKind::dendrite,
+ std::vector<double>{cable_radius, cable_radius, cable_radius, cable_radius},
+ std::vector<double>{cable_length, cable_length, cable_length}
+ );
+
+ EXPECT_EQ(c.num_segments(), 2);
+ }
+}
diff --git a/tests/test_point.cpp b/tests/test_point.cpp
index 9fbb62a2fd6819caf1e1b97a2c62013d98886f24..f795b72de29ff6bac8c931745f1d02a3d5866a5a 100644
--- a/tests/test_point.cpp
+++ b/tests/test_point.cpp
@@ -17,7 +17,7 @@ TEST(point, construction)
{
// initializer list
- point<float> p = {1, 2, 3};
+ point<float> p{1, 2, 3};
EXPECT_EQ(p.x, 1.);
EXPECT_EQ(p.y, 2.);
EXPECT_EQ(p.z, 3.);
diff --git a/tests/test_segment.cpp b/tests/test_segment.cpp
index 94dabca0b3f113124116dfa19b454aed98b47222..81f70aebe63cf0e97d3122f400c0d80e0ab6f593 100644
--- a/tests/test_segment.cpp
+++ b/tests/test_segment.cpp
@@ -4,29 +4,95 @@
#include "../src/segment.hpp"
-TEST(segments, sphere)
+
+TEST(segments, soma)
{
using namespace nestmc;
+ using nestmc::math::pi;
{
- auto s = make_segment<spherical_segment>(1.0);
+ auto s = make_segment<soma_segment>(1.0);
- EXPECT_EQ(s->volume(), nestmc::pi<double>()*4./3.);
- EXPECT_EQ(s->area(), nestmc::pi<double>()*4.);
- EXPECT_EQ(s->kind(), nestmc::segmentKind::soma);
+ EXPECT_EQ(s->volume(), pi<double>()*4./3.);
+ EXPECT_EQ(s->area(), pi<double>()*4.);
+ EXPECT_EQ(s->kind(), segmentKind::soma);
}
{
- auto s = make_segment<spherical_segment>(1.0, point<double>(0., 1., 2.));
+ auto s = make_segment<soma_segment>(1.0, point<double>(0., 1., 2.));
- EXPECT_EQ(s->volume(), nestmc::pi<double>()*4./3.);
- EXPECT_EQ(s->area(), nestmc::pi<double>()*4.);
- EXPECT_EQ(s->kind(), nestmc::segmentKind::soma);
+ EXPECT_EQ(s->volume(), pi<double>()*4./3.);
+ EXPECT_EQ(s->area(), pi<double>()*4.);
+ EXPECT_EQ(s->kind(), segmentKind::soma);
}
}
-TEST(segments, frustrum)
+TEST(segments, cable)
{
+ using namespace nestmc;
+ using nestmc::math::pi;
+
+ // take advantage of fact that a cable segment with constant radius 1 and
+ // length 1 has volume=1. and area=2
+ auto length = 1./pi<double>();
+ auto radius = 1.;
+
+ // single cylindrical frustrum
{
+ auto s = make_segment<cable_segment>(segmentKind::dendrite, radius, radius, length);
+
+ EXPECT_EQ(s->volume(), 1.0);
+ EXPECT_EQ(s->area(), 2.0);
+ EXPECT_EQ(s->kind(), segmentKind::dendrite);
+ }
+
+ // cable made up of three identical cylindrical frustrums
+ {
+ auto s =
+ make_segment<cable_segment>(
+ segmentKind::axon,
+ std::vector<double>{radius, radius, radius, radius},
+ std::vector<double>{length, length, length}
+ );
+
+ EXPECT_EQ(s->volume(), 3.0);
+ EXPECT_EQ(s->area(), 6.0);
+ EXPECT_EQ(s->kind(), segmentKind::axon);
+ }
+}
+
+TEST(segments, cable_positions)
+{
+ using namespace nestmc;
+ using nestmc::math::pi;
+
+ // single frustrum of length 1 and radii 1 and 2
+ // the centre of each end are at the origin (0,0,0) and (0,1,0)
+ {
+ auto s =
+ make_segment<cable_segment>(
+ segmentKind::dendrite,
+ 1, 2,
+ point<double>(0,0,0), point<double>(0,1,0)
+ );
+
+ EXPECT_EQ(s->volume(), math::volume_frustrum(1., 1., 2.));
+ EXPECT_EQ(s->area(), math::area_frustrum (1., 1., 2.));
+ EXPECT_EQ(s->kind(), segmentKind::dendrite);
+ }
+
+ // cable made up of three frustrums
+ // that emulate the frustrum from the previous single-frustrum case
+ {
+ auto s =
+ make_segment<cable_segment>(
+ segmentKind::axon,
+ std::vector<double>{1, 1.5, 2},
+ std::vector<point<double>>{ {0,0,0}, {0,0.5,0}, {0,1,0} }
+ );
+
+ EXPECT_EQ(s->volume(), math::volume_frustrum(1., 1., 2.));
+ EXPECT_EQ(s->area(), math::area_frustrum(1., 1., 2.));
+ EXPECT_EQ(s->kind(), segmentKind::axon);
}
}