diff --git a/miniapp/miniapp.cpp b/miniapp/miniapp.cpp
index f267066ccb35093057527d4da8165d41015d1452..dfff77fa544cb68be9f0454dec8380c99aa530f8 100644
--- a/miniapp/miniapp.cpp
+++ b/miniapp/miniapp.cpp
@@ -111,6 +111,7 @@ struct model {
double value;
};
std::string name;
+ std::string units;
index_type id;
std::vector<sample_type> samples;
};
@@ -139,9 +140,9 @@ struct model {
};
mc::sampler make_simple_sampler(
- index_type probe_gid, const std::string name, index_type id, float dt)
+ index_type probe_gid, const std::string& name, const std::string& units, index_type id, float dt)
{
- traces.push_back(trace_data{name, id});
+ traces.push_back(trace_data{name, units, id});
return {probe_gid, simple_sampler_functor(traces, traces.size()-1, dt)};
}
@@ -156,14 +157,20 @@ struct model {
auto path = "trace_" + std::to_string(trace.id)
+ "_" + trace.name + ".json";
- nlohmann::json json;
- json["name"] = trace.name;
+ nlohmann::json jrep;
+ jrep["name"] = trace.name;
+ jrep["units"] = trace.units;
+ jrep["id"] = trace.id;
+
+ auto& jt = jrep["data"]["time"];
+ auto& jy = jrep["data"][trace.name];
+
for (const auto& sample: trace.samples) {
- json["time"].push_back(sample.time);
- json["value"].push_back(sample.value);
+ jt.push_back(sample.time);
+ jy.push_back(sample.value);
}
std::ofstream file(path);
- file << std::setw(1) << json << std::endl;
+ file << std::setw(1) << jrep << std::endl;
}
}
};
@@ -344,8 +351,8 @@ void all_to_all_model(nest::mc::io::cl_options& options, model& m) {
// setup probes
//
- mc::util::profiler_leave();
- mc::util::profiler_enter("probes");
+ mc::util::profiler_leave();
+ mc::util::profiler_enter("probes");
// monitor soma and dendrite on a few cells
float sample_dt = 0.1;
@@ -358,12 +365,16 @@ void all_to_all_model(nest::mc::io::cl_options& options, model& m) {
auto lid = m.communicator.group_lid(gid);
auto probe_soma = m.cell_groups[lid].probe_gid_range().first;
auto probe_dend = probe_soma+1;
+ auto probe_dend_current = probe_soma+2;
m.cell_groups[lid].add_sampler(
- m.make_simple_sampler(probe_soma, "vsoma", gid, sample_dt)
+ m.make_simple_sampler(probe_soma, "vsoma", "mV", gid, sample_dt)
);
m.cell_groups[lid].add_sampler(
- m.make_simple_sampler(probe_dend, "vdend", gid, sample_dt)
+ m.make_simple_sampler(probe_dend, "vdend", "mV", gid, sample_dt)
+ );
+ m.cell_groups[lid].add_sampler(
+ m.make_simple_sampler(probe_dend_current, "idend", "mA/cm²", gid, sample_dt)
);
}
@@ -421,10 +432,12 @@ mc::cell make_cell(int compartments_per_segment, int num_synapses) {
// add probes:
auto probe_soma = cell.add_probe({0, 0}, mc::probeKind::membrane_voltage);
auto probe_dendrite = cell.add_probe({1, 0.5}, mc::probeKind::membrane_voltage);
+ auto probe_dendrite_current = cell.add_probe({1, 0.5}, mc::probeKind::membrane_current);
EXPECTS(probe_soma==0);
EXPECTS(probe_dendrite==1);
- (void)probe_soma, (void)probe_dendrite;
+ EXPECTS(probe_dendrite_current==2);
+ (void)probe_soma, (void)probe_dendrite, (void)probe_dendrite_current;
return cell;
}
diff --git a/scripts/README.md b/scripts/README.md
new file mode 100644
index 0000000000000000000000000000000000000000..e645a762eaf88c5f7196a995502589da44c2f89c
--- /dev/null
+++ b/scripts/README.md
@@ -0,0 +1,81 @@
+tsplot
+------
+
+The `tsplot` script is a wrapper around matplotlib for displaying a collection of
+time series plots.
+
+## Input data
+
+`tsplot` reads timeseries in JSON format, according to the following conventions.
+
+```
+{
+ "units": <units>
+ "name": <name>
+ <other key-value metadata>
+ "data": {
+ "time": [ <time values> ]
+ <trace name>: [ <trace values> ]
+ }
+}
+```
+
+The `data` object must contain numeric arrays, with at least one with the key `time`;
+other members of `data` correspond to traces sampled at the corresponding time values.
+
+The other members of the top level object are regarded as metadata, with some keys
+treated specially:
+ * `units` are used to distinguish different axes for plotting, and the labels for those
+ axes. It's value is either a string, where the specified unit is taken as applying to
+ all included traces, or an object representing a mapping of trace names to their
+ corresponding unit string.
+ * `name` is taken as the title of the corresponding plot, if it is unambiguous.
+ * `label` is ignored: the _label_ for a trace is its name in the data object.
+
+## Operation
+
+The basic usage is simply:
+```
+tsplot data.json ...
+```
+which will produce an interactive plot of the timeseries provided by the provided
+files, with one trace per subplot.
+
+### Grouping
+
+Traces can be gathered on to the same subplot by grouping by metadata with the
+`-g` or `--group` option. To collect all traces with the same value of the key
+'id' and the same units:
+```
+tsplot -g units,id data.json ...
+```
+A subplot can comprise data with to two differint units, and will be plotted
+with two differing vertical axes.
+
+Note that for the purposes of tsplot, the value of the key _label_ is the
+propertu name of the trace in its json representation.
+
+### Restricting data
+
+The `-t` or `--trange` option exlcudes any points that have a time range outside
+that specified. Ranges are given by two numbers separated by a comma, but one or
+the other can be omitted to indicate that there is no bound on that side. For
+example:
+```
+tsplot -t ,100 data.json ...
+```
+will display all points with a time value less than or equal to 100.
+
+Extreme values for data can be automatically excluded and marked on the plot
+with the `-x` or `--exclude` option, taking a parameter _N_. All values in a
+timeseries that lie outside the interval [ _m_ - _Nr_, _m_ + _Nr_ ] are omitted,
+where _m_ is the median of the finite values in the timeseries, and _r_ is
+the 90% interquantile gap, that is, the difference between the 5% and 95% quantile
+of the timeseries data.
+
+### Output to file
+
+Use the `-o` or `--output` option to save the plot as an image, instead of
+displaying it interactively.
+
+
diff --git a/scripts/tsplot b/scripts/tsplot
new file mode 100755
index 0000000000000000000000000000000000000000..81a5da78c06a24d71f727431ef237a18f2a13b0b
--- /dev/null
+++ b/scripts/tsplot
@@ -0,0 +1,363 @@
+#!env python
+#coding: utf-8
+
+import argparse
+import json
+import numpy as np
+import sys
+import re
+import logging
+import matplotlib as M
+import matplotlib.pyplot as P
+from distutils.version import LooseVersion
+from itertools import chain, islice, cycle
+
+# Run-time check for matplot lib version for line style functionality.
+if LooseVersion(M.__version__)<LooseVersion("1.5.0"):
+ logging.critical("require matplotlib version ≥ 1.5.0")
+ sys.exit(1)
+
+# Read timeseries data from multiple files, plot each in one planel, with common
+# time axis, and optionally sharing a vertical axis as governed by the configuration.
+
+def parse_clargs():
+ def float_or_none(s):
+ try: return float(s)
+ except ValueError: return None
+
+ def parse_range_spec(s):
+ l, r = (float_or_none(x) for x in s.split(','))
+ return (l,r)
+
+ P = argparse.ArgumentParser(description='Plot time series data on one or more graphs.')
+ 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,
+ help='restrict time axis to RANGE (see below)')
+ P.add_argument('-g', '--group', metavar='KEY,...', dest='groupby',
+ 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('-x', '--exclude', metavar='NUM', dest='exclude',
+ type=float,
+ help='remove extreme points outside NUM times the 0.9-interquantile range of the median')
+
+ P.epilog = 'A range is specifed by a pair of floating point numbers min,max where '
+ P.epilog += 'either may be omitted to indicate the minimum or maximum of the corresponding '
+ P.epilog += 'values in the data.'
+
+ # modify args to avoid argparse having a fit when it encounters an option
+ # argument of the form '<negative number>,...'
+
+ argsbis = [' '+a if re.match(r'-[\d.]',a) else a for a in sys.argv[1:]]
+ return P.parse_args(argsbis)
+
+def isstring(s):
+ return isinstance(s,str) or isinstance(s,unicode)
+
+def take(n, s):
+ return islice((i for i in s), 0, n)
+
+class TimeSeries:
+ def __init__(self, ts, ys, **kwargs):
+ self.t = np.array(ts)
+ n = self.t.shape[0]
+
+ self.y = np.full_like(self.t, np.nan)
+ ny = min(len(ys), len(self.y))
+ self.y[:ny] = ys[:ny]
+
+ self.meta = dict(kwargs)
+ self.ex_ts = None
+
+ def trestrict(self, bounds):
+ clip = range_meet(self.trange(), bounds)
+ self.t = np.ma.masked_outside(self.t, v1=clip[0], v2=clip[1])
+ self.y = np.ma.masked_array(self.y, mask=self.t.mask)
+
+ def exclude_outliers(self, iqr_factor):
+ yfinite = np.ma.masked_invalid(self.y).compressed()
+ 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)
+
+ tex = np.ma.masked_array(self.t, mask=yex.mask)
+ self.ex_ts = TimeSeries(tex.compressed(), yex.compressed())
+ self.ex_ts.meta = dict(self.meta)
+
+ self.y = np.ma.filled(np.ma.masked_array(self.y, mask=~yex.mask), np.nan)
+
+ def excluded(self):
+ return self.ex_ts
+
+ def name(self):
+ return self.meta.get('name',"") # value of 'name' attribute in source
+
+ def label(self):
+ return self.meta.get('label',"") # name of column in source
+
+ def units(self):
+ return self.meta.get('units',"")
+
+ def trange(self):
+ return (min(self.t), max(self.t))
+
+ def yrange(self):
+ return (min(self.y), max(self.y))
+
+
+def read_json_timeseries(source):
+ j = json.load(source)
+
+ # Convention:
+ #
+ # Time series data is represented by an object with a subobject 'data' and optionally
+ # other key/value entries comprising metadata for the time series.
+ #
+ # The 'data' object holds one array of numbers 'time' and zero or more other
+ # numeric arrays of sample values corresponding to the values in 'time'. The
+ # names of these other arrays are taken to be the labels for the plots.
+ #
+ # Units can be specified by a top level entry 'units' which is either a string
+ # (units are common for all data series in the object) or by a map that
+ # takes a label to a unit string.
+
+ ts_list = []
+ jdata = j['data']
+ ncol = len(jdata)
+
+ times = jdata['time']
+ nsample = len(times)
+
+ def units(label):
+ try:
+ unitmap = j['units']
+ if isstring(unitmap):
+ return unitmap
+ else:
+ return unitmap[label]
+ except:
+ return ""
+
+ i = 1
+ for key in jdata.keys():
+ if key=="time": continue
+
+ meta = dict(j.items() + {'label': key, 'data': None, 'units': units(key)}.items())
+ del meta['data']
+
+ 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]))
+
+def range_meet(r, s):
+ return (max(r[0], s[0]), min(r[1], s[1]))
+
+
+class PlotData:
+ def __init__(self, key_label=""):
+ self.series = []
+ self.group_key_label = key_label
+
+ def trange(self):
+ return reduce(range_join, [s.trange() for s in self.series])
+
+ def yrange(self):
+ return reduce(range_join, [s.yrange() for s in self.series])
+
+ def name(self):
+ return reduce(lambda n, s: n or s.name(), self.series, "")
+
+ def group_label(self):
+ return self.group_key_label
+
+ def unique_labels(self):
+ # attempt to create unique labels for plots in the group based on
+ # meta data
+ labels = [s.label() for s in self.series]
+ if len(labels)<2:
+ return labels
+
+ n = len(labels)
+ keyset = reduce(lambda k, s: k.union(s.meta.keys()), self.series, set())
+ keyi = iter(keyset)
+ try:
+ while len(set(labels)) != n:
+ k = next(keyi)
+ if k=='label':
+ continue
+
+ vs = [s.meta.get(k,None) for s in self.series]
+ if len(set(vs))==1:
+ continue
+
+ for i in xrange(n):
+ prefix = '' if k=='name' else k+'='
+ if vs[i] is not None:
+ labels[i] += u', '+k+u'='+unicode(vs[i])
+
+ except StopIteration:
+ pass
+
+ return labels
+
+# Input: list of TimeSeries objects; collection of metadata keys to group on
+# Return list of plot info (time series, data extents, metadata), one per plot.
+
+def gather_ts_plots(tss, groupby):
+ group_lookup = {}
+ plot_groups = []
+
+ for ts in tss:
+ key = tuple([ts.meta.get(g) for g in groupby])
+ if key is () or None in key or key not in group_lookup:
+ pretty_key=', '.join([unicode(k)+u'='+unicode(v) for k,v in zip(groupby, key) if v is not None])
+ pd = PlotData(pretty_key)
+ pd.series = [ts]
+ plot_groups.append(pd)
+ group_lookup[key] = len(plot_groups)-1
+ else:
+ plot_groups[group_lookup[key]].series.append(ts)
+
+ return plot_groups
+
+
+def make_palette(cm_name, n, cmin=0, cmax=1):
+ smap = M.cm.ScalarMappable(M.colors.Normalize(cmin/float(cmin-cmax),(cmin-1)/float(cmin-cmax)),
+ 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):
+ nplots = len(plot_groups)
+ plot_groups = sorted(plot_groups, key=lambda g: g.group_label())
+
+ # use same global time scale for all plots
+ trange = reduce(range_join, [g.trange() for g in plot_groups])
+
+ # use group names for titles?
+ group_titles = any((g.group_label() for g in plot_groups))
+
+ figure = P.figure()
+ for i in xrange(nplots):
+ group = plot_groups[i]
+ plot = figure.add_subplot(nplots, 1, i+1)
+
+ title = group.group_label() if group_titles else group.name()
+ plot.set_title(title)
+
+ # y-axis label: use timeseries label and units if only
+ # one series in group, otherwise use a legend with labels,
+ # and units alone on the axes. At most two different unit
+ # axes can be drawn.
+
+ def ylabel(unit):
+ if len(group.series)==1:
+ lab = group.series[0].label()
+ if unit:
+ lab += ' (' + unit + ')'
+ else:
+ lab = unit
+
+ return lab
+
+ uniq_units = list(set([s.units() for s in group.series]))
+ uniq_units.sort()
+ if len(uniq_units)>2:
+ logging.warning('more than two different units on the same plot')
+ uniq_units = uniq_units[:2]
+
+ # store each series in a slot corresponding to one of the units,
+ # together with a best-effort label
+
+ series_by_unit = [[] for i in xrange(len(uniq_units))]
+ unique_labels = group.unique_labels()
+
+ for si in xrange(len(group.series)):
+ s = group.series[si]
+ label = unique_labels[si]
+ try:
+ series_by_unit[uniq_units.index(s.units())].append((s,label))
+ except ValueError:
+ pass
+
+ # TODO: need to find a scheme of colour/line allocation that is
+ # double y-axis AND greyscale friendly.
+
+ palette = \
+ [make_palette(cm, n, 0, 0.5) for
+ 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:
+ plot = plot.twinx()
+
+ axis_color = palette[ui][0]
+ plot.set_ylabel(ylabel(uniq_units[ui]), color=axis_color)
+ for l in plot.get_yticklabels():
+ l.set_color(axis_color)
+
+ 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])
+ line = next(lines)
+ for s, l in series_by_unit[ui]:
+ c = next(colours)
+ plot.plot(s.t, s.y, color=c, ls=line, label=l)
+ # treat exluded points especially
+ ex = s.excluded()
+ if ex is not None:
+ ymin, ymax = s.yrange()
+ plot.plot(ex.t, np.clip(ex.y, ymin, ymax), marker='x', ls='', color=c)
+
+ if first_plot:
+ plot.legend(loc=2, fontsize='small')
+ plot.grid()
+ else:
+ plot.legend(loc=1, fontsize='small')
+
+ first_plot = False
+
+ # adapted from http://stackoverflow.com/questions/6963035
+ 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)
+ else:
+ P.show()
+
+args = parse_clargs()
+tss = []
+for filename in args.inputs:
+ with open(filename) as f:
+ tss.extend(read_json_timeseries(f))
+
+if args.trange:
+ for ts in tss:
+ ts.trestrict(args.trange)
+
+if args.exclude:
+ for ts in tss:
+ ts.exclude_outliers(args.exclude)
+
+groupby = args.groupby if args.groupby else []
+plots = gather_ts_plots(tss, groupby)
+
+if not args.outfile:
+ M.interactive(False)
+
+plot_plots(plots, save=args.outfile)