diff --git a/benchmarks/engines/busyring/arbor/ring.cpp b/benchmarks/engines/busyring/arbor/ring.cpp index bd7d97f..9e52086 100644 --- a/benchmarks/engines/busyring/arbor/ring.cpp +++ b/benchmarks/engines/busyring/arbor/ring.cpp @@ -32,7 +32,7 @@ using arb::cell_size_type; using arb::cell_member_type; using arb::cell_kind; using arb::time_type; -using arb::cell_probe_address; +using arb::cable_probe_membrane_voltage; // Writes voltage trace as a json file. void write_trace_json(std::string fname, const arb::trace_data& trace); @@ -127,12 +127,9 @@ class ring_recipe: public arb::recipe { } arb::probe_info get_probe(cell_member_type id) const override { - // Get the appropriate kind for measuring voltage. - cell_probe_address::probe_kind kind = cell_probe_address::membrane_voltage; // Measure at the soma. - arb::segment_location loc(0, 0.0); - - return arb::probe_info{id, kind, cell_probe_address{loc, kind}}; + arb::mlocation loc{0, 0.0}; + return arb::probe_info{id, 0, cable_probe_membrane_voltage{loc}}; } private: @@ -146,8 +143,7 @@ class ring_recipe: public arb::recipe { struct cell_stats { using size_type = unsigned; size_type ncells = 0; - size_type nsegs = 0; - size_type ncomp = 0; + size_type nbranch = 0; cell_stats(arb::recipe& r) { #ifdef ARB_MPI_ENABLED @@ -158,21 +154,17 @@ struct cell_stats { size_type cells_per_rank = ncells/nranks; size_type b = rank*cells_per_rank; size_type e = (rank==nranks-1)? ncells: (rank+1)*cells_per_rank; - size_type nsegs_tmp = 0; - size_type ncomp_tmp = 0; + size_type nbranch_tmp = 0; for (size_type i=b; i(r.get_cell_description(i)); - nsegs_tmp += c.num_segments(); - ncomp_tmp += c.num_compartments(); + nbranch_tmp += c.morphology().num_branches(); } - MPI_Allreduce(&nsegs_tmp, &nsegs, 1, MPI_UNSIGNED, MPI_SUM, MPI_COMM_WORLD); - MPI_Allreduce(&ncomp_tmp, &ncomp, 1, MPI_UNSIGNED, MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(&nbranch_tmp, &nbranch, 1, MPI_UNSIGNED, MPI_SUM, MPI_COMM_WORLD); #else ncells = r.num_cells(); for (size_type i=0; i(r.get_cell_description(i)); - nsegs += c.num_segments(); - ncomp += c.num_compartments(); + nbranch += c.morphology().num_branches(); } #endif } @@ -180,8 +172,8 @@ struct cell_stats { friend std::ostream& operator<<(std::ostream& o, const cell_stats& s) { return o << "cell stats: " << s.ncells << " cells; " - << s.nsegs << " segments; " - << s.ncomp << " compartments."; + << s.nbranch << " branches; " + << 0 << " compartments; "; } }; @@ -336,13 +328,11 @@ double interp(const std::array& r, unsigned i, unsigned n) { } arb::cable_cell branch_cell(arb::cell_gid_type gid, const cell_parameters& params) { - arb::cable_cell cell; - - cell.default_parameters.axial_resistivity = 100; + arb::sample_tree tree; // Add soma. - auto soma = cell.add_soma(12.6157/2.0); // For area of 500 μm². - soma->add_mechanism("hh"); + double soma_radius = 12.6157/2.0; + tree.append(arb::mnpos, {{0,0,0,soma_radius}, 1}); // For area of 500 μm². std::vector> levels; levels.push_back({0}); @@ -353,7 +343,7 @@ arb::cable_cell branch_cell(arb::cell_gid_type gid, const cell_parameters& param double dend_radius = 0.5; // Diameter of 1 μm for each cable. - unsigned nsec = 1; + double dist_from_soma = soma_radius; for (unsigned i=0; iset_compartments(nc); - dend->add_mechanism("pas"); + auto z = dist_from_soma; + auto p = tree.append(sec, {{0,0,z,dend_radius}, 3}); + if (nc>1) { + auto dz = l/nc; + for (unsigned k=1; k pos_dis(0, 1); - std::uniform_int_distribution seg_dis(1, nsec-1); + // Add additional synapses that will not be connected to anything. for (unsigned i=1u; i0 else 0 relerr_lb = abserr_lb/r_absmax if r_absmax>0 else 0 + relerr_rms = abserr_rms/r_absmax if r_absmax>0 else 0 + relerr_rms_lb = abserr_rms_lb/r_absmax if r_absmax>0 else 0 out[var+'.delta'] = delta out[var+'.interperr'] = interperr out[var+'.abserr'] = abserr out[var+'.abserr.lb'] = abserr_lb + out[var+'.abserr.rms'] = abserr_rms + out[var+'.abserr.rms.lb'] = abserr_rms_lb out[var+'.relerr'] = relerr out[var+'.relerr.lb'] = relerr_lb + out[var+'.relerr.rms'] = relerr_rms + out[var+'.relerr.rms.lb'] = relerr_rms_lb out.to_netcdf(opts.output) diff --git a/common/bin/tsplotx b/common/bin/tsplotx index 14ee26b..108e02d 100755 --- a/common/bin/tsplotx +++ b/common/bin/tsplotx @@ -302,7 +302,7 @@ def gather_ts_plots(tss, groupby): 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: + if key == () 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] diff --git a/docs/install.rst b/docs/install.rst index b0c6627..88f12a0 100644 --- a/docs/install.rst +++ b/docs/install.rst @@ -64,17 +64,18 @@ The simulation engines to install are provided as arguments. Further options for installing the simulation engines in a user-specified path and customising the build environment can be provided: -==================== ================= ====================================================== -Flag Default value Explanation -==================== ================= ====================================================== -simulator none Which simulation engines to download and install. - Any number of the following: {``arbor``, ``neuron``, ``coreneuron``, ``all``}. -``--prefix`` current path Path for downloading, compiling, installing simulation engines. - Also used to store inputs and outputs from benchmarks and validation tests. - Can be either a relative or absolute path. -``--env`` none Optional script for configuring the environment and build steps. - See :ref:`customenv`. -==================== ================= ====================================================== +========================= ================= ====================================================== +Flag Default value Explanation +========================= ================= ====================================================== +simulator none Which simulation engines to download and install. + Any number of the following: {``arbor``, ``neuron``, ``coreneuron``, ``all``}. +``--prefix`` current path Path for downloading, compiling, installing simulation engines. + Also used to store inputs and outputs from benchmarks and validation tests. + Can be either a relative or absolute path. +``--env`` none Optional script for configuring the environment and build steps. + See :ref:`customenv`. +``--disable-validation`` false Suppress building of validation. +======================== ================= ====================================================== .. _customenv: @@ -184,7 +185,7 @@ Variable Default value Explanat ``ns_arb_arch`` ``native`` `The CPU architecture target `_ for Arbor. Must be set when cross compiling. Default ``native`` targets the architecture used to configure NSuite. -``ns_arb_with_gpu`` ``OFF`` Whether to build Arbor with NVIDIA GPU support. +``ns_arb_with_gpu`` ``none`` Whether to build Arbor with GPU support, valid choices are (none|cuda|cuda-clang|hip) ``ns_arb_vectorize`` ``ON`` Whether to use explicit vectorization for Arbor. ======================== =========================================== ====================================================== diff --git a/install-local.sh b/install-local.sh index 9b0093c..d8dd757 100755 --- a/install-local.sh +++ b/install-local.sh @@ -6,16 +6,16 @@ Setup NSuite framework, build and install simulators, benchmarks, and validation tests. Options: - --pyvenv=VENVOPT Customize use of Python virtual environment (see below) - --env=SCRIPT Source SCRIPT before building. - --prefix=PATH Use PATH as base for install and other run-time - working directories. - -TARGET is one of: - arbor Build Arbor. - neuron Build NEURON. - coreneuron Build CoreNEURON. - all Build all simulators. + --pyvenv=VENVOPT Customize use of Python virtual environment (see below) + --env=SCRIPT Source SCRIPT before building. + --prefix=PATH Use PATH as base for install and other run-time + working directories. + --disable-validation Skip building of validation tests. +TARGET is one of: + arbor Build Arbor. + neuron Build NEURON. + coreneuron Build CoreNEURON. + all Build all simulators. Building a TARGET will also build any associated tests and benchmarks as required. @@ -23,9 +23,9 @@ benchmarks as required. By default, a Python virtual environment is created in which required modules are installed if they are missing in the system environment. VENVOPT is one of: - disable Do not use a Python virtualenv at all. - inherit Use a virtualenv, using site packages. - enable Use a virtualenv, ignoring site packages (default). + disable Do not use a Python virtualenv at all. + inherit Use a virtualenv, using site packages. + enable Use a virtualenv, ignoring site packages (default). _end_ } @@ -83,6 +83,9 @@ do shift ns_prefix=$1 ;; + --disable-validation ) + ns_disable_validation=true + ;; * ) echo "unknown option '$1'" usage @@ -201,11 +204,13 @@ cd "$ns_base_path" [ "$ns_build_coreneuron" = true ] && echo && source "$ns_base_path/scripts/build_coreneuron.sh" cd "$ns_base_path" -# Always attempt to build validation models/generators. -echo -msghi "Building validation tests and generators" -source "$ns_base_path/scripts/build_validation_models.sh" -cd "$ns_base_path" +# attempt to build validation models/generators. +if [ "$ns_disable_validation" != "true" ]; then + echo + msghi "Building validation tests and generators" + source "$ns_base_path/scripts/build_validation_models.sh" + cd "$ns_base_path" +fi echo msghi "Installation finished" diff --git a/scripts/build_arbor.sh b/scripts/build_arbor.sh index 1ebec2b..6cb7715 100644 --- a/scripts/build_arbor.sh +++ b/scripts/build_arbor.sh @@ -48,7 +48,9 @@ fi cd "$arb_build_path" cmake_args=-DCMAKE_INSTALL_PREFIX:PATH="$ns_install_path" cmake_args="$cmake_args -DARB_WITH_MPI=$ns_with_mpi" -cmake_args="$cmake_args -DARB_WITH_GPU=$ns_arb_with_gpu" +if [ "$ns_arb_with_gpu" != "none" ]; then + cmake_args="$cmake_args -DARB_GPU=$ns_arb_with_gpu" +fi cmake_args="$cmake_args -DARB_ARCH=$ns_arb_arch" cmake_args="$cmake_args -DARB_VECTORIZE=$ns_arb_vectorize" if [ "$ns_arb_xcompile_modcc" == "ON" ]; then diff --git a/scripts/environment.sh b/scripts/environment.sh index f0efbb2..e592d8b 100644 --- a/scripts/environment.sh +++ b/scripts/environment.sh @@ -67,10 +67,10 @@ default_environment() { # Arbor specific ns_arb_git_repo=https://github.com/arbor-sim/arbor.git - ns_arb_branch=v0.2.2 + ns_arb_branch=master ns_arb_arch=native - ns_arb_with_gpu=OFF + ns_arb_with_gpu=none ns_arb_vectorize=ON ns_arb_xcompile_modcc=OFF diff --git a/systems/daint-gpu.sh b/systems/daint-gpu.sh index 4fbe386..3028e06 100755 --- a/systems/daint-gpu.sh +++ b/systems/daint-gpu.sh @@ -28,7 +28,7 @@ ns_cc=$(which cc) ns_cxx=$(which CC) ns_with_mpi=ON -ns_arb_with_gpu=ON +ns_arb_with_gpu=cuda ns_arb_arch=haswell export CRAYPE_LINK_TYPE=dynamic diff --git a/systems/jureca-gpu.sh b/systems/jureca-gpu.sh index ec79e94..076d16f 100644 --- a/systems/jureca-gpu.sh +++ b/systems/jureca-gpu.sh @@ -30,7 +30,7 @@ ns_cc=$(which mpicc) ns_cxx=$(which mpicxx) ns_with_mpi=ON -ns_arb_with_gpu=ON +ns_arb_with_gpu=cuda ns_arb_arch=haswell ns_arb_branch=master diff --git a/systems/juwels-gpu.sh b/systems/juwels-gpu.sh index 458fab0..084e132 100644 --- a/systems/juwels-gpu.sh +++ b/systems/juwels-gpu.sh @@ -32,7 +32,7 @@ ns_cc=$(which mpicc) ns_cxx=$(which mpicxx) ns_with_mpi=ON -ns_arb_with_gpu=ON +ns_arb_with_gpu=cuda ns_arb_arch=skylake-avx512 ns_arb_branch=master diff --git a/validation/cable-steadystate/README.md b/validation/cable-steadystate/README.md new file mode 100644 index 0000000..860d338 --- /dev/null +++ b/validation/cable-steadystate/README.md @@ -0,0 +1,68 @@ +Model cable-steadystate +======================= + +Simulates a passive cable with the same electrical characteristics (excluding +reversal potential) as the Rallpack1 model (see [@bhalla1992]). The simulation +is run for 20·τ where τ is the electrical time constant, in order to reach +steady state. + +The electrical parameters are as follows: + +| Parameter | Value | +|--------------------------------|-----------| +| Cable diameter | 1.0 µm | +| Cable length | 1.0 mm | +| Membrane resistivity | 4.0 Ω m² | +| Membrane specific capacitiance | 0.01 F/m² | +| Membrane reversal potential | 0.0 mV | +| Injected current | 0.1 nA | + +The initial membrane voltage is the reversal potential. + +Test parameters +--------------- + + +| Parameter | Interpretation | +|-----------|----------------| +| `dt` | maximum simulation integration timestep [ms] | +| `d0` | diameter at left end of cable [µm] | +| `d1` | diameter at right end of cable [µm] | +| `n` | number of compartments/control volumes | + + +Acceptance critera +------------------ + +The computed membrane voltages at the end of the simulation period +are compared against the steady state analytic solution at +the centre of each control volume/compartment. + +The result is accepted if the maximum relative error is less than 0.1%. + +Implementation notes +-------------------- + +### NEURON + +Supported tags: +* `firstorder` + + Use first order integrator instead of default second order. + +--- +references: +- id: bhalla1992 + type: article-journal + author: + - { family: Bhalla, given: Upinder S. } + - { family: Bilitch, given: David H. } + - { family: Bower, given: James M. } + issued: { raw: "1992-11" } + title: 'Rallpacks: A set of benchmarks for neuronal simulators' + container-title: Trends in Neurosciences + page: 453–458 + volume: 15 + issue: 11 + DOI: 10.1016/0166-2236(92)90009-W +... diff --git a/validation/cable-steadystate/default.param b/validation/cable-steadystate/default.param new file mode 100644 index 0000000..e98b943 --- /dev/null +++ b/validation/cable-steadystate/default.param @@ -0,0 +1,4 @@ +dt=0.05 +n=1000 +d0=1.0 +d1=1.7 diff --git a/validation/cable-steadystate/generate-cable-steadystate b/validation/cable-steadystate/generate-cable-steadystate new file mode 100755 index 0000000..f4bd2a6 --- /dev/null +++ b/validation/cable-steadystate/generate-cable-steadystate @@ -0,0 +1,88 @@ +#!/usr/bin/env python + +from __future__ import print_function +from itertools import count +import math +import scipy.integrate as integrate +import scipy.special as special +import numpy as np +import xarray +import re +import sys + +import nsuite.stdarg as stdarg +import nsuite.stdattr as stdattr + +ra = 1.0 # axial resistivity [Ω m] +rm = 4.0 # membrane resistivity [Ω m²] +cm = 0.01 # memrane specific capacitance [F/m²] +Erev = 0 # reversal potential [mV] + +length = 1000 # cable length [µm] +iinj = 0.1 # current injection [nA] + +# Parameters +d0 = 1.0 # cable diameter (left) [µm] +d1 = 1.5 # cable diameter (right) [µm] +dt = 0 # (ignored) +n = 0 # used only to determine number of samples + +output, tags, params = stdarg.parse_run_stdarg() +for v in ['d0', 'd1', 'dt', 'n']: + if v in params: globals()[v] = params[v] + +n = int(max(1000, 3*n)) + +def k1(x): return special.kn(1, x) +def k2(x): return special.kn(2, x) +def i1(x): return special.iv(1, x) +def i2(x): return special.iv(2, x) + +def u(x, a, b): + if a==0: + return math.cosh(x)/math.sinh(b) + else: + ra = math.sqrt(a) + rb = math.sqrt(b) + rx = math.sqrt(x) + + return b/rx*(k2(2*ra)*i1(2*rx)+i2(2*ra)*k1(2*rx))/(k2(2*ra)*i2(2*rb)-i2(2*ra)*k2(2*rb)) + +def membrane_voltage(): + if (d0==d1): + r1 = d1/2 + lam = math.sqrt(rm*r1/(2*ra))*1e3 # [µm] + a = 0 + b = length/lam + E = -iinj*ra/(math.pi*r1*r1) + + else: + r0 = d0/2 + r1 = d1/2 + dr = r1-r0 + sintheta = math.sqrt(1-1/(1+(dr/length)**2)) + A = r0/dr * length + B = r1/dr * length + lam = rm/(2*ra)*sintheta*1e6 # [µm] + a = A/lam + b = B/lam + E = -iinj*ra/(math.pi*r1*r1) + + xs = np.linspace(a, b, n) + vs = np.zeros(xs.size) + for i in range(xs.size): + vs[i] = Erev - lam*E*u(xs[i], a, b) + + xs = (xs - a)*lam # rescale to physical coordinates + return xs, vs + + +xs, vs = membrane_voltage() + +out = xarray.Dataset({'v': (['x'], vs)}, coords={'x': xs}) +out.x.attrs['units'] = 'µm' +out.v.attrs['units'] = 'mV' + +stdattr.set_stdattr(out, model='cable-steadystate', simulator='reference', params=params) +out.to_netcdf(output) + diff --git a/validation/cable-steadystate/neuron-cable-steadystate.py b/validation/cable-steadystate/neuron-cable-steadystate.py new file mode 100755 index 0000000..872a819 --- /dev/null +++ b/validation/cable-steadystate/neuron-cable-steadystate.py @@ -0,0 +1,111 @@ +#!/usr/bin/env python + +from math import pi +import sys +import os + +from neuron import h +import nsuite.stdarg as stdarg +import nsuite.stdattr as stdattr +import numpy as np +import xarray + +ra = 1.0 # axial resistivity [Ω m] +rm = 4.0 # membrane resistivity [Ω m²] +cm = 0.01 # memrane specific capacitance [F/m²] +Erev = 0 # reversal potential [mV] + +length = 1000.0 # cable length [µm] +iinj = 0.1 # current injection [nA] + +# Parameters +d0 = 1.0 # cable diameter (left) [µm] +d1 = 1.5 # cable diameter (right) [µm] +dt = 0 # [ms] +n = 0 # number of compartments + +# One recognized tag: 'firstorder' + +output, tags, params = stdarg.parse_run_stdarg(tagset=['firstorder']) + +for v in ['d0', 'd1', 'dt', 'n']: + if v in params: globals()[v] = params[v] + +tau = rm*cm*1000 # time constant [ms] +tend = 20*tau # transient amplitude is proportional to exp(-t/tau). +tsample = tend-dt + +# TODO: resolve routines below with exisiting neuron support code. + +def hoc_execute_quiet(arg): + with open(os.devnull, 'wb') as null: + fd = sys.stdout.fileno() + keep = os.dup(fd) + sys.stdout.flush() + os.dup2(null.fileno(), fd) + h(arg) + sys.stdout.flush() + os.dup2(keep, fd) + +def hoc_setup(): + hoc_execute_quiet('load_file("stdrun.hoc")') + +# Make model + +hoc_setup() + +cable = h.Section(name='cable') +h.pt3dadd(0.0, 0.0, 0.0, d0, sec=cable) +h.pt3dadd(length, 0.0, 0.0, d1, sec=cable) + +cable.cm = 100*cm # [µF/cm² = 0.01 F/m²] +cable.Ra = 100*ra # [Ω cm = 0.01 Ω m] +cable.nseg = int(n) + +cable.insert('pas') +cable.g_pas = 0.0001/rm # [S/cm² = 10000 S/m²] +cable.e_pas = Erev + +stim = h.IClamp(cable(1)) +stim.delay = 0 +stim.dur = tend +stim.amp = iinj + +h.v_init = Erev + +# Run model + +# Take samples from the middle of each compartment. +xs = [(2.0*i+1)/(2*n) for i in range(int(n))] + +tvec = h.Vector([tsample]) +vrecs = [] +for x in xs: + vrecs.append(h.Vector()) + vrecs[-1].record(cable(x/length)._ref_v, tvec) + +t = h.Vector() +t.record(h._ref_t, tvec) + +h.dt = dt +h.steps_per_ms = 1/dt # or else NEURON might noisily fudge dt +if 'firstorder' in tags: + h.secondorder = 0 +else: + h.secondorder = 2 +h.tstop = tend +h.run() + +# Collect and save data + +vs = [rec.x[0] for rec in vrecs] + +out = xarray.Dataset({'v': (['x'], vs), 't_sample': ([], t.x[0])}, coords={'x': list(xs)}) +out.x.attrs['units'] = 'µm' +out.v.attrs['units'] = 'mV' +out.t_sample.attrs['units'] = 'ms' + +nrnver = h.nrnversion() +stdattr.set_stdattr(out, model='cable-steadystate', simulator='neuron', simulator_build=nrnver, tags=tags, params=params) + +out.to_netcdf(output) diff --git a/validation/cable-steadystate/run b/validation/cable-steadystate/run new file mode 100755 index 0000000..f3e8a20 --- /dev/null +++ b/validation/cable-steadystate/run @@ -0,0 +1,40 @@ +#!/usr/bin/env bash +# Invocation: +# run outputdir simname paramsetname +# run --list-tags simname + +# Change to script directory and attempt to find nsuite base directory. + +unset CDPATH +cd "${BASH_SOURCE[0]%/*}" +[ -n "$ns_base_path" ] || ns_base_path="$(cd ../..; pwd)" + +# Set up model paths, model_XXX variables, and functions +# exit_model_fail, exit_model_missing. + +outdir=$1 +shift + +source "$ns_base_path/scripts/model_common.sh" +model_setup cable-steadystate "$@" + +# Run sim-specific implementation with tags and parameter data. + +impl="./run-$model_sim" +[ -x "$impl" ] || exit_model_missing + +outfile="$outdir/run.nc" +"$impl" -o "$outfile" ${model_impl_stdargs[@]} || exit $? + +# Generate reference data if required. + +reffile=$(model_find_cacheable "ref-${model_name}-${model_param}.nc") || \ + ./generate-cable-steadystate -o "$reffile" $model_param_data || exit 1 + +# Run comparison. +# Pass if relative error is less than 0.1%. + +deltafile="$outdir/delta.nc" +comparex "$outfile" --warn --ref "$reffile" --interpolate x -o "$deltafile" || exit 1 + +thresholdx "$deltafile" -e "v.relerr.lb<0.001" || exit_model_fail diff --git a/validation/cable-steadystate/run-arbor b/validation/cable-steadystate/run-arbor new file mode 100755 index 0000000..84a2138 --- /dev/null +++ b/validation/cable-steadystate/run-arbor @@ -0,0 +1,2 @@ +#!/usr/bin/env bash +arbor-cable-steadystate "${@}" diff --git a/validation/cable-steadystate/run-neuron b/validation/cable-steadystate/run-neuron new file mode 100755 index 0000000..4951029 --- /dev/null +++ b/validation/cable-steadystate/run-neuron @@ -0,0 +1,2 @@ +#!/usr/bin/env bash +./neuron-cable-steadystate.py "${@}" diff --git a/validation/rallpack1/README.md b/validation/rallpack1/README.md new file mode 100644 index 0000000..bfb84d2 --- /dev/null +++ b/validation/rallpack1/README.md @@ -0,0 +1,72 @@ +Model rallpack1 +================ + +The model comprises the Rallpack 1 test from [@bhalla1992], viz. +a constant-radius cable with passive dynamics and a constant +current injection at one end. + +The electrical parameters are as follows: + +| Parameter | Value | +|--------------------------------|-----------| +| Cable diameter | 1.0 µm | +| Cable length | 1.0 mm | +| Membrane resistivity | 4.0 Ω m² | +| Membrane specific capacitiance | 0.01 F/m² | +| Membrane reversal potential | -65.0 mV | +| Injected current | 0.1 nA | + +The initial membrane voltage is the reversal potential. + +Rallpack 1 specifies that the discretization has 1000 compartments, +which is interpreted as 1000 CVs for Arbor. + + +Test parameters +--------------- + + +| Parameter | Interpretation | +|-----------|----------------| +| `dt` | maximum simulation integration timestep [ms] | +| `x0` | first measurement point (proportional distance) | +| `x1` | second measurement point (proportional distance) | +| `n` | number of compartments/control volumes | + + +Acceptance critera +------------------ + +The simulated membrane voltages at x0 and x1 are compared against the analytic +solution from t = 0 to t = 250 ms. The relative error is computed as +the RMS error over time, divided by the maximum absolute value of +the voltage at that point over the simulation time interval. + +The result is accepted if both relative RMS errors are less than 0.1%. + +Implementation notes +-------------------- + +### NEURON + +Supported tags: +* `firstorder` + + Use first order integrator instead of default second order. + +--- +references: +- id: bhalla1992 + type: article-journal + author: + - { family: Bhalla, given: Upinder S. } + - { family: Bilitch, given: David H. } + - { family: Bower, given: James M. } + issued: { raw: "1992-11" } + title: 'Rallpacks: A set of benchmarks for neuronal simulators' + container-title: Trends in Neurosciences + page: 453–458 + volume: 15 + issue: 11 + DOI: 10.1016/0166-2236(92)90009-W +... diff --git a/validation/rallpack1/default.param b/validation/rallpack1/default.param new file mode 100644 index 0000000..b928a08 --- /dev/null +++ b/validation/rallpack1/default.param @@ -0,0 +1,5 @@ +dt=0.05 +n=1000 +x0=0 +x1=1 +tend=250 diff --git a/validation/rallpack1/generate-rallpack1 b/validation/rallpack1/generate-rallpack1 new file mode 100755 index 0000000..5723f1b --- /dev/null +++ b/validation/rallpack1/generate-rallpack1 @@ -0,0 +1,93 @@ +#!/usr/bin/env python + +from __future__ import print_function +from itertools import count +import math +import scipy.integrate as integrate +import numpy as np +import xarray +import re +import sys + +import nsuite.stdarg as stdarg +import nsuite.stdattr as stdattr + +ra = 1.0 # axial resistivity [Ω m] +rm = 4.0 # membrane resistivity [Ω m²] +cm = 0.01 # memrane specific capacitance [F/m²] +Erev = -65 # reversal potential [mV] + +diam = 1.0 # cable diameter [µm] +length = 1000 # cable length [µm] +iinj = 0.1 # current injection [nA] + +max_dt = 0.01 # [ms] + +# Parameters +tend = 250 # [ms] +x0 = 0 +x1 = 1 +dt = max_dt # [ms] (only used for nsamp) + +output, tags, params = stdarg.parse_run_stdarg() +for v in ['dt', 'tend', 'x0', 'x1']: + if v in params: globals()[v] = params[v] + +dt = min(dt, max_dt) +nsamp = int(tend/dt) +tend_s = tend*0.001 # [s] + +def u(x, t, b, tol): + if t<=0: + return 0.0 + + uinf = math.cosh(x)/math.sinh(b) + + acc = math.exp(-t)/2 + rtol = tol*t*b + + sign = 1 + for k in count(1): + a = k*math.pi/b + l = 1+a*a + q = math.exp(-l*t)/l + sign = -sign + + if q +#include +#include +#include +#include + +#include + +#include + +#include +#include +#include +#include +#include +#include + +#include "common_args.h" +#include "common_attr.h" +#include "netcdf_wrap.h" + +using namespace arb; + +paramset default_parameters = { + {"dt", 0.0025}, // timestep [ms] + {"d0", 1.0}, // diameter at left end [µm] + {"d1", 1.5}, // diameter at right end [µm] + {"n", 1000} // number of CVs +}; + +struct rc_cable_recipe: public arb::recipe { + // Fixed parameters: + + static constexpr double length = 1000.0; // cable length [µm] + static constexpr double rl = 1.0; // bulk resistivity [Ωm] + static constexpr double rm = 4.0; // membrane resistivity [Ωm²] + static constexpr double cm = 0.01; // membrane specific capacitance [F/m²] + static constexpr double iinj = 0.1; // current injection at right end [nA] + + // Customizable parameters: + unsigned n = 0; // number of CV + double d0 = 0, d1 = 0; // diameters + + explicit rc_cable_recipe(const paramset& ps): + d0(ps.at("d0")), + d1(ps.at("d1")), + n(static_cast(ps.at("n"))) + {} + + cell_size_type num_cells() const override { return 1; } + cell_size_type num_targets(cell_gid_type) const override { return 0; } + cell_size_type num_probes(cell_gid_type) const override { return n; } + cell_kind get_cell_kind(cell_gid_type) const override { return cell_kind::cable; } + + util::any get_global_properties(cell_kind kind) const override { + cable_cell_global_properties prop; + + prop.default_parameters.init_membrane_potential = 0; + prop.default_parameters.axial_resistivity = 100*rl; // [Ω·cm] + prop.default_parameters.membrane_capacitance = (double)cm; // [F/m²] + prop.default_parameters.temperature_K = 0; + prop.ion_species.clear(); + return prop; + } + + probe_info get_probe(cell_member_type id) const override { + // n probes, centred over CVs. + double pos = probe_x(id.index)/length; + return probe_info{id, 0, cell_probe_membrane_voltage{{0, pos}}}; + } + + util::unique_any get_cell_description(cell_gid_type) const override { + sample_tree samples({msample{{0., 0., 0., d0/2}, 0}, msample{{0., 0., length, d1/2}, 0}}, {mnpos, 0u}); + + cable_cell c(samples); + c.default_parameters.discretization = cv_policy_fixed_per_branch(n); + + mechanism_desc pas("pas"); + pas["g"] = 1e-4/rm; // [S/cm^2] + pas["e"] = 0; // erev=0 + + c.paint(reg::all(), pas); + c.place(mlocation{0, 1.}, i_clamp{0, INFINITY, iinj}); + return c; + } + + // time constant in [ms] + static double tau() { return rm*cm*1000; } + + // probe position in [µm] + double probe_x(unsigned i) const { + return length*((2.0*i+1.0)/(2.0*n)); + } +}; + +domain_decomposition trivial_dd(const recipe& r) { + cell_size_type ncell = r.num_cells(); + + std::vector all_gids(ncell); + std::iota(all_gids.begin(), all_gids.end(), cell_gid_type(0)); + + return domain_decomposition{ + [](cell_gid_type) { return 0; }, // gid_domain map + 1, // num_domains + 0, // domain_id + ncell, // num_local_cells + ncell, // num_global_cells + {{r.get_cell_kind(0), all_gids, backend_kind::multicore}} // groups + }; +} + +int main(int argc, char** argv) { + common_args A; + A.params = default_parameters; + parse_common_args(A, argc, argv, {"binevents"}); + + auto ctx = make_context(); + rc_cable_recipe rec(A.params); + simulation sim(rec, trivial_dd(rec), ctx); + + time_type dt = A.params["dt"]; + time_type t_end = 20*rec.tau(); // transient amplitudes are proportional to exp(-t/tau). + + std::vector voltage(rec.num_probes(0)); + std::vector x(rec.num_probes(0)); + for (unsigned i = 0; i(); + t_sample = rec[0].time; + }); + + sim.run(t_end, dt); + + // Write to netcdf: + + std::size_t vlen = voltage.size(); + + int ncid; + nc_check(nc_create, A.output.c_str(), 0, &ncid); + + int x_dimid, xid, vid, tid; + nc_check(nc_def_dim, ncid, "x", vlen, &x_dimid); + nc_check(nc_def_var, ncid, "x", NC_DOUBLE, 1, &x_dimid, &xid); + nc_check(nc_def_var, ncid, "v", NC_DOUBLE, 1, &x_dimid, &vid); + nc_check(nc_def_var, ncid, "t_sample", NC_DOUBLE, 0, nullptr, &tid); + + auto nc_put_att_cstr = [](int ncid, int varid, const char* name, const char* value) { + nc_check(nc_put_att_text, ncid, varid, name, std::strlen(value), value); + }; + + nc_put_att_cstr(ncid, xid, "units", "µm"); + nc_put_att_cstr(ncid, vid, "units", "mV"); + nc_put_att_cstr(ncid, tid, "units", "ms"); + + common_attr attrs; + attrs.model = "cable-steadystate"; + attrs.simulator = "arbor"; + attrs.simulator_build = arb::version; + attrs.simulator_build += ' '; + attrs.simulator_build += arb::source_id; + + set_common_attr(ncid, attrs, A.tags, A.params); + nc_check(nc_enddef, ncid); + + nc_check(nc_put_var_double, ncid, tid, &t_sample) + nc_check(nc_put_var_double, ncid, vid, voltage.data()) + nc_check(nc_put_var_double, ncid, xid, x.data()) + nc_check(nc_close, ncid); +} diff --git a/validation/src/arbor-rallpack1/BUILDFOR b/validation/src/arbor-rallpack1/BUILDFOR new file mode 100644 index 0000000..1b0a0a6 --- /dev/null +++ b/validation/src/arbor-rallpack1/BUILDFOR @@ -0,0 +1 @@ +arbor diff --git a/validation/src/arbor-rallpack1/CMakeLists.txt b/validation/src/arbor-rallpack1/CMakeLists.txt new file mode 100644 index 0000000..9219984 --- /dev/null +++ b/validation/src/arbor-rallpack1/CMakeLists.txt @@ -0,0 +1,17 @@ +cmake_minimum_required(VERSION 3.9) +project(arbor-rallpack1 LANGUAGES CXX) + +set(CMAKE_CXX_STANDARD 14) + +set(CMAKE_MODULE_PATH "${CMAKE_SOURCE_DIR}/../cmake") +include(UseLibraryPath) + +find_package(arbor REQUIRED) +find_package(netcdf REQUIRED) + +add_executable(arbor-rallpack1 arbor-rallpack1.cpp) +target_link_libraries(arbor-rallpack1 arbor::arbor netcdf::netcdf) +target_include_directories(arbor-rallpack1 PRIVATE ${CMAKE_SOURCE_DIR}/../include) + +install(TARGETS arbor-rallpack1 DESTINATION bin) + diff --git a/validation/src/arbor-rallpack1/arbor-rallpack1.cpp b/validation/src/arbor-rallpack1/arbor-rallpack1.cpp new file mode 100644 index 0000000..a0a4cf3 --- /dev/null +++ b/validation/src/arbor-rallpack1/arbor-rallpack1.cpp @@ -0,0 +1,183 @@ +#include +#include +#include +#include +#include + +#include + +#include +#include +#include +#include +#include +#include + +#include "common_args.h" +#include "common_attr.h" +#include "netcdf_wrap.h" + +using namespace arb; + +paramset default_parameters = { + {"dt", 0.0025}, + {"tend", 250.0}, + {"x0", 0.0}, + {"x1", 1.0}, + {"n", 1000} +}; + +struct rc_rallpack1_recipe: public arb::recipe { + // Fixed parameters: + + static constexpr double d = 1.0; // cable diameter [µm] + static constexpr double length = 1000.0; // cable length [µm] + static constexpr double rl = 1.0; // bulk resistivity [Ωm] + static constexpr double rm = 4.0; // membrane resistivity [Ωm²] + static constexpr double cm = 0.01; // membrane specific capacitance [F/m²] + static constexpr double erev = -65; // reversal potential [mV] + static constexpr double iinj = 0.1; // current injection [nA] + + // Customizable parameters: + unsigned n = 0; // number of CV + double x0 = 0, x1 = 0; // (proportional) sample locations + + explicit rc_rallpack1_recipe(const paramset& ps): + x0(ps.at("x0")), + x1(ps.at("x1")), + n(static_cast(ps.at("n"))) + {} + + cell_size_type num_cells() const override { return 1; } + cell_size_type num_targets(cell_gid_type) const override { return 0; } + cell_size_type num_probes(cell_gid_type) const override { return 2; } + cell_kind get_cell_kind(cell_gid_type) const override { return cell_kind::cable; } + + util::any get_global_properties(cell_kind kind) const override { + cable_cell_global_properties prop; + + prop.default_parameters.init_membrane_potential = (double)erev; + prop.default_parameters.axial_resistivity = 100*rl; // [Ω·cm] + prop.default_parameters.membrane_capacitance = (double)cm; // [F/m²] + prop.default_parameters.temperature_K = 0; + prop.ion_species.clear(); + return prop; + } + + probe_info get_probe(cell_member_type id) const override { + double pos = id.index==0? x0: x1; + return probe_info{id, 0, cell_probe_membrane_voltage{{0, pos}}}; + } + + util::unique_any get_cell_description(cell_gid_type) const override { + sample_tree samples({msample{{0., 0., 0., d/2}, 0}, msample{{0., 0., length, d/2}, 0}}, {mnpos, 0u}); + + cable_cell c(samples); + c.default_parameters.discretization = cv_policy_fixed_per_branch(n); + + mechanism_desc pas("pas"); + pas["g"] = 1e-4/rm; // [S/cm^2] + pas["e"] = (double)erev; + + c.paint(reg::all(), pas); + c.place(mlocation{0, 0}, i_clamp{0, INFINITY, iinj}); + return c; + } +}; + +domain_decomposition trivial_dd(const recipe& r) { + cell_size_type ncell = r.num_cells(); + + std::vector all_gids(ncell); + std::iota(all_gids.begin(), all_gids.end(), cell_gid_type(0)); + + return domain_decomposition{ + [](cell_gid_type) { return 0; }, // gid_domain map + 1, // num_domains + 0, // domain_id + ncell, // num_local_cells + ncell, // num_global_cells + {{r.get_cell_kind(0), all_gids, backend_kind::multicore}} // groups + }; +} + +int main(int argc, char** argv) { + common_args A; + A.params = default_parameters; + parse_common_args(A, argc, argv, {"binevents"}); + + auto ctx = make_context(); + rc_rallpack1_recipe rec(A.params); + simulation sim(rec, trivial_dd(rec), ctx); + + time_type dt = A.params["dt"]; + time_type t_end = A.params["tend"]; + time_type sample_dt = dt>0.01? dt: 0.01; // [ms] + + trace_data vtrace0, vtrace1; + sim.add_sampler(one_probe({0, 0}), regular_schedule(sample_dt), make_simple_sampler(vtrace0)); + sim.add_sampler(one_probe({0, 1}), regular_schedule(sample_dt), make_simple_sampler(vtrace1)); + + sim.run(t_end, dt); + + // Split sample times from voltages; assert times align for both traces. + + std::vector times, v0, v1; + + if (vtrace0.size()!=vtrace1.size()) { + fputs("sample time mismatch", stderr); + return 1; + } + + for (std::size_t i = 0; i0 && vtrace0[i].t==vtrace0[i-1].t) { + // Multiple sample in same integration timestep; discard. + continue; + } + + times.push_back(vtrace0[i].t); + v0.push_back(vtrace0[i].v); + v1.push_back(vtrace1[i].v); + } + + // Write to netcdf: + + std::size_t vlen = times.size(); + + int ncid; + nc_check(nc_create, A.output.c_str(), 0, &ncid); + + int time_dimid, timeid, v0id, v1id; + nc_check(nc_def_dim, ncid, "time", vlen, &time_dimid); + nc_check(nc_def_var, ncid, "time", NC_DOUBLE, 1, &time_dimid, &timeid); + nc_check(nc_def_var, ncid, "v0", NC_DOUBLE, 1, &time_dimid, &v0id); + nc_check(nc_def_var, ncid, "v1", NC_DOUBLE, 1, &time_dimid, &v1id); + + auto nc_put_att_cstr = [](int ncid, int varid, const char* name, const char* value) { + nc_check(nc_put_att_text, ncid, varid, name, std::strlen(value), value); + }; + + nc_put_att_cstr(ncid, timeid, "units", "ms"); + nc_put_att_cstr(ncid, v0id, "units", "mV"); + nc_put_att_cstr(ncid, v1id, "units", "mV"); + + common_attr attrs; + attrs.model = "rallpack1"; + attrs.simulator = "arbor"; + attrs.simulator_build = arb::version; + attrs.simulator_build += ' '; + attrs.simulator_build += arb::source_id; + + set_common_attr(ncid, attrs, A.tags, A.params); + nc_check(nc_enddef, ncid); + + nc_check(nc_put_var_double, ncid, timeid, times.data()) + nc_check(nc_put_var_double, ncid, v0id, v0.data()) + nc_check(nc_put_var_double, ncid, v1id, v1.data()) + nc_check(nc_close, ncid); +} diff --git a/validation/src/arbor-rc-exp2syn-spike/arbor-rc-exp2syn-spike.cpp b/validation/src/arbor-rc-exp2syn-spike/arbor-rc-exp2syn-spike.cpp index bc26054..fd5d2e1 100644 --- a/validation/src/arbor-rc-exp2syn-spike/arbor-rc-exp2syn-spike.cpp +++ b/validation/src/arbor-rc-exp2syn-spike/arbor-rc-exp2syn-spike.cpp @@ -49,8 +49,8 @@ struct rc_exp2syn_spike_recipe: public arb::recipe { // Computed values: std::vector delay; // delay[i] is connection delay from gid 0 to gid i - static segment_location soma_centre() { - return segment_location(0u, 0.5); + static mlocation soma_centre() { + return {0u, 0.5}; } explicit rc_exp2syn_spike_recipe(const paramset& ps): @@ -88,7 +88,7 @@ struct rc_exp2syn_spike_recipe: public arb::recipe { } probe_info get_probe(cell_member_type id) const override { - return probe_info{id, 0, cell_probe_address{soma_centre(), cell_probe_address::membrane_voltage}}; + return probe_info{id, 0, cell_probe_membrane_voltage{soma_centre()}}; } std::vector event_generators(cell_gid_type gid) const override { @@ -103,23 +103,29 @@ struct rc_exp2syn_spike_recipe: public arb::recipe { } util::unique_any get_cell_description(cell_gid_type) const override { - cable_cell c; + sample_tree samples; + samples.append({{0, 0, 0, r*1e6}, 1}); mechanism_desc pas("pas"); pas["g"] = 1e-10/(rm*area); // [S/cm^2] pas["e"] = erev; - auto soma = c.add_soma(r*1e6); - soma->parameters.membrane_capacitance = cm*1e-9/area; // [F/m^2] - soma->add_mechanism(pas); - mechanism_desc expsyn("exp2syn"); expsyn["tau1"] = tau1; expsyn["tau2"] = tau2; expsyn["e"] = 0; - c.add_synapse(soma_centre(), expsyn); - c.add_detector(soma_centre(), threshold); + label_dict labels; + labels.set("soma", reg::tagged(1)); + labels.set("centre", soma_centre()); + + cable_cell c(morphology(samples), labels); + c.default_parameters.membrane_capacitance = cm*1e-9/area; // [F/m^2] + + c.paint("soma", pas); + c.place("centre", expsyn); + c.place("centre", threshold_detector{threshold}); + return c; } diff --git a/validation/src/arbor-rc-expsyn/arbor-rc-expsyn.cpp b/validation/src/arbor-rc-expsyn/arbor-rc-expsyn.cpp index 7bfc679..3d878ce 100644 --- a/validation/src/arbor-rc-expsyn/arbor-rc-expsyn.cpp +++ b/validation/src/arbor-rc-expsyn/arbor-rc-expsyn.cpp @@ -39,8 +39,8 @@ struct rc_expsyn_recipe: public arb::recipe { // Customizable parameters: double g0; // synaptic conductance at time 0 [µS] - static segment_location soma_centre() { - return segment_location(0u, 0.5); + static mlocation soma_centre() { + return mlocation{0u, 0.5}; } explicit rc_expsyn_recipe(const paramset& ps): g0(ps.at("g0")) {} @@ -63,7 +63,7 @@ struct rc_expsyn_recipe: public arb::recipe { } probe_info get_probe(cell_member_type id) const override { - return probe_info{id, 0, cell_probe_address{soma_centre(), cell_probe_address::membrane_voltage}}; + return probe_info{id, 0, cell_probe_membrane_voltage{soma_centre()}}; } std::vector event_generators(cell_gid_type) const override { @@ -76,20 +76,26 @@ struct rc_expsyn_recipe: public arb::recipe { } util::unique_any get_cell_description(cell_gid_type) const override { - cable_cell c; + sample_tree samples; + samples.append({{0, 0, 0, r*1e6}, 1}); mechanism_desc pas("pas"); pas["g"] = 1e-10/(rm*area); // [S/cm^2] pas["e"] = erev; - auto soma = c.add_soma(r*1e6); - soma->parameters.membrane_capacitance = cm*1e-9/area; // [F/m^2] - soma->add_mechanism(pas); - mechanism_desc expsyn("expsyn"); expsyn["tau"] = syntau; expsyn["e"] = 0; - c.add_synapse(soma_centre(), expsyn); + + label_dict labels; + labels.set("soma", reg::tagged(1)); + labels.set("centre", soma_centre()); + + cable_cell c(morphology(samples), labels); + c.default_parameters.membrane_capacitance = cm*1e-9/area; // [F/m^2] + + c.paint("soma", pas); + c.place("centre", expsyn); return c; }