diff --git a/src/cell.cpp b/src/cell.cpp index a124eb384837e35215823bc8812f399905f56b68..ab01007a206ea241befea21fbed83c531298a36d 100644 --- a/src/cell.cpp +++ b/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 diff --git a/src/cell.hpp b/src/cell.hpp index 0a254e1ec006dec417b89ef36fa8f28631a283b7..68bc4c5a21f99237b18557f690226c0441e580c6 100644 --- a/src/cell.hpp +++ b/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); } diff --git a/src/fvm.hpp b/src/fvm.hpp index 9e84119f526c531f3e80c082ea22d8dc7edd997b..1cab7b19e0ded99986648f2a0759e716b1112b2e 100644 --- a/src/fvm.hpp +++ b/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()) { diff --git a/src/matrix.hpp b/src/matrix.hpp index 3e5f4442373c24821bc9f6b8499db90a37a87835..16936fa47c26d3c632d08d008e59c76c2e4b0db0 100644 --- a/src/matrix.hpp +++ b/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 < diff --git a/tests/test_cell.cpp b/tests/test_cell.cpp index 09342b86e8d95a3185af68a903a7f59abdf77f2a..88c8dcaced167d4030529adb8a348a116158c17f 100644 --- a/tests/test_cell.cpp +++ b/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); diff --git a/tests/test_run.cpp b/tests/test_run.cpp index 712fca4346e2a1780a9c1c3dca080381eeb8218d..c4e3b1025e82ca257295fdcf8980fd82240d1956 100644 --- a/tests/test_run.cpp +++ b/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