diff --git a/external/modparser b/external/modparser
index b200bf6376a2dc30edea98fcc2375fc9be095135..ad8d926fd4a3e92cdcfb7b77c527550f859496d2 160000
--- a/external/modparser
+++ b/external/modparser
@@ -1 +1 @@
-Subproject commit b200bf6376a2dc30edea98fcc2375fc9be095135
+Subproject commit ad8d926fd4a3e92cdcfb7b77c527550f859496d2
diff --git a/external/vector b/external/vector
index c8678e80660cd63b293ade80ece8eed2249c1b06..6a9798d87f83246f36c999540d4a7810fb33d199 160000
--- a/external/vector
+++ b/external/vector
@@ -1 +1 @@
-Subproject commit c8678e80660cd63b293ade80ece8eed2249c1b06
+Subproject commit 6a9798d87f83246f36c999540d4a7810fb33d199
diff --git a/miniapp/io.cpp b/miniapp/io.cpp
index 344fd56135b232bef3563d6a3cc61b570868a360..348489ac6bb03491e3573d0c9d61b369325bba76 100644
--- a/miniapp/io.cpp
+++ b/miniapp/io.cpp
@@ -110,7 +110,7 @@ static void update_option(util::optional<T>& opt, const nlohmann::json& j, const
}
// Read options from (optional) json file and command line arguments.
-cl_options read_options(int argc, char** argv) {
+cl_options read_options(int argc, char** argv, bool allow_write) {
// Default options:
const cl_options defopts{
@@ -121,6 +121,8 @@ cl_options read_options(int argc, char** argv) {
100., // tfinal
0.025, // dt
false, // all_to_all
+ false, // ring
+ 1, // group_size
false, // probe_soma_only
0.0, // probe_ratio
"trace_", // trace_prefix
@@ -170,6 +172,11 @@ cl_options read_options(int argc, char** argv) {
false, defopts.dt, "time", cmd);
TCLAP::SwitchArg all_to_all_arg(
"m","alltoall","all to all network", cmd, false);
+ TCLAP::SwitchArg ring_arg(
+ "r","ring","ring network", cmd, false);
+ TCLAP::ValueArg<uint32_t> group_size_arg(
+ "g", "group-size", "number of cells per cell group",
+ false, defopts.compartments_per_segment, "integer", cmd);
TCLAP::ValueArg<double> probe_ratio_arg(
"p", "probe-ratio", "proportion between 0 and 1 of cells to probe",
false, defopts.probe_ratio, "proportion", cmd);
@@ -206,6 +213,8 @@ cl_options read_options(int argc, char** argv) {
update_option(options.dt, fopts, "dt");
update_option(options.tfinal, fopts, "tfinal");
update_option(options.all_to_all, fopts, "all_to_all");
+ update_option(options.ring, fopts, "ring");
+ update_option(options.group_size, fopts, "group_size");
update_option(options.probe_ratio, fopts, "probe_ratio");
update_option(options.probe_soma_only, fopts, "probe_soma_only");
update_option(options.trace_prefix, fopts, "trace_prefix");
@@ -239,12 +248,22 @@ cl_options read_options(int argc, char** argv) {
update_option(options.tfinal, tfinal_arg);
update_option(options.dt, dt_arg);
update_option(options.all_to_all, all_to_all_arg);
+ update_option(options.ring, ring_arg);
+ update_option(options.group_size, group_size_arg);
update_option(options.probe_ratio, probe_ratio_arg);
update_option(options.probe_soma_only, probe_soma_only_arg);
update_option(options.trace_prefix, trace_prefix_arg);
update_option(options.trace_max_gid, trace_max_gid_arg);
update_option(options.spike_file_output, spike_output_arg);
+ if (options.all_to_all && options.ring) {
+ throw usage_error("can specify at most one of --ring and --all-to-all");
+ }
+
+ if (options.group_size<1) {
+ throw usage_error("minimum of one cell per group");
+ }
+
save_file = ofile_arg.getValue();
}
catch (TCLAP::ArgException& e) {
@@ -252,7 +271,7 @@ cl_options read_options(int argc, char** argv) {
}
// Save option values if requested.
- if (save_file != "") {
+ if (save_file != "" && allow_write) {
std::ofstream fid(save_file);
if (fid) {
try {
@@ -265,6 +284,8 @@ cl_options read_options(int argc, char** argv) {
fopts["dt"] = options.dt;
fopts["tfinal"] = options.tfinal;
fopts["all_to_all"] = options.all_to_all;
+ fopts["ring"] = options.ring;
+ fopts["group_size"] = options.group_size;
fopts["probe_ratio"] = options.probe_ratio;
fopts["probe_soma_only"] = options.probe_soma_only;
fopts["trace_prefix"] = options.trace_prefix;
@@ -297,6 +318,8 @@ std::ostream& operator<<(std::ostream& o, const cl_options& options) {
o << " simulation time : " << options.tfinal << "\n";
o << " dt : " << options.dt << "\n";
o << " all to all network : " << (options.all_to_all ? "yes" : "no") << "\n";
+ o << " ring network : " << (options.ring ? "yes" : "no") << "\n";
+ o << " group size : " << options.group_size << "\n";
o << " probe ratio : " << options.probe_ratio << "\n";
o << " probe soma only : " << (options.probe_soma_only ? "yes" : "no") << "\n";
o << " trace prefix : " << options.trace_prefix << "\n";
diff --git a/miniapp/io.hpp b/miniapp/io.hpp
index c18e688c33ceff0c3966b5d66ae8eb1680bc48a1..ac769d436b6a36550afe2e64b33e0b35a69f3f80 100644
--- a/miniapp/io.hpp
+++ b/miniapp/io.hpp
@@ -21,6 +21,8 @@ struct cl_options {
double tfinal;
double dt;
bool all_to_all;
+ bool ring;
+ uint32_t group_size;
bool probe_soma_only;
double probe_ratio;
std::string trace_prefix;
@@ -49,7 +51,7 @@ public:
std::ostream& operator<<(std::ostream& o, const cl_options& opt);
-cl_options read_options(int argc, char** argv);
+cl_options read_options(int argc, char** argv, bool allow_write = true);
} // namespace io
diff --git a/miniapp/miniapp.cpp b/miniapp/miniapp.cpp
index bc4133cb15ceecc5c96b8b7649f286980acb232f..2b5df5d963305f82263d4cb3676d949989d90ea2 100644
--- a/miniapp/miniapp.cpp
+++ b/miniapp/miniapp.cpp
@@ -9,10 +9,9 @@
#include <common_types.hpp>
#include <cell.hpp>
-#include <cell_group.hpp>
#include <communication/communicator.hpp>
#include <communication/global_policy.hpp>
-#include <fvm_cell.hpp>
+#include <fvm_multicell.hpp>
#include <io/exporter_spike_file.hpp>
#include <mechanism_catalogue.hpp>
#include <model.hpp>
@@ -29,7 +28,8 @@
using namespace nest::mc;
using global_policy = communication::global_policy;
-using lowered_cell = fvm::fvm_cell<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>;
@@ -51,7 +51,7 @@ int main(int argc, char** argv) {
banner();
// read parameters
- io::cl_options options = io::read_options(argc, argv);
+ io::cl_options options = io::read_options(argc, argv, global_policy::id()==0);
std::cout << options << "\n";
std::cout << "\n";
std::cout << ":: simulation to " << options.tfinal << " ms in "
@@ -66,8 +66,16 @@ int main(int argc, char** argv) {
auto recipe = make_recipe(options, pdist);
auto cell_range = distribute_cells(recipe->num_cells());
+ std::vector<cell_gid_type> group_divisions;
+ for (auto i = cell_range.first; i<cell_range.second; i+=options.group_size) {
+ group_divisions.push_back(i);
+ }
+ group_divisions.push_back(cell_range.second);
+
+ EXPECTS(group_divisions.front() == cell_range.first);
+ EXPECTS(group_divisions.back() == cell_range.second);
- model_type m(*recipe, cell_range.first, cell_range.second);
+ model_type m(*recipe, util::partition_view(group_divisions));
auto register_exporter = [] (const io::cl_options& options) {
return
@@ -183,6 +191,9 @@ std::unique_ptr<recipe> make_recipe(const io::cl_options& options, const probe_d
if (options.all_to_all) {
return make_basic_kgraph_recipe(options.cells, p, pdist);
}
+ else if (options.ring) {
+ return make_basic_ring_recipe(options.cells, p, pdist);
+ }
else {
return make_basic_rgraph_recipe(options.cells, p, pdist);
}
diff --git a/scripts/tsplot b/scripts/tsplot
index 81a5da78c06a24d71f727431ef237a18f2a13b0b..15977538a3ce5d977c7559e536d2e03af845e4e4 100755
--- a/scripts/tsplot
+++ b/scripts/tsplot
@@ -33,13 +33,19 @@ def parse_clargs():
P.add_argument('inputs', metavar='FILE', nargs='+',
help='time series data in JSON format')
P.add_argument('-t', '--trange', metavar='RANGE', dest='trange',
- type=parse_range_spec,
+ type=parse_range_spec,
help='restrict time axis to RANGE (see below)')
P.add_argument('-g', '--group', metavar='KEY,...', dest='groupby',
- type=lambda s: s.split(','),
+ type=lambda s: s.split(','),
help='plot series with same KEYs on the same axes')
P.add_argument('-o', '--output', metavar='FILE', dest='outfile',
help='save plot to file FILE')
+ P.add_argument('--dpi', metavar='NUM', dest='dpi',
+ type=int,
+ help='set dpi for output image')
+ P.add_argument('--scale', metavar='NUM', dest='scale',
+ type=float,
+ help='scale size of output image by NUM')
P.add_argument('-x', '--exclude', metavar='NUM', dest='exclude',
type=float,
help='remove extreme points outside NUM times the 0.9-interquantile range of the median')
@@ -82,7 +88,7 @@ class TimeSeries:
l_, lq, median, uq, u_ = np.percentile(yfinite, [0, 5.0, 50.0, 95.0, 100])
lb = median - iqr_factor*(uq-lq)
ub = median + iqr_factor*(uq-lq)
-
+
np_err_save = np.seterr(all='ignore')
yex = np.ma.masked_where(np.isfinite(self.y)&(self.y<=ub)&(self.y>=lb), self.y)
np.seterr(**np_err_save)
@@ -104,7 +110,7 @@ class TimeSeries:
def units(self):
return self.meta.get('units',"")
-
+
def trange(self):
return (min(self.t), max(self.t))
@@ -155,7 +161,7 @@ def read_json_timeseries(source):
ts_list.append(TimeSeries(times, jdata[key], **meta))
return ts_list
-
+
def range_join(r, s):
return (min(r[0], s[0]), max(r[1], s[1]))
@@ -236,7 +242,7 @@ def make_palette(cm_name, n, cmin=0, cmax=1):
M.cm.get_cmap(cm_name))
return [smap.to_rgba((2*i+1)/float(2*n)) for i in xrange(n)]
-def plot_plots(plot_groups, save=None):
+def plot_plots(plot_groups, save=None, dpi=None, scale=None):
nplots = len(plot_groups)
plot_groups = sorted(plot_groups, key=lambda g: g.group_label())
@@ -268,7 +274,7 @@ def plot_plots(plot_groups, save=None):
lab = unit
return lab
-
+
uniq_units = list(set([s.units() for s in group.series]))
uniq_units.sort()
if len(uniq_units)>2:
@@ -297,7 +303,7 @@ def plot_plots(plot_groups, save=None):
cm, n in zip(['hot', 'winter'], [len(x) for x in series_by_unit])]
lines = cycle(["-",(0,(3,1))])
-
+
first_plot = True
for ui in xrange(len(uniq_units)):
if not first_plot:
@@ -310,7 +316,7 @@ def plot_plots(plot_groups, save=None):
plot.get_yaxis().get_major_formatter().set_useOffset(False)
plot.get_yaxis().set_major_locator(M.ticker.MaxNLocator(nbins=6))
-
+
plot.set_xlim(trange)
colours = cycle(palette[ui])
@@ -336,10 +342,14 @@ def plot_plots(plot_groups, save=None):
axis_ymin = min([ax.get_position().ymin for ax in figure.axes])
figure.text(0.5, axis_ymin - float(3)/figure.dpi, 'time', ha='center', va='center')
if save:
- figure.savefig(save)
+ if scale:
+ base = figure.get_size_inches()
+ figure.set_size_inches((base[0]*scale, base[1]*scale))
+
+ figure.savefig(save, dpi=dpi)
else:
P.show()
-
+
args = parse_clargs()
tss = []
for filename in args.inputs:
@@ -360,4 +370,4 @@ plots = gather_ts_plots(tss, groupby)
if not args.outfile:
M.interactive(False)
-plot_plots(plots, save=args.outfile)
+plot_plots(plots, save=args.outfile, dpi=args.dpi, scale=args.scale)
diff --git a/src/algorithms.hpp b/src/algorithms.hpp
index 37eeaf370acd54ef326f70cc4ae17e46f1c6f2b0..7be1406ac7581df6924c146e6b86aaa8201755bb 100644
--- a/src/algorithms.hpp
+++ b/src/algorithms.hpp
@@ -51,6 +51,17 @@ C make_index(C const& c)
return out;
}
+/// test for membership within half-open interval
+template <typename X, typename I, typename J>
+bool in_interval(const X& x, const I& lower, const J& upper) {
+ return x>=lower && x<upper;
+}
+
+template <typename X, typename I, typename J>
+bool in_interval(const X& x, const std::pair<I, J>& bounds) {
+ return x>=bounds.first && x<bounds.second;
+}
+
/// works like std::is_sorted(), but with stronger condition that succesive
/// elements must be greater than those before them
template <typename C>
@@ -243,7 +254,7 @@ std::vector<typename C::value_type> make_parent_index(
return {};
}
- EXPECTS(parent_index.size() == branch_index.back());
+ EXPECTS(parent_index.size()-branch_index.back() == 0);
EXPECTS(has_contiguous_compartments(parent_index));
EXPECTS(is_strictly_monotonic_increasing(branch_index));
diff --git a/src/cell.cpp b/src/cell.cpp
index e96cae69aca834d61d81c9c53c9488532a3b5b69..f247a2a8a83e327abce312ca31b31c5e560b8ce6 100644
--- a/src/cell.cpp
+++ b/src/cell.cpp
@@ -188,7 +188,7 @@ void cell::add_stimulus(segment_location loc, i_clamp stim)
)
);
}
- stimulii_.push_back({loc, std::move(stim)});
+ stimuli_.push_back({loc, std::move(stim)});
}
void cell::add_detector(segment_location loc, double threshold)
diff --git a/src/cell.hpp b/src/cell.hpp
index 746549289f1b718433388e92f5f1ca9d13d433c1..bd83b1b182238377cdd9e7642c60d0be31fe4892 100644
--- a/src/cell.hpp
+++ b/src/cell.hpp
@@ -129,18 +129,18 @@ public:
compartment_model model() const;
//////////////////
- // stimulii
+ // stimuli
//////////////////
void add_stimulus(segment_location loc, i_clamp stim);
std::vector<stimulus_instance>&
- stimulii() {
- return stimulii_;
+ stimuli() {
+ return stimuli_;
}
const std::vector<stimulus_instance>&
- stimulii() const {
- return stimulii_;
+ stimuli() const {
+ return stimuli_;
}
//////////////////
@@ -188,8 +188,8 @@ private:
// the segments
std::vector<segment_ptr> segments_;
- // the stimulii
- std::vector<stimulus_instance> stimulii_;
+ // the stimuli
+ std::vector<stimulus_instance> stimuli_;
// the synapses
std::vector<synapse_instance> synapses_;
diff --git a/src/cell_group.hpp b/src/cell_group.hpp
index 77276f49e52f12dc2b8205a7d24113d7409365f4..86a9059361110259d72a85389adf81dda55c8fdf 100644
--- a/src/cell_group.hpp
+++ b/src/cell_group.hpp
@@ -2,14 +2,18 @@
#include <cstdint>
#include <functional>
+#include <iterator>
#include <vector>
+#include <algorithms.hpp>
#include <cell.hpp>
#include <common_types.hpp>
#include <event_queue.hpp>
#include <spike.hpp>
#include <spike_source.hpp>
-#include <util/singleton.hpp>
+#include <util/debug.hpp>
+#include <util/partition.hpp>
+#include <util/range.hpp>
#include <profiling/profiler.hpp>
@@ -36,28 +40,37 @@ public:
cell_group() = default;
- cell_group(cell_gid_type gid, const cell& c) :
- gid_base_{gid}
+ template <typename Cells>
+ cell_group(cell_gid_type first_gid, const Cells& cells):
+ gid_base_{first_gid}
{
- detector_handles_.resize(c.detectors().size());
- target_handles_.resize(c.synapses().size());
- probe_handles_.resize(c.probes().size());
+ // Create lookup structure for probe and target ids.
+ build_handle_partitions(cells);
+ std::size_t n_probes = probe_handle_divisions_.back();
+ std::size_t n_targets = target_handle_divisions_.back();
+ std::size_t n_detectors =
+ algorithms::sum(util::transform_view(cells, [](const cell& c) { return c.detectors().size(); }));
- cell_.initialize(util::singleton_view(c), detector_handles_, target_handles_, probe_handles_);
+ // Allocate space to store handles.
+ detector_handles_.resize(n_detectors);
+ target_handles_.resize(n_targets);
+ probe_handles_.resize(n_probes);
- // Create spike detectors and associate them with globally unique source ids,
- // as specified by cell gid and cell-local zero-based index.
+ cell_.initialize(cells, detector_handles_, target_handles_, probe_handles_);
+ // Create spike detectors and associate them with globally unique source ids.
cell_gid_type source_gid = gid_base_;
- cell_lid_type source_lid = 0u;
-
unsigned i = 0;
- for (auto& d : c.detectors()) {
- cell_member_type source_id{source_gid, source_lid++};
-
- spike_sources_.push_back({
- source_id, spike_detector_type(cell_, detector_handles_[i], d.threshold, 0.f)
- });
+ for (const auto& cell: cells) {
+ cell_lid_type source_lid = 0u;
+ for (auto& d: cell.detectors()) {
+ cell_member_type source_id{source_gid, source_lid++};
+
+ spike_sources_.push_back({
+ source_id, spike_detector_type(cell_, detector_handles_[i++], d.threshold, 0.f)
+ });
+ }
+ ++source_gid;
}
}
@@ -65,7 +78,7 @@ public:
clear_spikes();
clear_events();
reset_samplers();
- //initialize_cells();
+ cell_.reset();
for (auto& spike_source: spike_sources_) {
spike_source.source.reset(cell_, 0.f);
}
@@ -112,13 +125,13 @@ public:
// apply events
if (next) {
- auto handle = target_handles_[next->target.index];
+ auto handle = get_target_handle(next->target);
cell_.deliver_event(handle, next->weight);
// apply events that are due within some epsilon of the current
// time step. This should be a parameter. e.g. with for variable
// order time stepping, use the minimum possible time step size.
while(auto e = events_.pop_if_before(cell_.time()+dt/10.)) {
- auto handle = target_handles_[e->target.index];
+ auto handle = get_target_handle(e->target);
cell_.deliver_event(handle, e->weight);
}
}
@@ -151,9 +164,10 @@ public:
}
void add_sampler(cell_member_type probe_id, sampler_function s, time_type start_time = 0) {
- EXPECTS(probe_id.gid==gid_base_);
+ auto handle = get_probe_handle(probe_id);
+
auto sampler_index = uint32_t(samplers_.size());
- samplers_.push_back({probe_handles_[probe_id.index], s});
+ samplers_.push_back({handle, s});
sampler_start_times_.push_back(start_time);
sample_events_.push({sampler_index, start_time});
}
@@ -173,7 +187,7 @@ public:
}
value_type probe(cell_member_type probe_id) const {
- return cell_.probe(probe_handles_[probe_id.index]);
+ return cell_.probe(get_probe_handle(probe_id));
}
private:
@@ -216,6 +230,50 @@ private:
/// collection of samplers to be run against probes in this group
std::vector<sampler_entry> samplers_;
+
+ /// lookup table for probe ids -> local probe handle indices
+ std::vector<std::size_t> probe_handle_divisions_;
+
+ /// lookup table for target ids -> local target handle indices
+ std::vector<std::size_t> target_handle_divisions_;
+
+ /// build handle index lookup tables
+ template <typename Cells>
+ void build_handle_partitions(const Cells& cells) {
+ auto probe_counts = util::transform_view(cells, [](const cell& c) { return c.probes().size(); });
+ auto target_counts = util::transform_view(cells, [](const cell& c) { return c.synapses().size(); });
+
+ make_partition(probe_handle_divisions_, probe_counts);
+ make_partition(target_handle_divisions_, target_counts);
+ }
+
+ /// use handle partition to get index from id
+ template <typename Divisions>
+ std::size_t handle_partition_lookup(const Divisions& divisions, cell_member_type id) const {
+ // NB: without any assertion checking, this would just be:
+ // return divisions[id.gid-gid_base_]+id.index;
+
+ EXPECTS(id.gid>=gid_base_);
+
+ auto handle_partition = util::partition_view(divisions);
+ EXPECTS(id.gid-gid_base_<handle_partition.size());
+
+ auto ival = handle_partition[id.gid-gid_base_];
+ std::size_t i = ival.first + id.index;
+ EXPECTS(i<ival.second);
+
+ return i;
+ }
+
+ /// get probe handle from probe id
+ probe_handle get_probe_handle(cell_member_type probe_id) const {
+ return probe_handles_[handle_partition_lookup(probe_handle_divisions_, probe_id)];
+ }
+
+ /// get target handle from target id
+ target_handle get_target_handle(cell_member_type target_id) const {
+ return target_handles_[handle_partition_lookup(target_handle_divisions_, target_id)];
+ }
};
} // namespace mc
diff --git a/src/communication/communicator.hpp b/src/communication/communicator.hpp
index 13d4a54303da74aefbcab0bc4d0ecc3636d7f472..764020dcfc4accda25c7074469fa515f7f755ab4 100644
--- a/src/communication/communicator.hpp
+++ b/src/communication/communicator.hpp
@@ -13,6 +13,7 @@
#include <spike.hpp>
#include <util/debug.hpp>
#include <util/double_buffer.hpp>
+#include <util/partition.hpp>
namespace nest {
namespace mc {
@@ -41,19 +42,18 @@ public:
using event_queue =
std::vector<postsynaptic_spike_event<time_type>>;
- communicator() = default;
+ using gid_partition_type =
+ util::partition_range<std::vector<cell_gid_type>::const_iterator>;
- // TODO
- // for now, still assuming one-to-one association cells <-> groups,
- // so that 'group' gids as represented by their first cell gid are
- // contiguous.
- communicator(id_type cell_from, id_type cell_to):
- cell_gid_from_(cell_from), cell_gid_to_(cell_to)
+ communicator() {}
+
+ explicit communicator(gid_partition_type cell_gid_partition):
+ cell_gid_partition_(cell_gid_partition)
{}
cell_local_size_type num_groups_local() const
{
- return cell_gid_to_-cell_gid_from_;
+ return cell_gid_partition_.size();
}
void add_connection(connection_type con) {
@@ -63,7 +63,7 @@ public:
/// returns true if the cell with gid is on the domain of the caller
bool is_local_cell(id_type gid) const {
- return gid>=cell_gid_from_ && gid<cell_gid_to_;
+ return algorithms::in_interval(gid, cell_gid_partition_.bounds());
}
/// builds the optimized data structure
@@ -135,9 +135,8 @@ public:
private:
std::size_t cell_group_index(cell_gid_type cell_gid) const {
- // this will be more elaborate when there is more than one cell per cell group
- EXPECTS(cell_gid>=cell_gid_from_ && cell_gid<cell_gid_to_);
- return cell_gid-cell_gid_from_;
+ EXPECTS(is_local_cell(cell_gid));
+ return cell_gid_partition_.index(cell_gid);
}
std::vector<connection_type> connections_;
@@ -145,8 +144,8 @@ private:
communication_policy_type communication_policy_;
uint64_t num_spikes_ = 0u;
- id_type cell_gid_from_;
- id_type cell_gid_to_;
+
+ gid_partition_type cell_gid_partition_;
};
} // namespace communication
diff --git a/src/communication/mpi.hpp b/src/communication/mpi.hpp
index 9ebb0ca412d83dbb600471417beb299fd1a55a8b..472c64039ba8b4b74e96adf97cc6f63782a69f72 100644
--- a/src/communication/mpi.hpp
+++ b/src/communication/mpi.hpp
@@ -11,6 +11,7 @@
#include <algorithms.hpp>
#include <communication/gathered_vector.hpp>
+#include <util/debug.hpp>
namespace nest {
namespace mc {
@@ -97,7 +98,7 @@ namespace mpi {
}
template <typename T>
- std::vector<T> gather_all(const std::vector<T> &values) {
+ std::vector<T> gather_all(const std::vector<T>& values) {
static_assert(
true,//std::is_trivially_copyable<T>::value,
"gather_all can only be performed on trivally copyable types");
@@ -149,6 +150,10 @@ namespace mpi {
MPI_COMM_WORLD
);
+ for (auto& d : displs) {
+ d /= traits::count();
+ }
+
return gathered_type(
std::move(buffer),
std::vector<count_type>(displs.begin(), displs.end())
diff --git a/src/compartment.hpp b/src/compartment.hpp
index da6bbcf57298593a60f7510822832060c7a89846..083155ffb7c5a3ea0d9d089740b6d7a3345ae234 100644
--- a/src/compartment.hpp
+++ b/src/compartment.hpp
@@ -3,7 +3,10 @@
#include <iterator>
#include <utility>
-#include "common_types.hpp"
+#include <common_types.hpp>
+#include <util/counter.hpp>
+#include <util/span.hpp>
+#include <util/transform.hpp>
namespace nest {
namespace mc {
@@ -34,132 +37,47 @@ 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_;
-};
+// (NB: auto type deduction and lambda in C++14 will simplify the following)
-class compartment_range {
+template <typename size_type, typename real_type>
+class compartment_maker {
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_maker(size_type n, real_type length, real_type rL, real_type rR):
+ r0_{rL},
+ dr_{(rR-rL)/n},
+ dx_{length/n}
{}
- compartment_iterator begin() const
- {
- return {0, compartment_length(), radius_L_, delta_radius()};
- }
-
- compartment_iterator cbegin() const
- {
- return begin();
+ compartment operator()(size_type i) const {
+ return compartment(i, dx_, r0_+i*dr_, r0_+(i+1)*dr_);
}
- /// With 0-based indexing compartment number "num_compartments_" is
- /// one past the end
- compartment_iterator end() const
- {
- return {num_compartments_, 0, 0, 0};
- }
+private:
+ real_type r0_;
+ real_type dr_;
+ real_type dx_;
+};
- compartment_iterator cend() const
- {
- return end();
- }
+template <typename size_type, typename real_type>
+using compartment_iterator =
+ util::transform_iterator<util::counter<size_type>, compartment_maker<size_type, real_type>>;
- real_type delta_radius() const
- {
- return (radius_R_ - radius_L_) / num_compartments_;
- }
+template <typename size_type, typename real_type>
+using compartment_range =
+ util::range<compartment_iterator<size_type, real_type>>;
- 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_;
-};
+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,
+ real_type length)
+{
+ return util::transform_view(
+ util::span<size_type>(0, num_compartments),
+ compartment_maker<size_type, real_type>(num_compartments, length, radius_L, radius_R));
+}
} // namespace mc
} // namespace nest
diff --git a/src/fvm_cell.hpp b/src/fvm_cell.hpp
index 83fd8b769974f8447b3930d8f2477a0732d139fe..1d9ff96568dc3c2817944145330e9a7178c19e05 100644
--- a/src/fvm_cell.hpp
+++ b/src/fvm_cell.hpp
@@ -288,7 +288,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_cell::*, uint32_t>> probes_;
@@ -531,10 +531,10 @@ void fvm_cell<T, I>::initialize(
ion_ca().internal_concentration()(all) = 5e-5; // mM
ion_ca().external_concentration()(all) = 2.0; // mM
- // add the stimulii
- for(const auto& stim : cell.stimulii()) {
+ // add the stimuli
+ for(const auto& stim : cell.stimuli()) {
auto idx = find_compartment_index(stim.location, graph);
- stimulii_.push_back( {idx, stim.clamp} );
+ stimuli_.push_back( {idx, stim.clamp} );
}
// record probe locations by index into corresponding state vector
@@ -633,8 +633,8 @@ void fvm_cell<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/fvm_multicell.hpp b/src/fvm_multicell.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..72d22b25270bb1eb71a7a3307e39c2ad17bec62d
--- /dev/null
+++ b/src/fvm_multicell.hpp
@@ -0,0 +1,625 @@
+#pragma once
+
+#include <algorithm>
+#include <iterator>
+#include <map>
+#include <set>
+#include <string>
+#include <vector>
+
+#include <algorithms.hpp>
+#include <cell.hpp>
+#include <event_queue.hpp>
+#include <ion.hpp>
+#include <math.hpp>
+#include <matrix.hpp>
+#include <mechanism.hpp>
+#include <mechanism_catalogue.hpp>
+#include <profiling/profiler.hpp>
+#include <segment.hpp>
+#include <stimulus.hpp>
+#include <util.hpp>
+#include <util/debug.hpp>
+#include <util/meta.hpp>
+#include <util/partition.hpp>
+#include <util/span.hpp>
+
+#include <vector/include/Vector.hpp>
+
+namespace nest {
+namespace mc {
+namespace fvm {
+
+template <typename Value, typename Index>
+class fvm_multicell {
+public:
+ fvm_multicell() = default;
+
+ /// the real number type
+ using value_type = Value;
+
+ /// the integral index type
+ using size_type = Index;
+
+ /// the container used for indexes
+ using index_type = memory::HostVector<size_type>;
+
+ /// the container used for values
+ using vector_type = memory::HostVector<value_type>;
+
+ /// API for cell_group (see above):
+
+ using detector_handle = size_type;
+ using target_handle = std::pair<size_type, size_type>;
+ using probe_handle = std::pair<const vector_type fvm_multicell::*, size_type>;
+
+ void resting_potential(value_type potential_mV) {
+ resting_potential_ = potential_mV;
+ }
+
+ template <typename Cells, typename Detectors, typename Targets, typename Probes>
+ void initialize(
+ const Cells& cells, // collection of nest::mc::cell descriptions
+ Detectors& detector_handles, // (write) where to store detector handles
+ Targets& target_handles, // (write) where to store target handles
+ Probes& probe_handles); // (write) where to store probe handles
+
+ void reset();
+
+ void deliver_event(target_handle h, value_type weight) {
+ mechanisms_[synapse_base_+h.first]->net_receive(h.second, weight);
+ }
+
+ value_type detector_voltage(detector_handle h) const {
+ return voltage_[h]; // detector_handle is just the compartment index
+ }
+
+ value_type probe(probe_handle h) const {
+ return (this->*h.first)[h.second];
+ }
+
+ void advance(value_type dt);
+
+ /// Following types and methods are public only for testing:
+
+ /// the type used to store matrix information
+ using matrix_type = matrix<value_type, size_type>;
+
+ /// mechanism type
+ using mechanism_type =
+ nest::mc::mechanisms::mechanism_ptr<value_type, size_type>;
+
+ /// ion species storage
+ using ion_type = mechanisms::ion<value_type, size_type>;
+
+ /// view into index container
+ using index_view = typename index_type::view_type;
+ using const_index_view = typename index_type::const_view_type;
+
+ /// view into value container
+ using vector_view = typename vector_type::view_type;
+ using const_vector_view = typename vector_type::const_view_type;
+
+ /// build the matrix for a given time step
+ void setup_matrix(value_type dt);
+
+ /// which requires const_view in the vector library
+ const matrix_type& jacobian() {
+ return matrix_;
+ }
+
+ /// return list of CV areas in :
+ /// um^2
+ /// 1e-6.mm^2
+ /// 1e-8.cm^2
+ const_vector_view cv_areas() const {
+ return cv_areas_;
+ }
+
+ /// return the capacitance of each CV surface
+ /// this is the total capacitance, not per unit area,
+ /// i.e. equivalent to sigma_i * c_m
+ const_vector_view cv_capacitance() const {
+ return cv_capacitance_;
+ }
+
+ /// return the voltage in each CV
+ vector_view voltage() { return voltage_; }
+ const_vector_view voltage() const { return voltage_; }
+
+ /// return the current in each CV
+ vector_view current() { return current_; }
+ const_vector_view current() const { return current_; }
+
+ std::size_t size() const { return matrix_.size(); }
+
+ /// return reference to in iterable container of the mechanisms
+ std::vector<mechanism_type>& mechanisms() { return mechanisms_; }
+
+ /// return reference to list of ions
+ std::map<mechanisms::ionKind, ion_type>& ions() { return ions_; }
+ std::map<mechanisms::ionKind, ion_type> const& ions() const { return ions_; }
+
+ /// return reference to sodium ion
+ ion_type& ion_na() { return ions_[mechanisms::ionKind::na]; }
+ ion_type const& ion_na() const { return ions_[mechanisms::ionKind::na]; }
+
+ /// return reference to calcium ion
+ ion_type& ion_ca() { return ions_[mechanisms::ionKind::ca]; }
+ ion_type const& ion_ca() const { return ions_[mechanisms::ionKind::ca]; }
+
+ /// return reference to pottasium ion
+ ion_type& ion_k() { return ions_[mechanisms::ionKind::k]; }
+ ion_type const& ion_k() const { return ions_[mechanisms::ionKind::k]; }
+
+ /// flags if solution is physically realistic.
+ /// here we define physically realistic as the voltage being within reasonable bounds.
+ /// use a simple test of the voltage at the soma is reasonable, i.e. in the range
+ /// v_soma \in (-1000mv, 1000mv)
+ bool is_physical_solution() const {
+ auto v = voltage_[0];
+ return (v>-1000.) && (v<1000.);
+ }
+
+ value_type time() const { return t_; }
+
+ std::size_t num_probes() const { return probes_.size(); }
+
+private:
+ /// current time [ms]
+ value_type t_ = value_type{0};
+
+ /// resting potential (initial voltage condition)
+ value_type resting_potential_ = -65;
+
+ /// the linear system for implicit time stepping of cell state
+ matrix_type matrix_;
+
+ /// index for fast lookup of compartment index ranges of segments
+ index_type segment_index_;
+
+ /// cv_areas_[i] is the surface area of CV i [µm^2]
+ vector_type cv_areas_;
+
+ /// alpha_[i] is the following value at the CV face between
+ /// CV i and its parent, required when constructing linear system
+ /// face_alpha_[i] = area_face / (c_m * r_L * delta_x);
+ vector_type face_alpha_; // [µm·m^2/cm/s ≡ 10^5 µm^2/ms]
+
+ /// cv_capacitance_[i] is the capacitance of CV i per unit area (i.e. c_m) [F/m^2]
+ vector_type cv_capacitance_;
+
+ /// the average current density over the surface of each CV [mA/cm^2]
+ /// current_ = i_m - i_e
+ vector_type current_;
+
+ /// the potential in each CV [mV]
+ vector_type voltage_;
+
+ /// Where point mechanisms start in the mechanisms_ list.
+ std::size_t synapse_base_;
+
+ /// the set of mechanisms present in the cell
+ std::vector<mechanism_type> mechanisms_;
+
+ /// the ion species
+ std::map<mechanisms::ionKind, ion_type> ions_;
+
+ std::vector<std::pair<uint32_t, i_clamp>> stimuli_;
+
+ std::vector<std::pair<const vector_type fvm_multicell::*, uint32_t>> probes_;
+
+ // mechanism factory
+ using mechanism_catalogue = nest::mc::mechanisms::catalogue<value_type, size_type>;
+
+ // perform area and capacitance calculation on initialization
+ void compute_cv_area_unnormalized_capacitance(
+ std::pair<size_type, size_type> comp_ival, const segment* seg, index_type &parent);
+};
+
+////////////////////////////////////////////////////////////////////////////////
+//////////////////////////////// Implementation ////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////
+
+template <typename T, typename I>
+void fvm_multicell<T, I>::compute_cv_area_unnormalized_capacitance(
+ std::pair<size_type, size_type> comp_ival,
+ const segment* seg,
+ index_type &parent)
+{
+ using util::left;
+ using util::right;
+ using util::size;
+
+ // precondition: group_parent_index[j] holds the correct value for
+ // j in [base_comp, base_comp+segment.num_compartments()].
+
+ auto ncomp = comp_ival.second-comp_ival.first;
+
+ if (auto soma = seg->as_soma()) {
+ // confirm assumption that there is one compartment in soma
+ if (ncomp!=1) {
+ throw std::logic_error("soma allocated more than one compartment");
+ }
+ auto i = comp_ival.first;
+ auto area = math::area_sphere(soma->radius());
+
+ cv_areas_[i] += area;
+ cv_capacitance_[i] += area * soma->mechanism("membrane").get("c_m").value;
+ }
+ 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
+ //
+ // __________________________________
+ // | ........ | .cvleft. | cv |
+ // | ........ L ........ C R
+ // |__________|__________|__________|
+ //
+ // The compartment has end points marked L and R (left and right).
+ // The left compartment is assumed to be closer to the soma
+ // (i.e. it follows the minimal degree ordering)
+ // The face is at the center, marked C.
+ // The full control volume to the left (marked with .)
+
+ auto c_m = cable->mechanism("membrane").get("c_m").value;
+ auto r_L = cable->mechanism("membrane").get("r_L").value;
+ const auto& compartments = cable->compartments();
+
+ EXPECTS(size(compartments)==ncomp);
+
+ for (auto i: util::make_span(comp_ival)) {
+ const auto& c = compartments[i-comp_ival.first];
+ auto j = parent[i];
+
+ auto radius_center = math::mean(c.radius);
+ auto area_face = math::area_circle(radius_center);
+ face_alpha_[i] = area_face / (c_m * r_L * c.length);
+
+ auto halflen = c.length/2;
+ auto al = math::area_frustrum(halflen, left(c.radius), radius_center);
+ auto ar = math::area_frustrum(halflen, right(c.radius), radius_center);
+
+ cv_areas_[j] += al;
+ cv_areas_[i] += ar;
+ cv_capacitance_[j] += al * c_m;
+ cv_capacitance_[i] += ar * c_m;
+ }
+ }
+ else {
+ throw std::domain_error("FVM lowering encountered unsuported segment type");
+ }
+}
+
+template <typename T, typename I>
+template <typename Cells, typename Detectors, typename Targets, typename Probes>
+void fvm_multicell<T, I>::initialize(
+ const Cells& cells,
+ Detectors& detector_handles,
+ Targets& target_handles,
+ Probes& probe_handles)
+{
+ using util::transform_view;
+ using util::make_partition;
+ using util::make_span;
+ using util::size;
+
+ // count total detectors, targets and probes for validation of handle container sizes
+ size_type detectors_count = 0u;
+ size_type targets_count = 0u;
+ size_type probes_count = 0u;
+ size_type detectors_size = size(detector_handles);
+ size_type targets_size = size(target_handles);
+ size_type probes_size = size(probe_handles);
+
+ auto ncell = size(cells);
+ auto cell_num_compartments =
+ transform_view(cells, [](const cell& c) { return c.num_compartments(); });
+
+ std::vector<cell_lid_type> cell_comp_bounds;
+ auto cell_comp_part = make_partition(cell_comp_bounds, cell_num_compartments);
+ auto ncomp = cell_comp_part.bounds().second;
+
+ // initialize storage from total compartment count
+ cv_areas_ = vector_type{ncomp, T{0}};
+ face_alpha_ = vector_type{ncomp, T{0}};
+ cv_capacitance_ = vector_type{ncomp, T{0}};
+ current_ = vector_type{ncomp, T{0}};
+ voltage_ = vector_type{ncomp, T{resting_potential_}};
+
+ // create maps for mechanism initialization.
+ 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_type group_parent_index{ncomp, 0};
+
+ // create each cell:
+ auto target_hi = target_handles.begin();
+ auto detector_hi = detector_handles.begin();
+ auto probe_hi = probe_handles.begin();
+
+ for (size_type i = 0; i<ncell; ++i) {
+ const auto& c = cells[i];
+ auto comp_ival = cell_comp_part[i];
+
+ auto graph = c.model();
+
+ for (auto k: make_span(comp_ival)) {
+ group_parent_index[k] = graph.parent_index[k-comp_ival.first]+comp_ival.first;
+ }
+
+ auto seg_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;
+ auto seg_comp_part = make_partition(seg_comp_bounds, seg_num_compartments, comp_ival.first);
+
+ for (size_type j = 0; j<nseg; ++j) {
+ const auto& seg = c.segment(j);
+ const auto& seg_comp_ival = seg_comp_part[j];
+
+ compute_cv_area_unnormalized_capacitance(seg_comp_ival, seg, group_parent_index);
+
+ for (const auto& mech: seg->mechanisms()) {
+ if (mech.name()!="membrane") {
+ mech_map[mech.name()].push_back(seg_comp_ival);
+ }
+ }
+ }
+
+ for (const auto& syn: c.synapses()) {
+ EXPECTS(targets_count < targets_size);
+
+ const auto& name = syn.mechanism.name();
+ std::size_t syn_mech_index = 0;
+ if (syn_mech_indices.count(name)==0) {
+ syn_mech_index = syn_mech_map.size();
+ syn_mech_indices[name] = syn_mech_index;
+ syn_mech_map.push_back(std::vector<size_type>{});
+ }
+ else {
+ syn_mech_index = syn_mech_indices[name];
+ }
+
+ auto& map_entry = syn_mech_map[syn_mech_index];
+
+ size_type syn_comp = comp_ival.first+find_compartment_index(syn.location, graph);
+ size_type syn_index = map_entry.size();
+ map_entry.push_back(syn_comp);
+
+ *target_hi++ = target_handle{syn_mech_index, syn_index};
+ ++targets_count;
+ }
+
+ // normalize capacitance across cell
+ for (auto k: make_span(comp_ival)) {
+ cv_capacitance_[k] /= cv_areas_[k];
+ }
+
+ // add the stimuli
+ for (const auto& stim: c.stimuli()) {
+ auto idx = comp_ival.first+find_compartment_index(stim.location, graph);
+ stimuli_.push_back({idx, stim.clamp});
+ }
+
+ // detector handles are just their corresponding compartment indices
+ for (const auto& detector: c.detectors()) {
+ EXPECTS(detectors_count < detectors_size);
+
+ auto comp = comp_ival.first+find_compartment_index(detector.location, graph);
+ *detector_hi++ = comp;
+ ++detectors_count;
+ }
+
+ // record probe locations by index into corresponding state vector
+ for (const auto& probe: c.probes()) {
+ EXPECTS(probes_count < probes_size);
+
+ auto comp = comp_ival.first+find_compartment_index(probe.location, graph);
+ switch (probe.kind) {
+ case probeKind::membrane_voltage:
+ *probe_hi++ = {&fvm_multicell::voltage_, comp};
+ break;
+ case probeKind::membrane_current:
+ *probe_hi++ = {&fvm_multicell::current_, comp};
+ break;
+ default:
+ throw std::logic_error("unrecognized probeKind");
+ }
+ ++probes_count;
+ }
+ }
+
+ // confirm write-parameters were appropriately sized
+ EXPECTS(detectors_size==detectors_count);
+ EXPECTS(targets_size==targets_count);
+ EXPECTS(probes_size==probes_count);
+
+ // initalize matrix
+ matrix_ = matrix_type(group_parent_index);
+
+ // create density mechanisms
+ std::vector<size_type> mech_comp_indices;
+ mech_comp_indices.reserve(ncomp);
+
+ for (auto& mech: mech_map) {
+ mech_comp_indices.clear();
+ for (auto comp_ival: mech.second) {
+ util::append(mech_comp_indices, make_span(comp_ival));
+ }
+
+ mechanisms_.push_back(
+ mechanism_catalogue::make(mech.first, voltage_, current_, mech_comp_indices)
+ );
+ }
+
+ // create point (synapse) mechanisms
+ synapse_base_ = mechanisms_.size();
+ for (const auto& syni: syn_mech_indices) {
+ const auto& mech_name = syni.first;
+
+ 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));
+ }
+
+ // build the ion species
+ for(auto ion : mechanisms::ion_kinds()) {
+ // find the compartment indexes of all compartments that have a
+ // mechanism that depends on/influences ion
+ std::set<int> index_set;
+ for(auto& mech : mechanisms_) {
+ if(mech->uses_ion(ion)) {
+ for(auto idx : mech->node_index()) {
+ index_set.insert(idx);
+ }
+ }
+ }
+ std::vector<cell_lid_type> indexes(index_set.begin(), index_set.end());
+
+ // create the ion state
+ if(indexes.size()) {
+ ions_.emplace(ion, index_type(indexes));
+ }
+
+ // join the ion reference in each mechanism into the cell-wide ion state
+ for(auto& mech : mechanisms_) {
+ if(mech->uses_ion(ion)) {
+ mech->set_ion(ion, ions_[ion]);
+ }
+ }
+ }
+
+ // FIXME: Hard code parameters for now.
+ // Take defaults for reversal potential of sodium and potassium from
+ // the default values in Neuron.
+ // Neuron's defaults are defined in the file
+ // nrn/src/nrnoc/membdef.h
+ auto all = memory::all;
+
+ constexpr value_type DEF_vrest = -65.0; // same name as #define in Neuron
+
+ ion_na().reversal_potential()(all) = 115+DEF_vrest; // mV
+ ion_na().internal_concentration()(all) = 10.0; // mM
+ ion_na().external_concentration()(all) = 140.0; // mM
+
+ ion_k().reversal_potential()(all) = -12.0+DEF_vrest;// mV
+ ion_k().internal_concentration()(all) = 54.4; // mM
+ ion_k().external_concentration()(all) = 2.5; // mM
+
+ ion_ca().reversal_potential()(all) = 12.5 * std::log(2.0/5e-5);// mV
+ ion_ca().internal_concentration()(all) = 5e-5; // mM
+ ion_ca().external_concentration()(all) = 2.0; // mM
+
+ // initialise mechanism and voltage state
+ reset();
+}
+
+template <typename T, typename I>
+void fvm_multicell<T, I>::setup_matrix(T dt) {
+ using memory::all;
+
+ // convenience accesors to matrix storage
+ auto l = matrix_.l();
+ auto d = matrix_.d();
+ auto u = matrix_.u();
+ auto p = matrix_.p();
+ auto rhs = matrix_.rhs();
+
+ // The matrix has the following layout in memory
+ // where j is the parent index of i, i.e. i<j
+ //
+ // d[i] is the diagonal entry at a_ii
+ // u[i] is the upper triangle entry at a_ji
+ // l[i] is the lower triangle entry at a_ij
+ //
+ // d[j] . . u[i]
+ // . . .
+ // . . .
+ // l[i] . . d[i]
+ //
+
+ d(all) = cv_areas_; // [µm^2]
+ for (auto i=1u; i<d.size(); ++i) {
+ auto a = 1e5*dt * face_alpha_[i];
+
+ d[i] += a;
+ l[i] = -a;
+ u[i] = -a;
+
+ // add contribution to the diagonal of parent
+ d[p[i]] += a;
+ }
+
+ // the RHS of the linear system is
+ // V[i] - dt/cm*(im - ie)
+ auto factor = 10.*dt; // units: 10·ms/(F/m^2)·(mA/cm^2) ≡ mV
+ for(auto i=0u; i<d.size(); ++i) {
+ rhs[i] = cv_areas_[i]*(voltage_[i] - factor/cv_capacitance_[i]*current_[i]);
+ }
+}
+
+template <typename T, typename I>
+void fvm_multicell<T, I>::reset() {
+ voltage_(memory::all) = resting_potential_;
+ t_ = 0.;
+ for (auto& m : mechanisms_) {
+ m->nrn_init();
+ }
+}
+
+template <typename T, typename I>
+void fvm_multicell<T, I>::advance(T dt) {
+ using memory::all;
+
+ PE("current");
+ current_(all) = 0.;
+
+ // update currents from ion channels
+ for(auto& m : mechanisms_) {
+ PE(m->name().c_str());
+ m->set_params(t_, dt);
+ m->nrn_current();
+ PL();
+ }
+
+ // add current contributions from stimuli
+ for (auto& stim : stimuli_) {
+ auto ie = stim.second.amplitude(t_); // [nA]
+ auto loc = stim.first;
+
+ // note: current_ in [mA/cm^2], ie in [nA], cv_areas_ in [µm^2].
+ // unit scale factor: [nA/µm^2]/[mA/cm^2] = 100
+ current_[loc] -= 100*ie/cv_areas_[loc];
+ }
+ PL();
+
+ // solve the linear system
+ PE("matrix", "setup");
+ setup_matrix(dt);
+ PL(); PE("solve");
+ matrix_.solve();
+ PL();
+ voltage_(all) = matrix_.rhs();
+ PL();
+
+ // integrate state of gating variables etc.
+ PE("state");
+ for(auto& m : mechanisms_) {
+ PE(m->name().c_str());
+ m->nrn_state();
+ PL();
+ }
+ PL();
+
+ t_ += dt;
+}
+
+} // namespace fvm
+} // namespace mc
+} // namespace nest
+
diff --git a/src/mechanism.hpp b/src/mechanism.hpp
index c9c801fd198c764d61f3acf653055c8957daf00e..d826b2994f4dac39ef15a53f720acf1bf7fae603 100644
--- a/src/mechanism.hpp
+++ b/src/mechanism.hpp
@@ -1,7 +1,9 @@
#pragma once
+#include <algorithm>
#include <memory>
#include <string>
+#include <util/meta.hpp>
#include "indexed_view.hpp"
#include "ion.hpp"
@@ -16,9 +18,7 @@ enum class mechanismKind {point, density};
template <typename T, typename I>
class mechanism {
-
public:
-
using value_type = T;
using size_type = I;
@@ -27,15 +27,13 @@ public:
using view_type = typename vector_type::view_type;
using index_type = memory::HostVector<size_type>;
using index_view = typename index_type::view_type;
+ using const_index_view = typename index_type::const_view_type;
using indexed_view_type = indexed_view<value_type, size_type>;
using ion_type = ion<value_type, size_type>;
- mechanism(view_type vec_v, view_type vec_i, index_view node_index)
- : vec_v_(vec_v)
- , vec_i_(vec_i)
- , node_index_(node_index)
- , vec_area_(nullptr, 0)
+ mechanism(view_type vec_v, view_type vec_i, const_index_view node_index):
+ vec_v_(vec_v), vec_i_(vec_i), node_index_(node_index), vec_area_(nullptr, 0)
{}
std::size_t size() const
@@ -88,7 +86,7 @@ mechanism_ptr<typename M::value_type, typename M::size_type>
make_mechanism(
typename M::view_type vec_v,
typename M::view_type vec_i,
- typename M::index_view node_indices
+ typename M::const_index_view node_indices
) {
return util::make_unique<M>(vec_v, vec_i, node_indices);
}
diff --git a/src/mechanism_catalogue.hpp b/src/mechanism_catalogue.hpp
index c9d3e6de2771a1b18fa05f104dd257b8c7fe6dd9..5e79cd9f489558717de1132b502b5c211eff09f9 100644
--- a/src/mechanism_catalogue.hpp
+++ b/src/mechanism_catalogue.hpp
@@ -18,19 +18,22 @@ template <typename T, typename I>
struct catalogue {
using view_type = typename mechanism<T, I>::view_type;
using index_view = typename mechanism<T, I>::index_view;
+ using const_index_view = typename mechanism<T, I>::const_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 const& 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);
+ const_index_view node_view(node_indices);
+ return entry->second(vec_v, vec_i, node_view);
}
static bool has(const std::string& name) {
@@ -38,11 +41,15 @@ struct catalogue {
}
private:
- using maker_type = mechanism_ptr<T, I> (*)(view_type, view_type, index_view);
+ using maker_type = mechanism_ptr<T, I> (*)(view_type, view_type, const_index_view);
static const std::map<std::string, maker_type> mech_map;
template <template <typename, typename> class mech>
- static mechanism_ptr<T, I> maker(view_type vec_v, view_type vec_i, index_view node_indices) {
+ static mechanism_ptr<T, I> maker(
+ view_type vec_v,
+ view_type vec_i,
+ const_index_view node_indices)
+ {
return make_mechanism<mech<T, I>>(vec_v, vec_i, node_indices);
}
};
diff --git a/src/model.hpp b/src/model.hpp
index 5d0ce4fd8bcbe010cec3629e4b5d1adde38bcd96..2298e5af8365151654b2ca8a266f0dce2f82041b 100644
--- a/src/model.hpp
+++ b/src/model.hpp
@@ -15,6 +15,8 @@
#include <recipe.hpp>
#include <thread_private_spike_store.hpp>
#include <util/nop.hpp>
+#include <util/partition.hpp>
+#include <util/range.hpp>
#include "trace_sampler.hpp"
@@ -37,27 +39,36 @@ public:
probe_spec probe;
};
- model(const recipe& rec, cell_gid_type cell_from, cell_gid_type cell_to):
- cell_from_(cell_from),
- cell_to_(cell_to),
- communicator_(cell_from, cell_to)
+ template <typename Iter>
+ model(const recipe& rec, const util::partition_range<Iter>& groups):
+ cell_group_divisions_(groups.divisions().begin(), groups.divisions().end())
{
+ // set up communicator based on partition
+ communicator_ = communicator_type{gid_partition()};
+
// generate the cell groups in parallel, with one task per cell group
- cell_groups_ = std::vector<cell_group_type>{cell_to_-cell_from_};
+ cell_groups_ = std::vector<cell_group_type>{gid_partition().size()};
threading::parallel_vector<probe_record> probes;
- threading::parallel_for::apply(cell_from_, cell_to_,
+ threading::parallel_for::apply(0, cell_groups_.size(),
[&](cell_gid_type i) {
PE("setup", "cells");
- auto cell = rec.get_cell(i);
- auto idx = i-cell_from_;
- cell_groups_[idx] = cell_group_type(i, cell);
-
- cell_lid_type j = 0;
- for (const auto& probe: cell.probes()) {
- cell_member_type probe_id{i,j++};
- probes.push_back({probe_id, probe});
+
+ auto gids = gid_partition()[i];
+ std::vector<cell> cells{gids.second-gids.first};
+
+ for (auto gid: util::make_span(gids)) {
+ auto i = gid-gids.first;
+ cells[i] = rec.get_cell(gid);
+
+ cell_lid_type j = 0;
+ for (const auto& probe: cells[i].probes()) {
+ cell_member_type probe_id{gid, j++};
+ probes.push_back({probe_id, probe});
+ }
}
+
+ cell_groups_[i] = cell_group_type(gids.first, cells);
PL(2);
});
@@ -65,7 +76,7 @@ public:
probes_.assign(probes.begin(), probes.end());
// generate the network connections
- for (cell_gid_type i=cell_from_; i<cell_to_; ++i) {
+ for (cell_gid_type i: util::make_span(gid_partition().bounds())) {
for (const auto& cc: rec.connections_on(i)) {
connection<time_type> conn{cc.source, cc.dest, cc.weight, cc.delay};
communicator_.add_connection(conn);
@@ -187,11 +198,11 @@ public:
}
void attach_sampler(cell_member_type probe_id, sampler_function f, time_type tfrom = 0) {
- // TODO: translate probe_id.gid to appropriate group, but for now 1-1.
- if (probe_id.gid<cell_from_ || probe_id.gid>=cell_to_) {
+ if (!algorithms::in_interval(probe_id.gid, gid_partition().bounds())) {
return;
}
- cell_groups_[probe_id.gid-cell_from_].add_sampler(probe_id, f, tfrom);
+
+ cell_groups_[gid_partition().index(probe_id.gid)].add_sampler(probe_id, f, tfrom);
}
const std::vector<probe_record>& probes() const { return probes_; }
@@ -199,12 +210,14 @@ public:
std::size_t num_spikes() const {
return communicator_.num_spikes();
}
+
std::size_t num_groups() const {
return cell_groups_.size();
}
+
std::size_t num_cells() const {
- // TODO: fix when the assumption that there is one cell per cell group is no longer valid
- return num_groups();
+ auto bounds = gid_partition().bounds();
+ return bounds.second-bounds.first;
}
// register a callback that will perform a export of the global
@@ -220,18 +233,22 @@ public:
}
private:
- cell_gid_type cell_from_;
- cell_gid_type cell_to_;
+ std::vector<cell_gid_type> cell_group_divisions_;
+
+ auto gid_partition() const -> decltype(util::partition_view(cell_group_divisions_)) {
+ return util::partition_view(cell_group_divisions_);
+ }
+
time_type t_ = 0.;
std::vector<cell_group_type> cell_groups_;
communicator_type communicator_;
std::vector<probe_record> probes_;
using event_queue_type = typename communicator_type::event_queue;
- util::double_buffer< std::vector<event_queue_type> > event_queues_;
+ util::double_buffer<std::vector<event_queue_type>> event_queues_;
using local_spike_store_type = thread_private_spike_store<time_type>;
- util::double_buffer< local_spike_store_type > local_spikes_;
+ util::double_buffer<local_spike_store_type> local_spikes_;
spike_export_function global_export_callback_ = util::nop_function;
spike_export_function local_export_callback_ = util::nop_function;
diff --git a/src/parameter_list.cpp b/src/parameter_list.cpp
index 86c641aa809e6efe1c8f03247cc64f9f1b3a01f0..2ca508175f58a402165fd0d56b744148cdb2b267 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 0c06988756a077fdbb860fc8f42043b178b1b7e7..eb3ee8fe2f3c3ce1018bb113c966cb556e12916e 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 248306a6c3a976d123b9d412b93c09156305086d..5a9782b9ce41d8748bb5af6234ab9640c349f360 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/debug.hpp b/src/util/debug.hpp
index be2b3d51e4d7bf9b5d4c5fb9b7c78855c5e4b4c5..f189852dc1bd49685920af2c0effdabf5a9ac132 100644
--- a/src/util/debug.hpp
+++ b/src/util/debug.hpp
@@ -78,5 +78,6 @@ void debug_emit_trace(const char* file, int line, const char* varlist, const Arg
(void)((condition) || \
nest::mc::util::global_failed_assertion_handler(#condition, __FILE__, __LINE__, DEBUG_FUNCTION_NAME))
#else
- #define EXPECTS(condition)
+ #define EXPECTS(condition) \
+ (void)(false && (condition))
#endif // def WITH_ASSERTIONS
diff --git a/src/util/partition.hpp b/src/util/partition.hpp
index b35f81de8fa1f40217e331dd96e09a11eba1c174..425d41f164c3d6ecc27df1fd6d844f9917a55d8e 100644
--- a/src/util/partition.hpp
+++ b/src/util/partition.hpp
@@ -29,12 +29,19 @@ class partition_range: public range<partition_iterator<I>> {
public:
using typename base::iterator;
using typename base::value_type;
+ using typename base::size_type;
using base::left;
using base::right;
using base::front;
using base::back;
using base::empty;
+ // `npos` is returned by the `index()` method if the search fails;
+ // analogous to `std::string::npos`.
+ static constexpr size_type npos = static_cast<size_type>(-1);
+
+ partition_range() = default;
+
template <typename Seq>
partition_range(const Seq& s): base{std::begin(s), upto(std::begin(s), std::end(s))} {
EXPECTS(is_valid());
@@ -62,6 +69,11 @@ public:
return iterator{std::prev(i)};
}
+ size_type index(const inner_value_type& x) const {
+ iterator i = find(x);
+ return i==right? npos: i-left;
+ }
+
// access to underlying divisions
range<I> divisions() const {
return {left.get(), std::next(right.get())};
diff --git a/src/util/partition_iterator.hpp b/src/util/partition_iterator.hpp
index f13b82f25d43c604680ff4e96b3efd6936e98aea..dae1475ce91382801d6d827fd3f66f30da9c7a8d 100644
--- a/src/util/partition_iterator.hpp
+++ b/src/util/partition_iterator.hpp
@@ -37,6 +37,8 @@ public:
using pointer = const value_type*;
using reference = const value_type&;
+ partition_iterator() = default;
+
template <
typename J,
typename = enable_if_t<!std::is_same<decay_t<J>, partition_iterator>::value>
diff --git a/src/util/range.hpp b/src/util/range.hpp
index 9e3bd62a620ca40ddad3b3d178be95c4738faef1..df5625d0e6dc1e315e18d1e2e1a0e4b4ccecbd94 100644
--- a/src/util/range.hpp
+++ b/src/util/range.hpp
@@ -172,6 +172,42 @@ range<const T*> singleton_view(const T& item) {
return {&item, &item+1};
}
+/*
+ * Range/container utility functions
+ */
+
+template <typename Container, typename Seq>
+Container& append(Container &c, const Seq& seq) {
+ auto canon = canonical_view(seq);
+ c.insert(c.end(), std::begin(canon), std::end(canon));
+ return c;
+}
+
+template <typename AssignableContainer, typename Seq>
+AssignableContainer& assign(AssignableContainer& c, const Seq& seq) {
+ auto canon = canonical_view(seq);
+ c.assign(std::begin(canon), std::end(canon));
+ return c;
+}
+
+template <typename Seq>
+range<typename sequence_traits<Seq>::iterator_type, typename sequence_traits<Seq>::sentinel_type>
+range_view(Seq& seq) {
+ return make_range(std::begin(seq), std::end(seq));
+}
+
+template <
+ typename Seq,
+ typename Iter = typename sequence_traits<Seq>::iterator_type,
+ typename Size = typename sequence_traits<Seq>::size_type
+>
+enable_if_t<is_forward_iterator<Iter>::value, range<Iter>>
+subrange_view(Seq& seq, Size bi, Size ei) {
+ Iter b = std::next(std::begin(seq), bi);
+ Iter e = std::next(b, ei-bi);
+ return make_range(b, e);
+}
+
} // namespace util
} // namespace mc
} // namespace nest
diff --git a/src/util/singleton.hpp b/src/util/singleton.hpp
deleted file mode 100644
index c718eef3624af2088daa9d84abb00401ed1bd6d4..0000000000000000000000000000000000000000
--- a/src/util/singleton.hpp
+++ /dev/null
@@ -1,63 +0,0 @@
-#pragma once
-
-/*
- * Present a single object as a (non-owning) container with one
- * element.
- *
- * (Will be subsumed by range/view code.)
- */
-
-#include <algorithm>
-
-namespace nest {
-namespace mc {
-namespace util {
-
-template <typename X>
-struct singleton_adaptor {
- using size_type = std::size_t;
- using difference_type = std::ptrdiff_t;
- using value_type = X;
- using reference = X&;
- using const_reference = const X&;
- using iterator = X*;
- using const_iterator = const X*;
-
- X* xp;
- singleton_adaptor(X& x): xp(&x) {}
-
- const X* cbegin() const { return xp; }
- const X* begin() const { return xp; }
- X* begin() { return xp; }
-
- const X* cend() const { return xp+1; }
- const X* end() const { return xp+1; }
- X* end() { return xp+1; }
-
- std::size_t size() const { return 1u; }
-
- bool empty() const { return false; }
-
- const X* front() const { return *xp; }
- X* front() { return *xp; }
-
- const X* back() const { return *xp; }
- X* back() { return *xp; }
-
- const X* operator[](difference_type) const { return *xp; }
- X* operator[](difference_type) { return *xp; }
-
- void swap(singleton_adaptor& s) { std::swap(xp, s.xp); }
- friend void swap(singleton_adaptor& r, singleton_adaptor& s) {
- r.swap(s);
- }
-};
-
-template <typename X>
-singleton_adaptor<X> singleton_view(X& x) {
- return singleton_adaptor<X>(x);
-}
-
-} // namespace util
-} // namespace mc
-} // namespace nest
diff --git a/src/util/transform.hpp b/src/util/transform.hpp
index d6c930dbdc2cb653e94433ff0424fdc5612cbe3f..a5ab0b348a61fb0f3a174058dc31fdb128776e24 100644
--- a/src/util/transform.hpp
+++ b/src/util/transform.hpp
@@ -33,10 +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&;
+ transform_iterator() = default;
+
template <typename J, typename G>
transform_iterator(J&& c, G&& g): inner_(std::forward<J>(c)), f_(std::forward<G>(g)) {}
diff --git a/tests/global_communication/CMakeLists.txt b/tests/global_communication/CMakeLists.txt
index cc1e1af8b6567a2fc887815fa4eaeb25c561841e..800e7d193dcb335d2ee8ae164ed38b947eaae313 100644
--- a/tests/global_communication/CMakeLists.txt
+++ b/tests/global_communication/CMakeLists.txt
@@ -4,6 +4,7 @@ set(HEADERS
set(COMMUNICATION_SOURCES
test_exporter_spike_file.cpp
test_communicator.cpp
+ test_mpi_gather_all.cpp
# unit test driver
test.cpp
diff --git a/tests/global_communication/test_mpi_gather_all.cpp b/tests/global_communication/test_mpi_gather_all.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..50fbf4052da14b057c6ed003d51cfa8b13436e93
--- /dev/null
+++ b/tests/global_communication/test_mpi_gather_all.cpp
@@ -0,0 +1,100 @@
+#include "gtest.h"
+
+#ifdef WITH_MPI
+
+#include <cstring>
+#include <vector>
+
+#include <communication/mpi_global_policy.hpp>
+#include <communication/mpi.hpp>
+#include <util/range.hpp>
+
+using namespace nest::mc;
+using namespace nest::mc::communication;
+
+struct big_thing {
+ big_thing() {}
+ big_thing(int i): value_(i) {}
+
+ bool operator==(const big_thing& other) const {
+ return value_==other.value_ && !std::memcmp(salt_, other.salt_, sizeof(salt_));
+ }
+
+ bool operator!=(const big_thing& other) const {
+ return !(*this==other);
+ }
+
+private:
+ int value_;
+ char salt_[32] = "it's a lovely day for a picnic";
+};
+
+TEST(mpi, gather_all) {
+ using policy = mpi_global_policy;
+
+ int id = policy::id();
+
+ std::vector<big_thing> data;
+ // odd ranks: three items; even ranks: one item.
+ if (id%2) {
+ data = { id, id+7, id+8 };
+ }
+ else {
+ data = { id };
+ }
+
+ std::vector<big_thing> expected;
+ for (int i = 0; i<policy::size(); ++i) {
+ if (i%2) {
+ int rank_data[] = { i, i+7, i+8 };
+ util::append(expected, rank_data);
+ }
+ else {
+ int rank_data[] = { i };
+ util::append(expected, rank_data);
+ }
+ }
+
+ auto gathered = mpi::gather_all(data);
+
+ EXPECT_EQ(expected, gathered);
+}
+
+TEST(mpi, gather_all_with_partition) {
+ using policy = mpi_global_policy;
+
+ int id = policy::id();
+
+ std::vector<big_thing> data;
+ // odd ranks: three items; even ranks: one item.
+ if (id%2) {
+ data = { id, id+7, id+8 };
+ }
+ else {
+ data = { id };
+ }
+
+ std::vector<big_thing> expected_values;
+ std::vector<unsigned> expected_divisions;
+
+ expected_divisions.push_back(0);
+ for (int i = 0; i<policy::size(); ++i) {
+ if (i%2) {
+ int rank_data[] = { i, i+7, i+8 };
+ util::append(expected_values, rank_data);
+ expected_divisions.push_back(expected_divisions.back()+util::size(rank_data));
+ }
+ else {
+ int rank_data[] = { i };
+ util::append(expected_values, rank_data);
+ expected_divisions.push_back(expected_divisions.back()+util::size(rank_data));
+ }
+ }
+
+ auto gathered = mpi::gather_all_with_partition(data);
+
+ EXPECT_EQ(expected_values, gathered.values());
+ EXPECT_EQ(expected_divisions, gathered.partition());
+}
+
+#endif // WITH_MPI
diff --git a/tests/test_common_cells.hpp b/tests/test_common_cells.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..c4641c879950fb6a5249cc25f38c3e6006df69ea
--- /dev/null
+++ b/tests/test_common_cells.hpp
@@ -0,0 +1,112 @@
+#include <cell.hpp>
+#include <parameter_list.hpp>
+
+namespace nest {
+namespace mc {
+
+/*
+ * Create cell with just a soma:
+ *
+ * Soma:
+ * diameter: 18.8 µm
+ * mechanisms: membrane, HH
+ * memrane resistance: 123 Ω·cm
+ *
+ * Stimuli:
+ * soma centre, t=[10 ms, 110 ms), 0.1 nA
+ */
+
+inline cell make_cell_soma_only() {
+ cell c;
+
+ auto soma = c.add_soma(18.8/2.0);
+ soma->mechanism("membrane").set("r_L", 123);
+ soma->add_mechanism(hh_parameters());
+
+ c.add_stimulus({0,0.5}, {10., 100., 0.1});
+
+ return c;
+}
+
+/*
+ * Create cell with a soma and unbranched dendrite:
+ *
+ * Soma:
+ * mechanisms: HH
+ * diameter: 12.6157 µm
+ *
+ * Dendrite:
+ * mechanisms: none
+ * diameter: 1 µm
+ * length: 200 µm
+ * membrane resistance: 100 Ω·cm
+ * compartments: 4
+ *
+ * Stimulus:
+ * end of dendrite, t=[5 ms, 85 ms), 0.3 nA
+ */
+
+inline cell make_cell_ball_and_stick() {
+ cell c;
+
+ auto soma = c.add_soma(12.6157/2.0);
+ soma->add_mechanism(hh_parameters());
+
+ auto dendrite = c.add_cable(0, segmentKind::dendrite, 1.0/2, 1.0/2, 200.0);
+ dendrite->add_mechanism(pas_parameters());
+ dendrite->mechanism("membrane").set("r_L", 100);
+ dendrite->set_compartments(4);
+
+ c.add_stimulus({1,1}, {5., 80., 0.3});
+
+ return c;
+}
+
+/*
+ * Create cell with a soma and three-segment dendrite with single branch point:
+ *
+ * O----======
+ *
+ * Soma:
+ * mechanisms: HH
+ * diameter: 12.6157 µm
+ *
+ * Dendrites:
+ * mechanisms: membrane
+ * diameter: 1 µm
+ * length: 100 µm
+ * membrane resistance: 100 Ω·cm
+ * compartments: 4
+ *
+ * Stimulus:
+ * end of first terminal branch, t=[5 ms, 85 ms), 0.45 nA
+ * end of second terminal branch, t=[40 ms, 50 ms), -0.2 nA
+ */
+
+inline cell make_cell_ball_and_3sticks() {
+ cell c;
+
+ auto soma = c.add_soma(12.6157/2.0);
+ soma->add_mechanism(hh_parameters());
+
+ // add dendrite of length 200 um and diameter 1 um with passive channel
+ c.add_cable(0, segmentKind::dendrite, 0.5, 0.5, 100);
+ c.add_cable(1, segmentKind::dendrite, 0.5, 0.5, 100);
+ c.add_cable(1, segmentKind::dendrite, 0.5, 0.5, 100);
+
+ for (auto& seg: c.segments()) {
+ if (seg->is_dendrite()) {
+ seg->add_mechanism(pas_parameters());
+ seg->mechanism("membrane").set("r_L", 100);
+ seg->set_compartments(4);
+ }
+ }
+
+ c.add_stimulus({2,1}, {5., 80., 0.45});
+ c.add_stimulus({3,1}, {40., 10.,-0.2});
+
+ return c;
+}
+
+} // namespace mc
+} // namespace nest
diff --git a/tests/unit/CMakeLists.txt b/tests/unit/CMakeLists.txt
index 6b34d70125bb4d9792f8abb6b68daf411dfe78d8..1c286863c3e92a9c8f2a11eba8679f66094b3c72 100644
--- a/tests/unit/CMakeLists.txt
+++ b/tests/unit/CMakeLists.txt
@@ -18,6 +18,7 @@ set(TEST_SOURCES
test_either.cpp
test_event_queue.cpp
test_fvm.cpp
+ test_fvm_multi.cpp
test_cell_group.cpp
test_lexcmp.cpp
test_mask_stream.cpp
diff --git a/tests/unit/test_cell_group.cpp b/tests/unit/test_cell_group.cpp
index 3aed222029f4535a85d9713849cff184a3684531..9b3f6e6a0b13f7338558889e40e47d7fe7ed99c8 100644
--- a/tests/unit/test_cell_group.cpp
+++ b/tests/unit/test_cell_group.cpp
@@ -1,27 +1,19 @@
#include "gtest.h"
+#include <cell_group.hpp>
#include <common_types.hpp>
#include <fvm_cell.hpp>
-#include <cell_group.hpp>
+#include <util/range.hpp>
+
+#include "../test_common_cells.hpp"
nest::mc::cell make_cell() {
using namespace nest::mc;
- nest::mc::cell cell;
-
- // Soma with diameter 12.6157 um and HH channel
- auto soma = cell.add_soma(12.6157/2.0);
- soma->add_mechanism(hh_parameters());
-
- // add dendrite of length 200 um and diameter 1 um with passive channel
- auto dendrite = cell.add_cable(0, segmentKind::dendrite, 0.5, 0.5, 200);
- dendrite->add_mechanism(pas_parameters());
- dendrite->set_compartments(101);
-
- dendrite->mechanism("membrane").set("r_L", 100);
+ nest::mc::cell cell = make_cell_ball_and_stick();
cell.add_detector({0, 0}, 0);
- cell.add_stimulus({1, 1}, {5., 80., 0.3});
+ cell.segment(1)->set_compartments(101);
return cell;
}
@@ -31,7 +23,7 @@ TEST(cell_group, test)
using namespace nest::mc;
using cell_group_type = cell_group<fvm::fvm_cell<double, cell_local_size_type>>;
- auto group = cell_group_type{0, make_cell()};
+ auto group = cell_group_type{0, util::singleton_view(make_cell())};
group.advance(50, 0.01);
@@ -53,7 +45,7 @@ TEST(cell_group, sources)
cell.add_detector({1, 0.3}, 2.3);
cell_gid_type first_gid = 37u;
- auto group = cell_group_type{first_gid, cell};
+ auto group = cell_group_type{first_gid, util::singleton_view(cell)};
// expect group sources to be lexicographically sorted by source id
// with gids in cell group's range and indices starting from zero
diff --git a/tests/unit/test_compartments.cpp b/tests/unit/test_compartments.cpp
index c167b16cce5585a43efc13fc9f136025b21f7811..8bc034ade553675ac8f85ce19bc8bfccf5d14d5d 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,98 +35,27 @@ TEST(compartments, compartment)
}
}
-TEST(compartments, compartment_iterator)
+TEST(compartments, make_compartment_range)
{
- // 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);
- }
+ using namespace nest::mc;
+ auto rng = make_compartment_range(10, 1.0, 2.0, 10.);
- // check postfix increment
+ EXPECT_EQ((*rng.begin()).index, 0u);
+ EXPECT_EQ((*rng.end()).index, 10u);
+ EXPECT_NE(rng.begin(), rng.end());
- // 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);
+ unsigned count = 0;
+ for (auto c : rng) {
+ EXPECT_EQ(c.index, count);
+ auto er = 1.0 + double(count)/10.;
+ EXPECT_DOUBLE_EQ(left(c.radius), er);
+ EXPECT_DOUBLE_EQ(right(c.radius), er+0.1);
+ 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());
- }
+ auto rng_empty = make_compartment_range(0, 1.0, 1.0, 0.);
+ EXPECT_EQ(rng_empty.begin(), rng_empty.end());
}
diff --git a/tests/unit/test_fvm.cpp b/tests/unit/test_fvm.cpp
index 39fdb6c837f71820c3fd6981f9be2bc38cee91c8..6d067182f972a0bc65147e8345f3043230d70d9a 100644
--- a/tests/unit/test_fvm.cpp
+++ b/tests/unit/test_fvm.cpp
@@ -5,42 +5,16 @@
#include <common_types.hpp>
#include <cell.hpp>
#include <fvm_cell.hpp>
-#include <util/singleton.hpp>
+#include <util/range.hpp>
+#include "../test_common_cells.hpp"
#include "../test_util.hpp"
TEST(fvm, cable)
{
using namespace nest::mc;
- nest::mc::cell cell;
-
- cell.add_soma(6e-4); // 6um in cm
-
- // 1um radius and 4mm long, all in cm
- cell.add_cable(0, segmentKind::dendrite, 1e-4, 1e-4, 4e-1);
- 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.model().tree.num_segments(), 3u);
-
- // add passive to all 3 segments in the cell
- for(auto& seg :cell.segments()) {
- seg->add_mechanism(pas_parameters());
- }
-
- cell.soma()->add_mechanism(hh_parameters());
- cell.segment(2)->add_mechanism(hh_parameters());
-
- auto& soma_hh = cell.soma()->mechanism("hh");
-
- soma_hh.set("gnabar", 0.12);
- soma_hh.set("gkbar", 0.036);
- soma_hh.set("gl", 0.0003);
- soma_hh.set("el", -54.387);
-
- cell.segment(1)->set_compartments(4);
- cell.segment(2)->set_compartments(4);
+ nest::mc::cell cell = make_cell_ball_and_3sticks();
using fvm_cell = fvm::fvm_cell<double, cell_lid_type>;
@@ -53,7 +27,8 @@ TEST(fvm, cable)
auto& J = fvcell.jacobian();
- EXPECT_EQ(cell.num_compartments(), 9u);
+ // 1 (soma) + 3 (dendritic segments) × 4 compartments
+ EXPECT_EQ(cell.num_compartments(), 13u);
// assert that the matrix has one row for each compartment
EXPECT_EQ(J.size(), cell.num_compartments());
@@ -69,21 +44,11 @@ TEST(fvm, init)
{
using namespace nest::mc;
- 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);
+ nest::mc::cell cell = make_cell_ball_and_stick();
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
- // - the only "state" is that used to define parameters
- cell.soma()->add_mechanism(hh_parameters());
-
auto& soma_hh = cell.soma()->mechanism("hh");
soma_hh.set("gnabar", 0.12);
@@ -112,4 +77,3 @@ TEST(fvm, init)
fvcell.setup_matrix(0.01);
}
-
diff --git a/tests/unit/test_fvm_multi.cpp b/tests/unit/test_fvm_multi.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..fbd4de99b91198128c6b3d0f0098f38f8440984d
--- /dev/null
+++ b/tests/unit/test_fvm_multi.cpp
@@ -0,0 +1,142 @@
+#include <fstream>
+
+#include "gtest.h"
+
+#include <common_types.hpp>
+#include <cell.hpp>
+#include <fvm_multicell.hpp>
+#include <util/range.hpp>
+
+#include "../test_util.hpp"
+#include "../test_common_cells.hpp"
+
+TEST(fvm_multi, cable)
+{
+ using namespace nest::mc;
+
+ nest::mc::cell cell=make_cell_ball_and_3sticks();
+
+ using fvm_cell = fvm::fvm_multicell<double, cell_lid_type>;
+
+ std::vector<fvm_cell::target_handle> targets;
+ std::vector<fvm_cell::detector_handle> detectors;
+ std::vector<fvm_cell::probe_handle> probes;
+
+ fvm_cell fvcell;
+ fvcell.initialize(util::singleton_view(cell), detectors, targets, probes);
+
+ auto& J = fvcell.jacobian();
+
+ // 1 (soma) + 3 (dendritic segments) × 4 compartments
+ EXPECT_EQ(cell.num_compartments(), 13u);
+
+ // assert that the matrix has one row for each compartment
+ EXPECT_EQ(J.size(), cell.num_compartments());
+
+ fvcell.setup_matrix(0.02);
+
+ // assert that the number of cv areas is the same as the matrix size
+ // i.e. both should equal the number of compartments
+ EXPECT_EQ(fvcell.cv_areas().size(), J.size());
+}
+
+TEST(fvm_multi, init)
+{
+ using namespace nest::mc;
+
+ nest::mc::cell cell = make_cell_ball_and_stick();
+
+ const auto m = cell.model();
+ EXPECT_EQ(m.tree.num_segments(), 2u);
+
+ auto& soma_hh = cell.soma()->mechanism("hh");
+
+ soma_hh.set("gnabar", 0.12);
+ soma_hh.set("gkbar", 0.036);
+ soma_hh.set("gl", 0.0003);
+ soma_hh.set("el", -54.3);
+
+ // check that parameter values were set correctly
+ EXPECT_EQ(cell.soma()->mechanism("hh").get("gnabar").value, 0.12);
+ EXPECT_EQ(cell.soma()->mechanism("hh").get("gkbar").value, 0.036);
+ 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(10);
+
+ using fvm_cell = fvm::fvm_multicell<double, cell_lid_type>;
+ std::vector<fvm_cell::target_handle> targets;
+ std::vector<fvm_cell::detector_handle> detectors;
+ std::vector<fvm_cell::probe_handle> probes;
+
+ fvm_cell fvcell;
+ fvcell.initialize(util::singleton_view(cell), detectors, targets, probes);
+
+ auto& J = fvcell.jacobian();
+ EXPECT_EQ(J.size(), 11u);
+
+ fvcell.setup_matrix(0.01);
+}
+
+TEST(fvm_multi, multi_init)
+{
+ using namespace nest::mc;
+
+ nest::mc::cell cells[] = {
+ make_cell_ball_and_stick(),
+ make_cell_ball_and_3sticks()
+ };
+
+ EXPECT_EQ(cells[0].num_segments(), 2u);
+ EXPECT_EQ(cells[0].segment(1)->num_compartments(), 4u);
+ EXPECT_EQ(cells[1].num_segments(), 4u);
+ EXPECT_EQ(cells[1].segment(1)->num_compartments(), 4u);
+ EXPECT_EQ(cells[1].segment(2)->num_compartments(), 4u);
+ EXPECT_EQ(cells[1].segment(3)->num_compartments(), 4u);
+
+ cells[0].add_synapse({1, 0.4}, parameter_list("expsyn"));
+ cells[0].add_synapse({1, 0.4}, parameter_list("expsyn"));
+ cells[1].add_synapse({2, 0.4}, parameter_list("exp2syn"));
+ cells[1].add_synapse({3, 0.4}, parameter_list("expsyn"));
+
+ cells[1].add_detector({0, 0}, 3.3);
+
+ using fvm_cell = fvm::fvm_multicell<double, cell_lid_type>;
+ std::vector<fvm_cell::target_handle> targets(4);
+ std::vector<fvm_cell::detector_handle> detectors(1);
+ std::vector<fvm_cell::probe_handle> probes;
+
+ fvm_cell fvcell;
+ fvcell.initialize(cells, detectors, targets, probes);
+
+ auto& J = fvcell.jacobian();
+ EXPECT_EQ(J.size(), 5u+13u);
+
+ // check indices in instantiated mechanisms
+ for (const auto& mech: fvcell.mechanisms()) {
+ if (mech->name()=="hh") {
+ // HH on somas of two cells, with group compartment indices
+ // 0 and 5.
+ ASSERT_EQ(mech->node_index().size(), 2u);
+ EXPECT_EQ(mech->node_index()[0], 0u);
+ EXPECT_EQ(mech->node_index()[1], 5u);
+ }
+ if (mech->name()=="expsyn") {
+ // Three expsyn synapses, two in second compartment
+ // of dendrite segment of first cell, one in second compartment
+ // of last segment of second cell.
+ ASSERT_EQ(mech->node_index().size(), 3u);
+ EXPECT_EQ(mech->node_index()[0], 2u);
+ EXPECT_EQ(mech->node_index()[1], 2u);
+ EXPECT_EQ(mech->node_index()[2], 15u);
+ }
+ if (mech->name()=="exp2syn") {
+ // One exp2syn synapse, in second compartment
+ // of penultimate segment of second cell.
+ ASSERT_EQ(mech->node_index().size(), 1u);
+ EXPECT_EQ(mech->node_index()[0], 11u);
+ }
+ }
+
+ fvcell.setup_matrix(0.01);
+}
diff --git a/tests/unit/test_probe.cpp b/tests/unit/test_probe.cpp
index 13a12ad6bd067ed91b37f2867f22b0f6405ba89a..cf7562c94fbf9ddcf33e16308d418dbdddd0ebf7 100644
--- a/tests/unit/test_probe.cpp
+++ b/tests/unit/test_probe.cpp
@@ -3,7 +3,7 @@
#include <common_types.hpp>
#include <cell.hpp>
#include <fvm_cell.hpp>
-#include <util/singleton.hpp>
+#include <util/range.hpp>
TEST(probe, instantiation)
{
diff --git a/tests/validation/validate_ball_and_stick.cpp b/tests/validation/validate_ball_and_stick.cpp
index 59e952de58de61ae32cdbfb721ac991ec2119f8b..2a5ec45ccd575dc5f12ad79f565996d584b41a22 100644
--- a/tests/validation/validate_ball_and_stick.cpp
+++ b/tests/validation/validate_ball_and_stick.cpp
@@ -3,10 +3,13 @@
#include <common_types.hpp>
#include <cell.hpp>
-#include <fvm_cell.hpp>
-#include <util/singleton.hpp>
+//#include <fvm_cell.hpp>
+#include <fvm_multicell.hpp>
+#include <util/range.hpp>
#include "gtest.h"
+
+#include "../test_common_cells.hpp"
#include "../test_util.hpp"
#include "validation_data.hpp"
@@ -16,20 +19,7 @@ TEST(ball_and_stick, neuron_baseline)
using namespace nest::mc;
using namespace nlohmann;
- nest::mc::cell cell;
-
- // Soma with diameter 12.6157 um and HH channel
- auto soma = cell.add_soma(12.6157/2.0);
- soma->add_mechanism(hh_parameters());
-
- // add dendrite of length 200 um and diameter 1 um with passive channel
- auto dendrite = cell.add_cable(0, segmentKind::dendrite, 0.5, 0.5, 200);
- dendrite->add_mechanism(pas_parameters());
-
- dendrite->mechanism("membrane").set("r_L", 100);
-
- // add stimulus
- cell.add_stimulus({1,1}, {5., 80., 0.3});
+ nest::mc::cell cell = make_cell_ball_and_stick();
// load data from file
auto cell_data = testing::g_validation_data.load("ball_and_stick.json");
@@ -83,16 +73,16 @@ TEST(ball_and_stick, neuron_baseline)
}
};
- using fvm_cell = fvm::fvm_cell<double, cell_local_size_type>;
+ using fvm_cell = fvm::fvm_multicell<double, cell_local_size_type>;
std::vector<fvm_cell::detector_handle> detectors(cell.detectors().size());
std::vector<fvm_cell::target_handle> targets(cell.synapses().size());
std::vector<fvm_cell::probe_handle> probes(cell.probes().size());
std::vector<result> results;
- for(auto run_index=0u; run_index<cell_data.size(); ++run_index) {
+ for (auto run_index=0u; run_index<cell_data.size(); ++run_index) {
auto& run = cell_data[run_index];
int num_compartments = run["nseg"];
- dendrite->set_compartments(num_compartments);
+ cell.segment(1)->set_compartments(num_compartments);
std::vector<std::vector<double>> v(3);
// make the lowered finite volume cell
@@ -164,26 +154,7 @@ TEST(ball_and_3stick, neuron_baseline)
using namespace nest::mc;
using namespace nlohmann;
- nest::mc::cell cell;
-
- // Soma with diameter 12.6157 um and HH channel
- auto soma = cell.add_soma(12.6157/2.0);
- soma->add_mechanism(hh_parameters());
-
- // add dendrite of length 200 um and diameter 1 um with passive channel
- std::vector<cable_segment*> dendrites;
- dendrites.push_back(cell.add_cable(0, segmentKind::dendrite, 0.5, 0.5, 100));
- dendrites.push_back(cell.add_cable(1, segmentKind::dendrite, 0.5, 0.5, 100));
- dendrites.push_back(cell.add_cable(1, segmentKind::dendrite, 0.5, 0.5, 100));
-
- for(auto dend : dendrites) {
- dend->add_mechanism(pas_parameters());
- dend->mechanism("membrane").set("r_L", 100);
- }
-
- // add stimulus
- cell.add_stimulus({2,1}, {5., 80., 0.45});
- cell.add_stimulus({3,1}, {40., 10.,-0.2});
+ nest::mc::cell cell = make_cell_ball_and_3sticks();
// load data from file
auto cell_data = testing::g_validation_data.load("ball_and_3stick.json");
@@ -237,7 +208,7 @@ TEST(ball_and_3stick, neuron_baseline)
}
};
- using fvm_cell = fvm::fvm_cell<double, cell_local_size_type>;
+ using fvm_cell = fvm::fvm_multicell<double, cell_local_size_type>;
std::vector<fvm_cell::detector_handle> detectors(cell.detectors().size());
std::vector<fvm_cell::target_handle> targets(cell.synapses().size());
std::vector<fvm_cell::probe_handle> probes(cell.probes().size());
@@ -247,8 +218,10 @@ TEST(ball_and_3stick, neuron_baseline)
for(auto run_index=0u; run_index<cell_data.size(); ++run_index) {
auto& run = cell_data[run_index];
int num_compartments = run["nseg"];
- for(auto dend : dendrites) {
- dend->set_compartments(num_compartments);
+ for (auto& seg: cell.segments()) {
+ if (seg->is_dendrite()) {
+ seg->set_compartments(num_compartments);
+ }
}
std::vector<std::vector<double>> v(3);
diff --git a/tests/validation/validate_soma.cpp b/tests/validation/validate_soma.cpp
index af4d32826acb12bcc522daf7aa3f84c80dcc6bb2..8952ba30c93d5c7121fa09c6092e281d39564eb4 100644
--- a/tests/validation/validate_soma.cpp
+++ b/tests/validation/validate_soma.cpp
@@ -4,10 +4,12 @@
#include <common_types.hpp>
#include <cell.hpp>
#include <fvm_cell.hpp>
-#include <util/singleton.hpp>
+#include <util/range.hpp>
#include "gtest.h"
+
#include "../test_util.hpp"
+#include "../test_common_cells.hpp"
#include "validation_data.hpp"
// compares results with those generated by nrn/soma.py
@@ -17,15 +19,7 @@ TEST(soma, neuron_baseline)
using namespace nest::mc;
using namespace nlohmann;
- nest::mc::cell cell;
-
- // Soma with diameter 18.8um and HH channel
- auto soma = cell.add_soma(18.8/2.0);
- soma->mechanism("membrane").set("r_L", 123); // no effect for single compartment cell
- soma->add_mechanism(hh_parameters());
-
- // add stimulus to the soma
- cell.add_stimulus({0,0.5}, {10., 100., 0.1});
+ nest::mc::cell cell = make_cell_soma_only();
// make the lowered finite volume cell
using fvm_cell = fvm::fvm_cell<double, cell_local_size_type>;
@@ -86,15 +80,7 @@ TEST(soma, convergence)
{
using namespace nest::mc;
- nest::mc::cell cell;
-
- // Soma with diameter 18.8um and HH channel
- auto soma = cell.add_soma(18.8/2.0);
- soma->mechanism("membrane").set("r_L", 123); // no effect for single compartment cell
- soma->add_mechanism(hh_parameters());
-
- // add stimulus to the soma
- cell.add_stimulus({0,0.5}, {10., 100., 0.1});
+ nest::mc::cell cell = make_cell_soma_only();
// make the lowered finite volume cell
using fvm_cell = fvm::fvm_cell<double, cell_local_size_type>;
diff --git a/tests/validation/validate_synapses.cpp b/tests/validation/validate_synapses.cpp
index 9aa90abef6f6178a6285f9a1d5db2e9d1cb1d4f1..d1de3df06945b33466ccbbc105d412196c0a5f0b 100644
--- a/tests/validation/validate_synapses.cpp
+++ b/tests/validation/validate_synapses.cpp
@@ -8,6 +8,7 @@
#include <cell_group.hpp>
#include <fvm_cell.hpp>
#include <mechanism_interface.hpp>
+#include <util/range.hpp>
#include "gtest.h"
#include "../test_util.hpp"
@@ -108,7 +109,7 @@ void run_neuron_baseline(const char* syn_type, const char* data_file)
std::vector<std::vector<double>> v(2);
// make the lowered finite volume cell
- cell_group<lowered_cell> group(0, cell);
+ cell_group<lowered_cell> group(0, util::singleton_view(cell));
// add the 3 spike events to the queue
group.enqueue_events(synthetic_events);