diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml new file mode 100644 index 0000000..393ddbd --- /dev/null +++ b/.github/workflows/main.yml @@ -0,0 +1,45 @@ +name: Test + +on: + push: + branches: [ master ] + pull_request: + branches: [ master ] + +jobs: + testallconfigs: + name: "Tests across OSes, versions, compilers, and build configs." + runs-on: ${{ matrix.config.os }} + strategy: + matrix: + config: + - { + name: "Linux Min GCC", + os: "ubuntu-20.04", + cc: "gcc-10", + cxx: "g++-10", + cmake: "3.22.x", + } + env: + CC: ${{ matrix.config.cc }} + CXX: ${{ matrix.config.cxx }} + steps: + - name: "Linux: get clang/gcc, libxml2" + if: ${{ startsWith(matrix.config.os, 'ubuntu') }} + run: | + sudo apt-get update + sudo apt-get install -y g++-10 + - name: Set up CMake + uses: jwlawson/actions-setup-cmake@v1.12 + with: + cmake-version: ${{ matrix.config.cmake }} + - name: Clone + uses: actions/checkout@v2 + with: + submodules: recursive + - name: Install + run: | + ./install-local.sh arbor --env=systems/gh.sh + - name: Run + run: | + ./run-bench.sh arbor --config=small diff --git a/benchmarks/engines/busyring/arbor/CMakeLists.txt b/benchmarks/engines/busyring/arbor/CMakeLists.txt index b90ee9b..49890c7 100644 --- a/benchmarks/engines/busyring/arbor/CMakeLists.txt +++ b/benchmarks/engines/busyring/arbor/CMakeLists.txt @@ -6,7 +6,9 @@ set (CMAKE_CXX_STANDARD 17) find_package(arbor REQUIRED) add_executable(ring ring.cpp) +target_link_libraries(ring PRIVATE ${CUDA_LIBRARIES}) target_link_libraries(ring PRIVATE arbor::arbor arbor::arborenv) + target_include_directories(ring PRIVATE ../../../../common/cpp/include) set_target_properties(ring PROPERTIES OUTPUT_NAME arbor-busyring) diff --git a/benchmarks/engines/busyring/arbor/parameters.hpp b/benchmarks/engines/busyring/arbor/parameters.hpp index 124ecf1..1f286a7 100644 --- a/benchmarks/engines/busyring/arbor/parameters.hpp +++ b/benchmarks/engines/busyring/arbor/parameters.hpp @@ -38,6 +38,7 @@ struct ring_params { double duration = 100; double dt = 0.025; bool record_voltage = false; + bool record_spikes = true; std::string odir = "."; cell_parameters cell; }; @@ -77,6 +78,7 @@ ring_params read_options(int argc, char** argv) { param_from_json(params.dt, "dt", json); param_from_json(params.min_delay, "min-delay", json); param_from_json(params.record_voltage, "record", json); + param_from_json(params.record_spikes, "spikes", json); param_from_json(params.cell.complex_cell, "complex", json); param_from_json(params.cell.max_depth, "depth", json); param_from_json(params.cell.branch_probs, "branch-probs", json); diff --git a/benchmarks/engines/busyring/arbor/ring.cpp b/benchmarks/engines/busyring/arbor/ring.cpp index d3a0bfa..10dedb4 100644 --- a/benchmarks/engines/busyring/arbor/ring.cpp +++ b/benchmarks/engines/busyring/arbor/ring.cpp @@ -1,6 +1,11 @@ #include #include #include +#include +#include +#include +#include +#include #include @@ -16,9 +21,11 @@ #include #include -#include +#include #include +#include + #include "parameters.hpp" #ifdef ARB_MPI_ENABLED @@ -34,6 +41,8 @@ using arb::cell_kind; using arb::time_type; using arb::cable_probe_membrane_voltage; +using namespace arborio::literals; + // Writes voltage trace as a json file. void write_trace_json(std::string fname, const arb::trace_data& trace); @@ -46,26 +55,24 @@ class ring_recipe: public arb::recipe { ring_recipe(ring_params params): num_cells_(params.num_cells), min_delay_(params.min_delay), - params_(params), - cat(arb::global_allen_catalogue()) + params_(params) { - cat.import(arb::global_default_catalogue(), ""); gprop.default_parameters = arb::neuron_parameter_defaults; - gprop.catalogue = &cat; + gprop.catalogue.import(arb::global_allen_catalogue(), ""); if (params.cell.complex_cell) { gprop.default_parameters.reversal_potential_method["ca"] = "nernst/ca"; gprop.default_parameters.axial_resistivity = 100; gprop.default_parameters.temperature_K = 34 + 273.15; gprop.default_parameters.init_membrane_potential = -90; - event_weight_*=5; + event_weight_ = 0.2; + event_freq_ = 0.2; } } - cell_size_type num_cells() const override { - return num_cells_; - } - + std::any get_global_properties(cell_kind kind) const override { return gprop; } + cell_size_type num_cells() const override { return num_cells_; } + cell_kind get_cell_kind(cell_gid_type gid) const override { return cell_kind::cable; } arb::util::unique_any get_cell_description(cell_gid_type gid) const override { if (params_.cell.complex_cell) { return complex_cell(gid, params_.cell); @@ -73,24 +80,6 @@ class ring_recipe: public arb::recipe { return branch_cell(gid, params_.cell); } - cell_kind get_cell_kind(cell_gid_type gid) const override { - return cell_kind::cable; - } - - // Each cell has one spike detector (at the soma). - cell_size_type num_sources(cell_gid_type gid) const override { - return 1; - } - - // The cell has one target synapse, which will be connected to cell gid-1. - cell_size_type num_targets(cell_gid_type gid) const override { - return params_.cell.synapses; - } - - std::any get_global_properties(cell_kind kind) const override { - return gprop; - } - // Each cell has one incoming connection, from cell with gid-1, // and fan_in-1 random connections with very low weight. std::vector connections_on(cell_gid_type gid) const override { @@ -103,7 +92,7 @@ class ring_recipe: public arb::recipe { const auto group_start = s*group; const auto group_end = std::min(group_start+s, num_cells_); cell_gid_type src = gid==group_start? group_end-1: gid-1; - cons.push_back(arb::cell_connection({src, 0}, {gid, 0}, event_weight_, min_delay_)); + cons.push_back(arb::cell_connection({src, "d"}, {"p"}, event_weight_, min_delay_)); // Used to pick source cell for a connection. std::uniform_int_distribution dist(0, num_cells_-2); @@ -116,24 +105,20 @@ class ring_recipe: public arb::recipe { src = dist(src_gen); if (src==gid) ++src; const float delay = min_delay_+delay_dist(src_gen); - //const float delay = min_delay_; cons.push_back( - arb::cell_connection({src, 0}, {gid, i}, 0.f, delay)); + arb::cell_connection({src, "d"}, {"p"}, 0.f, delay)); } - return cons; } // Return one event generator on the first cell of each ring. // This generates a single event that will kick start the spiking on the sub-ring. std::vector event_generators(cell_gid_type gid) const override { - std::vector gens; if (gid%params_.ring_size == 0) { - gens.push_back( - arb::explicit_generator( - arb::pse_vector{{{gid, 0}, 1.0, event_weight_}})); + return {arb::regular_generator({"p"}, event_weight_, 0.0, event_freq_)}; + } else { + return {}; } - return gens; } std::vector get_probes(cell_gid_type gid) const override { @@ -146,10 +131,10 @@ class ring_recipe: public arb::recipe { cell_size_type num_cells_; double min_delay_; ring_params params_; - float event_weight_ = 0.01; + float event_weight_ = 0.025; + float event_freq_ = 1.0; arb::cable_cell_global_properties gprop; - arb::mechanism_catalogue cat; }; struct cell_stats { @@ -203,12 +188,7 @@ int main(int argc, char** argv) { auto params = read_options(argc, argv); arb::proc_allocation resources; - if (auto nt = arbenv::get_env_num_threads()) { - resources.num_threads = nt; - } - else { - resources.num_threads = arbenv::thread_concurrency(); - } + resources.num_threads = arbenv::default_concurrency(); #ifdef ARB_MPI_ENABLED arbenv::with_mpi guard(argc, argv, false); @@ -240,13 +220,8 @@ int main(int argc, char** argv) { cell_stats stats(recipe); if (root) std::cout << stats << "\n"; - //arb::partition_hint_map hints; - //hints[cell_kind::cable1d_neuron].cpu_group_size = 4; - //auto decomp = arb::partition_load_balance(recipe, context, hints); - auto decomp = arb::partition_load_balance(recipe, context); - // Construct the model. - arb::simulation sim(recipe, decomp, context); + arb::simulation sim(recipe, context); // Set up the probe that will measure voltage in the cell. @@ -264,7 +239,7 @@ int main(int argc, char** argv) { // Set up recording of spikes to a vector on the root process. std::vector recorded_spikes; - if (root) { + if (root && params.record_spikes) { sim.set_global_spike_callback( [&recorded_spikes](const std::vector& spikes) { recorded_spikes.insert(recorded_spikes.end(), spikes.begin(), spikes.end()); @@ -286,17 +261,19 @@ int main(int argc, char** argv) { if (root) { std::cout << "\n" << ns << " spikes generated at rate of " << params.duration/ns << " ms between spikes\n"; - std::ofstream fid(params.odir + "/" + params.name + "_spikes.gdf"); - if (!fid.good()) { - std::cerr << "Warning: unable to open file spikes.gdf for spike output\n"; - } - else { - char linebuf[45]; - for (auto spike: recorded_spikes) { - auto n = std::snprintf( - linebuf, sizeof(linebuf), "%u %.4f\n", - unsigned{spike.source.gid}, float(spike.time)); - fid.write(linebuf, n); + if (!recorded_spikes.empty()) { + std::ofstream fid(params.odir + "/" + params.name + "_spikes.gdf"); + if (!fid.good()) { + std::cerr << "Warning: unable to open file spikes.gdf for spike output\n"; + } + else { + char linebuf[45]; + for (auto spike: recorded_spikes) { + auto n = std::snprintf( + linebuf, sizeof(linebuf), "%u %.4f\n", + unsigned{spike.source.gid}, float(spike.time)); + fid.write(linebuf, n); + } } } } @@ -400,53 +377,51 @@ arb::cable_cell complex_cell (arb::cell_gid_type gid, const cell_parameters& par dist_from_soma += l; } - arb::label_dict dict; - - dict.set("soma", tagged(1)); - dict.set("axon", tagged(2)); - dict.set("dend", tagged(3)); - dict.set("apic", tagged(4)); - dict.set("center", location(0, 0.5)); - if (params.synapses>1) { - dict.set("synapses", arb::ls::uniform(arb::reg::all(), 0, params.synapses-2, gid)); - } + using arb::reg::tagged; + auto rall = arb::reg::all(); + auto soma = tagged(1); + auto axon = tagged(2); + auto dend = tagged(3); + auto apic = tagged(4); + auto cntr = location(0, 0.5); + auto syns = arb::ls::uniform(rall, 0, params.synapses-2, gid); arb::decor decor; - decor.paint(all(), arb::init_reversal_potential{"k", -107.0}); - decor.paint(all(), arb::init_reversal_potential{"na", 53.0}); + decor.paint(rall, arb::init_reversal_potential{"k", -107.0}); + decor.paint(rall, arb::init_reversal_potential{"na", 53.0}); - decor.paint("\"soma\"", arb::axial_resistivity{133.577}); - decor.paint("\"soma\"", arb::membrane_capacitance{4.21567e-2}); + decor.paint(soma, arb::axial_resistivity{133.577}); + decor.paint(soma, arb::membrane_capacitance{4.21567e-2}); - decor.paint("\"dend\"", arb::axial_resistivity{68.355}); - decor.paint("\"dend\"", arb::membrane_capacitance{2.11248e-2}); + decor.paint(dend, arb::axial_resistivity{68.355}); + decor.paint(dend, arb::membrane_capacitance{2.11248e-2}); - decor.paint("\"soma\"", mech("pas").set("g", 0.000119174).set("e", -76.4024)); - decor.paint("\"soma\"", mech("NaV").set("gbar", 0.0499779)); - decor.paint("\"soma\"", mech("SK").set("gbar", 0.000733676)); - decor.paint("\"soma\"", mech("Kv3_1").set("gbar", 0.186718)); - decor.paint("\"soma\"", mech("Ca_HVA").set("gbar", 9.96973e-05)); - decor.paint("\"soma\"", mech("Ca_LVA").set("gbar", 0.00344495)); - decor.paint("\"soma\"", mech("CaDynamics").set("gamma", 0.0177038).set("decay", 42.2507)); - decor.paint("\"soma\"", mech("Ih").set("gbar", 1.07608e-07)); + decor.paint(soma, arb::density("pas/e=-76.4024", {{"g", 0.000119174}})); + decor.paint(soma, arb::density("NaV", {{"gbar", 0.0499779}})); + decor.paint(soma, arb::density("SK", {{"gbar", 0.000733676}})); + decor.paint(soma, arb::density("Kv3_1", {{"gbar", 0.186718}})); + decor.paint(soma, arb::density("Ca_HVA", {{"gbar", 9.96973e-05}})); + decor.paint(soma, arb::density("Ca_LVA", {{"gbar", 0.00344495}})); + decor.paint(soma, arb::density("CaDynamics", {{"gamma", 0.0177038}, {"decay", 42.2507}})); + decor.paint(soma, arb::density("Ih", {{"gbar", 1.07608e-07}})); - decor.paint("\"dend\"", mech("pas").set("g", 9.57001e-05).set("e", -88.2554)); - decor.paint("\"dend\"", mech("NaV").set("gbar", 0.0472215)); - decor.paint("\"dend\"", mech("Kv3_1").set("gbar", 0.186859)); - decor.paint("\"dend\"", mech("Im_v2").set("gbar", 0.00132163)); - decor.paint("\"dend\"", mech("Ih").set("gbar", 9.18815e-06)); + decor.paint(dend, arb::density("pas/e=-88.2554", {{"g", 9.57001e-05}})); + decor.paint(dend, arb::density("NaV", {{"gbar", 0.0472215}})); + decor.paint(dend, arb::density("Kv3_1", {{"gbar", 0.186859}})); + decor.paint(dend, arb::density("Im_v2", {{"gbar", 0.00132163}})); + decor.paint(dend, arb::density("Ih", {{"gbar", 9.18815e-06}})); - decor.place("\"center\"", mech("expsyn")); + decor.place(cntr, arb::synapse("expsyn"), "p"); if (params.synapses>1) { - decor.place("\"synapses\"", "expsyn"); + // decor.place(syns, arb::synapse("expsyn"), "s"); } - decor.place("\"center\"", arb::threshold_detector{-20.0}); + decor.place(cntr, arb::threshold_detector{-20.0}, "d"); decor.set_default(arb::cv_policy_every_segment()); - return {arb::morphology(tree), dict, decor}; + return {arb::morphology(tree), decor}; } arb::cable_cell branch_cell(arb::cell_gid_type gid, const cell_parameters& params) { @@ -498,34 +473,30 @@ arb::cable_cell branch_cell(arb::cell_gid_type gid, const cell_parameters& param dist_from_soma += l; } - arb::label_dict labels; - using arb::reg::tagged; - labels.set("soma", tagged(1)); - labels.set("dendrites", join(tagged(3), tagged(4))); - if (params.synapses>1) { - labels.set("synapses", arb::ls::uniform(arb::reg::all(), 0, params.synapses-2, gid)); - } + auto soma = tagged(1); + auto dnds = join(tagged(3), tagged(4)); + auto syns = arb::ls::uniform(arb::reg::all(), 0, params.synapses-2, gid); arb::decor decor; - decor.paint("\"soma\"", "hh"); - decor.paint("\"dendrites\"", "pas"); + decor.paint(soma, arb::density{"hh"}); + decor.paint(dnds, arb::density{"pas"}); decor.set_default(arb::axial_resistivity{100}); // [Ω·cm] // Add spike threshold detector at the soma. - decor.place(arb::mlocation{0,0}, arb::threshold_detector{10}); + decor.place(arb::mlocation{0,0}, arb::threshold_detector{10}, "d"); // Add a synapse to proximal end of first dendrite. - decor.place(arb::mlocation{1, 0}, "expsyn"); + decor.place(arb::mlocation{1, 0}, arb::synapse{"expsyn"}, "p"); // Add additional synapses that will not be connected to anything. if (params.synapses>1) { - decor.place("\"synapses\"", "expsyn"); + decor.place(syns, arb::synapse{"expsyn"}, "s"); } // Make a CV between every sample in the sample tree. decor.set_default(arb::cv_policy_every_segment()); - return {arb::morphology(tree), labels, decor}; + return {arb::morphology(tree), decor}; } diff --git a/common/bin/comparex b/common/bin/comparex index 5704e93..f0fe4db 100755 --- a/common/bin/comparex +++ b/common/bin/comparex @@ -8,7 +8,7 @@ import sys import numpy as np import xarray import scipy.interpolate as interpolate -import scipy.ndimage.filters as filters +from scipy.ndimage import maximum_filter1d def parse_clargs(): P = argparse.ArgumentParser() @@ -92,10 +92,10 @@ def interpolate_array(x, t, tnew): x4_est[0] *= 3 # end-interval fudge factors x4_est[-1] *= 3 - x4_windowmax = np.pad(filters.maximum_filter1d(x4_est, 4, origin=-2), (1, 0), mode='edge') + x4_windowmax = np.pad(maximum_filter1d(x4_est, 4, origin=-2), (1, 0), mode='edge') local_x4max = interpolate.interp1d(t, np.resize(x4_windowmax, x.size), assume_sorted=True, kind='previous', fill_value='extrapolate') - dt_windowmax = np.pad(filters.maximum_filter1d(np.ediff1d(t), 3, origin=-1), (1, 0), mode='edge') + dt_windowmax = np.pad(maximum_filter1d(np.ediff1d(t), 3, origin=-1), (1, 0), mode='edge') local_dtmax = interpolate.interp1d(t, np.resize(dt_windowmax, x.size), assume_sorted=True, kind='previous', fill_value='extrapolate') xi_err = 5./384.*pow(local_dtmax(tnew), 4)*local_x4max(tnew) diff --git a/common/bin/thresholdx b/common/bin/thresholdx index 1775abe..60bb5bc 100755 --- a/common/bin/thresholdx +++ b/common/bin/thresholdx @@ -63,17 +63,17 @@ for var, op, value in opts.exprs: else: v = opts.input[var] if op=='<': - status = np.asscalar(np.all(v': - status = np.asscalar(np.all(v>value)) + status = (np.all(v>value)).item() elif op=='>=': - status = np.asscalar(np.all(v>=value)) + status = (np.all(v>=value)).item() elif op=='==' or op=='=': - status = np.asscalar(np.all(v==value)) + status = (np.all(v==value)).item() elif op=='!=': - status = np.asscalar(np.all(v!=value)) + status = (np.all(v!=value)).item() else: status = 'unknown operation' diff --git a/scripts/build_arbor.sh b/scripts/build_arbor.sh index 74d01f1..4454ca6 100644 --- a/scripts/build_arbor.sh +++ b/scripts/build_arbor.sh @@ -48,7 +48,7 @@ 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_USE_BUNDLED_LIBS=ON" +cmake_args="$cmake_args -DARB_USE_BUNDLED_LIBS=ON -DCMAKE_INTERPROCEDURAL_OPTIMIZATION=OFF" cmake_args="$cmake_args -DARB_GPU=$ns_arb_gpu" cmake_args="$cmake_args -DARB_ARCH=$ns_arb_arch" cmake_args="$cmake_args -DARB_VECTORIZE=$ns_arb_vectorize" diff --git a/systems/gh.sh b/systems/gh.sh new file mode 100644 index 0000000..f3f483a --- /dev/null +++ b/systems/gh.sh @@ -0,0 +1,31 @@ +### environment ### + +# record system name +ns_sysname="GH" + +### compilation options ### +# -DPYTHON_EXECUTABLE=/usr/local/bin/python3.8 +ns_cc="gcc" +ns_cxx="g++" +ns_with_mpi=OFF +ns_arb_arch=native +ns_arb_gpu=none + +ns_arb_git_repo='https://github.com/arbor-sim/arbor.git' +ns_arb_branch='master' + +ns_makej=4 +ns_validate=disable +### benchmark execution options ### +ns_threads_per_core=2 +ns_cores_per_socket=12 +ns_sockets=1 +ns_threads_per_socket=2 + + +run_with_mpi() { + export ARB_NUM_THREADS=$ns_threads_per_socket + export OMP_NUM_THREADS=$ns_threads_per_socket + echo ${@} + ${@} +} diff --git a/systems/juwels-booster.sh b/systems/juwels-booster.sh new file mode 100644 index 0000000..28855e6 --- /dev/null +++ b/systems/juwels-booster.sh @@ -0,0 +1,36 @@ +### environment ### + +# As of March 2022 + +# record system name +ns_sysname="juwels-booster" + +# set up environment for building on the booster part of juwels +ml CMake/3.21.1 Python/3.9.6 NVHPC/22.1 ParaStationMPI/5.5.0-1 +ml Python/3.9.6 mpi4py/3.1.3 + +ns_validate=disable + +ns_python=$(which python3) +ns_cc=$(which mpicc) +ns_cxx=$(which mpicxx) +ns_with_mpi=ON + +ns_arb_gpu=cuda +ns_arb_vectorize=ON +ns_arb_arch=native + +ns_arb_branch=master + +ns_makej=20 + +ns_threads_per_core=2 +ns_cores_per_socket=24 +ns_sockets=2 +ns_threads_per_socket=24 + +# activate budget via jutil env activate -p -A before running the benchmark +run_with_mpi() { + echo srun "${@}" + srun "${@}" +} diff --git a/systems/juwels-gpu.sh b/systems/juwels-gpu.sh index ceb9c17..625144f 100644 --- a/systems/juwels-gpu.sh +++ b/systems/juwels-gpu.sh @@ -3,47 +3,29 @@ # record system name ns_sysname="juwels-gpu" -# set up environment for building on the multicore part of juwels -module purge -module use /gpfs/software/juwels/otherstages/ -module load Stages/2019a +# set up environment for building on the booster part of juwels +ml CMake/3.21.1 Python/3.9.6 NVHPC/22.1 ParaStationMPI/5.5.0-1 +ml Python/3.9.6 mpi4py/3.1.3 -module load GCC -module load MVAPICH2 -module load CMake +ns_validate=disable -module load CUDA -export MV2_ENABLE_AFFINITY=0 -export MV2_USE_GPUDIRECT_GDRCOPY=0 - -module load Python/3.6.8 -module load SciPy-Stack/2019a-Python-3.6.8 ns_python=$(which python3) - -# for (core)neuron -module load mpi4py/3.0.1-Python-3.6.8 -module load flex/2.6.4 - -# for validation tests -module load netCDF/4.6.3-serial - -### compilation options ### ns_cc=$(which mpicc) ns_cxx=$(which mpicxx) ns_with_mpi=ON ns_arb_gpu=cuda -ns_arb_arch=skylake-avx512 +ns_arb_vectorize=ON +ns_arb_arch=native ns_arb_branch=master ns_makej=20 -### benchmark execution options ### ns_threads_per_core=2 -ns_cores_per_socket=20 -ns_sockets=1 -ns_threads_per_socket=20 +ns_cores_per_socket=24 +ns_sockets=2 +ns_threads_per_socket=24 # activate budget via jutil env activate -p -A before running the benchmark run_with_mpi() { diff --git a/validation/rc-exp2syn-spike/generate-rc-exp2syn-spike b/validation/rc-exp2syn-spike/generate-rc-exp2syn-spike index 03ad389..287ec39 100755 --- a/validation/rc-exp2syn-spike/generate-rc-exp2syn-spike +++ b/validation/rc-exp2syn-spike/generate-rc-exp2syn-spike @@ -67,7 +67,7 @@ def run_integration(): crossing.direction = 1 r = integrate.solve_ivp(dv_dt, (0., tend), [Erev], method='LSODA', t_eval=ts, jac=jacobian, atol=1e-10/nsamp, rtol=1e-10/nsamp, events=crossing) - spike = np.asscalar(r.t_events[0]) if r.t_events[0].size>0 else np.NaN + spike = (r.t_events[0]).item() if r.t_events[0].size>0 else np.NaN return (r.y[0], spike) diff --git a/validation/src/arbor-cable-steadystate/arbor-cable-steadystate.cpp b/validation/src/arbor-cable-steadystate/arbor-cable-steadystate.cpp index d7144ad..c4be69d 100644 --- a/validation/src/arbor-cable-steadystate/arbor-cable-steadystate.cpp +++ b/validation/src/arbor-cable-steadystate/arbor-cable-steadystate.cpp @@ -14,6 +14,8 @@ #include #include #include +#include +#include #include "common_args.h" #include "common_attr.h" @@ -48,9 +50,18 @@ struct rc_cable_recipe: public arb::recipe { {} cell_size_type num_cells() const override { return 1; } - cell_size_type num_targets(cell_gid_type) const override { return 0; } cell_kind get_cell_kind(cell_gid_type) const override { return cell_kind::cable; } - cell_size_type num_probes(cell_gid_type) const { return n; } + cell_size_type num_cvs_per_branch(cell_gid_type) const { return n; } + std::vector get_probes(cell_gid_type gid) const override { + std::vector probes; + + // n probes, centred over CVs. + for (unsigned i = 0; i < num_probes(gid); ++i) { + arb::mlocation loc{0, probe_x(i)/length}; + probes.push_back(cable_probe_membrane_voltage{loc}); + } + return probes; + } std::any get_global_properties(cell_kind kind) const override { cable_cell_global_properties prop; @@ -58,87 +69,58 @@ struct rc_cable_recipe: public arb::recipe { 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.default_parameters.temperature_K = 300; prop.ion_species.clear(); return prop; } - std::vector get_probes(cell_gid_type gid) const override { - std::vector probes; - - // n probes, centred over CVs. - for (unsigned i = 0; i < num_probes(gid); ++i) { - arb::mlocation loc{0, probe_x(i)/length}; - probes.push_back(cable_probe_membrane_voltage{loc}); - } - return probes; - } util::unique_any get_cell_description(cell_gid_type) const override { segment_tree tree; tree.append(arb::mnpos, {0., 0., 0., d0/2}, {0., 0., length, d1/2}, 0); - mechanism_desc pas("pas"); + mechanism_desc pas("pas/e=0"); pas["g"] = 1e-4/rm; // [S/cm^2] - pas["e"] = 0; // erev=0 + // pas["e"] = 0; // erev=0 decor D; - D.paint(reg::all(), pas); - D.place(mlocation{0, 1.}, i_clamp{0, INFINITY, iinj}); + D.paint(reg::all(), arb::density{pas}); + D.place(mlocation{0, 1.}, i_clamp{iinj}, "ic"); D.set_default(cv_policy_fixed_per_branch(n)); - return cable_cell(tree, {}, D); + return cable_cell(tree, D); } // 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)); - } - + double probe_x(unsigned i) const { return length*((2.0*i+1.0)/(2.0*n)); } + cell_size_type num_probes(cell_gid_type) const { return 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); + simulation sim(rec); 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)); + auto n_entities = rec.num_cvs_per_branch(0); + std::vector voltage(n_entities); + std::vector x(n_entities); for (unsigned i = 0; i(); t_sample = rec[0].time; - }; - - sim.add_sampler(all_probes, explicit_schedule({t_end-dt}), sampler, sampling_policy::exact); + }; + sim.add_sampler(all_probes, explicit_schedule({t_end-dt}), sampler); sim.run(t_end, dt); diff --git a/validation/src/arbor-rallpack1/arbor-rallpack1.cpp b/validation/src/arbor-rallpack1/arbor-rallpack1.cpp index 7591251..48cef79 100644 --- a/validation/src/arbor-rallpack1/arbor-rallpack1.cpp +++ b/validation/src/arbor-rallpack1/arbor-rallpack1.cpp @@ -49,7 +49,6 @@ struct rc_rallpack1_recipe: public arb::recipe { {} cell_size_type num_cells() const override { return 1; } - cell_size_type num_targets(cell_gid_type) const override { return 0; } cell_kind get_cell_kind(cell_gid_type) const override { return cell_kind::cable; } cell_size_type num_probes(cell_gid_type) const { return 2; } @@ -79,35 +78,19 @@ struct rc_rallpack1_recipe: public arb::recipe { segment_tree tree; tree.append(arb::mnpos, {0., 0., 0., d/2}, {0., 0., length, d/2}, 0); - mechanism_desc pas("pas"); + mechanism_desc pas("pas/e=-65"); pas["g"] = 1e-4/rm; // [S/cm^2] - pas["e"] = (double)erev; + // pas["e"] = (double)erev; decor D; - D.paint(reg::all(), pas); - D.place(mlocation{0, 0}, i_clamp{0, INFINITY, iinj}); + D.paint(reg::all(), arb::density{pas}); + D.place(mlocation{0, 0}, i_clamp{0, INFINITY, iinj}, "ic"); D.set_default(cv_policy_fixed_per_branch(n)); - return cable_cell(tree, {}, D); + return cable_cell(tree, D); } }; -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; @@ -115,7 +98,7 @@ int main(int argc, char** argv) { auto ctx = make_context(); rc_rallpack1_recipe rec(A.params); - simulation sim(rec, trivial_dd(rec), ctx); + simulation sim(rec); time_type dt = A.params["dt"]; time_type t_end = A.params["tend"]; 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 e13c638..e80c4c8 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 @@ -12,12 +12,14 @@ #include #include #include +#include #include "common_args.h" #include "common_attr.h" #include "netcdf_wrap.h" using namespace arb; +using namespace arborio::literals; paramset default_parameters = { {"dt", 0.01}, @@ -64,11 +66,7 @@ struct rc_exp2syn_spike_recipe: public arb::recipe { } cell_size_type num_cells() const override { return ncell; } - cell_size_type num_sources(cell_gid_type) const override { return 1; } - cell_size_type num_targets(cell_gid_type) const override { return 1; } - cell_size_type num_probes(cell_gid_type gid) const { - return gid==0? 1: 0; - } + cell_size_type num_probes(cell_gid_type gid) const { return gid==0? 1: 0; } cell_kind get_cell_kind(cell_gid_type) const override { return cell_kind::cable; } std::any get_global_properties(cell_kind kind) const override { @@ -89,22 +87,18 @@ struct rc_exp2syn_spike_recipe: public arb::recipe { std::vector event_generators(cell_gid_type gid) const override { if (gid!=0) return {}; - - spike_event ev; - ev.target = {0u, 0u}; - ev.time = 0; - ev.weight = g0; - - return {explicit_generator(pse_vector{{ev}})}; + return {explicit_generator(arb::cell_local_label_type("syn"), + g0, + std::vector{0.0f})}; } util::unique_any get_cell_description(cell_gid_type) const override { segment_tree tree; tree.append(arb::mnpos, {0., 0., 0., r*1e6}, {0., 0., 2*r*1e6, r*1e6}, 1); - mechanism_desc pas("pas"); + mechanism_desc pas("pas/e=-65"); pas["g"] = 1e-10/(rm*area); // [S/cm^2] - pas["e"] = erev; + // pas["e"] = erev; mechanism_desc exp2syn("exp2syn"); exp2syn["tau1"] = tau1; @@ -117,44 +111,30 @@ struct rc_exp2syn_spike_recipe: public arb::recipe { decor D; D.set_default(membrane_capacitance{cm*1e-9/area}); // [F/m^2] - D.paint("\"soma\"", pas); - D.place("\"midpoint\"", exp2syn); - D.place("\"midpoint\"", threshold_detector{threshold}); + D.paint(arb::reg::named("soma"), arb::density{pas}); + D.place(arb::ls::named("midpoint"), arb::synapse{exp2syn}, "syn"); + D.place(arb::ls::named("midpoint"), threshold_detector{threshold}, "det"); - return cable_cell(tree, labels, D); + return cable_cell(tree, D, labels); } std::vector connections_on(cell_gid_type gid) const override { if (gid==0) return {}; - return {arb::cell_connection({0, 0}, {gid, 0}, g0, delay[gid])}; + return {arb::cell_connection(arb::cell_global_label_type{0, "det"}, + arb::cell_local_label_type{"syn"}, + g0, + delay[gid])}; } }; -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_exp2syn_spike_recipe rec(A.params); - simulation sim(rec, trivial_dd(rec), ctx); + simulation sim(rec); time_type t_end = 10., sample_dt = 0.05; // [ms] time_type dt = A.params["dt"]; diff --git a/validation/src/arbor-rc-expsyn/arbor-rc-expsyn.cpp b/validation/src/arbor-rc-expsyn/arbor-rc-expsyn.cpp index 72573f8..66ddd2f 100644 --- a/validation/src/arbor-rc-expsyn/arbor-rc-expsyn.cpp +++ b/validation/src/arbor-rc-expsyn/arbor-rc-expsyn.cpp @@ -42,7 +42,6 @@ struct rc_expsyn_recipe: public arb::recipe { explicit rc_expsyn_recipe(const paramset& ps): g0(ps.at("g0")) {} cell_size_type num_cells() const override { return 1; } - cell_size_type num_targets(cell_gid_type) const override { return 1; } cell_kind get_cell_kind(cell_gid_type) const override { return cell_kind::cable; } cell_size_type num_probes(cell_gid_type) const { return 1; } @@ -62,22 +61,19 @@ struct rc_expsyn_recipe: public arb::recipe { return {cable_probe_membrane_voltage{ls::named("midpoint")}}; } - std::vector event_generators(cell_gid_type) const override { - spike_event ev; - ev.target = {0u, 0u}; - ev.time = 0; - ev.weight = g0; - - return {explicit_generator(pse_vector{{ev}})}; + std::vector event_generators(cell_gid_type gid) const override { + return {explicit_generator(arb::cell_local_label_type("syn"), + g0, + std::vector{0.0f})}; } util::unique_any get_cell_description(cell_gid_type) const override { segment_tree tree; tree.append(arb::mnpos, {0., 0., 0., r*1e6}, {0., 0., 2*r*1e6, r*1e6}, 1); - mechanism_desc pas("pas"); + mechanism_desc pas("pas/e=-65"); pas["g"] = 1e-10/(rm*area); // [S/cm^2] - pas["e"] = erev; + // pas["e"] = erev; mechanism_desc expsyn("expsyn"); expsyn["tau"] = syntau; @@ -89,37 +85,20 @@ struct rc_expsyn_recipe: public arb::recipe { decor D; D.set_default(membrane_capacitance{cm*1e-9/area}); // [F/m^2] - D.paint("\"soma\"", pas); - D.place("\"midpoint\"", expsyn); + D.paint(arb::reg::named("soma"), arb::density{pas}); + D.place(arb::ls::named("midpoint"), arb::synapse{expsyn}, "syn"); - return cable_cell(tree, labels, D); + return cable_cell(tree, D, labels); } }; -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_expsyn_recipe rec(A.params); - simulation sim(rec, trivial_dd(rec), ctx); + simulation sim(rec); time_type t_end = 10., sample_dt = 0.05; // [ms] time_type dt = A.params["dt"];