diff --git a/docs/SelfArx.cls b/docs/model/SelfArx.cls
similarity index 100%
rename from docs/SelfArx.cls
rename to docs/model/SelfArx.cls
diff --git a/docs/appendix.tex b/docs/model/appendix.tex
similarity index 100%
rename from docs/appendix.tex
rename to docs/model/appendix.tex
diff --git a/docs/bibliography.bib b/docs/model/bibliography.bib
similarity index 100%
rename from docs/bibliography.bib
rename to docs/model/bibliography.bib
diff --git a/docs/formulation.tex b/docs/model/formulation.tex
similarity index 100%
rename from docs/formulation.tex
rename to docs/model/formulation.tex
diff --git a/docs/images/cable.tex b/docs/model/images/cable.tex
similarity index 100%
rename from docs/images/cable.tex
rename to docs/model/images/cable.tex
diff --git a/docs/images/makefile b/docs/model/images/makefile
similarity index 100%
rename from docs/images/makefile
rename to docs/model/images/makefile
diff --git a/docs/images/soma.tex b/docs/model/images/soma.tex
similarity index 100%
rename from docs/images/soma.tex
rename to docs/model/images/soma.tex
diff --git a/docs/makefile b/docs/model/makefile
similarity index 100%
rename from docs/makefile
rename to docs/model/makefile
diff --git a/docs/report.bbl b/docs/model/report.bbl
similarity index 100%
rename from docs/report.bbl
rename to docs/model/report.bbl
diff --git a/docs/report.tex b/docs/model/report.tex
similarity index 100%
rename from docs/report.tex
rename to docs/model/report.tex
diff --git a/docs/symbols.tex b/docs/model/symbols.tex
similarity index 100%
rename from docs/symbols.tex
rename to docs/model/symbols.tex
diff --git a/docs/pid/makefile b/docs/pid/makefile
new file mode 100644
index 0000000000000000000000000000000000000000..a529bb95b9e7e7935ab393fd34201de25f2b09bc
--- /dev/null
+++ b/docs/pid/makefile
@@ -0,0 +1,18 @@
+TEXFILES=*.tex
+
+report.pdf : $(TEXFILES)
+ pdflatex --shell-escape report
+
+force : report.pdf
+ pdflatex --shell-escape report
+
+open : report.pdf
+ open report.pdf
+
+clean:
+ rm -f report.pdf
+ rm -f *.aux
+ rm -f *.log
+ rm -f *.toc
+ rm -f *.out
+
diff --git a/docs/pid/report.tex b/docs/pid/report.tex
new file mode 100644
index 0000000000000000000000000000000000000000..bcf2b31e0fb227a9017cd576610995be16b97bb2
--- /dev/null
+++ b/docs/pid/report.tex
@@ -0,0 +1,196 @@
+\documentclass[11pt,a4paper]{article}
+
+\usepackage[colorlinks=true,linkcolor=blue!30!black]{hyperref}
+\usepackage[pdftex]{graphicx}
+\usepackage[cmex10]{amsmath}
+\usepackage[usenames,dvipsnames]{color}
+\usepackage{xspace}
+\usepackage{fullpage}
+\usepackage{tikz}
+\usetikzlibrary{positioning}
+\usetikzlibrary{shapes,arrows,backgrounds,fit,shapes.geometric,calc}
+\usetikzlibrary{pgfplots.groupplots}
+\usepackage{pgfplots}
+\usepackage{pgfplotstable}
+\usepackage{textcomp}
+\usepackage{multirow}
+\usepackage{url}
+\usepackage{listings}
+\usepackage{lstautogobble}
+\usepackage{framed}
+\usepackage{tcolorbox}
+\usepackage{enumitem}
+
+
+\newcommand{\hilight}[1]{\textit{\textcolor{Red}{#1}}}
+\newcommand{\julich}[0]{J\"ulich\xspace}
+\newcommand{\nestmc}[0]{Nest$\mathcal{MC}$\xspace}
+
+% redefine emph for more... emphasis
+\DeclareTextFontCommand{\emph}{\bfseries}
+
+\begin{document}
+
+% title and cover page
+\title{HBP WP 7.5.4 Project Initiation Document}
+\date{\today}
+\maketitle
+
+%-------------------------------------------------------------
+\section{Project Title}
+%-------------------------------------------------------------
+Application Services Support.
+
+%-------------------------------------------------------------
+\section{Project Background}
+%-------------------------------------------------------------
+
+Current multi-compartment (MC) simulators don't target or effectively utilize key contemporary HPC architectures.
+HPC resources that are a key part of the HBP platform will be hetergoneous, and some key systems will be very different to the systems currently supported by simulation software
+\begin{itemize}
+ \item The main HBP system to be installed at \julich will be either Power8+GPU or KNL based system
+ \item The other HPC resources available include GPU clusters, BG/Q and Intel x86 multicore systems.
+\end{itemize}
+
+There is an effort to develop a HPC-specific version of the widely-used NEURON simulator that supports GPU and Intel Xeon Phi processors, however:
+\begin{itemize}
+ \item NEURON is large, complicated and \hilight{poorly engineered}. There is a strong argument that it is not possible to refactor the code effectively to provide both support for new hardware and a good platform for developing new algorithms.
+ \item The developers are unwilling to share their work or progress, which makes it impossible to understand the risk of relying on their efforts.
+\end{itemize}
+
+
+%-------------------------------------------------------------
+\section{Project Objectives}
+%-------------------------------------------------------------
+
+\begin{itemize}
+ \item Facilitate multicompartment simulation ``at scale'' on all HPC systems available to researchers in the HBP.
+ \item Performance portability: high performance of the delivered software on all target systems, and a forward-looking design that will be adaptable to future architectures.
+ \item Promote open source solution with close collaboration inside and outside the HBP (i.e. the broader community).
+ \item Improve software engineering practice in HBP.
+\end{itemize}
+
+%-------------------------------------------------------------
+\section{Project Deliverables}
+%-------------------------------------------------------------
+The project will be considered a success if all of the following are delivered in full and on time.
+\hilight{might be a bit restrictive}.
+
+\subsection{Open source software product \nestmc}
+
+\begin{itemize}
+ \item Performance portable on all potential HBP systems and main systems at CSCS and \julich
+ \begin{itemize}
+ \item multi-core Intel x86
+ \item Hybrid NVIDIA GPU based nodes (both Intel x86 and IBM Power host processors)
+ \item Intel KNL
+ \end{itemize}
+ \item Source available openly on github (or similar platform). \hilight{not like CoreNeuron}
+ \item High quality documentation.
+ \item Comprehensive validation and unit testing framework.
+ \item Automated testing.
+\end{itemize}
+
+\subsection{Scientific Use Case}
+
+\begin{itemize}
+ \item Developed by a scientific collaborator.
+ \item Running at scale on HPC resources.
+\end{itemize}
+
+\subsection{Publication}
+
+\begin{itemize}
+ \item Paper for the SC18 conference
+\end{itemize}
+
+%-------------------------------------------------------------
+\section{Project Scope}
+%-------------------------------------------------------------
+\subsection{This Project \underline{Will} Include}
+
+\begin{itemize}
+ \item Development of the main software product \nestmc.
+ \item Work on adapting and improving existing algorithms and mathematical techniques in multicompartmental modelling.
+ \item Low-overhead promotion of product at conferences, workshops and hackathons.
+ \item Assistance of scientific collaborator to develop a scientific use case.
+ \item Writing one paper.
+\end{itemize}
+
+\subsection{This Project \underline{Will Not} Include}
+
+\begin{itemize}
+ \item Development of the scientific use case.
+ \begin{itemize}
+ \item we provide support for an external scientific collaborator to develop the use case.
+ \end{itemize}
+ \item Development of new methods (mathematical and algorithmic).
+ \begin{itemize}
+ \item we provide a flexible software framework that can be used in follow on projects to test new methods.
+ \end{itemize}
+\end{itemize}
+
+%-------------------------------------------------------------
+\section{Constraints}
+%-------------------------------------------------------------
+
+The key constraint for this project is \emph{scope}.
+The project is well resourced in terms of hardware and software, however manpower is very limited and the two year time frame does not give us license to experiment with radical new solutions.
+Given this, we aim to build a high-quality product with a modest feature set as the basis for future work.
+
+%-------------------------------------------------------------
+\section{Key Assumptions}
+%-------------------------------------------------------------
+
+\begin{itemize}
+ \item \emph{A scientific collaborator will be found}.
+ A collaborator is essential to the project:
+ \begin{itemize}
+ \item To gain credibility of the product.
+ \item To provide motivation for continuing the development of the software after the initial 2 year project completion.
+ \item To give direct user input into the development and selection of features, so as to deliver a useful product.
+ \end{itemize}
+ However the project was started without one. Finding a collaborator inside HBP might be challenging politically, so we are also looking outside the HBP in the context of the Brain initiative.
+ \item \emph{The NEST Initiative will endorse NestMC}. This is not a given, particularly because the NEST Initiative was not consulted widely before starting the project. A key step will be communicating our intentions, and gaining the support of the NEST Initiative. \hilight{motivate}
+ \item \emph{A user community capable of taking over development will be in place}.
+ Of the assumptions this is the weakest link.
+ It is possible to develop the software and a scientific use case with the given time and resources, however building a community that is capable and willing to carry on development is quite unlikely.
+ Plans should be made for evaluating the viability of the product toward the end of SGA1 to decide whether to continue funded development in subsequent phases of the HBP.
+ \hilight{motivate}
+\end{itemize}
+
+%-------------------------------------------------------------
+\section{Project Manager}
+%-------------------------------------------------------------
+Ben Cumming at CSCS.
+
+%-------------------------------------------------------------
+\section{Project Team}
+%-------------------------------------------------------------
+
+The project team is small.
+There will potentially be contributions from a large group of \emph{student contributors} (Masters/PhD, postdoctoral and internship).
+However, to effectively plan for concrete deliverables we only consider staff that are guarenteed to be devoted to the project.
+At best we can rely on two full time developers contributing code that passes code review to the software product.
+
+\begin{itemize}
+ \item \emph{Alex Peyser} (\julich) XX\%: community liason
+ \item \emph{Ben Cumming} (CSCS) 80\%: project management and developer
+ \item \emph{Sam Yates} (CSCS) 100\%: lead developer
+ \item \emph{Wouter Klijn} (\julich) XX\%: developer
+ \item Student contributions.
+\end{itemize}
+
+%-------------------------------------------------------------
+\section{Duration}
+%-------------------------------------------------------------
+Project duration is the two year period from April 1 2016 until April 1 2018.
+
+%-------------------------------------------------------------
+\section{Approval From Sponsor}
+%-------------------------------------------------------------
+\vspace{20pt}
+Date :\hspace{5cm} Signature :
+
+\end{document}
+
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/profstats b/scripts/profstats
new file mode 100755
index 0000000000000000000000000000000000000000..33f37356836104e57758f0407e1da5fd0625a111
--- /dev/null
+++ b/scripts/profstats
@@ -0,0 +1,93 @@
+#!env python
+#coding: utf-8
+
+import json
+import argparse
+import re
+import numpy as np
+from itertools import chain
+
+def parse_clargs():
+ P = argparse.ArgumentParser(description='Aggregate and analyse MPI profile output.')
+ P.add_argument('inputs', metavar='FILE', nargs='+',
+ help='MPI profile output in JSON format')
+ P.add_argument('-r', '--raw', action='store_true',
+ help='emit raw times in csv table')
+
+ return P.parse_args()
+
+def parse_profile_json(source):
+ j = json.load(source)
+ rank = j['rank']
+ if rank is None:
+ raise ValueError('missing rank information in profile')
+
+ tx = dict()
+
+ def collect_times(j, prefix):
+ t = j['time']
+ n = j['name']
+
+ if t is None or n is None:
+ return
+
+ prefix = prefix + n
+ tx[prefix] = t
+
+ try:
+ children = j['regions']
+ # special case for top level
+ if prefix == 'total':
+ prefix = ''
+ else:
+ prefix = prefix + '/'
+
+ for j in children:
+ collect_times(j, prefix)
+ except KeyError:
+ pass
+
+ collect_times(j['regions'], '')
+ return rank, tx
+
+def csv_escape(x):
+ s = re.sub('"','""',str(x))
+ if re.search('["\t\n,]',s):
+ s = '"'+s+'"'
+ return s
+
+def emit_csv(cols, rows):
+ print(",".join([csv_escape(c) for c in cols]))
+ for r in rows:
+ print(",".join([csv_escape(r[c]) if c in r else '' for c in cols]))
+
+args = parse_clargs()
+
+rank_times = dict()
+for filename in args.inputs:
+ with open(filename) as f:
+ rank, times = parse_profile_json(f)
+ rank_times[rank] = times
+
+if args.raw:
+ rows = [rank_times[rank] for rank in sorted(rank_times.keys())]
+ cols = sorted({col for tbl in rows for col in tbl.keys()})
+ emit_csv(cols, rows)
+else:
+ rank_entry = [rank_times[rank] for rank in sorted(rank_times.keys())]
+ bins = sorted({col for tbl in rank_entry for col in tbl.keys()})
+
+ rows = []
+ for b in bins:
+ qs = np.percentile([entry[b] for entry in rank_times.values() if b in entry],
+ [0., 0.25, 0.5, 0.75, 1.])
+ rows.append({
+ 'region': b,
+ 'min': qs[0],
+ 'q25': qs[1],
+ 'median': qs[2],
+ 'q75': qs[3],
+ 'max': qs[4]
+ })
+
+ emit_csv(['region','min','q25','median','q75','max'], rows)
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/gathered_vector.hpp b/src/communication/gathered_vector.hpp
index cf48192f628b665713a25b2b1a53f0df3af86936..ceb5c16cde44a872e9ce08d3b4c65c58f9f73c1a 100644
--- a/src/communication/gathered_vector.hpp
+++ b/src/communication/gathered_vector.hpp
@@ -20,7 +20,7 @@ public:
partition_(std::move(p))
{
EXPECTS(std::is_sorted(partition_.begin(), partition_.end()));
- EXPECTS(std::size_t(partition_.back()) == v.size());
+ EXPECTS(std::size_t(partition_.back()) == values_.size());
}
/// the partition of distribution
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/communication/serial_global_policy.hpp b/src/communication/serial_global_policy.hpp
index 73d7b96b378da6880ea337f702cc13995d2f05ec..486266cadcacde9b9db4c7393efad9540a5c74f8 100644
--- a/src/communication/serial_global_policy.hpp
+++ b/src/communication/serial_global_policy.hpp
@@ -15,7 +15,11 @@ struct serial_global_policy {
template <typename Spike>
static gathered_vector<Spike>
gather_spikes(const std::vector<Spike>& local_spikes) {
- return gathered_vector<Spike>(std::vector<Spike>(local_spikes), {0u, 1u});
+ using count_type = typename gathered_vector<Spike>::count_type;
+ return gathered_vector<Spike>(
+ std::vector<Spike>(local_spikes),
+ {0u, static_cast<count_type>(local_spikes.size())}
+ );
}
static int id() {
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/iterutil.hpp b/src/util/iterutil.hpp
index 6a9fd159ffe9f75eddf57f2adb5c4d42cb9e7861..685398579a14c32971911087d772182d17b3d12d 100644
--- a/src/util/iterutil.hpp
+++ b/src/util/iterutil.hpp
@@ -43,6 +43,30 @@ upto(I iter, E end) {
return iter==I{end}? iter: I{std::prev(end)};
}
+template <typename I, typename E,
+ typename C = common_random_access_iterator_t<I,E>>
+enable_if_t<std::is_same<I, E>::value ||
+ (has_common_random_access_iterator<I,E>::value &&
+ is_forward_iterator<I>::value),
+ typename std::iterator_traits<C>::difference_type>
+distance(I first, E last) {
+ return std::distance(static_cast<C>(first), static_cast<C>(last));
+}
+
+template <typename I, typename E>
+enable_if_t<!has_common_random_access_iterator<I, E>::value &&
+ is_forward_iterator<I>::value,
+ typename std::iterator_traits<I>::difference_type>
+distance(I first, E last) {
+ typename std::iterator_traits<I>::difference_type ret = 0;
+ while (first != last) {
+ ++first;
+ ++ret;
+ }
+
+ return ret;
+}
+
/*
* Provide a proxy object for operator->() for iterator adaptors that
* present rvalues on dereference.
diff --git a/src/util/meta.hpp b/src/util/meta.hpp
index 40f08d481f23199e8e06a47fc41acdd58fdab8c9..e3d7fc9566a221425bc0d326debaa412a50b16bc 100644
--- a/src/util/meta.hpp
+++ b/src/util/meta.hpp
@@ -151,6 +151,45 @@ struct is_forward_iterator<T, enable_if_t<
template <typename T>
using is_forward_iterator_t = typename is_forward_iterator<T>::type;
+
+template <typename I, typename E, typename = void, typename = void>
+struct common_random_access_iterator {};
+
+template <typename I, typename E>
+struct common_random_access_iterator<
+ I,
+ E,
+ void_t<decltype(false ? std::declval<I>() : std::declval<E>())>,
+ enable_if_t<
+ is_random_access_iterator<
+ decay_t<decltype(false ? std::declval<I>() : std::declval<E>())>
+ >::value
+ >
+> {
+ using type = decay_t<
+ decltype(false ? std::declval<I>() : std::declval<E>())
+ >;
+};
+
+template <typename I, typename E>
+using common_random_access_iterator_t = typename common_random_access_iterator<I, E>::type;
+
+// NB: using this version of SFINAE instead of a void default template
+// parameter because of a bug in gcc 4.9.3.
+template <typename I, typename E>
+struct has_common_random_access_iterator {
+private:
+ constexpr static bool test(...) { return false; }
+
+ template <typename U = I>
+ constexpr static bool test(typename common_random_access_iterator<U, E>::type*) { return true; }
+
+public:
+ constexpr static bool value = test(nullptr);
+};
+
+
+
} // namespace util
} // namespace mc
} // namespace nest
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 2439ca74ee785fac1f909e92b3189bbfcd6f5da2..df5625d0e6dc1e315e18d1e2e1a0e4b4ccecbd94 100644
--- a/src/util/range.hpp
+++ b/src/util/range.hpp
@@ -35,6 +35,7 @@
#include <util/either.hpp>
#include <util/iterutil.hpp>
#include <util/meta.hpp>
+#include <util/sentinel.hpp>
namespace nest {
namespace mc {
@@ -66,7 +67,7 @@ struct range {
range& operator=(const range&) = default;
range& operator=(range&&) = default;
- bool empty() const { return left==right; }
+ bool empty() const { return left == right; }
iterator begin() const { return left; }
const_iterator cbegin() const { return left; }
@@ -77,10 +78,12 @@ struct range {
template <typename V = iterator>
enable_if_t<is_forward_iterator<V>::value, size_type>
size() const {
- return std::distance(begin(), end());
+ return util::distance(begin(), end());
}
- constexpr size_type max_size() const { return std::numeric_limits<size_type>::max(); }
+ constexpr size_type max_size() const {
+ return std::numeric_limits<size_type>::max();
+ }
void swap(range& other) {
std::swap(left, other.left);
@@ -148,164 +151,6 @@ range<U, V> make_range(const U& left, const V& right) {
return range<U, V>(left, right);
}
-/*
- * Use a proxy iterator to present a range as having the same begin and
- * end types, for use with e.g. pre-C++17 ranged-for loops or STL
- * algorithms.
- */
-template <typename I, typename S>
-class sentinel_iterator {
- nest::mc::util::either<I, S> e_;
-
- bool is_sentinel() const { return e_.index()!=0; }
-
- I& iter() {
- EXPECTS(!is_sentinel());
- return e_.template unsafe_get<0>();
- }
-
- const I& iter() const {
- EXPECTS(!is_sentinel());
- return e_.template unsafe_get<0>();
- }
-
- S& sentinel() {
- EXPECTS(is_sentinel());
- return e_.template unsafe_get<1>();
- }
-
- const S& sentinel() const {
- EXPECTS(is_sentinel());
- return e_.template unsafe_get<1>();
- }
-
-public:
- using difference_type = typename std::iterator_traits<I>::difference_type;
- using value_type = typename std::iterator_traits<I>::value_type;
- using pointer = typename std::iterator_traits<I>::pointer;
- using reference = typename std::iterator_traits<I>::reference;
- using iterator_category = typename std::iterator_traits<I>::iterator_category;
-
- sentinel_iterator(I i): e_(i) {}
-
- template <typename V = S, typename = enable_if_t<!std::is_same<I, V>::value>>
- sentinel_iterator(S i): e_(i) {}
-
- sentinel_iterator() = default;
- sentinel_iterator(const sentinel_iterator&) = default;
- sentinel_iterator(sentinel_iterator&&) = default;
-
- sentinel_iterator& operator=(const sentinel_iterator&) = default;
- sentinel_iterator& operator=(sentinel_iterator&&) = default;
-
- // forward and input iterator requirements
-
- auto operator*() const -> decltype(*(this->iter())) { return *iter(); }
-
- I operator->() const { return e_.template ptr<0>(); }
-
- sentinel_iterator& operator++() {
- ++iter();
- return *this;
- }
-
- sentinel_iterator operator++(int) {
- sentinel_iterator c(*this);
- ++*this;
- return c;
- }
-
- bool operator==(const sentinel_iterator& x) const {
- if (is_sentinel()) {
- return x.is_sentinel() || x.iter()==sentinel();
- }
- else {
- return x.is_sentinel()? iter()==x.sentinel(): iter()==x.iter();
- }
- }
-
- bool operator!=(const sentinel_iterator& x) const {
- return !(*this==x);
- }
-
- // bidirectional iterator requirements
-
- sentinel_iterator& operator--() {
- --iter();
- return *this;
- }
-
- sentinel_iterator operator--(int) {
- sentinel_iterator c(*this);
- --*this;
- return c;
- }
-
- // random access iterator requirements
-
- sentinel_iterator &operator+=(difference_type n) {
- iter() += n;
- return *this;
- }
-
- sentinel_iterator operator+(difference_type n) const {
- sentinel_iterator c(*this);
- return c += n;
- }
-
- friend sentinel_iterator operator+(difference_type n, sentinel_iterator x) {
- return x+n;
- }
-
- sentinel_iterator& operator-=(difference_type n) {
- iter() -= n;
- return *this;
- }
-
- sentinel_iterator operator-(difference_type n) const {
- sentinel_iterator c(*this);
- return c -= n;
- }
-
- difference_type operator-(sentinel_iterator x) const {
- return iter()-x.iter();
- }
-
- auto operator[](difference_type n) const -> decltype(*(this->iter())) {
- return *(iter()+n);
- }
-
- bool operator<=(const sentinel_iterator& x) const {
- return x.is_sentinel() || (!is_sentinel() && iter()<=x.iter());
- }
-
- bool operator<(const sentinel_iterator& x) const {
- return !is_sentinel() && (x.is_sentinel() || iter()<=x.iter());
- }
-
- bool operator>=(const sentinel_iterator& x) const {
- return !(x<*this);
- }
-
- bool operator>(const sentinel_iterator& x) const {
- return !(x<=*this);
- }
-};
-
-template <typename I, typename S>
-using sentinel_iterator_t =
- typename std::conditional<std::is_same<I, S>::value, I, sentinel_iterator<I, S>>::type;
-
-template <typename I, typename S>
-sentinel_iterator_t<I, S> make_sentinel_iterator(const I& i, const S& s) {
- return sentinel_iterator_t<I, S>(i);
-}
-
-template <typename I, typename S>
-sentinel_iterator_t<I, S> make_sentinel_end(const I& i, const S& s) {
- return sentinel_iterator_t<I, S>(s);
-}
-
template <typename Seq>
auto canonical_view(const Seq& s) ->
range<sentinel_iterator_t<decltype(std::begin(s)), decltype(std::end(s))>>
@@ -327,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/sentinel.hpp b/src/util/sentinel.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..8d08ecd41fe77f9c11ed08b8b539536792459c36
--- /dev/null
+++ b/src/util/sentinel.hpp
@@ -0,0 +1,184 @@
+#pragma once
+
+#include <type_traits>
+
+#include <util/meta.hpp>
+
+/*
+ * Use a proxy iterator to present a range as having the same begin and
+ * end types, for use with e.g. pre-C++17 ranged-for loops or STL
+ * algorithms.
+ */
+
+namespace nest {
+namespace mc {
+namespace util {
+
+template<typename I, typename S, bool Same>
+struct iterator_category_select {
+ // default category
+ using iterator_category = std::forward_iterator_tag;
+};
+
+template<typename I, typename S>
+struct iterator_category_select<I,S,true> {
+ using iterator_category = typename std::iterator_traits<I>::iterator_category;
+};
+
+template <typename I, typename S>
+class sentinel_iterator {
+ nest::mc::util::either<I, S> e_;
+
+ bool is_sentinel() const { return e_.index()!=0; }
+
+ I& iter() {
+ EXPECTS(!is_sentinel());
+ return e_.template unsafe_get<0>();
+ }
+
+ const I& iter() const {
+ EXPECTS(!is_sentinel());
+ return e_.template unsafe_get<0>();
+ }
+
+ S& sentinel() {
+ EXPECTS(is_sentinel());
+ return e_.template unsafe_get<1>();
+ }
+
+ const S& sentinel() const {
+ EXPECTS(is_sentinel());
+ return e_.template unsafe_get<1>();
+ }
+
+public:
+ using difference_type = typename std::iterator_traits<I>::difference_type;
+ using value_type = typename std::iterator_traits<I>::value_type;
+ using pointer = typename std::iterator_traits<I>::pointer;
+ using reference = typename std::iterator_traits<I>::reference;
+ using iterator_category = typename iterator_category_select<
+ I,S,std::is_same<I,S>::value>::iterator_category;
+
+ sentinel_iterator(I i): e_(i) {}
+
+ template <typename V = S, typename = enable_if_t<!std::is_same<I, V>::value>>
+ sentinel_iterator(S i): e_(i) {}
+
+ sentinel_iterator() = default;
+ sentinel_iterator(const sentinel_iterator&) = default;
+ sentinel_iterator(sentinel_iterator&&) = default;
+
+ sentinel_iterator& operator=(const sentinel_iterator&) = default;
+ sentinel_iterator& operator=(sentinel_iterator&&) = default;
+
+ // forward and input iterator requirements
+
+ auto operator*() const -> decltype(*(this->iter())) { return *iter(); }
+
+ I operator->() const { return e_.template ptr<0>(); }
+
+ sentinel_iterator& operator++() {
+ ++iter();
+ return *this;
+ }
+
+ sentinel_iterator operator++(int) {
+ sentinel_iterator c(*this);
+ ++*this;
+ return c;
+ }
+
+ bool operator==(const sentinel_iterator& x) const {
+ if (is_sentinel()) {
+ return x.is_sentinel() || x.iter()==sentinel();
+ }
+ else {
+ return x.is_sentinel()? iter()==x.sentinel(): iter()==x.iter();
+ }
+ }
+
+ bool operator!=(const sentinel_iterator& x) const {
+ return !(*this==x);
+ }
+
+ // bidirectional iterator requirements
+
+ sentinel_iterator& operator--() {
+ --iter();
+ return *this;
+ }
+
+ sentinel_iterator operator--(int) {
+ sentinel_iterator c(*this);
+ --*this;
+ return c;
+ }
+
+ // random access iterator requirements
+
+ sentinel_iterator &operator+=(difference_type n) {
+ iter() += n;
+ return *this;
+ }
+
+ sentinel_iterator operator+(difference_type n) const {
+ sentinel_iterator c(*this);
+ return c += n;
+ }
+
+ friend sentinel_iterator operator+(difference_type n, sentinel_iterator x) {
+ return x+n;
+ }
+
+ sentinel_iterator& operator-=(difference_type n) {
+ iter() -= n;
+ return *this;
+ }
+
+ sentinel_iterator operator-(difference_type n) const {
+ sentinel_iterator c(*this);
+ return c -= n;
+ }
+
+ difference_type operator-(sentinel_iterator x) const {
+ return iter()-x.iter();
+ }
+
+ auto operator[](difference_type n) const -> decltype(*(this->iter())) {
+ return *(iter()+n);
+ }
+
+ bool operator<=(const sentinel_iterator& x) const {
+ return x.is_sentinel() || (!is_sentinel() && iter()<=x.iter());
+ }
+
+ bool operator<(const sentinel_iterator& x) const {
+ return !is_sentinel() && (x.is_sentinel() || iter()<=x.iter());
+ }
+
+ bool operator>=(const sentinel_iterator& x) const {
+ return !(x<*this);
+ }
+
+ bool operator>(const sentinel_iterator& x) const {
+ return !(x<=*this);
+ }
+};
+
+template <typename I, typename S>
+using sentinel_iterator_t =
+ typename std::conditional<std::is_same<I, S>::value, I, sentinel_iterator<I, S>>::type;
+
+template <typename I, typename S>
+sentinel_iterator_t<I, S> make_sentinel_iterator(const I& i, const S& s) {
+ return sentinel_iterator_t<I, S>(i);
+}
+
+template <typename I, typename S>
+sentinel_iterator_t<I, S> make_sentinel_end(const I& i, const S& s) {
+ return sentinel_iterator_t<I, S>(s);
+}
+
+}
+}
+}
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_counter.cpp b/tests/unit/test_counter.cpp
index 7c97f2277e08a865ebb0be0a7e7b9d4db18f7360..600d86ca171980ce6d17036d1d74936cafbd1ed0 100644
--- a/tests/unit/test_counter.cpp
+++ b/tests/unit/test_counter.cpp
@@ -128,4 +128,3 @@ REGISTER_TYPED_TEST_CASE_P(counter_test, value, compare, arithmetic, iterator_tr
using int_types = ::testing::Types<signed char, unsigned char, short, unsigned short, int, unsigned, std::size_t, std::ptrdiff_t>;
INSTANTIATE_TYPED_TEST_CASE_P(int_types, counter_test, int_types);
-
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/unit/test_range.cpp b/tests/unit/test_range.cpp
index 7847212e3f0b0546c1f974dbdebab3303a15ea1d..118a4dc12393e74b74a62f84eee645d0399a490a 100644
--- a/tests/unit/test_range.cpp
+++ b/tests/unit/test_range.cpp
@@ -11,14 +11,17 @@
#include <tbb/tbb_stddef.h>
#endif
+#include <util/counter.hpp>
+#include <util/meta.hpp>
#include <util/range.hpp>
+#include <util/sentinel.hpp>
using namespace nest::mc;
TEST(range, list_iterator) {
std::list<int> l = { 2, 4, 6, 8, 10 };
- auto s = util::make_range(l.begin(), l.end());
+ auto s = util::make_range(l.begin(), l.end());
EXPECT_EQ(s.left, l.begin());
EXPECT_EQ(s.right, l.end());
@@ -40,6 +43,10 @@ TEST(range, list_iterator) {
auto sum2 = std::accumulate(s.begin(), s.end(), 0);
EXPECT_EQ(check, sum2);
+
+ // Check that different begin/end iterators are treated correctly
+ auto sc = util::make_range(l.begin(), l.cend());
+ EXPECT_EQ(l.size(), sc.size());
}
TEST(range, pointer) {
@@ -70,10 +77,56 @@ TEST(range, pointer) {
EXPECT_TRUE(std::equal(s.begin(), s.end(), &xs[l]));
}
+TEST(range, empty) {
+ int xs[] = { 10, 11, 12, 13, 14, 15, 16 };
+ auto l = 2;
+ auto r = 5;
+
+ auto empty_range_ll = util::make_range((const int *) &xs[l], &xs[l]);
+ EXPECT_TRUE(empty_range_ll.empty());
+ EXPECT_EQ(empty_range_ll.begin() == empty_range_ll.end(),
+ empty_range_ll.empty());
+ EXPECT_EQ(0u, empty_range_ll.size());
+
+
+ auto empty_range_rr = util::make_range(&xs[r], (const int *) &xs[r]);
+ EXPECT_TRUE(empty_range_rr.empty());
+ EXPECT_EQ(empty_range_rr.begin() == empty_range_rr.end(),
+ empty_range_rr.empty());
+ EXPECT_EQ(0u, empty_range_rr.size());
+}
+
+// This is the same test for enabling the incrementing util::distance
+// implementation.
+template <typename I, typename E>
+constexpr bool uses_incrementing_distance() {
+ return !util::has_common_random_access_iterator<I, E>::value &&
+ util::is_forward_iterator<I>::value;
+}
+
+TEST(range, size_implementation) {
+ ASSERT_TRUE((uses_incrementing_distance<
+ std::list<int>::iterator,
+ std::list<int>::const_iterator
+ >()));
+
+ ASSERT_FALSE((uses_incrementing_distance<
+ std::vector<int>::iterator,
+ std::vector<int>::const_iterator
+ >()));
+
+ ASSERT_FALSE((uses_incrementing_distance<int*, int*>()));
+
+ ASSERT_FALSE((uses_incrementing_distance<const int*, int*>()));
+
+ ASSERT_FALSE((uses_incrementing_distance<int*, const int*>()));
+}
+
TEST(range, input_iterator) {
int nums[] = { 10, 9, 8, 7, 6 };
std::istringstream sin("10 9 8 7 6");
- auto s = util::make_range(std::istream_iterator<int>(sin), std::istream_iterator<int>());
+ auto s = util::make_range(std::istream_iterator<int>(sin),
+ std::istream_iterator<int>());
EXPECT_TRUE(std::equal(s.begin(), s.end(), &nums[0]));
}
@@ -115,6 +168,7 @@ TEST(range, sentinel) {
}
EXPECT_EQ(s, std::string(cstr));
+ EXPECT_EQ(s.size(), cstr_range.size());
s.clear();
for (auto c: canonical_view(cstr_range)) {
@@ -122,8 +176,100 @@ TEST(range, sentinel) {
}
EXPECT_EQ(s, std::string(cstr));
+
+ const char *empty_cstr = "";
+ auto empty_cstr_range = util::make_range(empty_cstr, null_terminated);
+ EXPECT_TRUE(empty_cstr_range.empty());
+ EXPECT_EQ(0u, empty_cstr_range.size());
}
+template <typename V>
+class counter_range: public ::testing::Test {};
+
+TYPED_TEST_CASE_P(counter_range);
+
+TYPED_TEST_P(counter_range, max_size) {
+ using int_type = TypeParam;
+ using unsigned_int_type = typename std::make_unsigned<int_type>::type;
+ using counter = util::counter<int_type>;
+
+ auto l = counter{int_type{1}};
+ auto r = counter{int_type{10}};
+ auto range = util::make_range(l, r);
+ EXPECT_EQ(std::numeric_limits<unsigned_int_type>::max(),
+ range.max_size());
+
+}
+
+TYPED_TEST_P(counter_range, extreme_size) {
+ using int_type = TypeParam;
+ using signed_int_type = typename std::make_signed<int_type>::type;
+ using unsigned_int_type = typename std::make_unsigned<int_type>::type;
+ using counter = util::counter<signed_int_type>;
+
+ auto l = counter{std::numeric_limits<signed_int_type>::min()};
+ auto r = counter{std::numeric_limits<signed_int_type>::max()};
+ auto range = util::make_range(l, r);
+ EXPECT_FALSE(range.empty());
+ EXPECT_EQ(std::numeric_limits<unsigned_int_type>::max(), range.size());
+
+}
+
+
+TYPED_TEST_P(counter_range, size) {
+ using int_type = TypeParam;
+ using signed_int_type = typename std::make_signed<int_type>::type;
+ using counter = util::counter<signed_int_type>;
+
+ auto l = counter{signed_int_type{-3}};
+ auto r = counter{signed_int_type{3}};
+ auto range = util::make_range(l, r);
+ EXPECT_FALSE(range.empty());
+ EXPECT_EQ(6u, range.size());
+
+}
+
+TYPED_TEST_P(counter_range, at) {
+ using int_type = TypeParam;
+ using signed_int_type = typename std::make_signed<int_type>::type;
+ using counter = util::counter<signed_int_type>;
+
+ auto l = counter{signed_int_type{-3}};
+ auto r = counter{signed_int_type{3}};
+ auto range = util::make_range(l, r);
+ EXPECT_EQ(range.front(), range.at(0));
+ EXPECT_EQ(0, range.at(3));
+ EXPECT_EQ(range.back(), range.at(range.size()-1));
+ EXPECT_THROW(range.at(range.size()), std::out_of_range);
+ EXPECT_THROW(range.at(30), std::out_of_range);
+}
+
+TYPED_TEST_P(counter_range, iteration) {
+ using int_type = TypeParam;
+ using signed_int_type = typename std::make_signed<int_type>::type;
+ using counter = util::counter<signed_int_type>;
+
+ auto j = signed_int_type{-3};
+ auto l = counter{signed_int_type{j}};
+ auto r = counter{signed_int_type{3}};
+
+ for (auto i : util::make_range(l, r)) {
+ EXPECT_EQ(j++, i);
+ }
+}
+
+REGISTER_TYPED_TEST_CASE_P(counter_range, max_size, extreme_size, size, at,
+ iteration);
+
+using int_types = ::testing::Types<signed char, unsigned char, short,
+ unsigned short, int, unsigned, long,
+ std::size_t, std::ptrdiff_t>;
+
+using signed_int_types = ::testing::Types<signed char, short,
+ int, long, std::ptrdiff_t>;
+
+INSTANTIATE_TYPED_TEST_CASE_P(int_types, counter_range, int_types);
+
#ifdef WITH_TBB
TEST(range, tbb_split) {
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);