Merge branch 'cell-construction' of github.com:eth-cscs/cell_algorithms into cell-construction

src/cell.cpp

@@ -9,7 +9,6 @@ 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
@@ -36,7 +35,6 @@ 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)
@@ -56,7 +54,6 @@ 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)
@@ -87,21 +84,14 @@ 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();
     }
@@ -145,38 +135,27 @@ std::vector<int> cell::compartment_counts() const
         comp_count.push_back(s->num_compartments());
     }
     return comp_count;
-
 }
 
-void cell::construct() const
+int cell::num_compartments() const
 {
-    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);
+    auto n = 0;
+    for(auto &s : segments_) {
+        n += s->num_compartments();
     }
-    stale_ = false;
+    return n;
 }
 
-std::vector<int> const& cell::parent_index() const
+compartment_model cell::model() const
 {
-    construct();
-    return parent_index_;
-}
+    compartment_model m;
 
-std::vector<int> const& cell::segment_index() const
-{
-    construct();
-    return segment_index_;
-}
+    m.tree = cell_tree(parents_);
+    auto counts = compartment_counts();
+    m.parent_index = make_parent_index(m.tree.graph(), counts);
+    m.segment_index = algorithms::make_index(counts);
 
-cell_tree const& cell::tree() const
-{
-    construct();
-    return tree_;
+    return m;
 }
 
 std::vector<int> const& cell::segment_parents() const

src/cell.hpp

@@ -11,7 +11,15 @@
 namespace nest {
 namespace mc {
 
-// high-level abstract representation of a cell and its segments
+/// wrapper around compartment layout information derived from a high level cell
+/// description
+struct compartment_model {
+    cell_tree tree;
+    std::vector<int> parent_index;
+    std::vector<int> segment_index;
+};
+
+/// high-level abstract representation of a cell and its segments
 class cell {
     public:
 
@@ -61,10 +69,10 @@ class cell {
     /// the surface area of the cell
     value_type area() const;
 
-    std::vector<segment_ptr> const& segments() const;
+    /// the total number of compartments over all segments
+    int num_compartments() const;
 
-    /// the connectivity graph for the cell segments
-    cell_tree const& tree() const;
+    std::vector<segment_ptr> const& segments() const;
 
     /// return reference to array that enumerates the index of the parent of
     /// each segment
@@ -73,46 +81,14 @@ class cell {
     /// 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;
+    compartment_model model() 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
@@ -125,7 +101,6 @@ void cell::add_cable(cell::index_type parent, Args ...args)
             "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);
 }

src/fvm.hpp

@@ -97,18 +97,19 @@ class fvm_cell {
 
 template <typename T, typename I>
 fvm_cell<T, I>::fvm_cell(nest::mc::cell const& cell)
-:   matrix_        {cell.parent_index()}
-,   cv_areas_      {size(), T(0)}
-,   face_alpha_    {size(), T(0)}
-,   cv_capacitance_{size(), T(0)}
-,   current_       {size(), T(0)}
-,   voltage_       {size(), T(0)}
+:   cv_areas_      {cell.num_compartments(), T(0)}
+,   face_alpha_    {cell.num_compartments(), T(0)}
+,   cv_capacitance_{cell.num_compartments(), T(0)}
+,   current_       {cell.num_compartments(), T(0)}
+,   voltage_       {cell.num_compartments(), T(0)}
 {
     using util::left;
     using util::right;
 
+    const auto graph = cell.model();
+    matrix_ = matrix_type(graph.parent_index);
     auto parent_index = matrix_.p();
-    auto const& segment_index = cell.segment_index();
+    auto const& segment_index = graph.segment_index;
 
     auto seg_idx = 0;
     for(auto const& s : cell.segments()) {

src/matrix.hpp

@@ -24,6 +24,8 @@ class matrix {
     using index_type  = memory::HostVector<size_type>;
     using index_view  = typename index_type::view_type;
 
+    matrix() = default;
+
     /// construct matrix for one or more cells, with combined parent index
     /// and a cell index
     template <

tests/test_cell.cpp

@@ -145,7 +145,8 @@ TEST(cell_type, multiple_cables)
         EXPECT_EQ(c.area(), 8. + math::area_sphere(soma_radius));
 
         // construct the graph
-        auto const& con = c.tree();
+        const auto model = c.model();
+        auto const& con = model.tree;
 
         EXPECT_EQ(con.num_segments(), 5u);
         EXPECT_EQ(con.parent(0), -1);

tests/test_run.cpp

@@ -15,7 +15,7 @@ TEST(run, cable)
     cell.add_cable(0, segmentKind::dendrite, 1e-4, 1e-4, 4e-1);
 
     std::cout << cell.segment(1)->area() << " is the area\n";
-    EXPECT_EQ(cell.tree().num_segments(), 2u);
+    EXPECT_EQ(cell.model().tree.num_segments(), 2u);
 
     cell.soma()->add_mechanism(hh_parameters());
 
@@ -49,8 +49,6 @@ TEST(run, cable)
     J.rhs()[0] = 10.;
 
     J.solve();
-
-    //std::cout << "x" << J.rhs() << "\n";
 }
 
 TEST(run, init)
@@ -64,7 +62,8 @@ TEST(run, init)
 
     cell.add_cable(0, segmentKind::dendrite, 0.5, 0.5, 200);
 
-    EXPECT_EQ(cell.tree().num_segments(), 2u);
+    const auto m = cell.model();
+    EXPECT_EQ(m.tree.num_segments(), 2u);
 
     // in this context (i.e. attached to a segment on a high-level cell)
     // a mechanism is essentially a set of parameters