diff --git a/src/fvm_multicell.hpp b/src/fvm_multicell.hpp
index 67324155ccb32157bec224d50638cf111bd3e265..709273adfac72a0f485162318757703f51d9a522 100644
--- a/src/fvm_multicell.hpp
+++ b/src/fvm_multicell.hpp
@@ -224,6 +224,9 @@ void fvm_multicell<T, I>::compute_cv_area_unnormalized_capacitance(
const segment& seg,
index_vector &parent)
{
+ using util::left;
+ using util::right;
+
// precondition: group_parent_index[j] holds the correct value for
// j in [base_comp, base_comp+segment.num_compartments()].
@@ -316,16 +319,23 @@ void fvm_multicell<T, I>::initialize(
current_ = vector_type{ncomp, T{0}};
voltage_ = vector_type{ncomp, T{0}};
+ // create maps for mechanism initialization.
+ std::map<std::string, 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};
+ auto target_hi = target_handles.begin();
for (size_type i = 0; i<ncell; ++i) {
const auto& c = cells[i];
auto comps = cell_comp_part[i];
- auto parent_index = c.model().parent_index;
+ auto graph = c.model();
for (auto k: make_span(comps)) {
- group_parent_index[k] = parent_index[k-comps.first]+comps.first;
+ group_parent_index[k] = graph.parent_index[k-comps.first]+comps.first;
}
auto seg_num_compartments =
@@ -336,7 +346,37 @@ void fvm_multicell<T, I>::initialize(
auto seg_comp_part = make_partition(seg_comp_bounds, seg_num_compartments, comps.first);
for (size_type j = 0; j<nseg; ++j) {
- compute_cv_area_unnormalized_capacitance(seg_comp_part[j], c.segment(j), group_parent_index);
+ const auto& seg = c.segment(j);
+ const auto& seg_comps = seg_comp_part[j];
+
+ compute_cv_area_unnormalized_capacitance(seg_comps, seg, group_parent_index);
+
+ for (const auto& mech: seg->mechanisms()) {
+ if (mech.name()!="membrane") {
+ mech_map[mech_name().push_back(seg_comps)];
+ }
+ }
+ }
+
+ for (const auto& syn: c.synapses()) {
+ const auto& name = syn.mechanism.name();
+ std::size_t index = 0;
+ if (syn_mech_indices.count(name)==0) {
+ index = syn_mech_map.size();
+ syn_mech_indices[name] = index;
+ syn_mech_map.push_back(std::vector<size_type>{});
+ }
+ else {
+ index = syn_mech_indices[name];
+ }
+
+ size_type syn_comp = comps.first+find_compartment_index(syn.location, graph);
+
+ EXPECTS(targets_total < c.targets.size());
+ *target_hi++ = target_handle{index, syn_mech_map[index].size()};
+ ++targets_total;
+
+ syn_mech_map[index].push_back(syn_comp);
}
// normalize capacitance across cell
@@ -344,12 +384,7 @@ void fvm_multicell<T, I>::initialize(
cv_capacitance_[k] /= cv_areas_[k];
}
- // add mechanisms to look up table for instantiation below
-
- // TODO: move mech_map, syn_map insertions here etc.
-
detectors_total += c.detectors().size();
- targets_total += c.synapses().size();
probess_total += c.probes.size();
}
@@ -361,90 +396,32 @@ void fvm_multicell<T, I>::initialize(
// initalize matrix
matrix_ = matrix_type(group_parent_index);
- using util::left;
- using util::right;
-
- // create mechanisms
-
- // for each mechanism in the cell record the indexes of the segments that
- // contain the mechanism
- std::map<std::string, std::vector<unsigned>> mech_map;
-
- for (unsigned i=0; i<cell.num_segments(); ++i) {
- for (const auto& mech : cell.segment(i)->mechanisms()) {
- // FIXME : Membrane has to be a proper mechanism,
- // because it is exposed via the public interface.
- // This if statement is bad
- if(mech.name() != "membrane") {
- mech_map[mech.name()].push_back(i);
- }
- }
- }
-
- // Create the mechanism implementations with the state for each mechanism
- // instance.
- // TODO : this works well for density mechanisms (e.g. ion channels), but
- // does it work for point processes (e.g. synapses)?
- for (auto& mech : mech_map) {
- // calculate the number of compartments that contain the mechanism
- auto num_comp = 0u;
- for (auto seg : mech.second) {
- num_comp += segment_index_[seg+1] - segment_index_[seg];
- }
+ // create density mechanisms
+ std::vector<size_type> mech_comp_indices;
+ mech_comp_indices.reserve(ncomp);
- // build a vector of the indexes of the compartments that contain
- // the mechanism
- index_type compartment_index(num_comp);
- auto pos = 0u;
- for (auto seg : mech.second) {
- auto seg_size = segment_index_[seg+1] - segment_index_[seg];
- std::iota(
- compartment_index.data() + pos,
- compartment_index.data() + pos + seg_size,
- segment_index_[seg]
- );
- pos += seg_size;
+ for (auto& mech: mech_map) {
+ mech_comp_indices.clear();
+ for (auto comp_ival: mech.second) {
+ util::append(mech_comp_indices, make_span(comp_ival));
}
- // instantiate the mechanism
mechanisms_.push_back(
- mechanism_catalogue::make(mech.first, voltage_, current_, compartment_index)
+ mechanism_catalogue::make(mech.first, voltage_, current_, mech_comp_indices)
);
}
+ // create point (synapse) mechanisms
synapse_base_ = mechanisms_.size();
-
- // Create the synapse mechanism implementations together with the target handles
- std::vector<std::vector<cell_lid_type>> syn_mech_map;
- std::map<std::string, std::size_t> syn_mech_indices;
- auto target_hi = target_handles.begin();
-
- for (const auto& syn : cell.synapses()) {
- const auto& name = syn.mechanism.name();
- std::size_t index = 0;
- if (syn_mech_indices.count(name)==0) {
- index = syn_mech_map.size();
- syn_mech_indices[name] = index;
- syn_mech_map.push_back(std::vector<cell_lid_type>{});
- }
- else {
- index = syn_mech_indices[name];
- }
-
- size_type comp = find_compartment_index(syn.location, graph);
- *target_hi++ = target_handle{index, syn_mech_map[index].size()};
- syn_mech_map[index].push_back(comp);
- }
-
- for (const auto& syni : syn_mech_indices) {
+ for (const auto& syni: syn_mech_indices) {
const auto& mech_name = syni.first;
- index_type compartment_index(syn_mech_map[syni.second]);
- auto mech = mechanism_catalogue::make(mech_name, voltage_, current_, compartment_index);
+ auto mech = mechanism_catalogue::make(mech_name, voltage_, current_, syn_mech_map[syni.second]);
mech->set_areas(cv_areas_);
mechanisms_.push_back(std::move(mech));
}
+ // TODO: FROM HERE...
/////////////////////////////////////////////
// build the ion species
diff --git a/src/util/range.hpp b/src/util/range.hpp
index 7777444d0f65d7feebc03f8f3e988a234462cdf6..590bdb45d14433cd5642e9042b8adaeb4feb76f2 100644
--- a/src/util/range.hpp
+++ b/src/util/range.hpp
@@ -327,6 +327,15 @@ range<const T*> singleton_view(const T& item) {
return {&item, &item+1};
}
+/*
+ * Range/container utility functions
+ */
+
+template <typename Container, typename Seq>
+void append(Container &c, const Seq& seq) {
+ c.insert(c.end(), seq.begin(), seq.end());
+}
+
} // namespace util
} // namespace mc
} // namespace nest