From 16c625a169b922c9256a846d33a78cee15e13dc2 Mon Sep 17 00:00:00 2001
From: Sam Yates <sam@quux.dropbear.id.au>
Date: Thu, 1 Sep 2016 03:12:45 +0200
Subject: [PATCH] First version of fvm_multicell that compiles.
Note: failed assertion at run time with both fvm_cell
and fvm_multicell.
---
miniapp/miniapp.cpp | 3 +-
src/compartment.hpp | 137 ++-----------------------------
src/fvm_multicell.hpp | 40 +++++----
src/mechanism_catalogue.hpp | 5 +-
src/parameter_list.cpp | 135 +++++++++++++++---------------
src/parameter_list.hpp | 1 +
src/segment.hpp | 4 +-
src/util/transform.hpp | 4 +-
tests/unit/test_compartments.cpp | 101 +----------------------
9 files changed, 110 insertions(+), 320 deletions(-)
diff --git a/miniapp/miniapp.cpp b/miniapp/miniapp.cpp
index 424098fc..f0983d05 100644
--- a/miniapp/miniapp.cpp
+++ b/miniapp/miniapp.cpp
@@ -28,7 +28,8 @@
using namespace nest::mc;
using global_policy = communication::global_policy;
-using lowered_cell = fvm::fvm_multicell<double, cell_local_size_type>;
+//using lowered_cell = fvm::fvm_multicell<double, cell_local_size_type>;
+using lowered_cell = fvm::fvm_cell<double, cell_local_size_type>;
using model_type = model<lowered_cell>;
using time_type = model_type::time_type;
using sample_trace_type = sample_trace<time_type, model_type::value_type>;
diff --git a/src/compartment.hpp b/src/compartment.hpp
index 326563ad..083155ff 100644
--- a/src/compartment.hpp
+++ b/src/compartment.hpp
@@ -37,134 +37,6 @@ struct compartment {
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;
- }
-
- compartment_iterator operator++(int)
- {
- compartment_iterator ret(*this);
- operator++();
- return ret;
- }
-
- compartment operator*() const
- {
- return
- compartment(
- index_, length_, radius_, radius_ + delta_radius_
- );
- }
-
- bool operator==(compartment_iterator const& other) const
- {
- return other.index_ == index_;
- }
-
- bool operator!=(compartment_iterator const& other) const
- {
- return !operator==(other);
- }
-
- private :
-
- 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
- {
- return {0, compartment_length(), radius_L_, delta_radius()};
- }
-
- compartment_iterator cbegin() const
- {
- return begin();
- }
-
- /// With 0-based indexing compartment number "num_compartments_" is
- /// one past the end
- compartment_iterator end() const
- {
- return {num_compartments_, 0, 0, 0};
- }
-
- compartment_iterator cend() const
- {
- return end();
- }
-
- real_type delta_radius() const
- {
- return (radius_R_ - radius_L_) / num_compartments_;
- }
-
- real_type compartment_length() const
- {
- return length_ / num_compartments_;
- }
-
- private:
-
- size_type num_compartments_;
- real_type radius_L_;
- real_type radius_R_;
- real_type length_;
-};
-
// (NB: auto type deduction and lambda in C++14 will simplify the following)
template <typename size_type, typename real_type>
@@ -187,11 +59,16 @@ private:
};
template <typename size_type, typename real_type>
-using compartment_iterator_bis =
+using compartment_iterator =
util::transform_iterator<util::counter<size_type>, compartment_maker<size_type, real_type>>;
template <typename size_type, typename real_type>
-util::range<compartment_iterator_bis<size_type, real_type>> make_compartment_range(
+using compartment_range =
+ util::range<compartment_iterator<size_type, real_type>>;
+
+
+template <typename size_type, typename real_type>
+compartment_range<size_type, real_type> make_compartment_range(
size_type num_compartments,
real_type radius_L,
real_type radius_R,
diff --git a/src/fvm_multicell.hpp b/src/fvm_multicell.hpp
index 8239c7e5..9100a8fb 100644
--- a/src/fvm_multicell.hpp
+++ b/src/fvm_multicell.hpp
@@ -20,6 +20,8 @@
#include <stimulus.hpp>
#include <util.hpp>
#include <util/meta.hpp>
+#include <util/partition.hpp>
+#include <util/span.hpp>
#include <vector/include/Vector.hpp>
@@ -202,7 +204,7 @@ private:
/// the ion species
std::map<mechanisms::ionKind, ion_type> ions_;
- std::vector<std::pair<uint32_t, i_clamp>> stimulii_;
+ std::vector<std::pair<uint32_t, i_clamp>> stimuli_;
std::vector<std::pair<const vector_type fvm_multicell::*, uint32_t>> probes_;
@@ -211,7 +213,7 @@ private:
// perform area and capacitance calculation on initialization
void compute_cv_area_unnormalized_capacitance(
- std::pair<size_type, size_type> comps, const segment& seg, index_type &parent);
+ std::pair<size_type, size_type> comps, const segment* seg, index_type &parent);
};
////////////////////////////////////////////////////////////////////////////////
@@ -221,7 +223,7 @@ private:
template <typename T, typename I>
void fvm_multicell<T, I>::compute_cv_area_unnormalized_capacitance(
std::pair<size_type, size_type> comps,
- const segment& seg,
+ const segment* seg,
index_type &parent)
{
using util::left;
@@ -230,9 +232,11 @@ void fvm_multicell<T, I>::compute_cv_area_unnormalized_capacitance(
// precondition: group_parent_index[j] holds the correct value for
// j in [base_comp, base_comp+segment.num_compartments()].
- if (auto soma = seg.as_soma()) {
+ auto ncomp = comps.second-comps.first;
+
+ if (auto soma = seg->as_soma()) {
// confirm assumption that there is one compartment in soma
- if (comps.size()!=1)
+ if (ncomp!=1) {
throw std::logic_error("soma allocated more than one compartment");
}
auto i = comps.first;
@@ -241,7 +245,7 @@ void fvm_multicell<T, I>::compute_cv_area_unnormalized_capacitance(
cv_areas_[i] += area;
cv_capacitance_[i] += area * soma->mechanism("membrane").get("c_m").value;
}
- else if (auto cable = s.as_cable()) {
+ else if (auto cable = seg->as_cable()) {
// loop over each compartment in the cable
// each compartment has the face between two CVs at its centre
// the centers of the CVs are the end points of the compartment
@@ -261,10 +265,10 @@ void fvm_multicell<T, I>::compute_cv_area_unnormalized_capacitance(
auto r_L = cable->mechanism("membrane").get("r_L").value;
const auto& compartments = cable->compartments();
- EXPECTS(util::size(compartments)==comps.second-comps.first);
+ EXPECTS(util::size(compartments)==ncomp);
for (auto i: util::make_span(comps)) {
- const auto c& = compartments[i-comps.first];
+ const auto& c = compartments[i-comps.first];
auto j = parent[i];
auto radius_center = math::mean(c.radius);
@@ -302,10 +306,10 @@ void fvm_multicell<T, I>::initialize(
// count total detectors, targets and probes for validation of handle container sizes
size_type detectors_count = 0u;
size_type targets_count = 0u;
- size_type probess_count = 0u;
+ size_type probes_count = 0u;
size_type detectors_size = util::size(detector_handles);
- size_type targets_size = util::size(targets_handles);
- size_type probes_size = util::size(probes_handles);
+ size_type targets_size = util::size(target_handles);
+ size_type probes_size = util::size(probe_handles);
auto ncell = util::size(cells);
auto cell_num_compartments =
@@ -323,12 +327,12 @@ void fvm_multicell<T, I>::initialize(
voltage_ = vector_type{ncomp, T{0}};
// create maps for mechanism initialization.
- std::map<std::string, std::pair<size_type, size_type>> mech_map;
+ std::map<std::string, std::vector<std::pair<size_type, size_type>>> mech_map;
std::vector<std::vector<cell_lid_type>> syn_mech_map;
std::map<std::string, std::size_t> syn_mech_indices;
// initialize vector used for matrix creation.
- index_vector group_parent_index{ncomp, 0};
+ index_type group_parent_index{ncomp, 0};
// create each cell:
auto target_hi = target_handles.begin();
@@ -346,7 +350,7 @@ void fvm_multicell<T, I>::initialize(
}
auto seg_num_compartments =
- transform_view(c.segments(), [](const segment& s) { return s.num_compartments(); });
+ transform_view(c.segments(), [](const segment_ptr& s) { return s->num_compartments(); });
auto nseg = seg_num_compartments.size();
std::vector<cell_lid_type> seg_comp_bounds;
@@ -360,7 +364,7 @@ void fvm_multicell<T, I>::initialize(
for (const auto& mech: seg->mechanisms()) {
if (mech.name()!="membrane") {
- mech_map[mech_name().push_back(seg_comps)];
+ mech_map[mech.name()].push_back(seg_comps);
}
}
}
@@ -391,7 +395,7 @@ void fvm_multicell<T, I>::initialize(
cv_capacitance_[k] /= cv_areas_[k];
}
- // add the stimulii
+ // add the stimuli
for (const auto& stim: c.stimuli()) {
auto idx = comps.first+find_compartment_index(stim.location, graph);
stimuli_.push_back({idx, stim.clamp});
@@ -578,8 +582,8 @@ void fvm_multicell<T, I>::advance(T dt) {
PL();
}
- // add current contributions from stimulii
- for (auto& stim : stimulii_) {
+ // add current contributions from stimuli
+ for (auto& stim : stimuli_) {
auto ie = stim.second.amplitude(t_); // [nA]
auto loc = stim.first;
diff --git a/src/mechanism_catalogue.hpp b/src/mechanism_catalogue.hpp
index c9d3e6de..aecea5ce 100644
--- a/src/mechanism_catalogue.hpp
+++ b/src/mechanism_catalogue.hpp
@@ -19,18 +19,19 @@ struct catalogue {
using view_type = typename mechanism<T, I>::view_type;
using index_view = typename mechanism<T, I>::index_view;
+ template <typename Indices>
static mechanism_ptr<T, I> make(
const std::string& name,
view_type vec_v,
view_type vec_i,
- index_view node_indices)
+ Indices node_indices)
{
auto entry = mech_map.find(name);
if (entry==mech_map.end()) {
throw std::out_of_range("no such mechanism");
}
- return entry->second(vec_v, vec_i, node_indices);
+ return entry->second(vec_v, vec_i, index_view{node_indices});
}
static bool has(const std::string& name) {
diff --git a/src/parameter_list.cpp b/src/parameter_list.cpp
index 86c641aa..2ca50817 100644
--- a/src/parameter_list.cpp
+++ b/src/parameter_list.cpp
@@ -6,88 +6,92 @@
namespace nest {
namespace mc {
- bool parameter_list::add_parameter(parameter p)
- {
- if(has_parameter(p.name)) {
- return false;
- }
+bool parameter_list::add_parameter(parameter p) {
+ if (has_parameter(p.name)) {
+ return false;
+ }
- parameters_.push_back(std::move(p));
+ parameters_.push_back(std::move(p));
- return true;
- }
+ return true;
+}
- bool parameter_list::has_parameter(std::string const& n) const
- {
- return find_by_name(n) != parameters_.end();
- }
+bool parameter_list::has_parameter(std::string const& n) const {
+ return find_by_name(n) != parameters_.end();
+}
- int parameter_list::num_parameters() const
- {
- return parameters_.size();
- }
+int parameter_list::num_parameters() const {
+ return parameters_.size();
+}
- // returns true if parameter was succesfully updated
- // returns false if parameter was not updated, i.e. if
- // - no parameter with name n exists
- // - value is not in the valid range
- bool parameter_list::set(std::string const& n, parameter_list::value_type v)
- {
- auto p = find_by_name(n);
- if(p!=parameters_.end()) {
- if(p->is_in_range(v)) {
- p->value = v;
- return true;
- }
+// returns true if parameter was succesfully updated
+// returns false if parameter was not updated, i.e. if
+// - no parameter with name n exists
+// - value is not in the valid range
+bool parameter_list::set(std::string const& n, parameter_list::value_type v) {
+ auto p = find_by_name(n);
+ if(p!=parameters_.end()) {
+ if(p->is_in_range(v)) {
+ p->value = v;
+ return true;
}
- return false;
}
+ return false;
+}
- parameter& parameter_list::get(std::string const& n)
- {
- auto it = find_by_name(n);
- if(it==parameters_.end()) {
- throw std::domain_error(
- "parameter list does not contain parameter"
- );
- }
- return *it;
+parameter& parameter_list::get(std::string const& n) {
+ auto it = find_by_name(n);
+ if (it==parameters_.end()) {
+ throw std::domain_error(
+ "parameter list does not contain parameter"
+ );
}
+ return *it;
+}
- std::string const& parameter_list::name() const {
- return mechanism_name_;
+const parameter& parameter_list::get(std::string const& n) const {
+ auto it = find_by_name(n);
+ if (it==parameters_.end()) {
+ throw std::domain_error(
+ "parameter list does not contain parameter"
+ );
}
+ return *it;
+}
- std::vector<parameter> const& parameter_list::parameters() const
- {
- return parameters_;
- }
+std::string const& parameter_list::name() const {
+ return mechanism_name_;
+}
- auto parameter_list::find_by_name(std::string const& n)
- -> decltype(parameters_.begin())
- {
- return
- std::find_if(
- parameters_.begin(), parameters_.end(),
- [&n](parameter const& p) {return p.name == n;}
- );
- }
+std::vector<parameter> const& parameter_list::parameters() const {
+ return parameters_;
+}
+
+auto parameter_list::find_by_name(std::string const& n)
+ -> decltype(parameters_.begin())
+{
+ return
+ std::find_if(
+ parameters_.begin(), parameters_.end(),
+ [&n](parameter const& p) {return p.name == n;}
+ );
+}
+
+auto parameter_list::find_by_name(std::string const& n) const
+ -> decltype(parameters_.begin())
+{
+ return
+ std::find_if(
+ parameters_.begin(), parameters_.end(),
+ [&n](parameter const& p) {return p.name == n;}
+ );
+}
- auto parameter_list::find_by_name(std::string const& n) const
- -> decltype(parameters_.begin())
- {
- return
- std::find_if(
- parameters_.begin(), parameters_.end(),
- [&n](parameter const& p) {return p.name == n;}
- );
- }
} // namespace mc
} // namespace nest
std::ostream&
-operator<<(std::ostream& o, nest::mc::parameter const& p)
-{
+operator<<(std::ostream& o, nest::mc::parameter const& p) {
return o
<< "parameter("
<< "name " << p.name
@@ -97,8 +101,7 @@ operator<<(std::ostream& o, nest::mc::parameter const& p)
}
std::ostream&
-operator<<(std::ostream& o, nest::mc::parameter_list const& l)
-{
+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";
diff --git a/src/parameter_list.hpp b/src/parameter_list.hpp
index 0c069887..eb3ee8fe 100644
--- a/src/parameter_list.hpp
+++ b/src/parameter_list.hpp
@@ -102,6 +102,7 @@ namespace mc {
// - value is not in the valid range
bool set(std::string const& n, value_type v);
parameter& get(std::string const& n);
+ const parameter& get(std::string const& n) const;
std::string const& name() const;
diff --git a/src/segment.hpp b/src/segment.hpp
index 248306a6..5a9782b9 100644
--- a/src/segment.hpp
+++ b/src/segment.hpp
@@ -418,9 +418,9 @@ public:
}
/// iterable range type for simple compartment representation
- compartment_range compartments() const
+ compartment_range<size_type, value_type> compartments() const
{
- return {num_compartments(), radii_.front(), radii_.back(), length()};
+ return make_compartment_range(num_compartments(), radii_.front(), radii_.back(), length());
}
private:
diff --git a/src/util/transform.hpp b/src/util/transform.hpp
index f7d9fed0..d20b0061 100644
--- a/src/util/transform.hpp
+++ b/src/util/transform.hpp
@@ -33,12 +33,12 @@ class transform_iterator: public iterator_adaptor<transform_iterator<I, F>, I> {
public:
using typename base::difference_type;
- using value_type = decltype(f_(*inner_));
+ using value_type = typename std::result_of<F (inner_value_type)>::type;
using pointer = const value_type*;
using reference = const value_type&;
template <typename J, typename G>
- transform_iterator(J&& c, G&& g): inner_{std::forward<J>(c)}, f_{std::forward<G>(g)} {}
+ transform_iterator(J&& c, G&& g): inner_(std::forward<J>(c)), f_(std::forward<G>(g)) {}
transform_iterator(const transform_iterator&) = default;
transform_iterator(transform_iterator&&) = default;
diff --git a/tests/unit/test_compartments.cpp b/tests/unit/test_compartments.cpp
index 176bf181..8bc034ad 100644
--- a/tests/unit/test_compartments.cpp
+++ b/tests/unit/test_compartments.cpp
@@ -2,8 +2,8 @@
#include "gtest.h"
-#include "../src/compartment.hpp"
-#include "../src/util.hpp"
+#include <compartment.hpp>
+#include <util.hpp>
using nest::mc::util::left;
using nest::mc::util::right;
@@ -12,7 +12,6 @@ using nest::mc::util::right;
// are correctly stored in members
TEST(compartments, compartment)
{
-
{
nest::mc::compartment c(100, 1.2, 2.1, 2.2);
EXPECT_EQ(c.index, 100u);
@@ -36,102 +35,6 @@ TEST(compartments, compartment)
}
}
-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, 1u);
- 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, 1u);
- 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, 2u);
- 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, 2u);
- 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, 3u);
- 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, 0u);
- EXPECT_EQ((*rng.end()).index, 10u);
- EXPECT_NE(rng.begin(), rng.end());
-
- unsigned 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, 10u);
- }
-
- // test case of zero length range
- {
- nest::mc::compartment_range rng(0, 1.0, 1.0, 0.);
-
- EXPECT_EQ(rng.begin(), rng.end());
- }
-}
-
TEST(compartments, make_compartment_range)
{
using namespace nest::mc;
--
GitLab