From 85ee79c8a0299f8a9ba918607a34836ac0ca17c4 Mon Sep 17 00:00:00 2001
From: bcumming <bcumming@cscs.ch>
Date: Fri, 18 Mar 2016 16:46:38 +0100
Subject: [PATCH] start working on compartment model
---
src/compartment.hpp | 127 +++++++++++++++++++++++++++++++++
src/math.hpp | 26 +++++++
src/point.hpp | 68 ++++++++++++++++++
src/segment.hpp | 168 ++++++++++++++++++++++++++++++++++++++++++++
src/util.hpp | 10 +++
5 files changed, 399 insertions(+)
create mode 100644 src/compartment.hpp
create mode 100644 src/math.hpp
create mode 100644 src/point.hpp
create mode 100644 src/segment.hpp
diff --git a/src/compartment.hpp b/src/compartment.hpp
new file mode 100644
index 00000000..350caed3
--- /dev/null
+++ b/src/compartment.hpp
@@ -0,0 +1,127 @@
+#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
+
+
diff --git a/src/math.hpp b/src/math.hpp
new file mode 100644
index 00000000..d9368249
--- /dev/null
+++ b/src/math.hpp
@@ -0,0 +1,26 @@
+#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
diff --git a/src/point.hpp b/src/point.hpp
new file mode 100644
index 00000000..fc3f0dae
--- /dev/null
+++ b/src/point.hpp
@@ -0,0 +1,68 @@
+#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 << "]";
+}
+
+
diff --git a/src/segment.hpp b/src/segment.hpp
new file mode 100644
index 00000000..6cefbd78
--- /dev/null
+++ b/src/segment.hpp
@@ -0,0 +1,168 @@
+#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
+
diff --git a/src/util.hpp b/src/util.hpp
index 17a3a49b..580981be 100644
--- a/src/util.hpp
+++ b/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) ...));
+ }
+}
+
--
GitLab