diff --git a/src/algorithms.hpp b/src/algorithms.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..db8b33988f9a12b8cb49d5f510c5661c5904f18f
--- /dev/null
+++ b/src/algorithms.hpp
@@ -0,0 +1,103 @@
+#pragma once
+
+#include <algorithm>
+#include <type_traits>
+
+/*
+ * Some simple wrappers around stl algorithms to improve readability of code
+ * that uses them.
+ *
+ * For example, a simple sum() wrapper for finding the sum of a container,
+ * is much more readable than using std::accumulate()
+ *
+ */
+
+namespace nest {
+namespace mc {
+namespace algorithms{
+
+ template <typename C>
+ typename C::value_type
+ sum(C const& c)
+ {
+ using value_type = typename C::value_type;
+ return std::accumulate(c.begin(), c.end(), value_type{0});
+ }
+
+ template <typename C>
+ typename C::value_type
+ mean(C const& c)
+ {
+ return sum(c)/c.size();
+ }
+
+ template <typename C>
+ C make_index(C const& c)
+ {
+ C out(c.size()+1);
+ out[0] = 0;
+ std::partial_sum(c.begin(), c.end(), out.begin()+1);
+ return out;
+ }
+
+ /// works like std::is_sorted(), but with stronger condition that succesive
+ /// elements must be greater than those before them
+ template <typename C>
+ bool is_strictly_monotonic_increasing(C const& c)
+ {
+ using value_type = typename C::value_type;
+ return std::is_sorted(
+ c.begin(),
+ c.end(),
+ [] (value_type const& lhs, value_type const& rhs) {
+ return lhs <= rhs;
+ }
+ );
+ }
+
+ template <typename C>
+ bool is_strictly_monotonic_decreasing(C const& c)
+ {
+ using value_type = typename C::value_type;
+ return std::is_sorted(
+ c.begin(),
+ c.end(),
+ [] (value_type const& lhs, value_type const& rhs) {
+ return lhs >= rhs;
+ }
+ );
+ }
+
+ template <typename C>
+ bool is_minimal_degree(C const& c)
+ {
+ static_assert(
+ std::is_integral<typename C::value_type>::value,
+ "is_minimal_degree only applies to integral types"
+ );
+ for(auto i=0; i<c.size(); ++i) {
+ if(i<c[i]) {
+ return false;
+ }
+ }
+ return true;
+ }
+
+ template <typename C>
+ bool is_positive(C const& c)
+ {
+ static_assert(
+ std::is_integral<typename C::value_type>::value,
+ "is_positive only applies to integral types"
+ );
+ for(auto i=0; i<c.size(); ++i) {
+ if(c[i]<1) {
+ return false;
+ }
+ }
+ return true;
+ }
+
+} // namespace algorithms
+} // namespace mc
+} // namespace nest
diff --git a/src/cell.cpp b/src/cell.cpp
index c0d51bc32741c17b73da8215d7076b5325b21cb9..a124eb384837e35215823bc8812f399905f56b68 100644
--- a/src/cell.cpp
+++ b/src/cell.cpp
@@ -1,12 +1,15 @@
#include "cell.hpp"
+#include "tree.hpp"
-namespace nestmc {
+namespace nest {
+namespace mc {
cell::cell()
{
// insert a placeholder segment for the soma
segments_.push_back(make_segment<placeholder_segment>());
parents_.push_back(0);
+ stale_ = true;
}
int cell::num_segments() const
@@ -33,6 +36,7 @@ void cell::add_soma(value_type radius, point_type center)
else {
segments_[0] = make_segment<soma_segment>(radius);
}
+ stale_ = true;
}
void cell::add_cable(cell::index_type parent, segment_ptr&& cable)
@@ -52,6 +56,7 @@ void cell::add_cable(cell::index_type parent, segment_ptr&& cable)
}
segments_.push_back(std::move(cable));
parents_.push_back(parent);
+ stale_ = true;
}
segment* cell::segment(int index)
@@ -82,14 +87,21 @@ bool cell::has_soma() const
soma_segment* cell::soma()
{
+ stale_ = true;
if(has_soma()) {
return segment(0)->as_soma();
}
return nullptr;
}
+// this breaks the use of stale_, because the user could modify a cable obtained from
+// this call, without the stale_ tag being set.
+//
+// set stale_ to true, then expect the user to make any modifications before calling
+// graph() etc.
cable_segment* cell::cable(int index)
{
+ stale_ = true;
if(index>0 && index<num_segments()) {
return segment(index)->as_cable();
}
@@ -125,14 +137,43 @@ std::vector<segment_ptr> const& cell::segments() const
return segments_;
}
+std::vector<int> cell::compartment_counts() const
+{
+ std::vector<int> comp_count;
+ comp_count.reserve(num_segments());
+ for(auto const& s : segments()) {
+ comp_count.push_back(s->num_compartments());
+ }
+ return comp_count;
+
+}
+
void cell::construct() const
{
- if(num_segments()) {
+ std::lock_guard<std::mutex> lock(mutex_);
+
+ if(stale_) {
tree_ = cell_tree(parents_);
+ auto counts = compartment_counts();
+ parent_index_ = make_parent_index(tree_.graph(), counts);
+ segment_index_ = algorithms::make_index(counts);
}
+ stale_ = false;
+}
+
+std::vector<int> const& cell::parent_index() const
+{
+ construct();
+ return parent_index_;
+}
+
+std::vector<int> const& cell::segment_index() const
+{
+ construct();
+ return segment_index_;
}
-cell_tree const& cell::graph() const
+cell_tree const& cell::tree() const
{
construct();
return tree_;
@@ -143,4 +184,5 @@ std::vector<int> const& cell::segment_parents() const
return parents_;
}
-} // namespace nestmc
+} // namespace mc
+} // namespace nest
diff --git a/src/cell.hpp b/src/cell.hpp
index b3970ce0bbb2aa78fab1fcc1b3e9c6e375292e08..0a254e1ec006dec417b89ef36fa8f28631a283b7 100644
--- a/src/cell.hpp
+++ b/src/cell.hpp
@@ -8,108 +8,128 @@
#include "segment.hpp"
#include "cell_tree.hpp"
-namespace nestmc {
+namespace nest {
+namespace mc {
- // high-level abstract representation of a cell and its segments
- class cell {
- public:
+// high-level abstract representation of a cell and its segments
+class cell {
+ public:
- // types
- using index_type = int;
- using value_type = double;
- using point_type = point<value_type>;
+ // types
+ using index_type = int;
+ using value_type = double;
+ using point_type = point<value_type>;
- // constructor
- cell();
+ // constructor
+ cell();
- /// add a soma to the cell
- /// radius must be specified
- void add_soma(value_type radius, point_type center=point_type());
+ /// add a soma to the cell
+ /// radius must be specified
+ void add_soma(value_type radius, point_type center=point_type());
- /// 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
+ /// 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);
-
- /// the number of segments in the cell
- int num_segments() const;
-
- bool has_soma() const;
-
- class segment* segment(int index);
- class segment const* segment(int index) const;
-
- /// access pointer to the soma
- /// returns nullptr if the cell has no soma
- soma_segment* soma();
-
- /// access pointer to a cable segment
- /// will throw an std::out_of_range exception if
- /// the cable index is not valid
- cable_segment* cable(int index);
-
- /// 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;
-
- /// 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:
-
- /// generate the internal representation of the connectivity
- /// graph for the cell segments
- void construct() const;
-
-
- //
- // 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_;
-
- //
- // fixed cell description, which is computed from the layout description
- // this computed whenever a call to the graph() method is made
- // the graph method is const, so tree_ is mutable
- //
- // note: this isn't remotely thread safe!
-
- mutable cell_tree tree_;
- };
-
- // create a cable by forwarding cable construction parameters provided by the user
+ /// 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 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);
+ void add_cable(index_type parent, Args ...args);
+
+ /// the number of segments in the cell
+ int num_segments() const;
+
+ bool has_soma() const;
+
+ class segment* segment(int index);
+ class segment const* segment(int index) const;
+
+ /// access pointer to the soma
+ /// returns nullptr if the cell has no soma
+ soma_segment* soma();
+
+ /// access pointer to a cable segment
+ /// will throw an std::out_of_range exception if
+ /// the cable index is not valid
+ cable_segment* cable(int index);
+
+ /// 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;
+
+ /// the connectivity graph for the cell segments
+ cell_tree const& tree() const;
+
+ /// return reference to array that enumerates the index of the parent of
+ /// each segment
+ std::vector<int> const& segment_parents() const;
+
+ /// return a vector with the compartment count for each segment in the cell
+ std::vector<int> compartment_counts() const;
+
+ /// return the parent index for the compartments
+ std::vector<int> const& parent_index() const;
+
+ /// Return the segment index for the compartments
+ /// the segment index is an index into parent_index for looking
+ /// up the set of compartments associated with a segment.
+ /// i.e. the compartments for segment i are in the half open range
+ /// [segment_index()[i], segmend_index()[i+1])
+ std::vector<int> const& segment_index() const;
+
+ private:
+
+ /// generate the internal representation of the connectivity
+ /// graph for the cell segments
+ void construct() const;
+
+ //
+ // 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_;
+
+ // used internally to mark whether derived data (tree_, parent_index_, etc)
+ // are out of date
+ mutable bool stale_ = true;
+
+ //
+ // fixed cell description, which is computed from the layout description
+ // this computed whenever a call to the graph() method is made
+ // the graph method is const, so tree_ is mutable
+ //
+
+ mutable std::mutex mutex_;
+ mutable cell_tree tree_;
+ mutable std::vector<int> parent_index_;
+ mutable std::vector<int> segment_index_;
+};
+
+// 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"
+ );
}
+ stale_ = true;
+ segments_.push_back(make_segment<cable_segment>(std::forward<Args>(args)...));
+ parents_.push_back(parent);
+}
-
-} // namespace nestmc
+} // namespace mc
+} // namespace nest
diff --git a/src/cell_tree.hpp b/src/cell_tree.hpp
index 16f7e43f8141375e64855b8bf99d8dd63eea7b59..43fac2429a657545251159a3fa97002b7243fe86 100644
--- a/src/cell_tree.hpp
+++ b/src/cell_tree.hpp
@@ -13,6 +13,9 @@
#include "tree.hpp"
#include "util.hpp"
+namespace nest {
+namespace mc {
+
/// The tree data structure that describes the segments of a cell tree.
/// A cell is represented as a tree where each node may have any number of
/// children. Typically in a cell only the soma has more than two segments,
@@ -25,8 +28,8 @@
/// sets it appears that the soma was always index 0, however we need more
/// flexibility in choosing the root.
class cell_tree {
- public :
-
+ using range = memory::Range;
+public :
// use a signed 16-bit integer for storage of indexes, which is reasonable given
// that typical cells have at most 1000-2000 segments
using int_type = int16_t;
@@ -40,15 +43,12 @@ class cell_tree {
cell_tree(std::vector<int> const& parent_index)
{
// handle the case of an empty parent list, which implies a single-compartment model
- std::vector<int> segment_index;
if(parent_index.size()>0) {
- segment_index = tree_.init_from_parent_index(parent_index);
+ tree_ = tree(parent_index);
}
else {
- segment_index = tree_.init_from_parent_index(std::vector<int>({0}));
+ tree_ = tree(std::vector<int>({0}));
}
-
- // if needed, calculate meta-data like length[] and end[] arrays for data
}
/// construct from a tree
@@ -93,6 +93,10 @@ class cell_tree {
*this = std::move(other);
}
+ tree const& graph() const {
+ return tree_;
+ }
+
int_type soma() const {
return soma_;
}
@@ -169,8 +173,7 @@ class cell_tree {
return depth;
}
- private :
-
+private :
/// helper type for sub-tree computation
/// use in balance()
@@ -292,3 +295,6 @@ class cell_tree {
int_type soma_ = 0;
};
+} // namespace mc
+} // namespace nest
+
diff --git a/src/compartment.hpp b/src/compartment.hpp
index 350caed3277fadfe996922c2925649cf56e9c7d4..b92360a560ba5de30a036505d50d457b7b5a5a9c 100644
--- a/src/compartment.hpp
+++ b/src/compartment.hpp
@@ -1,127 +1,166 @@
#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()}
- {}
+#include <iterator>
+#include <utility>
+
+namespace nest {
+namespace mc {
+
+/// Defines the simplest type of compartment
+/// The compartment is a conic frustrum
+struct compartment {
+ using value_type = double;
+ using size_type = int;
+ using value_pair = std::pair<value_type, value_type>;
- constexpr compartment(
- point_type const& p1,
- point_type const& p2,
+ compartment() = delete;
+
+ compartment(
+ size_type idx,
+ value_type len,
value_type r1,
value_type r2
)
- : p1_{p1},
- p2_{p2},
- radius1_{r1},
- radius2_{r2}
+ : index{idx},
+ radius{r1, r2},
+ length{len}
{}
- value_type length() const {
- return norm(p1_-p2_);
+
+ size_type index;
+ std::pair<value_type, value_type> radius;
+ value_type length;
+};
+
+/// The simplest type of compartment iterator :
+/// - divide a segment into n compartments of equal length
+/// - assume that the radius varies linearly from one end of the segment
+/// to the other
+class compartment_iterator :
+ public std::iterator<std::forward_iterator_tag, compartment>
+{
+
+ public:
+
+ using base = std::iterator<std::forward_iterator_tag, compartment>;
+ using size_type = base::value_type::size_type;
+ using real_type = base::value_type::value_type;
+
+ compartment_iterator() = delete;
+
+ compartment_iterator(
+ size_type idx,
+ real_type len,
+ real_type rad,
+ real_type delta_rad
+ )
+ : index_(idx),
+ radius_(rad),
+ delta_radius_(delta_rad),
+ length_(len)
+ { }
+
+ compartment_iterator(compartment_iterator const& other) = default;
+
+ compartment_iterator& operator++()
+ {
+ index_++;
+ radius_ += delta_radius_;
+ return *this;
}
- value_type area() const {
- return volume_frustrum(length(), radius1_, radius2_);
+ compartment_iterator operator++(int)
+ {
+ compartment_iterator ret(*this);
+ operator++();
+ return ret;
}
- value_type volume() const {
- return volume_frustrum(length(), radius1_, radius2_);
+ compartment operator*() const
+ {
+ return
+ compartment(
+ index_, length_, radius_, radius_ + delta_radius_
+ );
}
- constexpr point_type midpoint() const {
- return 0.5*(p1_+p2_);
+ bool operator==(compartment_iterator const& other) const
+ {
+ return other.index_ == index_;
}
- constexpr point_type midradius() const {
- return 0.5*(radius1_+radius2_);
+ bool operator!=(compartment_iterator const& other) const
+ {
+ return !operator==(other);
}
private :
- value_type area_frustrum(
- value_type L,
- value_type r1,
- value_type r2
- ) const
+ size_type index_;
+ real_type radius_;
+ const real_type delta_radius_;
+ const real_type length_;
+};
+
+class compartment_range {
+ public:
+
+ using size_type = compartment_iterator::size_type;
+ using real_type = compartment_iterator::real_type;
+
+ compartment_range(
+ size_type num_compartments,
+ real_type radius_L,
+ real_type radius_R,
+ real_type length
+ )
+ : num_compartments_(num_compartments),
+ radius_L_(radius_L),
+ radius_R_(radius_R),
+ length_(length)
+ {}
+
+ compartment_iterator begin() const
{
- auto dr = r1 - r2;
- return pi<double>() * (r1+r2) * std::sqrt(L*L + dr*dr);
+ return {0, compartment_length(), radius_L_, delta_radius()};
}
- value_type volume_frustrum(
- value_type L,
- value_type r1,
- value_type r2
- ) const
+ compartment_iterator cbegin() const
{
- auto meanr = (r1+r2) / 2.;
- return pi<double>() * meanr * meanr * L;
+ return begin();
}
- point_type p1_;
- point_type p2_;
- value_type radius1_;
- value_type radius2_;
-};
+ /// With 0-based indexing compartment number "num_compartments_" is
+ /// one past the end
+ compartment_iterator end() const
+ {
+ return {num_compartments_, 0, 0, 0};
+ }
-class segment {
+ compartment_iterator cend() const
+ {
+ return end();
+ }
- private :
+ real_type delta_radius() const
+ {
+ return (radius_R_ - radius_L_) / num_compartments_;
+ }
- double length_;
- double radius_start_;
- double radius_end_;
- std::vector<compartment> compartments_;
- segmentKind kind_;
-};
+ real_type compartment_length() const
+ {
+ return length_ / num_compartments_;
+ }
-class abstract_cell {
- abstract_cell() = default;
+ private:
- private :
+ size_type num_compartments_;
+ real_type radius_L_;
+ real_type radius_R_;
+ real_type length_;
};
-} // namespace nestmc
+} // namespace mc
+} // namespace nest
diff --git a/src/fvm.hpp b/src/fvm.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..8ddb2938f45b3cf7740becfcc94356b938a99ca9
--- /dev/null
+++ b/src/fvm.hpp
@@ -0,0 +1,95 @@
+#pragma once
+
+#include <algorithm>
+
+#include "cell.hpp"
+#include "segment.hpp"
+
+#include "math.hpp"
+#include "util.hpp"
+#include "algorithms.hpp"
+
+#include "../vector/include/Vector.hpp"
+
+namespace nest {
+namespace mc {
+namespace fvm {
+
+template <typename T, typename I>
+class fvm_cell {
+ public :
+
+ /// the real number type
+ using value_type = T;
+ /// the integral index type
+ using size_type = I;
+
+ /// the container used for indexes
+ using index_type = memory::HostVector<size_type>;
+ /// view into index container
+ using index_view = typename index_type::view_type;
+
+ fvm_cell(nest::mc::cell const& cell);
+};
+
+////////////////////////////////////////////////////////////
+// Implementation
+////////////////////////////////////////////////////////////
+
+template <typename T, typename I>
+fvm_cell<T, I>::fvm_cell(nest::mc::cell const& cell)
+{
+ auto const& parent_index = cell.parent_index();
+ auto const& segment_index = cell.segment_index();
+ auto num_fv = segment_index.back();
+
+ // storage for the volume and surface area of the finite volumes
+ std::vector<value_type> fv_areas(num_fv);
+ std::vector<value_type> fv_volumes(num_fv);
+
+ auto seg_idx = 0;
+ for(auto const& s : cell.segments()) {
+ if(auto soma = s->as_soma()) {
+ // make the assumption that the soma is at 0
+ if(seg_idx!=0) {
+ throw std::domain_error(
+ "FVM lowering encountered soma with non-zero index"
+ );
+ }
+ fv_volumes[0] += math::volume_sphere(soma->radius());
+ fv_areas[0] += math::area_sphere (soma->radius());
+ }
+ else if(auto cable = s->as_cable()) {
+ using util::left;
+ using util::right;
+
+ for(auto c : cable->compartments()) {
+ auto rhs = segment_index[seg_idx] + c.index;
+ auto lhs = parent_index[rhs];
+
+ auto rad_C = math::mean(left(c.radius), right(c.radius));
+ auto len = c.length/2;
+
+ fv_volumes[lhs] += math::volume_frustrum(len, left(c.radius), rad_C);
+ fv_areas[lhs] += math::area_frustrum (len, left(c.radius), rad_C);
+
+ fv_volumes[rhs] += math::volume_frustrum(len, right(c.radius), rad_C);
+ fv_areas[rhs] += math::area_frustrum (len, right(c.radius), rad_C);
+ }
+ }
+ else {
+ throw std::domain_error(
+ "FVM lowering encountered unsuported segment type"
+ );
+ }
+ ++seg_idx;
+ }
+
+ //std::cout << "volumes " << fv_volumes << " : " << sum(fv_volumes) << "\n";
+ std::cout << "areas " << sum(fv_areas) << ", " << cell.volume() << "\n";
+ std::cout << "volumes " << sum(fv_volumes) << ", " << cell.area() << "\n";
+}
+
+} // namespace fvm
+} // namespace mc
+} // namespace nest
diff --git a/src/helpers.hpp b/src/helpers.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..cae1e30a5b7b43835cced5292d7f06e11b1c5d7f
--- /dev/null
+++ b/src/helpers.hpp
@@ -0,0 +1,15 @@
+#include <algorithm>
+
+namespace nestmc {
+namespace range{
+
+ template <typename C>
+ typename C::value_type
+ sum(C const& c)
+ {
+ using value_type = typename C::value_type;
+ return std::accumulate(c.begin(), c.end(), value_type{0});
+ }
+
+}
+}
diff --git a/src/math.hpp b/src/math.hpp
index c32c029eb5f03c9f1978b83caa836a8623bb9dc6..af0c402c40b603c8397fbeeb631695e68fa8a58b 100644
--- a/src/math.hpp
+++ b/src/math.hpp
@@ -2,7 +2,8 @@
#include <cmath>
-namespace nestmc {
+namespace nest {
+namespace mc {
namespace math {
template <typename T>
T constexpr pi()
@@ -61,6 +62,7 @@ namespace math {
return T(4) * pi<T>() * square(r);
}
-} // namespace nestmc::math
-} // namespace nestmc
+} // namespace math
+} // namespace mc
+} // namespace nest
diff --git a/src/parameter_list.cpp b/src/parameter_list.cpp
index 4394e800b32137501f3fc48b57ee573898926cad..a6d401b588c6eb8f543dd5f005201410cbaab827 100644
--- a/src/parameter_list.cpp
+++ b/src/parameter_list.cpp
@@ -3,7 +3,8 @@
#include "parameter_list.hpp"
-namespace nestmc {
+namespace nest {
+namespace mc {
bool parameter_list::add_parameter(parameter p)
{
@@ -75,25 +76,28 @@ namespace nestmc {
[&n](parameter const& p) {return p.name == n;}
);
}
-
-
- std::ostream& operator<<(std::ostream& o, parameter const& p)
- {
- return o
- << "parameter("
- << "name " << p.name
- << " : value " << p.value
- << " : range " << p.range
- << ")";
- }
-
- std::ostream& operator<<(std::ostream& o, parameter_list const& l)
- {
- o << "parameters \"" << l.name() << "\" :\n";
- for(auto const& p : l.parameters()) {
- o << " " << p << "\n";
- }
- return o;
+} // namespace mc
+} // namespace nest
+
+/*
+static std::ostream&
+operator<<(std::ostream& o, nest::mc::parameter const& p)
+{
+ return o
+ << "parameter("
+ << "name " << p.name
+ << " : value " << p.value
+ << " : range " << p.range
+ << ")";
+}
+
+static std::ostream&
+operator<<(std::ostream& o, nest::mc::parameter_list const& l)
+{
+ o << "parameters \"" << l.name() << "\" :\n";
+ for(nest::mc::parameter const& p : l.parameters()) {
+ o << " " << p << "\n";
}
-
-} // namespace nestmc
+ return o;
+}
+*/
diff --git a/src/parameter_list.hpp b/src/parameter_list.hpp
index ef790a640d40264585c1dcdb25c9d917a72624bb..114a18d07198aaa2b8230206e517ff4817bda445 100644
--- a/src/parameter_list.hpp
+++ b/src/parameter_list.hpp
@@ -6,7 +6,8 @@
#include <string>
#include <vector>
-namespace nestmc {
+namespace nest {
+namespace mc {
template <typename T>
struct value_range {
@@ -159,5 +160,6 @@ namespace nestmc {
}
};
-} // namespace nestmc
+} // namespace mc
+} // namespace nest
diff --git a/src/point.hpp b/src/point.hpp
index 00ded16ca2b70e3179cc18bb26a004344972c02b..cf74f3238d705e4e28c9430bb42cc85aa2b1ee14 100644
--- a/src/point.hpp
+++ b/src/point.hpp
@@ -3,6 +3,9 @@
#include <limits>
#include <ostream>
+namespace nest {
+namespace mc {
+
template <typename T>
struct point {
using value_type = T;
@@ -21,51 +24,61 @@ struct point {
z(std::numeric_limits<T>::quiet_NaN())
{}
- constexpr bool is_set() const {
+ constexpr bool is_set() const
+ {
return (x==x && y==y && z==z);
}
};
template <typename T>
-constexpr point<T> operator+ (
+constexpr point<T>
+operator+ (
point<T> const& lhs,
- point<T> const& rhs)
-{
+ 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- (
+constexpr point<T>
+operator- (
point<T> const& lhs,
- point<T> const& rhs)
-{
+ 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* (
+constexpr point<T>
+operator* (
T alpha,
- point<T> const& p)
-{
+ 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)
+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(
+constexpr T
+dot(
point<T> const& lhs,
- point<T> const& rhs)
-{
+ point<T> const& rhs
+) {
return lhs.x*rhs.x + lhs.y*rhs.y + lhs.z*rhs.z;
}
+} // namespace mc
+} // namespace nest
+
template <typename T>
-std::ostream& operator << (std::ostream& o, point<T> const& p) {
+std::ostream& operator << (std::ostream& o, nest::mc::point<T> const& p)
+{
return o << "[" << p.x << ", " << p.y << ", " << p.z << "]";
}
diff --git a/src/segment.hpp b/src/segment.hpp
index ca7251b898ec8800ec71711777bcd44863e618f7..a73dd9219cba674f90fb03d4e11ea44e0ab29dfd 100644
--- a/src/segment.hpp
+++ b/src/segment.hpp
@@ -4,12 +4,15 @@
#include <vector>
+#include "compartment.hpp"
#include "math.hpp"
#include "parameter_list.hpp"
#include "point.hpp"
+#include "algorithms.hpp"
#include "util.hpp"
-namespace nestmc {
+namespace nest {
+namespace mc {
template <typename T,
typename valid = typename std::is_floating_point<T>::type>
@@ -55,6 +58,9 @@ class segment {
return kind_==segmentKind::axon;
}
+ virtual int num_compartments() const = 0;
+ virtual void set_compartments(int) = 0;
+
virtual value_type volume() const = 0;
virtual value_type area() const = 0;
@@ -140,6 +146,14 @@ class placeholder_segment : public segment
{
return true;
}
+
+ int num_compartments() const override
+ {
+ return 0;
+ }
+
+ virtual void set_compartments(int) override
+ { }
};
class soma_segment : public segment
@@ -191,6 +205,15 @@ class soma_segment : public segment
return this;
}
+ /// soma has one and one only compartments
+ int num_compartments() const override
+ {
+ return 1;
+ }
+
+ void set_compartments(int n) override
+ { }
+
private :
// store the center and radius of the soma
@@ -230,7 +253,7 @@ class cable_segment : public segment
value_type r2,
value_type len
)
- : cable_segment{k, {r1, r2}, std::vector<value_type>{len}}
+ : cable_segment{k, {r1, r2}, {len}}
{ }
// constructor that lets the user describe the cable as a
@@ -268,7 +291,6 @@ class cable_segment : public segment
sum += math::volume_frustrum(lengths_[i], radii_[i], radii_[i+1]);
}
return sum;
-
}
value_type area() const override
@@ -280,12 +302,19 @@ class cable_segment : public segment
return sum;
}
- bool has_locations() const {
+ value_type length() const
+ {
+ return algorithms::sum(lengths_);
+ }
+
+ bool has_locations() const
+ {
return locations_.size() > 0;
}
// the number sub-segments that define the cable segment
- int num_sub_segments() const {
+ int num_sub_segments() const
+ {
return radii_.size()-1;
}
@@ -299,6 +328,53 @@ class cable_segment : public segment
return this;
}
+ int num_compartments() const override
+ {
+ return num_compartments_;
+ }
+
+ void set_compartments(int n) override
+ {
+ if(n<1) {
+ throw std::out_of_range(
+ "number of compartments in a segment must be one or more"
+ );
+ }
+ num_compartments_ = n;
+ }
+
+ value_type radius(value_type loc) const
+ {
+ if(loc>=1.) return radii_.back();
+ if(loc<=0.) return radii_.front();
+
+ auto len = length();
+ value_type pos = loc*len;
+
+ // This could be cached using a partial sum.
+ // In fact a lot of this stuff can be cached if
+ // we find ourselves having to do it over and over again.
+ // The time to cache it might be when update_lengths() is called.
+ auto sum = value_type(0);
+ auto i=0;
+ for(i=0; i<num_sub_segments(); ++i) {
+ if(sum+lengths_[i]>pos) {
+ break;
+ }
+ sum += lengths_[i];
+ }
+
+ auto rel = (len - sum)/lengths_[i];
+
+ return rel*radii_[i] + (1.-rel)*radii_[i+1];
+ }
+
+ /// iterable range type for simple compartment representation
+ compartment_range compartments() const
+ {
+ return {num_compartments(), radii_.front(), radii_.back(), length()};
+ }
+
private :
void update_lengths()
@@ -311,17 +387,24 @@ class cable_segment : public segment
}
}
+ int num_compartments_ = 1;
std::vector<value_type> lengths_;
std::vector<value_type> radii_;
std::vector<point_type> locations_;
};
+/// Unique pointer wrapper for abstract segment base class
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)...));
+/// Helper for constructing segments in a segment_ptr unique pointer wrapper.
+/// Forwards the supplied arguments to construct a segment of type SegmentType.
+/// e.g. auto my_cable = make_segment<cable>(segmentKind::dendrite, ... );
+template <typename SegmentType, typename... Args>
+segment_ptr make_segment(Args&&... args)
+{
+ return segment_ptr(new SegmentType(std::forward<Args>(args)...));
}
-} // namespace nestmc
+} // namespace mc
+} // namespace nest
diff --git a/src/tree.hpp b/src/tree.hpp
index eca2fca6c7e8dd198974128784a6cc2b03b5db3b..27822daad87a1ed8480223e784bbe8ee6b0e8cf1 100644
--- a/src/tree.hpp
+++ b/src/tree.hpp
@@ -7,11 +7,15 @@
#include <cassert>
#include "vector/include/Vector.hpp"
+#include "algorithms.hpp"
#include "util.hpp"
-using range = memory::Range;
+namespace nest {
+namespace mc {
class tree {
+ using range = memory::Range;
+
public :
using int_type = int16_t;
@@ -49,14 +53,21 @@ class tree {
/// create the tree from a parent_index
template <typename I>
- std::vector<I>
- init_from_parent_index(std::vector<I> const& parent_index)
+ tree(std::vector<I> const& parent_index)
{
+ // validate the inputs
+ if(!algorithms::is_minimal_degree(parent_index)) {
+ throw std::domain_error(
+ "parent index used to build a tree did not satisfy minimal degree ordering"
+ );
+ }
+
// n = number of compartment in cell
auto n = parent_index.size();
- // On completion of this loop child_count[i] is the number of children of compartment i
- // compensate count for compartment 0, which has itself as its own parent
+ // On completion of this loop child_count[i] is the number of children
+ // of compartment i compensate count for compartment 0, which has itself
+ // as its own parent
index_type child_count(n, 0);
child_count[0] = -1;
for(auto i : parent_index) {
@@ -137,9 +148,6 @@ class tree {
child_index_[i+1] = child_index_[i];
}
child_index_[0] = 0;
-
- // return the branch index to the caller for later use
- return branch_index;
}
size_t num_children() const {
@@ -344,8 +352,46 @@ class tree {
// 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);
+ index_view children_ = data_(0, 0);
+ index_view child_index_= data_(0, 0);
+ index_view parents_ = data_(0, 0);
};
+template <typename C>
+std::vector<int> make_parent_index(tree const& t, C const& counts)
+{
+ using range = memory::Range;
+
+ if( !algorithms::is_positive(counts)
+ || counts.size() != t.num_nodes() )
+ {
+ throw std::domain_error(
+ "make_parent_index requires one non-zero count per segment"
+ );
+ }
+ auto index = algorithms::make_index(counts);
+ auto num_compartments = index.back();
+ std::vector<int> parent_index(num_compartments);
+ auto pos = 0;
+ for(auto i : range(0, t.num_nodes())) {
+ // get the parent of this segment
+ // taking care for the case where the root node has -1 as its parent
+ auto parent = t.parent(i);
+ parent = parent>=0 ? parent : 0;
+ parent_index[pos++] = index[parent];
+ while(pos<index[i]) {
+ parent_index[pos] = pos-1;
+ pos++;
+ }
+ }
+
+ // if one of these assertions is tripped, we have to improve
+ // the input validation above
+ assert(pos==num_compartments);
+ assert(algorithms::is_minimal_degree(parent_index));
+
+ return parent_index;
+}
+
+} // namespace mc
+} // namespace nest
diff --git a/src/util.hpp b/src/util.hpp
index 580981be6e10f4f676057423d1533da5def51a68..738c3d0f6302cdc467a28a467427c50b8a96a593 100644
--- a/src/util.hpp
+++ b/src/util.hpp
@@ -2,6 +2,7 @@
#include "vector/include/Vector.hpp"
+
/*
using memory::util::red;
using memory::util::yellow;
@@ -13,6 +14,7 @@ using memory::util::cyan;
#include <memory>
#include <ostream>
+#include <utility>
#include <vector>
template <typename T>
@@ -27,12 +29,32 @@ operator << (std::ostream &o, std::vector<T>const& v)
return o;
}
+namespace nest {
+namespace mc {
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) ...));
}
-}
+
+ /// helper for taking the first value in a std::pair
+ template <typename L, typename R>
+ L const& left(std::pair<L, R> const& p)
+ {
+ return p.first;
+ }
+
+ /// helper for taking the second value in a std::pair
+ template <typename L, typename R>
+ R const& right(std::pair<L, R> const& p)
+ {
+ return p.second;
+ }
+
+} // namespace util
+} // namespace mc
+} // namespace nest
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 53db3ff016a2765cecef53d807b20d09f5c8eaa1..15d3ec160e9f9a701564c802ee188ebd8845b48e 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -14,12 +14,14 @@ set(TEST_SOURCES
# unit tests
test_cell.cpp
+ test_compartments.cpp
test_point.cpp
+ test_algorithms.cpp
+ test_run.cpp
test_segment.cpp
test_stimulus.cpp
test_swcio.cpp
test_tree.cpp
- test_run.cpp
# unit test driver
main.cpp
diff --git a/tests/test_algorithms.cpp b/tests/test_algorithms.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..342c094ae50a31baf7083f1db89c60e056828a6d
--- /dev/null
+++ b/tests/test_algorithms.cpp
@@ -0,0 +1,169 @@
+#include "gtest.h"
+
+#include <vector>
+
+#include "../src/algorithms.hpp"
+#include "../src/util.hpp"
+
+TEST(algorithms, sum)
+{
+ // sum of 10 times 2 is 20
+ std::vector<int> v1(10, 2);
+ EXPECT_EQ(10*2, nest::mc::algorithms::sum(v1));
+
+ // make an array 1:20 and sum it up using formula for arithmetic sequence
+ std::vector<int> v2(20);
+ std::iota(v2.begin(), v2.end(), 1);
+ auto n = 20;
+ EXPECT_EQ((n+1)*n/2, nest::mc::algorithms::sum(v2));
+}
+
+TEST(algorithms, make_index)
+{
+ {
+ std::vector<int> v(10, 1);
+ auto index = nest::mc::algorithms::make_index(v);
+
+ EXPECT_EQ(index.size(), 11u);
+ EXPECT_EQ(index.back(), nest::mc::algorithms::sum(v));
+ }
+
+ {
+ std::vector<int> v(10);
+ std::iota(v.begin(), v.end(), 1);
+ auto index = nest::mc::algorithms::make_index(v);
+
+ EXPECT_EQ(index.size(), 11u);
+ EXPECT_EQ(index.back(), nest::mc::algorithms::sum(v));
+ }
+}
+
+TEST(algorithms, minimal_degree)
+{
+ {
+ std::vector<int> v = {0};
+ EXPECT_TRUE(nest::mc::algorithms::is_minimal_degree(v));
+ }
+
+ {
+ std::vector<int> v = {0, 1, 2, 3, 4};
+ EXPECT_TRUE(nest::mc::algorithms::is_minimal_degree(v));
+ }
+
+ {
+ std::vector<int> v = {0, 0, 1, 2, 3, 4};
+ EXPECT_TRUE(nest::mc::algorithms::is_minimal_degree(v));
+ }
+
+ {
+ std::vector<int> v = {0, 0, 1, 2, 0, 4};
+ EXPECT_TRUE(nest::mc::algorithms::is_minimal_degree(v));
+ }
+
+ {
+ std::vector<int> v = {0, 0, 1, 2, 0, 4, 5, 4};
+ EXPECT_TRUE(nest::mc::algorithms::is_minimal_degree(v));
+ }
+
+ {
+ std::vector<int> v = {1};
+ EXPECT_FALSE(nest::mc::algorithms::is_minimal_degree(v));
+ }
+
+ {
+ std::vector<int> v = {0, 2};
+ EXPECT_FALSE(nest::mc::algorithms::is_minimal_degree(v));
+ }
+}
+
+TEST(algorithms, is_strictly_monotonic_increasing)
+{
+ EXPECT_TRUE(
+ nest::mc::algorithms::is_strictly_monotonic_increasing(
+ std::vector<int>{0}
+ )
+ );
+ EXPECT_TRUE(
+ nest::mc::algorithms::is_strictly_monotonic_increasing(
+ std::vector<int>{0, 1, 2, 3}
+ )
+ );
+ EXPECT_TRUE(
+ nest::mc::algorithms::is_strictly_monotonic_increasing(
+ std::vector<int>{8, 20, 42, 89}
+ )
+ );
+ EXPECT_FALSE(
+ nest::mc::algorithms::is_strictly_monotonic_increasing(
+ std::vector<int>{0, 0}
+ )
+ );
+ EXPECT_FALSE(
+ nest::mc::algorithms::is_strictly_monotonic_increasing(
+ std::vector<int>{8, 20, 20, 89}
+ )
+ );
+ EXPECT_FALSE(
+ nest::mc::algorithms::is_strictly_monotonic_increasing(
+ std::vector<int>{3, 2, 1, 0}
+ )
+ );
+}
+
+TEST(algorithms, is_strictly_monotonic_decreasing)
+{
+ EXPECT_TRUE(
+ nest::mc::algorithms::is_strictly_monotonic_decreasing(
+ std::vector<int>{0}
+ )
+ );
+ EXPECT_TRUE(
+ nest::mc::algorithms::is_strictly_monotonic_decreasing(
+ std::vector<int>{3, 2, 1, 0}
+ )
+ );
+ EXPECT_FALSE(
+ nest::mc::algorithms::is_strictly_monotonic_decreasing(
+ std::vector<int>{0, 1, 2, 3}
+ )
+ );
+ EXPECT_FALSE(
+ nest::mc::algorithms::is_strictly_monotonic_decreasing(
+ std::vector<int>{8, 20, 42, 89}
+ )
+ );
+ EXPECT_FALSE(
+ nest::mc::algorithms::is_strictly_monotonic_decreasing(
+ std::vector<int>{0, 0}
+ )
+ );
+ EXPECT_FALSE(
+ nest::mc::algorithms::is_strictly_monotonic_decreasing(
+ std::vector<int>{8, 20, 20, 89}
+ )
+ );
+}
+
+TEST(algorithms, is_positive)
+{
+ EXPECT_TRUE(
+ nest::mc::algorithms::is_positive(
+ std::vector<int>{}
+ )
+ );
+ EXPECT_TRUE(
+ nest::mc::algorithms::is_positive(
+ std::vector<int>{3, 2, 1}
+ )
+ );
+ EXPECT_FALSE(
+ nest::mc::algorithms::is_positive(
+ std::vector<int>{3, 2, 1, 0}
+ )
+ );
+ EXPECT_FALSE(
+ nest::mc::algorithms::is_positive(
+ std::vector<int>{-1}
+ )
+ );
+}
diff --git a/tests/test_cell.cpp b/tests/test_cell.cpp
index 68748512f594702d84a68d89c3eb72ec8010af25..09342b86e8d95a3185af68a903a7f59abdf77f2a 100644
--- a/tests/test_cell.cpp
+++ b/tests/test_cell.cpp
@@ -7,7 +7,7 @@ TEST(cell_type, soma)
// test that insertion of a soma works
// define with no centre point
{
- nestmc::cell c;
+ nest::mc::cell c;
auto soma_radius = 2.1;
EXPECT_EQ(c.has_soma(), false);
@@ -22,7 +22,7 @@ TEST(cell_type, soma)
// test that insertion of a soma works
// define with centre point @ (0,0,1)
{
- nestmc::cell c;
+ nest::mc::cell c;
auto soma_radius = 3.2;
EXPECT_EQ(c.has_soma(), false);
@@ -40,7 +40,7 @@ TEST(cell_type, soma)
TEST(cell_type, add_segment)
{
- using namespace nestmc;
+ using namespace nest::mc;
// add a pre-defined segment
{
cell c;
@@ -103,7 +103,7 @@ TEST(cell_type, add_segment)
TEST(cell_type, multiple_cables)
{
- using namespace nestmc;
+ using namespace nest::mc;
// generate a cylindrical cable segment of length 1/pi and radius 1
// volume = 1
@@ -145,7 +145,7 @@ TEST(cell_type, multiple_cables)
EXPECT_EQ(c.area(), 8. + math::area_sphere(soma_radius));
// construct the graph
- auto const& con = c.graph();
+ auto const& con = c.tree();
EXPECT_EQ(con.num_segments(), 5u);
EXPECT_EQ(con.parent(0), -1);
diff --git a/tests/test_compartments.cpp b/tests/test_compartments.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..4a40f1f8bf69be50d7d17413d117dc4fc0b979c6
--- /dev/null
+++ b/tests/test_compartments.cpp
@@ -0,0 +1,131 @@
+#include "gtest.h"
+
+#include "../src/compartment.hpp"
+#include "../src/util.hpp"
+
+using nest::mc::util::left;
+using nest::mc::util::right;
+
+// not much to test here: just test that values passed into the constructor
+// are correctly stored in members
+TEST(compartments, compartment)
+{
+
+ {
+ nest::mc::compartment c(100, 1.2, 2.1, 2.2);
+ EXPECT_EQ(c.index, 100);
+ EXPECT_EQ(c.length, 1.2);
+ EXPECT_EQ(left(c.radius), 2.1);
+ EXPECT_EQ(right(c.radius), 2.2);
+
+ auto c2 = c;
+ EXPECT_EQ(c2.index, 100);
+ EXPECT_EQ(c2.length, 1.2);
+ EXPECT_EQ(left(c2.radius), 2.1);
+ EXPECT_EQ(right(c2.radius), 2.2);
+ }
+
+ {
+ nest::mc::compartment c{100, 1, 2, 3};
+ EXPECT_EQ(c.index, 100);
+ EXPECT_EQ(c.length, 1.);
+ EXPECT_EQ(left(c.radius), 2.);
+ EXPECT_EQ(right(c.radius), 3.);
+ }
+}
+
+TEST(compartments, compartment_iterator)
+{
+ // iterator with arguments
+ // idx = 0
+ // len = 2.5
+ // rad = 1
+ // delta_rad = 2
+ nest::mc::compartment_iterator it(0, 2.5, 1, 2);
+
+ // so each time the iterator is incremented
+ // idx is incremented by 1
+ // len is unchanged
+ // rad is incremented by 2
+
+ // check the prefix increment
+ ++it;
+ {
+ auto c = *it;
+ EXPECT_EQ(c.index, 1);
+ EXPECT_EQ(left(c.radius), 3.0);
+ EXPECT_EQ(right(c.radius), 5.0);
+ EXPECT_EQ(c.length, 2.5);
+ }
+
+ // check postfix increment
+
+ // returned iterator should be unchanged
+ {
+ auto c = *(it++);
+ EXPECT_EQ(c.index, 1);
+ EXPECT_EQ(left(c.radius), 3.0);
+ EXPECT_EQ(right(c.radius), 5.0);
+ EXPECT_EQ(c.length, 2.5);
+ }
+ // while the iterator itself was updated
+ {
+ auto c = *it;
+ EXPECT_EQ(c.index, 2);
+ EXPECT_EQ(left(c.radius), 5.0);
+ EXPECT_EQ(right(c.radius), 7.0);
+ EXPECT_EQ(c.length, 2.5);
+ }
+
+ // check that it can be copied and compared
+ {
+ // copy iterator
+ auto it2 = it;
+ auto c = *it2;
+ EXPECT_EQ(c.index, 2);
+ EXPECT_EQ(left(c.radius), 5.0);
+ EXPECT_EQ(right(c.radius), 7.0);
+ EXPECT_EQ(c.length, 2.5);
+
+ // comparison
+ EXPECT_EQ(it2, it);
+ it2++;
+ EXPECT_NE(it2, it);
+
+ // check the copy has updated correctly when incremented
+ c= *it2;
+ EXPECT_EQ(c.index, 3);
+ EXPECT_EQ(left(c.radius), 7.0);
+ EXPECT_EQ(right(c.radius), 9.0);
+ EXPECT_EQ(c.length, 2.5);
+ }
+}
+
+TEST(compartments, compartment_range)
+{
+ {
+ nest::mc::compartment_range rng(10, 1.0, 2.0, 10.);
+
+ EXPECT_EQ((*rng.begin()).index, 0);
+ EXPECT_EQ((*rng.end()).index, 10);
+ EXPECT_NE(rng.begin(), rng.end());
+
+ auto count = 0;
+ for(auto c : rng) {
+ EXPECT_EQ(c.index, count);
+ auto er = 1.0 + double(count)/10.;
+ EXPECT_TRUE(std::fabs(left(c.radius)-er)<1e-15);
+ EXPECT_TRUE(std::fabs(right(c.radius)-(er+0.1))<1e-15);
+ EXPECT_EQ(c.length, 1.0);
+ ++count;
+ }
+ EXPECT_EQ(count, 10);
+ }
+
+ // test case of zero length range
+ {
+ nest::mc::compartment_range rng(0, 1.0, 1.0, 0.);
+
+ EXPECT_EQ(rng.begin(), rng.end());
+ }
+}
diff --git a/tests/test_fvm.cpp b/tests/test_fvm.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/tests/test_point.cpp b/tests/test_point.cpp
index 57a8567b94e4a690ca40eb246987120c5721a8d3..9666a1387f895d408778f40afc6b07db5363bec6 100644
--- a/tests/test_point.cpp
+++ b/tests/test_point.cpp
@@ -5,6 +5,8 @@
#include "../src/point.hpp"
+using namespace nest::mc;
+
TEST(point, construction)
{
{
diff --git a/tests/test_run.cpp b/tests/test_run.cpp
index 618d9e79560ea0764d7c22ec844749dc8ec51ffa..74b505c6e71c66bc32ffca346caae1ec57262b3c 100644
--- a/tests/test_run.cpp
+++ b/tests/test_run.cpp
@@ -1,19 +1,20 @@
#include "gtest.h"
#include "../src/cell.hpp"
+#include "../src/fvm.hpp"
TEST(run, init)
{
- using namespace nestmc;
+ using namespace nest::mc;
- nestmc::cell cell;
+ nest::mc::cell cell;
cell.add_soma(12.6157/2.0);
//auto& props = cell.soma()->properties;
cell.add_cable(0, segmentKind::dendrite, 0.5, 0.5, 200);
- EXPECT_EQ(cell.graph().num_segments(), 2u);
+ EXPECT_EQ(cell.tree().num_segments(), 2u);
/*
for(auto &s : cell.segments()) {
@@ -41,19 +42,23 @@ TEST(run, init)
EXPECT_EQ(cell.soma()->mechanism("hh").get("gl").value, 0.0003);
EXPECT_EQ(cell.soma()->mechanism("hh").get("el").value, -54.3);
+
+ cell.segment(1)->set_compartments(2);
+
+ //using fvm_cell = fvm::fvm_cell<double, int>;
+ //fvm_cell fvcell(cell);
// print out the parameters if you want...
//std::cout << soma_hh << "\n";
-
}
// test out the parameter infrastructure
TEST(run, parameters)
{
- nestmc::parameter_list list("test");
+ nest::mc::parameter_list list("test");
EXPECT_EQ(list.name(), "test");
EXPECT_EQ(list.num_parameters(), 0);
- nestmc::parameter p("a", 0.12, {0, 10});
+ nest::mc::parameter p("a", 0.12, {0, 10});
// add_parameter() returns a bool that indicates whether
// it was able to successfull add the parameter
diff --git a/tests/test_segment.cpp b/tests/test_segment.cpp
index 81f70aebe63cf0e97d3122f400c0d80e0ab6f593..9c3df625db9ef87110474d35a5ce88f3f4e07f09 100644
--- a/tests/test_segment.cpp
+++ b/tests/test_segment.cpp
@@ -7,8 +7,8 @@
TEST(segments, soma)
{
- using namespace nestmc;
- using nestmc::math::pi;
+ using namespace nest::mc;
+ using nest::mc::math::pi;
{
auto s = make_segment<soma_segment>(1.0);
@@ -29,8 +29,8 @@ TEST(segments, soma)
TEST(segments, cable)
{
- using namespace nestmc;
- using nestmc::math::pi;
+ using namespace nest::mc;
+ using nest::mc::math::pi;
// take advantage of fact that a cable segment with constant radius 1 and
// length 1 has volume=1. and area=2
@@ -63,8 +63,8 @@ TEST(segments, cable)
TEST(segments, cable_positions)
{
- using namespace nestmc;
- using nestmc::math::pi;
+ using namespace nest::mc;
+ using nest::mc::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)
diff --git a/tests/test_tree.cpp b/tests/test_tree.cpp
index 58435e15240344665921799dbf07ff8eca211f5c..36b273a064a9d3fe15102a745f732b05f7f71f99 100644
--- a/tests/test_tree.cpp
+++ b/tests/test_tree.cpp
@@ -11,6 +11,8 @@
using json = nlohmann::json;
using range = memory::Range;
+using namespace nest::mc;
+
TEST(cell_tree, from_parent_index) {
// tree with single branch corresponding to the root node
// this is equivalent to a single compartment model
@@ -137,8 +139,7 @@ TEST(tree, change_root) {
// |
// 2
std::vector<int> parent_index = {0,0,0};
- tree t;
- t.init_from_parent_index(parent_index);
+ tree t(parent_index);
t.change_root(1);
EXPECT_EQ(t.num_nodes(), 3u);
@@ -156,8 +157,7 @@ TEST(tree, change_root) {
// / \ |
// 3 4 4
std::vector<int> parent_index = {0,0,0,1,1};
- tree t;
- t.init_from_parent_index(parent_index);
+ tree t(parent_index);
t.change_root(1u);
EXPECT_EQ(t.num_nodes(), 5u);
@@ -181,8 +181,7 @@ TEST(tree, change_root) {
// / \.
// 5 6
std::vector<int> parent_index = {0,0,0,1,1,4,4};
- tree t;
- t.init_from_parent_index(parent_index);
+ tree t(parent_index);
t.change_root(1);
@@ -239,7 +238,8 @@ TEST(cell_tree, balance) {
// this test doesn't test anything yet... it just loads each cell in turn
// from a json file and creates a .dot file for it
-TEST(cell_tree, json_load) {
+TEST(cell_tree, json_load)
+{
json cell_data;
std::ifstream("../data/cells_small.json") >> cell_data;
@@ -251,3 +251,13 @@ TEST(cell_tree, json_load) {
}
}
+TEST(tree, make_parent_index)
+{
+ {
+ std::vector<int> parent_index = {0,0,1,2,3,4,3,6,0,8};
+ std::vector<int> counts = {1,3,2,2,2};
+ tree t(parent_index);
+ auto new_pid = make_parent_index(t, counts);
+ EXPECT_EQ(parent_index.size(), new_pid.size());
+ }
+}