diff --git a/.travis.yml b/.travis.yml
index 9f9d3effab76d1da68a2cfc5e0100016c4949ab4..31ca2dbd28de07594d9d42c5b0a92fce5c6d5d04 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -1,8 +1,7 @@
######## Testing minimal compiler requirements ########
# GCC 6.4.0
-# Clang 4.0
-# Apple Clang 900.0.39.2
-# Python 3.6
+# Clang 7.0
+# Apple Clang 1100.0.33.16
#######################################################
language: cpp
@@ -14,10 +13,10 @@ matrix:
## test gcc6 - single node/rank with threading backend ##
- name: "osx, gcc, serial, py"
os: osx
- osx_image: xcode9.2
+ osx_image: xcode11.3
python: 3.6
env:
- - MATRIX_EVAL="brew install gcc@6 && brew link --force --overwrite gcc@6 && CC=gcc-6 && CXX=g++-6"
+ - MATRIX_EVAL="brew install gcc@6 && brew link --force --overwrite gcc@6 && brew install cmake && CC=gcc-6 && CXX=g++-6"
- BUILD_NAME=cthread-osx-gcc-py
- WITH_DISTRIBUTED=serial WITH_PYTHON=true PY=3
compiler: gcc-6
@@ -25,10 +24,10 @@ matrix:
## test gcc6 - mpi with threading backend ##
- name: "osx, gcc, mpi, py"
os: osx
- osx_image: xcode9.2
+ osx_image: xcode11.3
python: 3.6
env:
- - MATRIX_EVAL="brew install gcc@6 && brew link --force --overwrite gcc@6 && CC=gcc-6 && CXX=g++-6"
+ - MATRIX_EVAL="brew install gcc@6 && brew link --force --overwrite gcc@6 && brew install cmake && CC=gcc-6 && CXX=g++-6"
- BUILD_NAME=mpi-osx-gcc-py
- WITH_DISTRIBUTED=mpi WITH_PYTHON=true PY=3
compiler: gcc-6
@@ -36,7 +35,7 @@ matrix:
## test clang9 - single node/rank with threading backend ##
- name: "osx, apple clang, serial, py"
os: osx
- osx_image: xcode9.2
+ osx_image: xcode11.3
python: 3.6
env:
- MATRIX_EVAL="CC=clang && CXX=clang++"
@@ -47,7 +46,7 @@ matrix:
## test clang9 - mpi with threading backend ##
- name: "osx, apple clang, mpi, py"
os: osx
- osx_image: xcode9.2
+ osx_image: xcode11.3
python: 3.6
env:
- MATRIX_EVAL="CC=clang && CXX=clang++"
@@ -59,12 +58,10 @@ matrix:
## test gcc6 - single node/rank with threading backend ##
- name: "linux, gcc, serial, py"
os: linux
- dist: trusty
- python: 3.6
+ dist: bionic
addons:
apt:
sources:
- - ubuntu-toolchain-r-test
packages:
- g++-6
- openmpi-bin
@@ -78,12 +75,10 @@ matrix:
## test gcc6 - mpi with threading backend ##
- name: "linux, gcc, mpi, py"
os: linux
- dist: trusty
- python: 3.6
+ dist: bionic
addons:
apt:
sources:
- - ubuntu-toolchain-r-test
packages:
- g++-6
- openmpi-bin
@@ -97,20 +92,15 @@ matrix:
## test clang4 - single node/rank with threading backend ##
- name: "linux, clang, serial, py"
os: linux
- dist: trusty
- python: 3.6
+ dist: bionic
addons:
apt:
sources:
- - ubuntu-toolchain-r-test
- - llvm-toolchain-trusty-4.0
packages:
- - clang-4.0
- - g++-6
- openmpi-bin
- libopenmpi-dev
env:
- - MATRIX_EVAL="CC=clang-4.0 && CXX=clang++-4.0"
+ - MATRIX_EVAL="CC=clang && CXX=clang++"
- BUILD_NAME=cthread-linux-clang-py
- WITH_DISTRIBUTED=serial WITH_PYTHON=true PY=3
compiler: clang-4.0
@@ -118,36 +108,21 @@ matrix:
## test clang4 - mpi with threading backend ##
- name: "linux, clang, mpi, py"
os: linux
- dist: trusty
- python: 3.6
+ dist: bionic
addons:
apt:
sources:
- - ubuntu-toolchain-r-test
- - llvm-toolchain-trusty-4.0
packages:
- - clang-4.0
- - g++-6
- openmpi-bin
- libopenmpi-dev
env:
- - MATRIX_EVAL="CC=clang-4.0 && CXX=clang++-4.0"
+ - MATRIX_EVAL="CC=clang && CXX=clang++"
- BUILD_NAME=mpi-linux-clang-py
- WITH_DISTRIBUTED=mpi WITH_PYTHON=true PY=3
compiler: clang-4.0
before_install:
- if [[ "$TRAVIS_OS_NAME" == "osx" ]]; then export HOMEBREW_NO_AUTO_UPDATE=1; brew cask uninstall --force oclint; fi
- - |
- if [[ "$WITH_PYTHON" == "true" ]]; then
- if [[ "$TRAVIS_OS_NAME" == "linux" ]]; then
- source ~/virtualenv/python3.6/bin/activate
- elif [[ "$TRAVIS_OS_NAME" == "osx" ]]; then
- brew unlink python
- brew install https://raw.githubusercontent.com/Homebrew/homebrew-core/e128fa1bce3377de32cbf11bd8e46f7334dfd7a6/Formula/python.rb
- brew switch python 3.6.5
- fi
- fi
install:
- |
diff --git a/.ycm_extra_conf.py b/.ycm_extra_conf.py
index 8fb825ba3b1abe8630ae830f74807e7897bfda8e..a61e461ac159eaeefa9debeb1790b91c436434ac 100644
--- a/.ycm_extra_conf.py
+++ b/.ycm_extra_conf.py
@@ -36,19 +36,12 @@ import ycm_core
# CHANGE THIS LIST OF FLAGS. YES, THIS IS THE DROID YOU HAVE BEEN LOOKING FOR.
flags = [
'-DNDEBUG',
- '-DARB_HAVE_CTHREAD',
'-std=c++14',
'-x',
'c++',
'-I',
- 'src',
- '-I',
- '.',
- '-I',
'modcc',
'-I',
- 'include',
- '-I',
'arbor',
'-I',
'arbor/include',
@@ -57,17 +50,15 @@ flags = [
'-I',
'ext/json/single_include',
'-I',
- 'build/include',
- '-I',
- '/cm/shared/apps/cuda/8.0.44/include', # TODO: run a command to find this on "any" system
- '-I',
'python/pybind11/include',
'-I',
- '/usr/include/python3.6m', # TODO: run a command to find this on "any" system
- '-I',
'sup/include',
'-I',
- '/usr/include/python3.7m',
+ 'test',
+ '-I',
+ '/usr/include/python3.8',
+ '-I',
+ 'python',
]
# Set this to the absolute path to the folder (NOT the file!) containing the
@@ -153,13 +144,6 @@ def FlagsForFile( filename, **kwargs ):
compilation_info.compiler_flags_,
compilation_info.compiler_working_dir_ )
- # NOTE: This is just for YouCompleteMe; it's highly likely that your project
- # does NOT need to remove the stdlib flag. DO NOT USE THIS IN YOUR
- # ycm_extra_conf IF YOU'RE NOT 100% SURE YOU NEED IT.
- try:
- final_flags.remove( '-stdlib=libc++' )
- except ValueError:
- pass
else:
relative_to = DirectoryOfThisScript()
final_flags = MakeRelativePathsInFlagsAbsolute( flags, relative_to )
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 29a2c9066886554fc561f0a645f15459fbea99d9..fc86e4465e89791c936c4cead21aed11ceaed038 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -1,5 +1,4 @@
-# 3.9 requirement for CUDA language support and MPI support.
-cmake_minimum_required(VERSION 3.9)
+cmake_minimum_required(VERSION 3.12)
file(READ VERSION FULL_VERSION_STRING)
string(STRIP "${FULL_VERSION_STRING}" FULL_VERSION_STRING)
diff --git a/cmake/CompilerOptions.cmake b/cmake/CompilerOptions.cmake
index dac4325ea64f169877d72b27e2e09b5f4e02539e..8ff6415acc9bbc11b0d6fe6f1afc2e10ccf4a116 100644
--- a/cmake/CompilerOptions.cmake
+++ b/cmake/CompilerOptions.cmake
@@ -96,55 +96,6 @@ function(set_arch_target optvar arch)
else()
set(arch_opt "-mcpu=${arch}")
endif()
-
- elseif(CMAKE_CXX_COMPILER_ID MATCHES "Intel")
- # Translate target architecture names to Intel-compatible names.
- # icc 17 recognizes the following specific microarchitecture names for -mtune:
- # broadwell, haswell, ivybridge, knl, sandybridge, skylake
-
- if(arch MATCHES "sandybridge")
- set(tune "${arch}")
- set(arch "AVX")
- elseif(arch MATCHES "ivybridge")
- set(tune "${arch}")
- set(arch "CORE-AVX-I")
- elseif(arch MATCHES "broadwell|haswell|skylake")
- set(tune "${arch}")
- set(arch "CORE-AVX2")
- elseif(arch MATCHES "knl")
- set(tune "${arch}")
- set(arch "MIC-AVX512")
- elseif(arch MATCHES "nehalem|westmere")
- set(tune "corei7")
- set(arch "SSE4.2")
- elseif(arch MATCHES "core2")
- set(tune "core2")
- set(arch "SSSE3")
- elseif(arch MATCHES "native")
- unset(tune)
- set(arch "Host")
- else()
- set(tune "generic")
- set(arch "SSE2") # default for icc
- endif()
-
- if(tune)
- set(arch_opt "-x${arch};-mtune=${tune}")
- else()
- set(arch_opt "-x${arch}")
- endif()
-
- elseif(CMAKE_CXX_COMPILER_ID MATCHES "XL")
- # xlC 13 for Linux uses -mcpu. Not even attempting to get xlC 12 for BG/Q right
- # at this point: use CXXFLAGS as required!
- #
- # xlC, gcc, and clang all recognize power8 and power9 as architecture keywords.
-
- if(arch MATCHES "native")
- set(arch_opt "-qarch=auto")
- else()
- set(arch_opt "-mcpu=${arch}")
- endif()
endif()
get_property(enabled_languages GLOBAL PROPERTY ENABLED_LANGUAGES)
diff --git a/doc/index.rst b/doc/index.rst
index b87dc5b1aa19b83fb6e97ee72d0396125c88de82..6f4e7339f16aeb443cf7071f9fa3fe8003d05b05 100644
--- a/doc/index.rst
+++ b/doc/index.rst
@@ -72,6 +72,8 @@ Alternative citation formats for the paper can be `downloaded here <https://ieee
:caption: Getting Stared:
install
+ python
+ single_cell
.. toctree::
:caption: Arbor Models:
diff --git a/doc/python.rst b/doc/python.rst
new file mode 100644
index 0000000000000000000000000000000000000000..d3238b346766ac0f09635615a182cb9cfbbf9d8e
--- /dev/null
+++ b/doc/python.rst
@@ -0,0 +1,48 @@
+.. _getstarted_python:
+
+Python
+======
+
+Arbor provides access to all of the C++ library's functionality in Python,
+which is the only interface for many users.
+The getting started guides will introduce Arbor via the Python interface.
+
+Before starting Arbor needs to be installed with the Python interface enabled,
+as per the `installation guide <_installarbor>`_.
+To test that Arbor is available, open `Python 3 <python2_>`_, and try the following:
+
+.. container:: example-code
+
+ .. code-block:: python
+
+ import arbor
+ arbor.config()
+
+ {'mpi': True, 'mpi4py': True, 'gpu': False, 'version': '0.3'}
+
+Calling ``arbor.config()`` returns a dictionary with information about the Arbor installation.
+This can be used to check that Arbor supports features that you require to run your model,
+or even to dynamically decide how to run a model.
+Single cell models like the one introduced here we do not require parallelism like
+that provided by MPI or GPUs, so the ``'mpi'`` and ``'gpu'`` fields can be ``False``.
+
+Performance
+--------------
+
+The Python interface can incur significant performance overheads relative to C++
+during the *model building* phase, however simulation performance will be the same
+for both interfaces.
+
+.. _python2:
+
+Python 2
+----------
+
+Python 2.7 will reach `end of life <https://pythonclock.org/>`_ in January 2020.
+Arbor should be installed using Python 3, and all examples in the documentation are in
+Python 3. It might be possible to install and run Arbor using Python 2.7, however Arbor support
+will only be provided for Python 3.
+
+NEURON users that use Python 2 can take the opportunity of trying Arbor to make
+the move to Python 3.
+
diff --git a/doc/single_cell.rst b/doc/single_cell.rst
new file mode 100644
index 0000000000000000000000000000000000000000..94de6c7306a6eb8d8d33335ba1a3ad263861fd64
--- /dev/null
+++ b/doc/single_cell.rst
@@ -0,0 +1,90 @@
+.. _single_cell:
+
+Single Cell Models
+==================
+
+Arbor breaks down building *single cell models* into the following steps:
+
+1. Defining the `morphology <single_morpho_>`_ of the cell.
+2. Labeling regions and locations on the morphology.
+3. Defining the mechanisms that will be applied to the cell.
+4. Adding ion channels and synapses (mechanisms) to labeled regions and locations.
+5. Attaching stimuli, spike detectors, event generators and probes to locations (inputs & outputs).
+
+This guide will provide a simple example workflow of this kind that the reader
+can follow along in Python, with links to separate documentation that covers
+relevant topics in more detail.
+
+.. note::
+ Most readers will be familiar with NEURON. Boxes like this
+ will be used to highlight differences between NEURON and Arbor
+ throughout the guide.
+
+ NEURON users will recognise that Arbor uses many similar concepts, and
+ an effort has been made to use the same nomenclature wherever possible.
+
+ Arbor takes a more structured approach to model building,
+ from morphology descriptions up to network connectivity, to allow model
+ descriptions that are more scalable and portable.
+
+.. _single_morpho:
+
+Morphology
+----------
+
+The first step in building a cell model is to define the cell's *morphology*.
+Conceptually, Arbor treats morphologies as a tree of truncated frustums, with
+an optional spherical segment at the root of the tree.
+These are represented as a tree of sample points, where each sample has a 3D location,
+a radius, and a tag, and a parent sample.
+
+.. note::
+ NEURON represents morphologies as a tree of cylindrical *segments*, whereas
+ in Arbor the radius can vary linearly between two sample locations.
+
+ A cylinder with equal diameter and length is used to model spherical somata
+ in NEURON, which coincidently has the same surface area as a sphere of the same diameter.
+ Arbor allows the user to optionally use a spherical section at the root
+ of the tree to represent spherical somata.
+
+.. note::
+ In NEURON cell morphologies are constructed by creating individual sections,
+ then connecting them together. In Arbor we start with an "empty"
+ sample tree, to which samples are appended to build a connected morphology.
+
+Let's start with a simple "ball and stick" model cell.
+
+.. container:: example-code
+
+ .. code-block:: python
+
+ import arbor
+ builder = arbor.flat_cell_builder()
+
+ # Start with a spherical segment with radius 10 μm.
+ # Label this segment 'soma'.
+ p = builder.add_sphere(radius=10, name='soma')
+
+ # Attach a cable to the soma with length 100 μm and constant radius 4 μm.
+ q = builder.add_cable(parent=p, length=100, radius=4, name='dend')
+
+ # Attach two dendrites to the first of length 50 μm, that taper from 4 μm to 2 μm.
+ p = builder.add_cable(parent=q, length=50, radius=(4,2), name='dend')
+ p = builder.add_cable(parent=q, length=50, radius=(4,2), name='dend')
+
+
+Building the morphology there are two approaches: construct it manually using
+``sample_tree`` or ``flat_cell_builder``, or load from swc file.
+
+TODO: cover all methods here?
+ - we could just ``flat_cell_builder`` because it is most comfortable for
+ users coming over from NEURON.
+ - have links to another page that goes into detail on all the different
+ methods for morphology building. That page could take a moderately
+ complicated, well-defined, morphology, and illustrate how to build
+ it using all of the different methods.
+
+Labeling Regions and Locations
+------------------------------
+
+
diff --git a/example/ring/ring.cpp b/example/ring/ring.cpp
index 611ca77c9d5f54034a4ee697f017263896eaef5c..91fc043baf48c1d46af4bf14c6a702ade317d7f5 100644
--- a/example/ring/ring.cpp
+++ b/example/ring/ring.cpp
@@ -226,7 +226,7 @@ int main(int argc, char** argv) {
// The id of the only probe on the cell: the cell_member type points to (cell 0, probe 0)
auto probe_id = cell_member_type{0, 0};
// The schedule for sampling is 10 samples every 1 ms.
- auto sched = arb::regular_schedule(0.1);
+ auto sched = arb::regular_schedule(1);
// This is where the voltage samples will be stored as (time, value) pairs
arb::trace_data<double> voltage;
// Now attach the sampler at probe_id, with sampling schedule sched, writing to voltage
diff --git a/python/CMakeLists.txt b/python/CMakeLists.txt
index ad34fce7044a8855f9b11a882fdac6fb52623b0b..bf4d5f006fcda97c60658f4cccac595b5eb0180f 100644
--- a/python/CMakeLists.txt
+++ b/python/CMakeLists.txt
@@ -14,16 +14,17 @@ endif()
set(PYBIND11_CPP_STANDARD -std=c++14)
add_subdirectory(pybind11)
-# The Python library. MODULE will make a Python-exclusive model.
-add_library(pyarb MODULE
+set(pyarb_source
cells.cpp
config.cpp
context.cpp
domain_decomposition.cpp
error.cpp
event_generator.cpp
+ flat_cell_builder.cpp
identifiers.cpp
morphology.cpp
+ morph_parse.cpp
mpi.cpp
profiler.cpp
pyarb.cpp
@@ -31,10 +32,19 @@ add_library(pyarb MODULE
sampling.cpp
schedule.cpp
simulation.cpp
+ single_cell_model.cpp
spikes.cpp
+ s_expr.cpp
)
-target_link_libraries(pyarb PRIVATE arbor pybind11::module)
+# compile the pyarb sources into an object library that will be
+# use by both the Python wrapper target (pyarb) and for the
+# unit tests of the C++ components in the Python wrapper.
+add_library(pyarb_obj OBJECT ${pyarb_source})
+target_link_libraries(pyarb_obj PRIVATE arbor pybind11::module)
+
+# The Python library. MODULE will make a Python-exclusive model.
+add_library(pyarb MODULE $<TARGET_OBJECTS:pyarb_obj>)
# The output name of the pyarb .so file is "arbor", to facilitate "import arbor"
set_target_properties(pyarb PROPERTIES OUTPUT_NAME arbor)
@@ -42,6 +52,10 @@ set_target_properties(pyarb PROPERTIES OUTPUT_NAME arbor)
# arbor.cpython-37m-x86_64-linux-gnu.so
set_target_properties(pyarb PROPERTIES PREFIX "${PYTHON_MODULE_PREFIX}" SUFFIX "${PYTHON_MODULE_EXTENSION}")
+# this dependency has to be spelt out again, despite being added to
+# pyarb_obj because CMake.
+target_link_libraries(pyarb PRIVATE arbor pybind11::module)
+
# Add support for mpi4py if available.
if (ARB_WITH_MPI)
find_python_module(mpi4py)
@@ -55,3 +69,6 @@ endif()
find_package(PythonInterp REQUIRED)
set(ARB_PYEXECDIR "${CMAKE_INSTALL_LIBDIR}/python${PYTHON_VERSION_MAJOR}.${PYTHON_VERSION_MINOR}/site-packages")
install(TARGETS pyarb LIBRARY DESTINATION ${ARB_PYEXECDIR})
+
+# for unit tests on C++ side of Python wrappers
+add_subdirectory(test)
diff --git a/python/cells.cpp b/python/cells.cpp
index 6eaa3cb09f5caf1133c56913b0ff75998360c1f4..f9807142e457c19e57b862cf82fdd381514f4f8a 100644
--- a/python/cells.cpp
+++ b/python/cells.cpp
@@ -4,11 +4,17 @@
#include <arbor/benchmark_cell.hpp>
#include <arbor/cable_cell.hpp>
#include <arbor/lif_cell.hpp>
+#include <arbor/morph/label_dict.hpp>
+#include <arbor/morph/locset.hpp>
+#include <arbor/morph/region.hpp>
#include <arbor/schedule.hpp>
#include <arbor/spike_source_cell.hpp>
+#include <arbor/util/any.hpp>
#include <arbor/util/unique_any.hpp>
+#include "conversion.hpp"
#include "error.hpp"
+#include "morph_parse.hpp"
#include "schedule.hpp"
#include "strprintf.hpp"
@@ -44,119 +50,107 @@ arb::util::unique_any convert_cell(pybind11::object o) {
}
//
-// Somewhat hacky bit of code for generating cells with random morphologies.
+// proxies
//
-// Parameters used to generate the random cell morphologies.
-struct cell_parameters {
- cell_parameters() = default;
-
- // Maximum number of levels in the cell (not including the soma)
- unsigned max_depth = 5;
-
- // The following parameters are described as ranges.
- // The first value is at the soma, and the last value is used on the last level.
- // Values at levels in between are found by linear interpolation.
- std::array<double,2> branch_probs = {1.0, 0.5}; // Probability of a branch occuring.
- std::array<unsigned,2> compartments = {20, 2}; // Compartment count on a branch.
- std::array<double,2> lengths = {200, 20}; // Length of branch in μm.
-
- // The number of synapses per cell.
- unsigned synapses = 1;
-
- friend std::ostream& operator<<(std::ostream& o, const cell_parameters& p) {
- return
- o << "<cell_parameters: depth " << p.max_depth
- << ", synapses " << p.synapses
- << ", branch_probs [" << p.branch_probs[0] << ":" << p.branch_probs[1] << "]"
- << ", compartments [" << p.compartments[0] << ":" << p.compartments[1] << "]"
- << ", lengths [" << p.lengths[0] << ":" << p.lengths[1] << "]>";
- }
-};
+struct label_dict_proxy {
+ using str_map = std::unordered_map<std::string, std::string>;
+ arb::label_dict dict;
+ str_map cache;
+ std::vector<std::string> locsets;
+ std::vector<std::string> regions;
-// Helper used to interpolate in branch_cell.
-template <typename T>
-double interp(const std::array<T,2>& r, unsigned i, unsigned n) {
- double p = i * 1./(n-1);
- double r0 = r[0];
- double r1 = r[1];
- return r[0] + p*(r1-r0);
-}
+ label_dict_proxy() = default;
-arb::cable_cell make_cable_cell(arb::cell_gid_type gid, const cell_parameters& params) {
- arb::sample_tree tree;
-
- // Add soma.
- double soma_radius = 12.6157/2.0;
- tree.append(arb::mnpos, {{0,0,0,soma_radius}, 1}); // For area of 500 μm².
-
- std::vector<std::vector<unsigned>> levels;
- levels.push_back({0});
-
- // Standard mersenne_twister_engine seeded with gid.
- std::mt19937 gen(gid);
- std::uniform_real_distribution<double> dis(0, 1);
-
- double dend_radius = 0.5; // Diameter of 1 μm for each cable.
-
- double dist_from_soma = soma_radius;
- for (unsigned i=0; i<params.max_depth; ++i) {
- // Branch prob at this level.
- double bp = interp(params.branch_probs, i, params.max_depth);
- // Length at this level.
- double l = interp(params.lengths, i, params.max_depth);
- // Number of compartments at this level.
- unsigned nc = std::round(interp(params.compartments, i, params.max_depth));
-
- std::vector<unsigned> sec_ids;
- for (unsigned sec: levels[i]) {
- for (unsigned j=0; j<2; ++j) {
- if (dis(gen)<bp) {
- 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<nc; ++k) {
- p = tree.append(p, {{0,0,z+k*dz, dend_radius}, 3});
- }
- }
- sec_ids.push_back(tree.append(p, {{0,0,z+l,dend_radius}, 3}));
- }
- }
+ label_dict_proxy(const str_map& in) {
+ for (auto& i: in) {
+ set(i.first.c_str(), i.second.c_str());
}
- if (sec_ids.empty()) {
- break;
- }
- levels.push_back(sec_ids);
-
- dist_from_soma += l;
}
- arb::label_dict d;
-
- using arb::reg::tagged;
- d.set("soma", tagged(1));
- d.set("dendrites", join(tagged(3), tagged(4)));
-
- arb::cable_cell cell(arb::morphology(tree, true), d);
-
- cell.paint("soma", "hh");
- cell.paint("dendrites", "pas");
- cell.default_parameters.axial_resistivity = 100; // [Ω·cm]
-
- // Add spike threshold detector at the soma.
- cell.place(arb::mlocation{0,0}, arb::threshold_detector{10});
+ std::size_t size() const {
+ return locsets.size() + regions.size();
+ }
- // Add a synapse to the mid point of the first dendrite.
- cell.place(arb::mlocation{1, 0.5}, "expsyn");
+ void set(const char* name, const char* desc) {
+ using namespace std::string_literals;
+ // The following code takes an input name and a region or locset
+ // description, e.g.:
+ // name='reg', desc='(tag 4)'
+ // name='loc', desc='(terminal)'
+ // name='foo', desc='(join (tag 2) (tag 3))'
+ // Then it parses the description, and tests whether the description
+ // is a region or locset, and updates the label dictionary appropriately.
+ // Errors occur when:
+ // * the name is not a valid name.
+ // * a region is described with a name that matches an existing locset
+ // (and vice versa.)
+ // * the description is not well formed, e.g. it contains a syntax error.
+ // * the description is well-formed, but describes neither a region or locset.
+ try{
+ // Test that the identifier is valid, i.e.
+ // * only numbers, letters and underscore.
+ // * no leading number or underscore.
+ if (!test_identifier(name)) {
+ throw std::string(util::pprintf("'{}' is not a valid label name.", name));
+ }
+ // Parse the input string into an s-expression.
+ auto parsed = parse(desc);
+ // Evaluate the s-expression to build a region/locset.
+ auto result = eval(parsed);
+ if (!result) { // an error parsing / evaluating description.
+ throw std::string(result.error().message);
+ }
+ else if (result->type()==typeid(arb::region)) { // describes a region.
+ dict.set(name, std::move(arb::util::any_cast<arb::region&>(*result)));
+ auto it = std::lower_bound(regions.begin(), regions.end(), name);
+ if (it==regions.end() || *it!=name) regions.insert(it, name);
+ }
+ else if (result->type()==typeid(arb::locset)) { // describes a locset.
+ dict.set(name, std::move(arb::util::any_cast<arb::locset&>(*result)));
+ auto it = std::lower_bound(locsets.begin(), locsets.end(), name);
+ if (it==locsets.end() || *it!=name) locsets.insert(it, name);
+ }
+ else {
+ // Successfully parsed an expression that is neither region nor locset.
+ throw util::pprintf("The defninition of '{} = {}' does not define a valid region or locset.", name, desc);
+ }
+ // The entry was added succesfully: store it in the cache.
+ cache[name] = desc;
+ }
+ catch (std::string msg) {
+ const char* base = "\nError adding the label '{}' = '{}'\n{}\n";
- // Add additional synapses that will not be connected to anything.
- for (unsigned i=1u; i<params.synapses; ++i) {
- cell.place(arb::mlocation{1, 0.5}, "expsyn");
+ throw std::runtime_error(util::pprintf(base, name, desc, msg));
+ }
+ // Exceptions are thrown in parse or eval if an unexpected error occured.
+ catch (std::exception& e) {
+ const char* msg =
+ "\n----- internal error -------------------------------------------"
+ "\nError parsing the label: '{}' = '{}'"
+ "\n"
+ "\n{}"
+ "\n"
+ "\nPlease file a bug report with this full error message at:"
+ "\n github.com/arbor-sim/arbor/issues"
+ "\n----------------------------------------------------------------";
+ throw arb::arbor_internal_error(util::pprintf(msg, name, desc, e.what()));
+ }
}
- return cell;
-}
+ std::string to_string() const {
+ std::string s;
+ s += "(label_dict";
+ for (auto& x: dict.regions()) {
+ s += util::pprintf(" (region '{}' {})", x.first, x.second);
+ }
+ for (auto& x: dict.locsets()) {
+ s += util::pprintf(" (locset '{}' {})", x.first, x.second);
+ }
+ s += ")";
+ return s;
+ }
+};
//
// string printers
@@ -164,12 +158,21 @@ arb::cable_cell make_cable_cell(arb::cell_gid_type gid, const cell_parameters& p
std::string lif_str(const arb::lif_cell& c){
return util::pprintf(
- "<arbor.lif_cell: tau_m {}, V_th {}, C_m {}, E_L {}, V_m {}, t_ref {}, V_reset {}>",
- c.tau_m, c.V_th, c.C_m, c.E_L, c.V_m, c.t_ref, c.V_reset);
+ "<arbor.lif_cell: tau_m {}, V_th {}, C_m {}, E_L {}, V_m {}, t_ref {}, V_reset {}>",
+ c.tau_m, c.V_th, c.C_m, c.E_L, c.V_m, c.t_ref, c.V_reset);
+}
+
+
+std::string mechanism_desc_str(const arb::mechanism_desc& md) {
+ return util::pprintf("<arbor.mechanism: name '{}', parameters {}",
+ md.name(), util::dictionary_csv(md.values()));
}
void register_cells(pybind11::module& m) {
using namespace pybind11::literals;
+ using arb::util::optional;
+
+ // arb::spike_source_cell
pybind11::class_<arb::spike_source_cell> spike_source_cell(m, "spike_source_cell",
"A spike source cell, that generates a user-defined sequence of spikes that act as inputs for other cells in the network.");
@@ -190,6 +193,8 @@ void register_cells(pybind11::module& m) {
.def("__repr__", [](const arb::spike_source_cell&){return "<arbor.spike_source_cell>";})
.def("__str__", [](const arb::spike_source_cell&){return "<arbor.spike_source_cell>";});
+ // arb::benchmark_cell
+
pybind11::class_<arb::benchmark_cell> benchmark_cell(m, "benchmark_cell",
"A benchmarking cell, used by Arbor developers to test communication performance.\n"
"A benchmark cell generates spikes at a user-defined sequence of time points, and\n"
@@ -216,47 +221,325 @@ void register_cells(pybind11::module& m) {
.def("__repr__", [](const arb::benchmark_cell&){return "<arbor.benchmark_cell>";})
.def("__str__", [](const arb::benchmark_cell&){return "<arbor.benchmark_cell>";});
+ // arb::lif_cell
+
pybind11::class_<arb::lif_cell> lif_cell(m, "lif_cell",
"A benchmarking cell, used by Arbor developers to test communication performance.");
lif_cell
.def(pybind11::init<>())
- .def_readwrite("tau_m", &arb::lif_cell::tau_m, "Membrane potential decaying constant [ms].")
- .def_readwrite("V_th", &arb::lif_cell::V_th, "Firing threshold [mV].")
- .def_readwrite("C_m", &arb::lif_cell::C_m, "Membrane capacitance [pF].")
- .def_readwrite("E_L", &arb::lif_cell::E_L, "Resting potential [mV].")
- .def_readwrite("V_m", &arb::lif_cell::V_m, "Initial value of the Membrane potential [mV].")
- .def_readwrite("t_ref", &arb::lif_cell::t_ref, "Refractory period [ms].")
- .def_readwrite("V_reset", &arb::lif_cell::V_reset, "Reset potential [mV].")
+ .def_readwrite("tau_m", &arb::lif_cell::tau_m,
+ "Membrane potential decaying constant [ms].")
+ .def_readwrite("V_th", &arb::lif_cell::V_th,
+ "Firing threshold [mV].")
+ .def_readwrite("C_m", &arb::lif_cell::C_m,
+ " Membrane capacitance [pF].")
+ .def_readwrite("E_L", &arb::lif_cell::E_L,
+ "Resting potential [mV].")
+ .def_readwrite("V_m", &arb::lif_cell::V_m,
+ "Initial value of the Membrane potential [mV].")
+ .def_readwrite("t_ref", &arb::lif_cell::t_ref,
+ "Refractory period [ms].")
+ .def_readwrite("V_reset", &arb::lif_cell::V_reset,
+ "Reset potential [mV].")
.def("__repr__", &lif_str)
.def("__str__", &lif_str);
- pybind11::class_<cell_parameters> cell_params(m, "cell_parameters", "Parameters used to generate the random cell morphologies.");
- cell_params
+ // arb::label_dict
+
+ pybind11::class_<label_dict_proxy> label_dict(m, "label_dict",
+ "A dictionary of labelled region and locset definitions, with a\n"
+ "unique label is assigned to each definition.");
+ label_dict
+ .def(pybind11::init<>(), "Create an empty label dictionary.")
+ .def(pybind11::init<const std::unordered_map<std::string, std::string>&>(),
+ "Initialize a label dictionary from a dictionary with string labels as keys,"
+ " and corresponding definitions as strings.")
+ .def("__setitem__",
+ [](label_dict_proxy& l, const char* name, const char* desc) {
+ l.set(name, desc);})
+ .def("__getitem__",
+ [](label_dict_proxy& l, const char* name) {
+ if (!l.cache.count(name)) {
+ throw std::runtime_error(util::pprintf("\nKeyError: '{}'", name));
+ }
+ return l.cache.at(name);
+ })
+ .def("__len__", &label_dict_proxy::size)
+ .def("__iter__",
+ [](const label_dict_proxy &ld) {
+ return pybind11::make_key_iterator(ld.cache.begin(), ld.cache.end());},
+ pybind11::keep_alive<0, 1>())
+ .def_readonly("regions", &label_dict_proxy::regions,
+ "The region definitions.")
+ .def_readonly("locsets", &label_dict_proxy::locsets,
+ "The locset definitions.")
+ .def("__repr__", [](const label_dict_proxy& d){return d.to_string();})
+ .def("__str__", [](const label_dict_proxy& d){return d.to_string();});
+
+ // Data structures used to describe mechanisms, electrical properties,
+ // gap junction properties, etc.
+
+ // arb::cable_cell_ion_data
+
+ pybind11::class_<arb::initial_ion_data> ion_data(m, "ion",
+ "For setting ion properties (internal and external concentration and reversal potential) on cells and regions.");
+ ion_data
+ .def(pybind11::init(
+ [](const char* name,
+ optional<double> int_con, optional<double> ext_con,
+ optional<double> rev_pot)
+ {
+ arb::initial_ion_data x;
+ x.ion = name;
+ if (int_con) x.initial.init_int_concentration = *int_con;
+ if (ext_con) x.initial.init_int_concentration = *ext_con;
+ if (rev_pot) x.initial.init_reversal_potential = *rev_pot;
+ return x;
+ }
+ ),
+ "ion_name"_a,
+ pybind11::arg_v("int_con", pybind11::none(), "Intial internal concentration [mM]"),
+ pybind11::arg_v("ext_con", pybind11::none(), "Intial external concentration [mM]"),
+ pybind11::arg_v("rev_pot", pybind11::none(), "Intial reversal potential [mV]"),
+ "If concentrations or reversal potential are specified as 'None',"
+ " cell default or global default value will be used, in that order if set.");
+
+ // arb::mechanism_desc
+
+ pybind11::class_<arb::mechanism_desc> mechanism_desc(m, "mechanism");
+ mechanism_desc
+ .def(pybind11::init([](const char* n) {return arb::mechanism_desc{n};}))
+ // allow construction of a description with parameters provided in a dictionary:
+ // mech = arbor.mechanism('mech_name', {'param1': 1.2, 'param2': 3.14})
+ .def(pybind11::init(
+ [](const char* name, std::unordered_map<std::string, double> params) {
+ arb::mechanism_desc md(name);
+ for (const auto& p: params) md.set(p.first, p.second);
+ return md;
+ }),
+ "name"_a, "The name of the mechanism",
+ "params"_a, "A dictionary of parameter values, with parameter name as key.",
+ "Example usage setting pararmeters:\n"
+ " m = arbor.mechanism('expsyn', {'tau': 1.4})\n"
+ "will create parameters for the 'expsyn' mechanism, with the provided value\n"
+ "for 'tau' overrides the default. If a parameter is not set, the default\n"
+ "(as defined in NMODL) is used.\n\n"
+ "Example overriding a global parameter:\n"
+ " m = arbor.mechanism('nernst/R=8.3145,F=96485')")
+ .def("set",
+ [](arb::mechanism_desc& md, std::string name, double value) {
+ md.set(name, value);
+ },
+ "name"_a, "value"_a, "Set parameter value.")
+ .def_property_readonly("name",
+ [](const arb::mechanism_desc& md) {
+ return md.name();
+ },
+ "The name of the mechanism.")
+ .def_property_readonly("values",
+ [](const arb::mechanism_desc& md) {
+ return md.values();
+ }, "A dictionary of parameter values with parameter name as key.")
+ .def("__repr__",
+ [](const arb::mechanism_desc& md) {
+ return util::pprintf("<arbor.mechanism: name '{}', parameters {}", md.name(), util::dictionary_csv(md.values())); })
+ .def("__str__",
+ [](const arb::mechanism_desc& md) {
+ return util::pprintf("('{}' {})", md.name(), util::dictionary_csv(md.values())); });
+
+ // arb::gap_junction_site
+
+ pybind11::class_<arb::gap_junction_site> gjsite(m, "gap_junction",
+ "For marking a location on a cell morphology as a gap junction site.");
+ gjsite
.def(pybind11::init<>())
- .def_readwrite("depth", &cell_parameters::max_depth, "The maximum depth of the branch structure.")
- .def_readwrite("lengths", &cell_parameters::lengths, "Length of branch [μm], given as range.")
- .def_readwrite("synapses", &cell_parameters::synapses, "The number of randomly generated synapses on the cell.")
- .def_readwrite("branch_probs", &cell_parameters::branch_probs, "Probability of a branch occuring, given as range.")
- .def_readwrite("compartments", &cell_parameters::compartments, "Compartment count on a branch, given as range.")
- .def("__repr__", util::to_string<cell_parameters>)
- .def("__str__", util::to_string<cell_parameters>);
-
- // Wrap cable cell description type.
- // Cable cells are going to be replaced with a saner API, so we don't go
- // adding much in the way of interface. Instead we just provide a helper
- // that will generate random cell morphologies for benchmarking.
- pybind11::class_<arb::cable_cell> cable_cell(m, "cable_cell");
-
+ .def("__repr__", [](const arb::gap_junction_site&){return "<arbor.gap_junction>";})
+ .def("__str__", [](const arb::gap_junction_site&){return "<arbor.gap_junction>";});
+
+ // arb::i_clamp
+
+ pybind11::class_<arb::i_clamp> i_clamp(m, "iclamp",
+ "A current clamp, for injecting a single pulse of current with fixed duration and current.");
+ i_clamp
+ .def(pybind11::init(
+ [](double ts, double dur, double cur) {
+ return arb::i_clamp{ts, dur, cur};
+ }), "tstart"_a=0, "duration"_a=0, "current"_a=0)
+ .def_readonly("tstart", &arb::i_clamp::delay, "Time at which current starts [ms]")
+ .def_readonly("duration", &arb::i_clamp::duration, "Duration of the current injection [ms]")
+ .def_readonly("current", &arb::i_clamp::amplitude, "Amplitude of the injected current [nA]")
+ .def("__repr__", [](const arb::i_clamp& c){
+ return util::pprintf("(iclamp {} {} {})", c.delay, c.duration, c.amplitude);})
+ .def("__str__", [](const arb::i_clamp& c){
+ return util::pprintf("<arbor.iclamp: tstart {} ms, duration {} ms, current {} nA>",
+ c.delay, c.duration, c.amplitude);});
+
+ // arb::threshold_detector
+
+ pybind11::class_<arb::threshold_detector> detector(m, "spike_detector",
+ "A spike detector, generates a spike when voltage crosses a threshold.");
+ detector
+ .def(pybind11::init(
+ [](double thresh) {
+ return arb::threshold_detector{thresh};
+ }), "threshold"_a)
+ .def_readonly("threshold", &arb::threshold_detector::threshold, "Voltage threshold of spike detector [ms]")
+ .def("__repr__", [](const arb::threshold_detector& d){
+ return util::pprintf("<arbor.threshold_detector: threshold {} mV>", d.threshold);})
+ .def("__str__", [](const arb::threshold_detector& d){
+ return util::pprintf("(threshold_detector {})", d.threshold);});
+
+ // arb::cable_cell
+
+ pybind11::class_<arb::cable_cell> cable_cell(m, "cable_cell",
+ "Represents morphologically-detailed cell models, with morphology represented as a\n"
+ "tree of one-dimensional cable segments.");
cable_cell
+ .def(pybind11::init(
+ [](const arb::morphology& m, const label_dict_proxy& labels) {
+ return arb::cable_cell(m, labels.dict);
+ }), "morphology"_a, "labels"_a)
+ .def(pybind11::init(
+ [](const arb::sample_tree& t, const label_dict_proxy& labels) {
+ return arb::cable_cell(arb::morphology(t), labels.dict);
+ }),
+ "morphology"_a, "labels"_a,
+ "Construct with a morphology derived from a sample_tree, with automatic detection of whether\n"
+ "the morphology has a spherical root/soma.")
+ .def_property_readonly("num_branches",
+ [](const arb::cable_cell& c) {return c.morphology().num_branches();},
+ "The number of unbranched cable sections in the morphology.")
+ // Set cell-wide properties
+ .def("set_properties",
+ [](arb::cable_cell& c,
+ optional<double> Vm, optional<double> cm,
+ optional<double> rL, optional<double> tempK)
+ {
+ if (Vm) c.default_parameters.init_membrane_potential = Vm;
+ if (cm) c.default_parameters.membrane_capacitance=cm;
+ if (rL) c.default_parameters.axial_resistivity=rL;
+ if (tempK) c.default_parameters.temperature_K=tempK;
+ },
+ "Vm"_a=pybind11::none(), "cm"_a=pybind11::none(), "rL"_a=pybind11::none(), "tempK"_a=pybind11::none(),
+ "Set default values for cable and cell properties. These values can be overridden on specific regions using the paint interface.\n"
+ " Vm: initial membrane voltage [mV].\n"
+ " cm: membrane capacitance [F/m²].\n"
+ " rL: axial resistivity [Ω·cm].\n"
+ " tempK: temperature [Kelvin].")
+ .def("set_ion",
+ [](arb::cable_cell& c, const char* ion,
+ optional<double> int_con, optional<double> ext_con,
+ optional<double> rev_pot, optional<arb::mechanism_desc> method)
+ {
+ auto& data = c.default_parameters.ion_data[ion];
+ if (int_con) data.init_int_concentration = *int_con;
+ if (ext_con) data.init_ext_concentration = *ext_con;
+ if (rev_pot) data.init_reversal_potential = *rev_pot;
+ if (method) c.default_parameters.reversal_potential_method[ion] = *method;
+ },
+ "ion"_a,
+ "int_con"_a=pybind11::none(),
+ "ext_con"_a=pybind11::none(),
+ "rev_pot"_a=pybind11::none(),
+ "method"_a=pybind11::none(),
+ "Set the propoerties of ion species named 'ion' that will be applied\n"
+ "by default everywhere on the cell. Species concentrations and reversal\n"
+ "potential can be overridden on specific regions using the paint interface, \n"
+ "while the method for calculating reversal potential is global for all\n"
+ "compartments in the cell, and can't be overriden locally.\n"
+ " ion: name of ion species.\n"
+ " int_con: initial internal concentration [mM].\n"
+ " ext_con: initial external concentration [mM].\n"
+ " rev_pot: reversal potential [mV].\n"
+ " method: method for calculating reversal potential.")
+ // Paint mechanisms.
+ .def("paint",
+ [](arb::cable_cell& c, const char* region, const arb::mechanism_desc& d) {
+ c.paint(region, d);
+ },
+ "region"_a, "mechanism"_a,
+ "Associate a mechanism with a region.")
+ .def("paint",
+ [](arb::cable_cell& c, const char* region, const char* mech_name) {
+ c.paint(region, mech_name);
+ },
+ "region"_a, "mechanism"_a,
+ "Associate a mechanism with a region.")
+ // Paint membrane/static properties.
+ .def("paint",
+ [](arb::cable_cell& c,
+ const char* region,
+ optional<double> Vm, optional<double> cm,
+ optional<double> rL, optional<double> tempK)
+ {
+ if (Vm) c.paint(region, arb::init_membrane_potential{*Vm});
+ if (cm) c.paint(region, arb::membrane_capacitance{*cm});
+ if (rL) c.paint(region, arb::axial_resistivity{*rL});
+ if (tempK) c.paint(region, arb::temperature_K{*tempK});
+ },
+ "region"_a, "Vm"_a=pybind11::none(), "cm"_a=pybind11::none(), "rL"_a=pybind11::none(), "tempK"_a=pybind11::none(),
+ "Set cable properties on a region.\n"
+ " region: initial membrane voltage [mV].\n"
+ " cm: membrane capacitance [F/m²].\n"
+ " rL: axial resistivity [Ω·cm].\n"
+ " tempK: temperature [Kelvin].")
+
+ // Paint ion species initial conditions on a region.
+ .def("paint",
+ [](arb::cable_cell& c, const char* region, const arb::initial_ion_data& d) {
+ c.paint(region, d);
+ },
+ "region"_a, "ion_data"_a,
+ "Set ion species properties conditions on a region.")
+ // Place synapses
+ .def("place",
+ [](arb::cable_cell& c, const char* locset, const arb::mechanism_desc& d) {
+ c.place(locset, d); },
+ "locations"_a, "mechanism"_a,
+ "Place one instance of synapse described by 'mechanism' to each location in 'locations'.")
+ .def("place",
+ [](arb::cable_cell& c, const char* locset, const char* mech_name) {
+ c.place(locset, mech_name);
+ },
+ "locations"_a, "mechanism"_a,
+ "Place one instance of synapse described by 'mechanism' to each location in 'locations'.")
+ // Place gap junctions.
+ .def("place",
+ [](arb::cable_cell& c, const char* locset, const arb::gap_junction_site& site) {
+ c.place(locset, site);
+ },
+ "locations"_a, "gapjunction"_a,
+ "Place one gap junction site at each location in 'locations'.")
+ // Place current clamp stimulus.
+ .def("place",
+ [](arb::cable_cell& c, const char* locset, const arb::i_clamp& stim) {
+ c.place(locset, stim);
+ },
+ "locations"_a, "iclamp"_a,
+ "Add a current stimulus at each location in locations.")
+ // Place spike detector.
+ .def("place",
+ [](arb::cable_cell& c, const char* locset, const arb::threshold_detector& d) {
+ c.place(locset, d);
+ },
+ "locations"_a, "detector"_a,
+ "Add a voltage spike detector at each location in locations.")
+ // Discretization control.
+ .def("compartments_on_samples",
+ [](arb::cable_cell& c) {c.default_parameters.discretization = arb::cv_policy_every_sample{};},
+ "Decompose each branch into compartments defined by sample locations.")
+ .def("compartments_length",
+ [](arb::cable_cell& c, double len) {
+ c.default_parameters.discretization = arb::cv_policy_max_extent{len};
+ },
+ "maxlen"_a, "Decompose each branch into compartments of equal length, not exceeding maxlen.")
+ .def("compartments_per_branch",
+ [](arb::cable_cell& c, unsigned n) {c.default_parameters.discretization = arb::cv_policy_fixed_per_branch{n};},
+ "n"_a, "Decompose each branch into n compartments of equal length.")
+ // Stringification
.def("__repr__", [](const arb::cable_cell&){return "<arbor.cable_cell>";})
.def("__str__", [](const arb::cable_cell&){return "<arbor.cable_cell>";});
-
- m.def("make_cable_cell", &make_cable_cell,
- "Construct a branching cell with a random morphology and synapse end points locations described by params.\n"
- "seed is an integral value used to seed the random number generator, for which the gid of the cell is a good default.",
- "seed"_a,
- "params"_a=cell_parameters());
}
-}
+} // namespace pyarb
diff --git a/python/conversion.hpp b/python/conversion.hpp
index 6db53d8ec4287be78cdfc9112d108e5809dfb649..5d58748adcc479fc4e5529390a3593e474e7448d 100644
--- a/python/conversion.hpp
+++ b/python/conversion.hpp
@@ -17,10 +17,17 @@ namespace pybind11 { namespace detail {
namespace pyarb {
struct is_nonneg {
+ template<typename T>
+ constexpr bool operator()(const T& v) {
+ return v>=T(0);
+ }
+};
+
+struct is_positive {
template<typename T>
constexpr
bool operator()(const T& v) {
- return v>=T(0);
+ return v>T(0);
}
};
@@ -65,4 +72,19 @@ arb::util::optional<T> py2optional(pybind11::object o, const char* msg) {
return o.is_none()? arb::util::nullopt: arb::util::optional<T>(std::move(value));
}
+// Attempt to cast a Python object to a C++ type T.
+// Returns an optional that is set if o could be cast
+// to T, otherwise it is empty. Hence not being able
+// to cast is not an error.
+template <typename T>
+arb::util::optional<T> try_cast(pybind11::object o) {
+ try {
+ return o.cast<T>();
+ }
+ // Ignore cast_error: if unable to perform cast.
+ catch (pybind11::cast_error& e) {}
+
+ return arb::util::nullopt;
}
+
+} // namespace arb
diff --git a/python/error.hpp b/python/error.hpp
index d475f2fa64eb6d81316dc80e5e576cee0e1c04e0..c97dad4331812268066aaa43676d51f0f0fdd631 100644
--- a/python/error.hpp
+++ b/python/error.hpp
@@ -4,6 +4,13 @@
#include <stdexcept>
#include <string>
+#include <pybind11/pybind11.h>
+
+#include <arbor/arbexcept.hpp>
+#include <arbor/util/either.hpp>
+
+#include "strprintf.hpp"
+
namespace pyarb {
extern std::exception_ptr py_exception;
@@ -45,4 +52,67 @@ auto try_catch_pyexception(L func, const char* msg){
throw;
}
}
+
+template <typename T, typename E>
+struct hopefully {
+ using value_type = T;
+ using error_type = E;
+ arb::util::either<value_type, error_type> state;
+
+ hopefully(const hopefully&) = default;
+
+ hopefully(value_type x): state(std::move(x)) {}
+ hopefully(error_type x): state(std::move(x)) {}
+
+ const value_type& operator*() const {
+ return try_get();
+ }
+ value_type& operator*() {
+ return try_get();
+ }
+ const value_type* operator->() const {
+ return &try_get();
+ }
+ value_type* operator->() {
+ return &try_get();
+ }
+
+ operator bool() const {
+ return (bool)state;
+ }
+
+ const error_type& error() const {
+ try {
+ return state.template get<1>();
+ }
+ catch(arb::util::either_invalid_access& e) {
+ throw arb::arbor_internal_error(
+ "Attempt to get an error from a valid hopefully wrapper.");
+ }
+ }
+
+private:
+
+ const value_type& try_get() const {
+ try {
+ return state.template get<0>();
+ }
+ catch(arb::util::either_invalid_access& e) {
+ throw arbor_internal_error(util::pprintf(
+ "Attempt to unwrap a hopefully with error state '{}'",
+ error().message));
+ }
+ }
+ value_type& try_get() {
+ try {
+ return state.template get<0>();
+ }
+ catch(arb::util::either_invalid_access& e) {
+ throw arb::arbor_internal_error(util::pprintf(
+ "Attempt to unwrap a hopefully with error state '{}'",
+ error().message));
+ }
+ }
+};
+
} // namespace pyarb
diff --git a/python/example/ring.py b/python/example/ring.py
index aec4f5cbaccd6eb405e569599f1b30a2f9dd326b..94136094a93e0acc72a124383ce20cfe75439462 100644
--- a/python/example/ring.py
+++ b/python/example/ring.py
@@ -2,14 +2,53 @@ import sys
import arbor
import matplotlib.pyplot as plt
+# Construct a cell with the following morphology.
+# The soma (at the root of the tree) is marked 's', and
+# the end of each branch i is marked 'bi'.
+#
+# b2
+# /
+# s----b1
+# \
+# b3
+
+def make_cable_cell(gid):
+ b = arbor.flat_cell_builder()
+
+ # Soma with radius 6 μm.
+ s = b.add_sphere(6, "soma")
+ # Single dendrite of length 100 μm and radius 2 μm attached to soma.
+ b1 = b.add_cable(parent=s, length=100, radius=2, name="dend", ncomp=1)
+ # Attach two dendrites of length 50 μm to the end of the first dendrite.
+ # Radius tapers from 2 to 0.5 μm over the length of the dendrite.
+ b2 = b.add_cable(parent=b1, length=50, radius=(2,0.5), name="dend", ncomp=1)
+ # Constant radius of 1 μm over the length of the dendrite.
+ b3 = b.add_cable(parent=b1, length=50, radius=1, name="dend", ncomp=1)
+
+ # Mark location for synapse at the midpoint of branch 1 (the first dendrite).
+ b.add_label('synapse_site', '(location 1 0.5)')
+ # Mark the root of the tree.
+ b.add_label('root', '(root)')
+
+ cell = b.build()
+
+ # Put hh dynamics on soma, and passive properties on the dendrites.
+ cell.paint('soma', 'hh')
+ cell.paint('dend', 'pas')
+ # Attach a single synapse.
+ cell.place('synapse_site', 'expsyn')
+ # Attach a spike detector with threshold of -10 mV.
+ cell.place('root', arbor.spike_detector(-10))
+
+ return cell
+
class ring_recipe (arbor.recipe):
- def __init__(self, n=4):
+ def __init__(self, n=10):
# The base C++ class constructor must be called first, to ensure that
# all memory in the C++ class is initialized correctly.
arbor.recipe.__init__(self)
self.ncells = n
- self.params = arbor.cell_parameters()
# The num_cells method that returns the total number of cells in the model
# must be implemented.
@@ -18,7 +57,7 @@ class ring_recipe (arbor.recipe):
# The cell_description method returns a cell
def cell_description(self, gid):
- return arbor.make_cable_cell(gid, self.params)
+ return make_cable_cell(gid)
def num_targets(self, gid):
return 1
@@ -35,7 +74,7 @@ class ring_recipe (arbor.recipe):
def connections_on(self, gid):
src = (gid-1)%self.ncells
w = 0.01
- d = 10
+ d = 5
return [arbor.connection(arbor.cell_member(src,0), arbor.cell_member(gid,0), w, d)]
# Attach a generator to the first cell in the ring.
@@ -59,7 +98,8 @@ print(context)
meters = arbor.meter_manager()
meters.start(context)
-recipe = ring_recipe(10)
+ncells = 4
+recipe = ring_recipe(ncells)
print(f'{recipe}')
meters.checkpoint('recipe-create', context)
@@ -84,37 +124,42 @@ meters.checkpoint('simulation-init', context)
spike_recorder = arbor.attach_spike_recorder(sim)
-pid = arbor.cell_member(0,0) # cell 0, probe 0
# Attach a sampler to the voltage probe on cell 0.
-# Sample rate of 1 sample every ms.
-sampler = arbor.attach_sampler(sim, 1, pid)
+# Sample rate of 10 sample every ms.
+samplers = [arbor.attach_sampler(sim, 0.1, arbor.cell_member(gid,0)) for gid in range(ncells)]
-sim.run(100)
+tfinal=1
+sim.run(tfinal)
print(f'{sim} finished')
meters.checkpoint('simulation-run', context)
+# Print profiling information
print(f'{arbor.meter_report(meters, context)}')
+# Print spike times
+print('spikes:')
for sp in spike_recorder.spikes:
- print(sp)
+ print(' ', sp)
-print('voltage samples for probe id ', end = '')
-print(pid, end = '')
-print(':')
+# Plot the voltage trace at the soma of each cell.
+fig, ax = plt.subplots()
+for gid in range(ncells):
+ times = [s.time for s in samplers[gid].samples(arbor.cell_member(gid,0))]
+ volts = [s.value for s in samplers[gid].samples(arbor.cell_member(gid,0))]
+ ax.plot(times, volts, '.')
-time = []
-value = []
-for sa in sampler.samples(pid):
- print(sa)
- time.append(sa.time)
- value.append(sa.value)
+legends = ['cell {}'.format(gid) for gid in range(ncells)]
+ax.legend(legends)
-# plot the recorded voltages over time
-fig, ax = plt.subplots()
-ax.plot(time, value)
ax.set(xlabel='time (ms)', ylabel='voltage (mV)', title='ring demo')
-ax.legend(['voltage'])
-plt.xlim(0,100)
+plt.xlim(0,tfinal)
+plt.ylim(-80,40)
ax.grid()
-fig.savefig("voltages.png", dpi=300)
+
+plot_to_file=False
+if plot_to_file:
+ fig.savefig("voltages.png", dpi=300)
+ print('voltage samples saved to voltages.png')
+else:
+ plt.show()
diff --git a/python/example/single_cell_builder.py b/python/example/single_cell_builder.py
new file mode 100644
index 0000000000000000000000000000000000000000..7c7cd7a2d03ec79d87cba51870a5f5b62076af67
--- /dev/null
+++ b/python/example/single_cell_builder.py
@@ -0,0 +1,103 @@
+import arbor
+from arbor import mechanism as mech
+from arbor import location as loc
+import matplotlib.pyplot as plt
+
+# Make a ball and stick cell model
+b = arbor.flat_cell_builder()
+
+# Construct a cell with the following morphology.
+# The soma (at the root of the tree) is marked 's', and
+# the end of each branch i is marked 'bi'.
+#
+# b5
+# /
+# /
+# b2---b4
+# /
+# /
+# s-------b1
+# \
+# \
+# b3
+
+# Add a spherical soma with radius 6 um
+s = b.add_sphere(6, "soma")
+
+# Add the dendrite cables, labelling those closest to the soma "dendn",
+# and those furthest with "dendx" because we will set different electrical
+# properties for the two regions.
+b1 = b.add_cable(parent=s, length=100, radius=2, name="dendn", ncomp=100)
+# Radius tapers from 2 to 0.5 over the length of the branch.
+b2 = b.add_cable(parent=b1, length= 50, radius=(2,0.5), name="dendn", ncomp=50)
+b3 = b.add_cable(parent=b1, length= 50, radius=1, name="dendn", ncomp=50)
+b4 = b.add_cable(parent=b2, length= 50, radius=1, name="dendx", ncomp=50)
+b5 = b.add_cable(parent=b2, length= 50, radius=1, name="dendx", ncomp=50)
+
+# Combine the "dendn" and "dendx" regions into a single "dend" region.
+# The dendrites were labelled as such so that we can set different
+# properties on each sub-region, and then combined so that we can
+# set other properties on the whole dendrites.
+b.add_label('dend', '(join (region "dendn") (region "dendx"))')
+# Location of stimuli, in the middle of branch 2.
+b.add_label('stim_site', '(location 2 0.5)')
+# The root of the tree (equivalent to '(location 0 0)')
+b.add_label('root', '(root)')
+
+# Extract the cable cell from the builder.
+cell = b.build()
+
+# Set initial membrane potential everywhere on the cell to -40 mV.
+cell.set_properties(Vm=-40)
+# Put hh dynamics on soma, and passive properties on the dendrites.
+cell.paint('soma', 'hh')
+cell.paint('dend', 'pas')
+# Set axial resistivity in dendrite regions (Ohm.cm)
+cell.paint('dendn', rL=500)
+cell.paint('dendx', rL=10000)
+# Attach stimuli with duration of 2 ms and current of 0.8 nA.
+# There are three stimuli, which activate at 10 ms, 50 ms and 80 ms.
+cell.place('stim_site', arbor.iclamp( 10, 2, 0.8))
+cell.place('stim_site', arbor.iclamp( 50, 2, 0.8))
+cell.place('stim_site', arbor.iclamp( 80, 2, 0.8))
+# Add a spike detector with threshold of -10 mV.
+cell.place('root', arbor.spike_detector(-10))
+
+# make single cell model
+m = arbor.single_cell_model(cell)
+
+# Attach voltage probes, sampling at 10 kHz.
+m.probe('voltage', loc(0,0), 10000) # at the soma
+m.probe('voltage', loc(3,1), 10000) # at the end of branch 3
+m.probe('voltage', loc(4,1), 10000) # at the end of branch 4
+
+# Run simulation for 100 ms of simulated activity.
+tfinal=100
+m.run(tfinal)
+
+# Print spike times.
+if len(m.spikes)>0:
+ print('{} spikes:'.format(len(m.spikes)))
+ for s in m.spikes:
+ print(' {:7.4f}'.format(s))
+else:
+ print('no spikes')
+
+# Plot the recorded voltages over time.
+fig, ax = plt.subplots()
+for t in m.traces:
+ ax.plot(t.time, t.value)
+
+legend_labels = ['{}: {}'.format(s.variable, s.location) for s in m.traces]
+ax.legend(legend_labels)
+ax.set(xlabel='time (ms)', ylabel='voltage (mV)', title='cell builder demo')
+plt.xlim(0,tfinal)
+plt.ylim(-80,50)
+ax.grid()
+
+# Set to True to save the image to file instead of opening a plot window.
+plot_to_file=False
+if plot_to_file:
+ fig.savefig("voltages.png", dpi=300)
+else:
+ plt.show()
diff --git a/python/example/single_cell_swc.py b/python/example/single_cell_swc.py
new file mode 100644
index 0000000000000000000000000000000000000000..951f502f311a90bcd5f9db3f35a79277c69770db
--- /dev/null
+++ b/python/example/single_cell_swc.py
@@ -0,0 +1,78 @@
+import arbor
+from arbor import mechanism as mech
+from arbor import location as loc
+import matplotlib.pyplot as plt
+
+# Load a cell morphology from an swc file.
+# The model has 31 branches, including soma, dendrites and axon.
+tree = arbor.load_swc('../../test/unit/swc/example.swc')
+
+# Define the regions and locsets in the model.
+defs = {'soma': '(tag 1)', # soma has tag 1 in swc files.
+ 'axon': '(tag 2)', # axon has tag 1 in swc files.
+ 'dend': '(tag 3)', # dendrites have tag 1 in swc files.
+ 'root': '(root)', # the start of the soma in this morphology is at the root of the cell.
+ 'stim_site': '(location 1 0.5)'} # site for the stimulus, in the middle of branch 1.
+labels = arbor.label_dict(defs)
+
+# Combine morphology with region and locset definitions to make a cable cell.
+cell = arbor.cable_cell(tree, labels)
+
+# Set initial membrane potential to -55 mV
+cell.set_properties(Vm=-55)
+# Use Nernst to calculate reversal potential for calcium.
+cell.set_ion('ca', method=mech('nernst/x=ca'))
+# hh mechanism on the soma and axon.
+cell.paint('soma', 'hh')
+cell.paint('axon', 'hh')
+# pas mechanism the dendrites.
+cell.paint('dend', 'pas')
+# Increase resistivity on dendrites.
+cell.paint('dend', rL=500)
+# Attach stimuli that inject 0.8 nA currents for 1 ms, starting at 3 and 8 ms.
+cell.place('stim_site', arbor.iclamp(3, 1, current=0.8))
+cell.place('stim_site', arbor.iclamp(8, 1, 0.8))
+# Detect spikes at the soma with a voltage threshold of -10 mV.
+cell.place('root', arbor.spike_detector(-10))
+
+# Have one compartment between each sample point.
+cell.compartments_on_samples()
+
+# Make single cell model.
+m = arbor.single_cell_model(cell)
+
+# Attach voltage probes that sample at 50 kHz.
+m.probe('voltage', where=loc(0,0), frequency=50000)
+m.probe('voltage', where=loc(2,1), frequency=50000)
+m.probe('voltage', where=loc(4,1), frequency=50000)
+m.probe('voltage', where=loc(30,1), frequency=50000)
+
+# Simulate the cell for 15 ms.
+tfinal=15
+m.run(tfinal)
+
+# Print spike times.
+if len(m.spikes)>0:
+ print('{} spikes:'.format(len(m.spikes)))
+ for s in m.spikes:
+ print(' {:7.4f}'.format(s))
+else:
+ print('no spikes')
+
+# Plot the recorded voltages over time.
+fig, ax = plt.subplots()
+for t in m.traces:
+ ax.plot(t.time, t.value)
+
+legend_labels = ['{}: {}'.format(s.variable, s.location) for s in m.traces]
+ax.legend(legend_labels)
+ax.set(xlabel='time (ms)', ylabel='voltage (mV)', title='swc morphology demo')
+plt.xlim(0,tfinal)
+plt.ylim(-80,50)
+ax.grid()
+
+plot_to_file=False
+if plot_to_file:
+ fig.savefig("voltages.png", dpi=300)
+else:
+ plt.show()
diff --git a/python/flat_cell_builder.cpp b/python/flat_cell_builder.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..53662f86b25a948062f5548d4d1380cff52c1a99
--- /dev/null
+++ b/python/flat_cell_builder.cpp
@@ -0,0 +1,268 @@
+#include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+
+#include <arbor/cable_cell.hpp>
+#include <arbor/morph/morphology.hpp>
+#include <arbor/morph/primitives.hpp>
+#include <arbor/morph/sample_tree.hpp>
+
+#include "conversion.hpp"
+#include "error.hpp"
+#include "morph_parse.hpp"
+#include "s_expr.hpp"
+#include "strprintf.hpp"
+
+namespace pyarb {
+
+class flat_cell_builder {
+ // The sample tree describing the morphology: constructed on the fly as
+ // the cables/spheres are added with add_cable/add_sphere.
+ arb::sample_tree tree_;
+
+ // The distal sample id of each cable.
+ std::vector<arb::msize_t> cable_distal_id_;
+
+ // The number of unique region names used to label cables as they
+ // are added to the cell.
+ int tag_count_ = 0;
+ // Map from region names to the tag used to identify them.
+ std::unordered_map<std::string, int> tag_map_;
+
+ arb::label_dict dict_;
+
+ // The morphology is cached, and only updated on request when it is out of date.
+ mutable bool cached_morpho_ = true;
+ mutable arb::morphology morpho_;
+ mutable std::mutex mutex_;
+
+ // Set on construction and unchanged thereafter.
+ // Indicates whether
+ bool spherical_ = false;
+
+public:
+
+ flat_cell_builder() = default;
+
+ arb::msize_t add_sphere(double radius, const char* name) {
+ cached_morpho_ = false;
+ spherical_ = true;
+ if (size()) {
+ throw pyarb_error("Add soma to non-empty cell.");
+ }
+ tree_.append({{0,0,0,radius}, get_tag(name)});
+ cable_distal_id_.push_back(0);
+ return 0;
+ }
+
+ // Add a new branch that is attached to parent.
+ // Returns the id of the new cable.
+ arb::msize_t add_cable(arb::msize_t parent, double len,
+ double r1, double r2, const char* region, int ncomp)
+ {
+ using arb::mnpos;
+
+ cached_morpho_ = false;
+
+ if (!test_identifier(region)) {
+ throw pyarb_error(util::pprintf("'{}' is not a valid label name.", region));
+ }
+
+ // Get tag id of region (add a new tag if region does not already exist).
+ int tag = get_tag(region);
+ const bool at_root = parent==mnpos;
+
+ // Can't attach a cable to the root on a cell with spherical soma.
+ if (at_root && spherical_) {
+ throw pyarb_error("Invalid parent id.");
+ }
+ // Parent id must be in the range [0, size())
+ if (!at_root && parent>=size()) {
+ throw pyarb_error("Invalid parent id.");
+ }
+
+ // Calculating the sample that is the parent of this branch is
+ // neccesarily complicated to handle cable segments that are attached
+ // to the root.
+ arb::msize_t p = at_root? (size()? 0: mnpos): cable_distal_id_[parent];
+
+ double z = at_root? 0: // attach to root of non-spherical cell
+ spherical_&&!parent? soma_rad(): // attach to spherical root
+ tree_.samples()[p].loc.z; // attach to end of a cable
+
+ // Only add a first point at the very beginning of the cable if
+ // the cable is not attached to another
+ const bool add_first_point = p==arb::mnpos // attached to the "root"
+ || (!p && spherical_) // attached to a spherical root
+ || (r1!=tree_.samples()[p].loc.radius);
+ // proximal radius does not match r1
+ if (add_first_point) {
+ p = tree_.append(p, {{0,0,z,r1}, tag});
+ }
+ if (ncomp>1) {
+ double dz = len/ncomp;
+ double dr = (r2-r1)/ncomp;
+ for (auto i=1; i<ncomp; ++i) {
+ p = tree_.append(p, {{0,0,z+i*dz, r1+i*dr}, tag});
+ }
+ }
+ p = tree_.append(p, {{0,0,z+len,r2}, tag});
+ cable_distal_id_.push_back(p);
+
+ return size()-1;
+ }
+
+ void add_label(const char* name, const char* description) {
+ if (!test_identifier(name)) {
+ throw pyarb_error(util::pprintf("'{}' is not a valid label name.", name));
+ }
+
+ if (auto result = eval(parse(description)) ) {
+ // The description is a region.
+ if (result->type()==typeid(arb::region)) {
+ if (dict_.locset(name)) {
+ throw pyarb_error("Region name clashes with a locset.");
+ }
+ auto& reg = arb::util::any_cast<arb::region&>(*result);
+ if (auto r = dict_.region(name)) {
+ dict_.set(name, join(std::move(reg), std::move(*r)));
+ }
+ else {
+ dict_.set(name, std::move(reg));
+ }
+ }
+ else if (result->type()==typeid(arb::locset)) {
+ if (dict_.region(name)) {
+ throw pyarb_error("Locset name clashes with a region.");
+ }
+ auto& loc = arb::util::any_cast<arb::locset&>(*result);
+ if (auto l = dict_.locset(name)) {
+ dict_.set(name, sum(std::move(loc), std::move(*l)));
+ }
+ else {
+ dict_.set(name, std::move(loc));
+ }
+ }
+ else {
+ throw pyarb_error("Label describes neither a region nor a locset.");
+ }
+ }
+ else {
+ throw pyarb_error(result.error().message);
+ }
+ }
+
+ const arb::sample_tree& samples() const {
+ return tree_;
+ }
+
+ std::unordered_map<std::string, std::string> labels() const {
+ std::unordered_map<std::string, std::string> map;
+ for (auto& r: dict_.regions()) {
+ map[r.first] = util::pprintf("{}", r.second);
+ }
+ for (auto& l: dict_.locsets()) {
+ map[l.first] = util::pprintf("{}", l.second);
+ }
+
+ return map;
+ }
+
+ const arb::morphology& morphology() const {
+ const std::lock_guard<std::mutex> guard(mutex_);
+ if (!cached_morpho_) {
+ morpho_ = arb::morphology(tree_, spherical_);
+ cached_morpho_ = true;
+ }
+ return morpho_;
+ }
+
+ arb::cable_cell build() const {
+ // Make cable_cell from sample tree and dictionary.
+ auto c = arb::cable_cell(morphology(), dict_);
+ c.default_parameters.discretization = arb::cv_policy_every_sample{};
+ return c;
+ }
+
+ private:
+
+ // Get tag id of region with name.
+ // Add a new tag if region with that name has not already had a tag associated with it.
+ int get_tag(const std::string& name) {
+ using arb::reg::tagged;
+
+ // Name is in the map: return the tag.
+ auto it = tag_map_.find(name);
+ if (it!=tag_map_.end()) {
+ return it->second;
+ }
+ // If the name is not in the map, the next step depends on
+ // whether the name is used for a locst, or a region that is
+ // not in the map, or is not used at all.
+
+ // Name is a locset: error.
+ if (dict_.locset(name)) {
+ throw pyarb_error(util::pprintf("'{}' is a label for a locset."));
+ }
+ // Name is a region: add tag to region definition.
+ else if(auto reg = dict_.region(name)) {
+ tag_map_[name] = ++tag_count_;
+ dict_.set(name, join(*reg, tagged(tag_count_)));
+ return tag_count_;
+ }
+ // Name has not been registerd: make a unique tag and new region.
+ else {
+ tag_map_[name] = ++tag_count_;
+ dict_.set(name, tagged(tag_count_));
+ return tag_count_;
+ }
+ }
+
+ // Only valid if called on a non-empty tree with spherical soma.
+ double soma_rad() const {
+ return tree_.samples()[0].loc.radius;
+ }
+
+ // The number of cable segements (plus one optional soma) used to construct the cell.
+ std::size_t size() const {
+ return cable_distal_id_.size();
+ }
+};
+
+void register_flat_builder(pybind11::module& m) {
+ using namespace pybind11::literals;
+
+ pybind11::class_<flat_cell_builder> builder(m, "flat_cell_builder");
+ builder
+ .def(pybind11::init<>())
+ .def("add_sphere", &flat_cell_builder::add_sphere,
+ "radius"_a, "name"_a)
+ .def("add_cable",
+ [](flat_cell_builder& b, arb::msize_t p, double len, pybind11::object rad, const char* name, int ncomp) {
+ using pybind11::isinstance;
+ using pybind11::cast;
+ if (auto radius = try_cast<double>(rad) ) {
+ return b.add_cable(p, len, *radius, *radius, name, ncomp);
+ }
+
+ if (auto radii = try_cast<std::pair<double, double>>(rad)) {
+ return b.add_cable(p, len, radii->first, radii->second, name, ncomp);
+ }
+ else {
+ throw pyarb_error(
+ "Radius parameter is not a scalar (constant branch radius) or "
+ "a tuple (radius at proximal and distal ends respectively).");
+ }
+ },
+ "parent"_a, "length"_a, "radius"_a, "name"_a, "ncomp"_a=1)
+ .def("add_label", &flat_cell_builder::add_label,
+ "name"_a, "description"_a)
+ .def_property_readonly("samples",
+ [](const flat_cell_builder& b) { return b.samples(); })
+ .def_property_readonly("labels",
+ [](const flat_cell_builder& b) { return b.labels(); })
+ .def_property_readonly("morphology",
+ [](const flat_cell_builder& b) { return b.morphology(); })
+ .def("build", &flat_cell_builder::build);
+}
+
+} // namespace pyarb
diff --git a/python/morph_parse.cpp b/python/morph_parse.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..bef6fef2b5c650d86cdbb5209deb0b41f810f70d
--- /dev/null
+++ b/python/morph_parse.cpp
@@ -0,0 +1,349 @@
+#include <arbor/util/any.hpp>
+#include <arbor/morph/region.hpp>
+#include <arbor/morph/locset.hpp>
+
+#include "error.hpp"
+#include "s_expr.hpp"
+#include "morph_parse.hpp"
+
+namespace pyarb {
+
+struct nil_tag {};
+
+template <typename T>
+bool match(const std::type_info& info) {
+ return info == typeid(T);
+}
+
+template <>
+bool match<double>(const std::type_info& info) {
+ return info == typeid(double) || info == typeid(int);
+}
+
+template <>
+bool match<arb::region>(const std::type_info& info) {
+ return info == typeid(arb::region) || info == typeid(nil_tag);
+}
+
+template <>
+bool match<arb::locset>(const std::type_info& info) {
+ return info == typeid(arb::locset) || info == typeid(nil_tag);
+}
+
+template <typename T>
+T eval_cast(arb::util::any arg) {
+ return std::move(arb::util::any_cast<T&>(arg));
+}
+
+template <>
+double eval_cast<double>(arb::util::any arg) {
+ if (arg.type()==typeid(int)) return arb::util::any_cast<int>(arg);
+ return arb::util::any_cast<double>(arg);
+}
+
+template <>
+arb::region eval_cast<arb::region>(arb::util::any arg) {
+ if (arg.type()==typeid(arb::region)) return arb::util::any_cast<arb::region>(arg);
+ return arb::reg::nil();
+}
+
+template <>
+arb::locset eval_cast<arb::locset>(arb::util::any arg) {
+ if (arg.type()==typeid(arb::locset)) return arb::util::any_cast<arb::locset>(arg);
+ return arb::ls::nil();
+}
+
+template <typename... Args>
+struct call_eval {
+ using ftype = std::function<arb::util::any(Args...)>;
+ ftype f;
+ call_eval(ftype f): f(std::move(f)) {}
+
+ template<std::size_t... I>
+ arb::util::any expand_args_then_eval(std::vector<arb::util::any> args, std::index_sequence<I...>) {
+ return f(eval_cast<Args>(std::move(args[I]))...);
+ }
+
+ arb::util::any operator()(std::vector<arb::util::any> args) {
+ return expand_args_then_eval(std::move(args), std::make_index_sequence<sizeof...(Args)>());
+ }
+};
+
+template <typename... Args>
+struct call_match {
+ template <std::size_t I, typename T, typename Q, typename... Rest>
+ bool match_args_impl(const std::vector<arb::util::any>& args) const {
+ return match<T>(args[I].type()) && match_args_impl<I+1, Q, Rest...>(args);
+ }
+
+ template <std::size_t I, typename T>
+ bool match_args_impl(const std::vector<arb::util::any>& args) const {
+ return match<T>(args[I].type());
+ }
+
+ template <std::size_t I>
+ bool match_args_impl(const std::vector<arb::util::any>& args) const {
+ return true;
+ }
+
+ bool operator()(const std::vector<arb::util::any>& args) const {
+ const auto nargs_in = args.size();
+ const auto nargs_ex = sizeof...(Args);
+ return nargs_in==nargs_ex? match_args_impl<0, Args...>(args): false;
+ }
+};
+
+template <typename T>
+struct fold_eval {
+ using fold_fn = std::function<T(T, T)>;
+ fold_fn f;
+
+ using anyvec = std::vector<arb::util::any>;
+ using iterator = anyvec::iterator;
+
+ fold_eval(fold_fn f): f(std::move(f)) {}
+
+ T fold_impl(iterator left, iterator right) {
+ if (std::distance(left,right)==1u) {
+ return eval_cast<T>(std::move(*left));
+ }
+ return f(eval_cast<T>(std::move(*left)), fold_impl(left+1, right));
+ }
+
+ arb::util::any operator()(anyvec args) {
+ return fold_impl(args.begin(), args.end());
+ }
+};
+
+template <typename T>
+struct fold_match {
+ using anyvec = std::vector<arb::util::any>;
+ bool operator()(const anyvec& args) const {
+ if (args.size()<2u) return false;
+ for (auto& a: args) {
+ if (!match<T>(a.type())) return false;
+ }
+ return true;
+ }
+};
+
+struct evaluator {
+ using any_vec = std::vector<arb::util::any>;
+ using eval_fn = std::function<arb::util::any(any_vec)>;
+ using args_fn = std::function<bool(const any_vec&)>;
+
+ eval_fn eval;
+ args_fn match_args;
+ const char* message;
+
+ evaluator(eval_fn f, args_fn a, const char* m):
+ eval(std::move(f)),
+ match_args(std::move(a)),
+ message(m)
+ {}
+};
+
+template <typename... Args>
+struct make_call {
+ evaluator state;
+
+ template <typename F>
+ make_call(F&& f, const char* msg="call"):
+ state(call_eval<Args...>(std::forward<F>(f)), call_match<Args...>(), msg)
+ {}
+
+ operator evaluator() const {
+ return state;
+ }
+};
+
+template <typename T>
+struct make_fold {
+ evaluator state;
+
+ template <typename F>
+ make_fold(F&& f, const char* msg="fold"):
+ state(fold_eval<T>(std::forward<F>(f)), fold_match<T>(), msg)
+ {}
+
+ operator evaluator() const {
+ return state;
+ }
+};
+
+std::unordered_multimap<std::string, evaluator> eval_map {
+ // Functions that return regions
+ {"nil", make_call<>(arb::reg::nil,
+ "'nil' with 0 arguments")},
+ {"all", make_call<>(arb::reg::all,
+ "'all' with 0 arguments")},
+ {"tag", make_call<int>(arb::reg::tagged,
+ "'tag' with 1 argment: (tag_id:integer)")},
+ {"branch", make_call<int>(arb::reg::branch,
+ "'branch' with 1 argument: (branch_id:integer)")},
+ {"cable", make_call<int, double, double>(arb::reg::cable,
+ "'cable' with 3 arguments: (branch_id:integer prox:real dist:real)")},
+ {"region", make_call<std::string>(arb::reg::named,
+ "'region' with 1 argument: (name:string)")},
+ {"distal_interval", make_call<arb::locset, double>(arb::reg::distal_interval,
+ "'distal_interval' with 2 arguments: (start:locset extent:real)")},
+ {"proximal_interval",make_call<arb::locset, double>(arb::reg::proximal_interval,
+ "'proximal_interval' with 2 arguments: (start:locset extent:real)")},
+ {"radius_lt",make_call<arb::region, double>(arb::reg::radius_lt,
+ "'radius_lt' with 2 arguments: (reg:region radius:real)")},
+ {"radius_le",make_call<arb::region, double>(arb::reg::radius_le,
+ "'radius_le' with 2 arguments: (reg:region radius:real)")},
+ {"radius_gt",make_call<arb::region, double>(arb::reg::radius_gt,
+ "'radius_gt' with 2 arguments: (reg:region radius:real)")},
+ {"radius_ge",make_call<arb::region, double>(arb::reg::radius_ge,
+ "'radius_ge' with 2 arguments: (reg:region radius:real)")},
+ {"join", make_fold<arb::region>(static_cast<arb::region(*)(arb::region, arb::region)>(arb::join),
+ "'join' with at least 2 arguments: (region region [...region])")},
+ {"intersect",make_fold<arb::region>(static_cast<arb::region(*)(arb::region, arb::region)>(arb::intersect),
+ "'intersect' with at least 2 arguments: (region region [...region])")},
+
+ // Functions that return locsets
+ {"root", make_call<>(arb::ls::root,
+ "'root' with 0 arguments")},
+
+ {"location", make_call<int, double>([](int bid, double pos){return arb::ls::location(arb::msize_t(bid), pos);},
+ "'location' with 2 arguments: (branch_id:integer position:real)")},
+ {"terminal", make_call<>(arb::ls::terminal,
+ "'terminal' with 0 arguments")},
+ {"sample", make_call<int>(arb::ls::sample,
+ "'sample' with 1 argument: (sample_id:integer)")},
+ {"distal", make_call<arb::region>(arb::ls::most_distal,
+ "'distal' with 1 argument: (reg:region)")},
+ {"proximal",make_call<arb::region>(arb::ls::most_proximal,
+ "'proximal' with 1 argument: (reg:region)")},
+ {"uniform",make_call<arb::region, int, int, int>(arb::ls::uniform,
+ "'uniform' with 4 arguments: (reg:region, first:int, last:int, seed:int)")},
+ {"on_branches",make_call<double>(arb::ls::on_branches,
+ "'on_branches' with 1 argument: (pos:double)")},
+ {"locset", make_call<std::string>(arb::ls::named,
+ "'locset' with 1 argument: (name:string)")},
+ {"join", make_fold<arb::locset>(static_cast<arb::locset(*)(arb::locset, arb::locset)>(arb::join),
+ "'join' with at least 2 arguments: (locset locset [...locset])")},
+ {"sum", make_fold<arb::locset>(static_cast<arb::locset(*)(arb::locset, arb::locset)>(arb::sum),
+ "'sum' with at least 2 arguments: (locset locset [...locset])")},
+};
+
+parse_hopefully<arb::util::any> eval(const s_expr& e);
+
+parse_hopefully<std::vector<arb::util::any>> eval_args(const s_expr& e) {
+ if (!e) return {std::vector<arb::util::any>{}}; // empty argument list
+ const s_expr* h = &e;
+ std::vector<arb::util::any> args;
+ while (*h) {
+ auto arg = eval(h->head());
+ if (!arg) return std::move(arg.error());
+ args.push_back(std::move(*arg));
+ h = &h->tail();
+ }
+ return args;
+}
+
+// Generate a string description of a function evaluation of the form:
+// Example output:
+// 'foo' with 1 argument: (real)
+// 'bar' with 0 arguments
+// 'cat' with 3 arguments: (locset region integer)
+// Where 'foo', 'bar' and 'cat' are the name of the function, and the
+// types (integer, real, region, locset) are inferred from the arguments.
+std::string eval_description(const char* name, const std::vector<arb::util::any>& args) {
+ auto type_string = [](const std::type_info& t) -> const char* {
+ if (t==typeid(int)) return "integer";
+ if (t==typeid(double)) return "real";
+ if (t==typeid(arb::region)) return "region";
+ if (t==typeid(arb::locset)) return "locset";
+ if (t==typeid(nil_tag)) return "()";
+ return "unknown";
+ };
+
+ const auto nargs = args.size();
+ std::string msg =
+ util::pprintf("'{}' with {} argument{}",
+ name, nargs,
+ nargs==0?"s": nargs==1u?":": "s:");
+ if (nargs) {
+ msg += " (";
+ bool first = true;
+ for (auto& a: args) {
+ msg += util::pprintf("{}{}", first?"":" ", type_string(a.type()));
+ first = false;
+ }
+ msg += ")";
+ }
+ return msg;
+}
+
+// Evaluate an s expression.
+// On success the result is wrapped in util::any, where the result is one of:
+// int : an integer atom
+// double : a real atom
+// arb::region : a region
+// arb::locset : a locset
+//
+// If there invalid input is detected, hopefully return value contains
+// a parse_error_state with an error string and location.
+//
+// If there was an unexpected/fatal error, an exception will be thrown.
+parse_hopefully<arb::util::any> eval(const s_expr& e) {
+ if (e.is_atom()) {
+ auto& t = e.atom();
+ switch (t.kind) {
+ case tok::integer:
+ return {std::stoi(t.spelling)};
+ case tok::real:
+ return {std::stod(t.spelling)};
+ case tok::nil:
+ return {nil_tag()};
+ case tok::string:
+ return {std::string(t.spelling)};
+ case tok::error:
+ return parse_error_state{e.atom().spelling, location(e)};
+ default:
+ return parse_error_state{
+ util::pprintf("Unexpected term: {}", e), location(e)};
+ }
+ }
+ if (e.head().is_atom()) {
+ // This must be a function evaluation, where head is the function name, and
+ // tail is a list of arguments.
+
+ // Evaluate the arguments, and return error state if an error ocurred.
+ auto args = eval_args(e.tail());
+ if (!args) {
+ return args.error();
+ }
+
+ // Find all candidate functions that match the name of the function.
+ auto& name = e.head().atom().spelling;
+ auto matches = eval_map.equal_range(name);
+
+ // Search for a candidate that matches the argument list.
+ for (auto i=matches.first; i!=matches.second; ++i) {
+ if (i->second.match_args(*args)) { // found a match: evaluate and return.
+ return i->second.eval(*args);
+ }
+ }
+
+ // Unable to find a match: try to return a helpful error message.
+ const auto nc = std::distance(matches.first, matches.second);
+ auto msg = util::pprintf("No matches for {}", eval_description(name.c_str(), *args));
+ msg += util::pprintf("\n There are {} potential candiates{}", nc, nc?":":".");
+ int count = 0;
+ for (auto i=matches.first; i!=matches.second; ++i) {
+ msg += util::pprintf("\n Candidate {} {}", ++count, i->second.message);
+ }
+ return parse_error_state{std::move(msg), location(e)};
+ }
+
+ return parse_error_state{
+ util::pprintf("Unable to evaluate '{}': expression must be either integer, real expression of the form (op <args>)", e),
+ location(e)};
+}
+
+} // namespace pyarb
+
+
diff --git a/python/morph_parse.hpp b/python/morph_parse.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..cb78cb57ee20f9f688e460201b70c44f3899d8f9
--- /dev/null
+++ b/python/morph_parse.hpp
@@ -0,0 +1,20 @@
+#pragma once
+
+#include "error.hpp"
+#include "s_expr.hpp"
+
+#include <arbor/util/any.hpp>
+
+namespace pyarb {
+
+struct parse_error_state {
+ std::string message;
+ int location;
+};
+
+template <typename T>
+using parse_hopefully = hopefully<T, parse_error_state>;
+
+parse_hopefully<arb::util::any> eval(const s_expr&);
+
+} // namespace pyarb
diff --git a/python/morphology.cpp b/python/morphology.cpp
index 6b995db640c0bdd79d71bc519cb1993c0fc1bb0a..0e503bf52701f225b0530d2ec1312b1a2f7d0a84 100644
--- a/python/morphology.cpp
+++ b/python/morphology.cpp
@@ -1,6 +1,12 @@
#include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+#include <fstream>
+
+#include <arbor/morph/morphology.hpp>
#include <arbor/morph/primitives.hpp>
+#include <arbor/morph/sample_tree.hpp>
+#include <arbor/swcio.hpp>
#include "error.hpp"
#include "strprintf.hpp"
@@ -10,7 +16,13 @@ namespace pyarb {
void register_morphology(pybind11::module& m) {
using namespace pybind11::literals;
- // segment location
+ //
+ // primitives: points, samples, locations, cables... etc.
+ //
+
+ m.attr("mnpos") = arb::mnpos;
+
+ // arb::mlocation
pybind11::class_<arb::mlocation> location(m, "location",
"A location on a cable cell.");
location
@@ -24,10 +36,178 @@ void register_morphology(pybind11::module& m) {
"Construct a location specification holding:\n"
" branch: The id of the branch.\n"
" position: The relative position (from 0., proximal, to 1., distal) on the branch.\n")
- .def_readonly("branch", &arb::mlocation::branch, "The id of the branch.")
- .def_readonly("position", &arb::mlocation::pos, "The relative position on the branch (∈ [0.,1.], where 0. means proximal and 1. distal).")
- .def("__str__", [](arb::mlocation l) {return util::pprintf("<arbor.location: branch {}, position {}>", l.branch, l.pos);})
- .def("__repr__", [](arb::mlocation l) {return util::pprintf("<arbor.location: branch {}, position {}>", l.branch, l.pos);});
+ .def_readonly("branch", &arb::mlocation::branch,
+ "The id of the branch.")
+ .def_readonly("position", &arb::mlocation::pos,
+ "The relative position on the branch (∈ [0.,1.], where 0. means proximal and 1. distal).")
+ .def("__str__",
+ [](arb::mlocation l) {
+ return util::pprintf("(location {} {})", l.branch, l.pos);
+ })
+ .def("__repr__",
+ [](arb::mlocation l) {
+ return util::pprintf("<arbor.location: branch {}, position {}>", l.branch, l.pos);
+ });
+
+ // arb::mpoint
+ pybind11::class_<arb::mpoint> mpoint(m, "mpoint");
+ mpoint
+ .def(pybind11::init(
+ [](double x, double y, double z, double r) {
+ return arb::mpoint{x,y,z,r};
+ }),
+ "x"_a, "y"_a, "z"_a, "radius"_a, "All values in μm.")
+ .def_readonly("x", &arb::mpoint::x, "X coordinate [μm].")
+ .def_readonly("y", &arb::mpoint::y, "Y coordinate [μm].")
+ .def_readonly("z", &arb::mpoint::z, "Z coordinate [μm].")
+ .def_readonly("radius", &arb::mpoint::radius,
+ "Radius of cable at sample location centered at coordinates [μm].")
+ .def("__str__",
+ [](const arb::mpoint& p) {
+ return util::pprintf("<arbor.mpoint: x {}, y {}, z {}, radius {}>", p.x, p.y, p.z, p.radius);
+ })
+ .def("__repr__",
+ [](const arb::mpoint& p) {return util::pprintf("{}>", p);});
+ // arb::msample
+ pybind11::class_<arb::msample> msample(m, "msample");
+ msample
+ .def(pybind11::init(
+ [](arb::mpoint loc, int tag){
+ return arb::msample{loc, tag};}),
+ "location"_a, "tag"_a)
+ .def(pybind11::init(
+ [](double x, double y, double z, double r, int tag){
+ return arb::msample{{x,y,z,r}, tag};}),
+ "x"_a, "y"_a, "z"_a, "radius"_a, "tag"_a,
+ "spatial values {x, y, z, radius} are in μm.")
+ .def_readonly("loc", &arb::msample::loc,
+ "Location of sample.")
+ .def_readonly("tag", &arb::msample::tag,
+ "Property tag of sample. "
+ "If loaded from standard SWC file the following tags are used: soma=1, axon=2, dendrite=3, apical dendrite=4, however arbitrary tags can be used.")
+ .def("__str__",
+ [](const arb::msample& s) {
+ return util::pprintf("<arbor.mpoint: x {}, y {}, z {}, radius {}, tag {}>",
+ s.loc.x, s.loc.y, s.loc.z, s.loc.radius, s.tag); })
+ .def("__repr__",
+ [](const arb::msample& s) {
+ return util::pprintf("{}>", s);});
+
+ // arb::mcable
+ pybind11::class_<arb::mcable> cable(m, "cable");
+ cable
+ .def(pybind11::init(
+ [](arb::msize_t bid, double prox, double dist) {
+ arb::mcable c{bid, prox, dist};
+ if (!test_invariants(c)) {
+ throw pyarb_error("Invalid cable description. Cable segments must have proximal and distal end points in the range [0,1].");
+ }
+ return c;
+ }),
+ "branch_id"_a, "prox"_a, "dist"_a)
+ .def_readonly("prox_pos", &arb::mcable::prox_pos,
+ "The relative position of the proximal end of the cable on its branch ∈ [0,1].")
+ .def_readonly("dist_pos", &arb::mcable::dist_pos,
+ "The relative position of the distal end of the cable on its branch ∈ [0,1].")
+ .def_readonly("branch", &arb::mcable::branch,
+ "The id of the branch on which the cable lies.")
+ .def("__str__", [](const arb::mcable& c) {
+ return util::pprintf("<arbor.cable: branch {}, prox {}, dist {}", c.branch, c.prox_pos, c.dist_pos); })
+ .def("__repr__", [](const arb::mcable& c) { return util::pprintf("{}", c); });
+
+ //
+ // Higher-level data structures (sample_tree, morphology)
+ //
+
+ // arb::sample_tree
+ pybind11::class_<arb::sample_tree> sample_tree(m, "sample_tree");
+ sample_tree
+ // constructors
+ .def(pybind11::init<>())
+ // modifiers
+ .def("reserve", &arb::sample_tree::reserve)
+ .def("append", [](arb::sample_tree& t, arb::msample s){return t.append(s);})
+ .def("append", [](arb::sample_tree& t, arb::msize_t p, arb::msample s){return t.append(p, s);})
+ // properties
+ .def_property_readonly("empty", [](const arb::sample_tree& st){return st.empty();},
+ "Indicates whether the sample tree is empty (i.e. whether it has size 0)")
+ .def_property_readonly("size", [](const arb::sample_tree& st){return st.size();},
+ "The number of samples in the sample tree.")
+ .def_property_readonly("parents", [](const arb::sample_tree& st){return st.parents();},
+ "A list with the parent index of each sample.")
+ .def("__str__", [](const arb::sample_tree& s) {
+ return util::pprintf("<arbor.sample_tree:\n{}>", s);});
+
+ // Function that creates a sample_tree from an swc file.
+ // Wraps calls to C++ functions arb::parse_swc_file() and arb::swc_as_sample_tree().
+ m.def("load_swc",
+ [](std::string fname) {
+ std::ifstream fid{fname};
+ if (!fid.good()) {
+ throw pyarb_error(util::pprintf("can't open file '{}'", fname));
+ }
+ try {
+ auto records = arb::parse_swc_file(fid);
+ arb::swc_canonicalize(records);
+ return arb::swc_as_sample_tree(records);
+ }
+ catch (arb::swc_error& e) {
+ // Try to produce helpful error messages for SWC parsing errors.
+ throw pyarb_error(
+ util::pprintf("error parsing line {} of '{}': {}.",
+ e.line_number, fname, e.what()));
+ }
+ },
+ "Load an swc file and convert to a sample_tree.");
+
+ // arb::morphology
+
+ pybind11::class_<arb::morphology> morph(m, "morphology");
+ morph
+ // constructors
+ .def(pybind11::init(
+ [](arb::sample_tree t){
+ return arb::morphology(std::move(t));
+ }))
+ .def(pybind11::init(
+ [](arb::sample_tree t, bool spherical_root) {
+ return arb::morphology(std::move(t), spherical_root);
+ }), "sample_tree"_a, "spherical_root"_a)
+ // morphology's interface is read-only by design, so most of it can
+ // be implemented as read-only properties.
+ .def_property_readonly("empty",
+ [](const arb::morphology& m){return m.empty();},
+ "A list with the parent index of each sample.")
+ .def_property_readonly("spherical_root",
+ [](const arb::morphology& m){return m.spherical_root();},
+ "Whether the root of the morphology is spherical.")
+ .def_property_readonly("num_branches",
+ [](const arb::morphology& m){return m.num_branches();},
+ "The number of branches in the morphology.")
+ .def_property_readonly("num_samples",
+ [](const arb::morphology& m){return m.num_samples();},
+ "The number of samples in the morphology.")
+ .def_property_readonly("samples",
+ [](const arb::morphology& m){return m.samples();},
+ "All of the samples in the morphology.")
+ .def_property_readonly("sample_parents",
+ [](const arb::morphology& m){return m.sample_parents();},
+ "The parent indexes of each sample.")
+ .def("branch_parent", &arb::morphology::branch_parent,
+ "i"_a, "The parent branch of branch i.")
+ .def("branch_children", &arb::morphology::branch_children,
+ "i"_a, "The child branches of branch i.")
+ .def("branch_indexes",
+ [](const arb::morphology& m, arb::msize_t i) {
+ auto p = m.branch_indexes(i);
+ return std::vector<arb::msize_t>(p.first, p.second);
+ },
+ "i"_a, "Range of indexes into the sample points in branch i.")
+ .def("__str__",
+ [](const arb::morphology& m) {
+ return util::pprintf("<arbor.morphology:\n{}>", m);
+ });
}
+
} // namespace pyarb
diff --git a/python/pyarb.cpp b/python/pyarb.cpp
index e17cbdfd919c6c68e948d6a4ad5b4c1cbc191703..bab7005fdc227c6188e3ba6593873a13cc539b23 100644
--- a/python/pyarb.cpp
+++ b/python/pyarb.cpp
@@ -11,6 +11,7 @@ void register_config(pybind11::module& m);
void register_contexts(pybind11::module& m);
void register_domain_decomposition(pybind11::module& m);
void register_event_generators(pybind11::module& m);
+void register_flat_builder(pybind11::module& m);
void register_identifiers(pybind11::module& m);
void register_morphology(pybind11::module& m);
void register_profiler(pybind11::module& m);
@@ -18,6 +19,7 @@ void register_recipe(pybind11::module& m);
void register_sampling(pybind11::module& m);
void register_schedules(pybind11::module& m);
void register_simulation(pybind11::module& m);
+void register_single_cell(pybind11::module& m);
void register_spike_handling(pybind11::module& m);
#ifdef ARB_MPI_ENABLED
@@ -26,7 +28,7 @@ void register_mpi(pybind11::module& m);
}
PYBIND11_MODULE(arbor, m) {
- m.doc() = "arbor: Python bindings for Arbor.";
+ m.doc() = "arbor: multicompartment neural network models.";
m.attr("__version__") = ARB_VERSION;
pyarb::register_cells(m);
@@ -34,6 +36,7 @@ PYBIND11_MODULE(arbor, m) {
pyarb::register_contexts(m);
pyarb::register_domain_decomposition(m);
pyarb::register_event_generators(m);
+ pyarb::register_flat_builder(m);
pyarb::register_identifiers(m);
pyarb::register_morphology(m);
pyarb::register_profiler(m);
@@ -41,8 +44,9 @@ PYBIND11_MODULE(arbor, m) {
pyarb::register_sampling(m);
pyarb::register_schedules(m);
pyarb::register_simulation(m);
+ pyarb::register_single_cell(m);
pyarb::register_spike_handling(m);
-
+
#ifdef ARB_MPI_ENABLED
pyarb::register_mpi(m);
#endif
diff --git a/python/recipe.cpp b/python/recipe.cpp
index 35283549326e2ea1d4e16ba02f93fb16d4d8732e..1a4c0f7a602161397ac01ea74ceade9fcd952874 100644
--- a/python/recipe.cpp
+++ b/python/recipe.cpp
@@ -29,18 +29,18 @@ arb::util::unique_any py_recipe_shim::get_cell_description(arb::cell_gid_type gi
"Python error already thrown");
}
-arb::probe_info cable_probe(std::string kind, arb::cell_member_type id, arb::mlocation loc) {
- arb::cell_probe_address::probe_kind pkind;
- if (kind == "voltage") {
- pkind = arb::cell_probe_address::probe_kind::membrane_voltage;
+arb::cell_probe_address::probe_kind probe_kind_from_string(const std::string& name) {
+ if (name == "voltage") {
+ return arb::cell_probe_address::probe_kind::membrane_voltage;
}
- else if (kind == "current") {
- pkind = arb::cell_probe_address::probe_kind::membrane_current;
+ else if (name == "current") {
+ return arb::cell_probe_address::probe_kind::membrane_current;
}
- else throw pyarb_error(
- util::pprintf(
- "invalid probe kind: {}, neither voltage nor current", kind));
+ else throw pyarb_error(util::pprintf("invalid probe kind: {}, neither voltage nor current", name));
+}
+arb::probe_info cable_probe(std::string kind, arb::cell_member_type id, arb::mlocation loc) {
+ auto pkind = probe_kind_from_string(kind);
arb::cell_probe_address probe{loc, pkind};
return arb::probe_info{id, pkind, probe};
};
diff --git a/python/recipe.hpp b/python/recipe.hpp
index ce1ef326f26d2952a64e2ffc67c9679d0b94bd2a..d2e177688d94dbaeb9d2c76acc577447d4ca4ec1 100644
--- a/python/recipe.hpp
+++ b/python/recipe.hpp
@@ -15,6 +15,8 @@
namespace pyarb {
+arb::probe_info cable_probe(std::string kind, arb::cell_member_type id, arb::mlocation loc);
+
// pyarb::recipe is the recipe interface used by Python.
// Calls that return generic types return pybind11::object, to avoid
// having to wrap some C++ types used by the C++ interface (specifically
@@ -166,7 +168,7 @@ public:
return try_catch_pyexception([&](){ return impl_->get_probe(id); }, msg);
}
- // TODO: wrap and make thread safe
+ // TODO: make thread safe
arb::util::any get_global_properties(arb::cell_kind kind) const override {
if (kind==arb::cell_kind::cable) {
arb::cable_cell_global_properties gprop;
diff --git a/python/s_expr.cpp b/python/s_expr.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..a7e40e8aa8298c1f583882fd86bcd41c45d4fea2
--- /dev/null
+++ b/python/s_expr.cpp
@@ -0,0 +1,404 @@
+#include <iostream>
+
+#include <cctype>
+#include <cstring>
+#include <string>
+#include <memory>
+#include <ostream>
+#include <vector>
+
+#include <arbor/util/either.hpp>
+#include <arbor/arbexcept.hpp>
+
+#include "s_expr.hpp"
+#include "strprintf.hpp"
+
+namespace pyarb {
+
+inline bool is_alphanumeric(char c) {
+ return std::isdigit(c) || std::isalpha(c);
+}
+inline bool is_plusminus(char c) {
+ return (c=='-' || c=='+');
+}
+
+std::ostream& operator<<(std::ostream& o, const tok& t) {
+ switch (t) {
+ case tok::nil: return o << "nil";
+ case tok::lparen: return o << "lparen";
+ case tok::rparen: return o << "rparen";
+ case tok::real: return o << "real";
+ case tok::integer:return o << "integer";
+ case tok::name: return o << "name";
+ case tok::string: return o << "string";
+ case tok::eof: return o << "eof";
+ case tok::error: return o << "error";
+ }
+ return o << "<unknown>";
+}
+
+std::ostream& operator<<(std::ostream& o, const token& t) {
+ return o << util::pprintf("{}", t.spelling);
+}
+
+//
+// lexer
+//
+
+struct parser_error: public arb::arbor_exception {
+ int loc;
+ parser_error(const std::string& msg, int l):
+ arbor_exception(msg), loc(l)
+ {}
+};
+
+static std::unordered_map<tok, const char*> tok_to_keyword = {
+ {tok::nil, "nil"},
+};
+
+static std::unordered_map<std::string, tok> keyword_to_tok = {
+ {"nil", tok::nil},
+};
+
+class lexer {
+ const char* data_;;
+ const char* end_;;
+ const char* current_;
+ int loc_;
+ token token_;
+
+public:
+
+ lexer(const char* s):
+ data_(s),
+ end_(data_ + std::strlen(data_)),
+ current_(data_)
+ {
+ // Prime the first token.
+ parse();
+ }
+
+ // Return the current token in the stream.
+ const token& current() {
+ return token_;
+ }
+
+ const token& next() {
+ parse();
+ return token_;
+ }
+
+private:
+
+ // Consume the and return the next token in the stream.
+ void parse() {
+ using namespace std::string_literals;
+
+ while (current_!=end_) {
+ loc_ = current_-data_;
+ switch (*current_) {
+ // end of file
+ case 0 : // end of string
+ token_ = {loc_, tok::eof, "eof"s};
+ return;
+
+ // white space
+ case ' ' :
+ case '\t' :
+ case '\v' :
+ case '\f' :
+ case '\n' :
+ character();
+ continue; // skip to next character
+
+ case '(':
+ token_ = {loc_, tok::lparen, {character()}};
+ return;
+ case ')':
+ token_ = {loc_, tok::rparen, {character()}};
+ return;
+ case 'a' ... 'z':
+ case 'A' ... 'Z':
+ token_ = name();
+ return;
+ case '0' ... '9':
+ token_ = number();
+ return;
+ case '"':
+ token_ = string();
+ return;
+ case '-':
+ case '+':
+ {
+ char c = current_[1];
+ if (std::isdigit(c) or c=='.') {
+ token_ = number();
+ return;
+ }
+ }
+ token_ = {loc_, tok::error,
+ util::pprintf("Unexpected character '{}'.", character())};
+ return;
+
+ default:
+ token_ = {loc_, tok::error,
+ util::pprintf("Unexpected character '{}'.", character())};
+ return;
+ }
+ }
+
+ if (current_!=end_) {
+ // todo: handle error: should never hit this
+ }
+ token_ = {loc_, tok::eof, "eof"s};
+ return;
+ }
+
+ // Parse alphanumeric sequence that starts with an alphabet character,
+ // and my contain alphabet, numeric or underscor '_' characters.
+ //
+ // Valid names:
+ // sub_dendrite
+ // temp_
+ // branch3
+ // A
+ // Invalid names:
+ // _cat ; can't start with underscore
+ // 2ndvar ; can't start with numeric character
+ //
+ // Returns the appropriate token kind if name is a keyword.
+ token name() {
+ std::string name;
+ char c = *current_;
+
+ // Assert that current position is at the start of an identifier
+ if( !(std::isalpha(c)) ) {
+ throw parser_error(
+ "Lexer attempting to read identifier when none is available", loc_);
+ }
+
+ name += c;
+ ++current_;
+ while(1) {
+ c = *current_;
+
+ if(is_alphanumeric(c) || c=='_') {
+ name += character();
+ }
+ else {
+ break;
+ }
+ }
+
+ // test if the name matches a keyword
+ auto it = keyword_to_tok.find(name.c_str());
+ if (it!=keyword_to_tok.end()) {
+ return {loc_, it->second, std::move(name)};
+ }
+ return {loc_, tok::name, std::move(name)};
+ }
+
+ token string() {
+ using namespace std::string_literals;
+
+ ++current_;
+ const char* begin = current_;
+ while (current_!=end_ && character()!='"');
+
+ if (current_==end_) return {loc_, tok::error, "string missing closing \""};
+
+ return {loc_, tok::string, std::string(begin, current_-1)};
+ }
+
+ token number() {
+ using namespace std::string_literals;
+
+ std::string str;
+ char c = *current_;
+
+ // Start counting the number of points in the number.
+ auto num_point = (c=='.' ? 1 : 0);
+ auto uses_scientific_notation = 0;
+
+ str += c;
+ current_++;
+ while(1) {
+ c = *current_;
+ if(std::isdigit(c)) {
+ str += c;
+ current_++;
+ }
+ else if(c=='.') {
+ if (++num_point>1) {
+ // Can't have more than one '.' in a number
+ return {int(current_-data_), tok::error, "unexpected '.'"s};
+ }
+ str += c;
+ current_++;
+ if(uses_scientific_notation) {
+ // Can't have a '.' in the mantissa
+ return {int(current_-data_), tok::error, "unexpected '.'"s};
+ }
+ }
+ else if(!uses_scientific_notation && (c=='e' || c=='E')) {
+ if(std::isdigit(current_[1]) ||
+ (is_plusminus(current_[1]) && std::isdigit(current_[2])))
+ {
+ uses_scientific_notation++;
+ str += c;
+ current_++;
+ // Consume the next char if +/-
+ if (is_plusminus(*current_)) {
+ str += *current_++;
+ }
+ }
+ else {
+ // the 'e' or 'E' is the beginning of a new token
+ break;
+ }
+ }
+ else {
+ break;
+ }
+ }
+
+ const bool is_real = uses_scientific_notation || num_point>0;
+ return {loc_, (is_real? tok::real: tok::integer), std::move(str)};
+ }
+
+ char character() {
+ return *current_++;
+ }
+};
+
+bool test_identifier(const char* in) {
+ lexer L(in);
+ auto x = L.current();
+ return x.kind==tok::name && x.spelling==in;
+}
+
+//
+// s expression members
+//
+
+bool s_expr::is_atom() const {
+ return (bool)state;
+}
+
+const token& s_expr::atom() const {
+ return state.get<0>();
+}
+
+const s_expr& s_expr::head() const {
+ return state.get<1>().head.get();
+}
+
+const s_expr& s_expr::tail() const {
+ return state.get<1>().tail.get();
+}
+
+s_expr& s_expr::head() {
+ return state.get<1>().head.get();
+}
+
+s_expr& s_expr::tail() {
+ return state.get<1>().tail.get();
+}
+
+s_expr::operator bool() const {
+ return !(is_atom() && atom().kind==tok::nil);
+}
+
+std::ostream& operator<<(std::ostream& o, const s_expr& x) {
+ if (x.is_atom()) return o << x.atom();
+#if 0 // print full tree with terminating 'nil'
+ return o << "(" << x.head() << " . " << x.tail() << ")";
+#else // print '(a . nil)' as 'a'
+ return x.tail()? o << "(" << x.head() << " . " << x.tail() << ")"
+ : o << x.head();
+#endif
+}
+
+std::size_t length(const s_expr& l) {
+ if (l.is_atom() && l) {
+ throw arb::arbor_internal_error(
+ util::pprintf("Internal error: can't take length of an atom in '{}'.", l));
+ }
+ if (!l) { // nil
+ return 0u;
+ }
+ return 1+length(l.tail());
+}
+
+int location(const s_expr& l) {
+ if (l.is_atom()) return l.atom().column;
+ return location(l.head());
+}
+
+//
+// parsing s expressions
+//
+
+// If there is a parsing error, then an atom with kind==tok::error is returned
+// with the error string in its spelling.
+s_expr parse(lexer& L) {
+ using namespace std::string_literals;
+
+ s_expr node;
+ auto t = L.current();
+
+ if (t.kind==tok::lparen) {
+ t = L.next();
+ s_expr* n = &node;
+ while (true) {
+ if (t.kind == tok::eof) {
+ return token{t.column, tok::error,
+ "Unexpected end of input. Missing a closing parenthesis ')'."};;
+ }
+ if (t.kind == tok::error) {
+ return t;
+ }
+ else if (t.kind == tok::rparen) {
+ *n = token{t.column, tok::nil, "nil"};
+ break;
+ }
+ else if (t.kind == tok::lparen) {
+ auto e = parse(L);
+ if (e.is_atom() && e.atom().kind==tok::error) return e;
+ *n = {std::move(e), {}};
+ }
+ else {
+ *n = {s_expr(t), {}};
+ }
+
+ n = &n->tail();
+ t = L.next();
+ }
+ }
+ else {
+ return token{t.column, tok::error, "Missing opening parenthesis'('."};;
+ }
+
+ return node;
+}
+
+s_expr parse(const char* in) {
+ lexer l(in);
+ s_expr result = parse(l);
+ const bool err = result.is_atom()? result.atom().kind==tok::error: false;
+ if (!err) {
+ auto t = l.next(); // pop the last rparen token.
+ if (t.kind!=tok::eof) {
+ return token{t.column, tok::error,
+ util::pprintf("Unexpected '{}' at the end of input.", t)};
+ }
+ }
+ return result;
+}
+
+s_expr parse(const std::string& in) {
+ return parse(in.c_str());
+}
+
+} // namespace pyarb
+
diff --git a/python/s_expr.hpp b/python/s_expr.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..f24e5f7cdd1f8b3fadb93bcc008356d69e3220fa
--- /dev/null
+++ b/python/s_expr.hpp
@@ -0,0 +1,125 @@
+#pragma once
+
+#include <string>
+#include <memory>
+#include <vector>
+
+#include <arbor/util/either.hpp>
+#include <arbor/util/optional.hpp>
+
+namespace pyarb {
+
+enum class tok {
+ nil,
+ real, // real number
+ integer, // integer
+ name, // name
+ lparen, // left parenthesis '('
+ rparen, // right parenthesis ')'
+ string, // string, written as "spelling"
+ eof, // end of file/input
+ error // special error state marker
+};
+
+std::ostream& operator<<(std::ostream&, const tok&);
+
+struct token {
+ int column;
+ tok kind;
+ std::string spelling;
+};
+
+std::ostream& operator<<(std::ostream&, const token&);
+
+std::vector<token> tokenize(const char* line);
+
+struct s_expr {
+ template <typename U>
+ struct s_pair {
+ U head = U();
+ U tail = U();
+ s_pair(U l, U r): head(std::move(l)), tail(std::move(r)) {}
+ };
+
+ // This value_wrapper is used to wrap the shared pointer
+ template <typename T>
+ struct value_wrapper{
+ using state_t = std::unique_ptr<T>;
+ state_t state;
+
+ value_wrapper() = default;
+
+ value_wrapper(const T& v):
+ state(std::make_unique<T>(v)) {}
+
+ value_wrapper(T&& v):
+ state(std::make_unique<T>(std::move(v))) {}
+
+ value_wrapper(const value_wrapper& other):
+ state(std::make_unique<T>(other.get())) {}
+
+ value_wrapper& operator=(const value_wrapper& other) {
+ state = std::make_unique<T>(other.get());
+ return *this;
+ }
+
+ value_wrapper(value_wrapper&& other) = default;
+
+ friend std::ostream& operator<<(std::ostream& o, const value_wrapper& w) {
+ return o << *w.state;
+ }
+
+ operator T() const {
+ return *state;
+ }
+
+ const T& get() const {
+ return *state;
+ }
+
+ T& get() {
+ return *state;
+ }
+ };
+
+ // An s_expr can be one of
+ // 1. an atom
+ // 2. a pair of s_expr (head and tail)
+ // The s_expr uses a util::either to represent these two possible states,
+ // which requires using an incomplete definition of s_expr, requiring
+ // with a std::shared_ptr.
+
+ using pair_type = s_pair<value_wrapper<s_expr>>;
+ arb::util::either<token, pair_type> state;
+
+ s_expr(const s_expr& s): state(s.state) {}
+ s_expr() = default;
+ s_expr(token t): state(std::move(t)) {}
+ s_expr(s_expr l, s_expr r):
+ state(pair_type(std::move(l), std::move(r)))
+ {}
+
+ bool is_atom() const;
+
+ const token& atom() const;
+
+ operator bool() const;
+
+ const s_expr& head() const;
+ const s_expr& tail() const;
+ s_expr& head();
+ s_expr& tail();
+
+ friend std::ostream& operator<<(std::ostream& o, const s_expr& x);
+};
+
+std::size_t length(const s_expr& l);
+int location(const s_expr& l);
+
+s_expr parse(const char* line);
+s_expr parse(const std::string& line);
+
+bool test_identifier(const char* in);
+
+} // namespace pyarb
+
diff --git a/python/single_cell_model.cpp b/python/single_cell_model.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..6d5084c47c558ee5e266eaeb9ad269dd0a408c51
--- /dev/null
+++ b/python/single_cell_model.cpp
@@ -0,0 +1,266 @@
+#include <string>
+#include <vector>
+
+#include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+
+#include <arbor/cable_cell.hpp>
+#include <arbor/load_balance.hpp>
+#include <arbor/recipe.hpp>
+#include <arbor/simulation.hpp>
+
+#include "error.hpp"
+
+namespace pyarb {
+
+//
+// Implementation of sampling infrastructure for Python single cell models.
+//
+// Note that only voltage sampling is currently supported, which simplifies the
+// code somewhat.
+//
+
+// Stores the location and sampling frequency for a probe in a single cell model.
+struct probe_site {
+ arb::mlocation site; // Location of sample on morphology.
+ double frequency; // Sampling frequency [Hz].
+};
+
+// Stores a single trace, which can be queried and viewed by the user at the end
+// of a simulation run.
+struct trace {
+ std::string variable; // Name of the variable being recorded.
+ arb::mlocation loc; // Sample site on morphology.
+ std::vector<arb::time_type> t; // Sample times [ms].
+ std::vector<double> v; // Sample values [units specific to sample variable].
+};
+
+// Callback provided to sampling API that records into a trace variable.
+struct trace_callback {
+ trace& trace_;
+
+ trace_callback(trace& t): trace_(t) {}
+
+ void operator()(arb::cell_member_type probe_id, arb::probe_tag tag, std::size_t n, const arb::sample_record* recs) {
+ // Push each (time, value) pair from the last epoch into trace_.
+ for (std::size_t i=0; i<n; ++i) {
+ if (auto p = arb::util::any_cast<const double*>(recs[i].data)) {
+ trace_.t.push_back(recs[i].time);
+ trace_.v.push_back(*p);
+ }
+ else {
+ throw std::runtime_error("unexpected sample type in simple_sampler");
+ }
+ }
+ }
+};
+
+// Used internally by the single cell model to expose model information to the
+// arb::simulation API when a model is instantiated.
+// Model descriptors, i.e. the cable_cell and probes, are instantiated
+// in the single_cell_model by user calls. The recipe is generated lazily, just
+// before simulation construction, so the recipe can use const references to all
+// of the model descriptors.
+struct single_cell_recipe: arb::recipe {
+ const arb::cable_cell& cell_;
+ const std::vector<probe_site>& probes_;
+ const arb::cable_cell_global_properties& gprop_;
+
+ single_cell_recipe(
+ const arb::cable_cell& c,
+ const std::vector<probe_site>& probes,
+ const arb::cable_cell_global_properties& props):
+ cell_(c), probes_(probes), gprop_(props)
+ {}
+
+ virtual arb::cell_size_type num_cells() const override {
+ return 1;
+ }
+
+ virtual arb::util::unique_any get_cell_description(arb::cell_gid_type gid) const override {
+ return cell_;
+ }
+
+ virtual arb::cell_kind get_cell_kind(arb::cell_gid_type) const override {
+ return arb::cell_kind::cable;
+ }
+
+ virtual arb::cell_size_type num_sources(arb::cell_gid_type) const override {
+ return cell_.detectors().size();
+ }
+
+ // synapses, connections and event generators
+
+ virtual arb::cell_size_type num_targets(arb::cell_gid_type) const override {
+ return cell_.synapses().size();
+ }
+
+ virtual std::vector<arb::cell_connection> connections_on(arb::cell_gid_type) const override {
+ return {}; // no connections on a single cell model
+ }
+
+ virtual std::vector<arb::event_generator> event_generators(arb::cell_gid_type) const override {
+ return {};
+ }
+
+ // probes
+
+ virtual arb::cell_size_type num_probes(arb::cell_gid_type) const override {
+ return probes_.size();
+ }
+
+ virtual arb::probe_info get_probe(arb::cell_member_type probe_id) const override {
+ // Test that a valid probe site is requested.
+ if (probe_id.gid || probe_id.index>=probes_.size()) {
+ throw arb::bad_probe_id(probe_id);
+ }
+
+ // For now only voltage can be selected for measurement.
+ auto kind = arb::cell_probe_address::membrane_voltage;
+ const auto& loc = probes_[probe_id.index].site;
+ return arb::probe_info{probe_id, kind, arb::cell_probe_address{loc, kind}};
+ }
+
+ // gap junctions
+
+ virtual arb::cell_size_type num_gap_junction_sites(arb::cell_gid_type gid) const override {
+ return 0; // No gap junctions on a single cell model.
+ }
+
+ virtual std::vector<arb::gap_junction_connection> gap_junctions_on(arb::cell_gid_type) const override {
+ return {}; // No gap junctions on a single cell model.
+ }
+
+ virtual arb::util::any get_global_properties(arb::cell_kind) const override {
+ return gprop_;
+ }
+};
+
+class single_cell_model {
+ arb::cable_cell cell_;
+ arb::context ctx_;
+ bool run_ = false;
+ arb::cable_cell_global_properties gprop_;
+
+ std::vector<probe_site> probes_;
+ std::unique_ptr<arb::simulation> sim_;
+ std::vector<double> spike_times_;
+ // Create one trace for each probe.
+ std::vector<trace> traces_;
+
+public:
+ single_cell_model(arb::cable_cell c):
+ cell_(std::move(c)), ctx_(arb::make_context())
+ {
+ gprop_.default_parameters = arb::neuron_parameter_defaults;
+ }
+
+ // example use:
+ // m.probe('voltage', arbor.location(2,0.5))
+ void probe(const std::string& what, const arb::mlocation& where, double frequency) {
+ if (what != "voltage") {
+ throw pyarb_error(
+ util::pprintf("{} does not name a valid variable to trace (currently only 'voltage' is supported)", what));
+ }
+ if (frequency<=0) {
+ throw pyarb_error(
+ util::pprintf("sampling frequency is not greater than zero", what));
+ }
+ if (where.branch>=cell_.morphology().num_branches()) {
+ throw pyarb_error(
+ util::pprintf("invalid location", what));
+ }
+ probes_.push_back({where, frequency});
+ }
+
+ void add_ion(const std::string& ion, double valence, double int_con, double ext_con, double rev_pot) {
+ gprop_.add_ion(ion, valence, int_con, ext_con, rev_pot);
+ }
+
+ void run(double tfinal) {
+ single_cell_recipe rec(cell_, probes_, gprop_);
+
+ auto domdec = arb::partition_load_balance(rec, ctx_);
+
+ sim_ = std::make_unique<arb::simulation>(rec, domdec, ctx_);
+
+ // Create one trace for each probe.
+ traces_.reserve(probes_.size());
+
+ // Add probes
+ for (arb::cell_lid_type i=0; i<probes_.size(); ++i) {
+ const auto& p = probes_[i];
+
+ traces_.push_back({"voltage", p.site, {}, {}});
+
+ auto sched = arb::regular_schedule(1000./p.frequency);
+
+ // Now attach the sampler at probe site, with sampling schedule sched, writing to voltage
+ sim_->add_sampler(arb::one_probe({0,i}), sched, trace_callback(traces_[i]));
+ }
+
+ // Set callback that records spike times.
+ sim_->set_global_spike_callback(
+ [this](const std::vector<arb::spike>& spikes) {
+ for (auto& s: spikes) {
+ spike_times_.push_back(s.time);
+ }
+ });
+
+ sim_->run(tfinal, 0.025);
+
+ run_ = true;
+ }
+
+ const std::vector<double>& spike_times() const {
+ return spike_times_;
+ }
+
+ const std::vector<trace>& traces() const {
+ return traces_;
+ }
+};
+
+void register_single_cell(pybind11::module& m) {
+ using namespace pybind11::literals;
+
+ pybind11::class_<trace> tr(m, "trace", "Values and meta-data for a sample-trace on a single cell model.");
+ tr
+ .def_readonly("variable", &trace::variable, "Name of the variable being recorded.")
+ .def_readonly("location", &trace::loc, "Location on cell morphology.")
+ .def_readonly("time", &trace::t, "Time stamps of samples [ms].")
+ .def_readonly("value", &trace::v, "Sample values.");
+
+ pybind11::class_<single_cell_model> model(m, "single_cell_model",
+ "Wrapper for simplified description, and execution, of single cell models.");
+
+ model
+ .def(pybind11::init<arb::cable_cell>(),
+ "cell"_a, "Initialise a single cell model for a cable cell.")
+ .def("run", &single_cell_model::run, "tfinal"_a, "Run model from t=0 to t=tfinal ms.")
+ .def("probe", &single_cell_model::probe,
+ "what"_a, "where"_a, "frequency"_a,
+ "Sample a variable on the cell.\n"
+ " what: Name of the variable to record (currently only 'voltage').\n"
+ " where: Location on cell morphology at which to sample the variable.\n"
+ " frequency: The target frequency at which to sample [Hz].")
+ .def("add_ion", &single_cell_model::add_ion,
+ "ion"_a, "valence"_a, "int_con"_a, "ext_con"_a, "rev_pot"_a,
+ "Add a new ion species to the model.\n"
+ " ion: name of the ion species.\n"
+ " valence: valence of the ion species.\n"
+ " int_con: initial internal concentration [mM].\n"
+ " ext_con: initial external concentration [mM].\n"
+ " rev_pot: reversal potential [mV].")
+ .def_property_readonly("spikes",
+ [](const single_cell_model& m) {
+ return m.spike_times();}, "Holds spike times [ms] after a call to run().")
+ .def_property_readonly("traces",
+ [](const single_cell_model& m) {
+ return m.traces();}, "Holds sample traces after a call to run().")
+ .def("__repr__", [](const single_cell_model&){return "<arbor.single_cell_model>";})
+ .def("__str__", [](const single_cell_model&){return "<arbor.single_cell_model>";});
+}
+
+} // namespace pyarb
+
diff --git a/python/strprintf.hpp b/python/strprintf.hpp
index e3b9dd49914318627f96c7f07c301df7f9b3fbbf..e8f96100fbfff0f74320eec7dc5eb751bdb645c1 100644
--- a/python/strprintf.hpp
+++ b/python/strprintf.hpp
@@ -7,6 +7,7 @@
#include <string>
#include <sstream>
#include <system_error>
+#include <unordered_map>
#include <utility>
#include <vector>
@@ -132,6 +133,25 @@ namespace impl {
}
};
+ template <typename Seq, typename F>
+ struct sepval_transform {
+ const Seq& seq_;
+ const char* sep_;
+ const F f_;
+
+ sepval_transform(const Seq& seq, const char* sep, F&& f): seq_(seq), sep_(sep), f_(std::forward(f)) {}
+
+ friend std::ostream& operator<<(std::ostream& o, const sepval_transform& s) {
+ bool first = true;
+ for (auto& x: s.seq_) {
+ if (!first) o << s.sep_;
+ first = false;
+ o << s.f(x);
+ }
+ return o;
+ }
+ };
+
template <typename Seq>
struct sepval_lim {
const Seq& seq_;
@@ -179,6 +199,21 @@ impl::sepval_lim<Seq> csv(const Seq& seq, unsigned n) {
return impl::sepval_lim<Seq>(seq, ", ", n);
}
+// Print dictionary: this could be done easily with range adaptors in C++17
+template <typename Key, typename T>
+std::string dictionary_csv(const std::unordered_map<Key, T>& dict) {
+ constexpr bool string_key = std::is_same<std::string, std::decay_t<Key>>::value;
+ std::string s = "{";
+ bool first = true;
+ for (auto& p: dict) {
+ if (!first) s += ", ";
+ first = false;
+ s += pprintf(string_key? "'{}': {}": "{}: {}", p.first, p.second);
+ }
+ s += "}";
+ return s;
+}
+
} // namespace util
template <typename T>
diff --git a/python/test/CMakeLists.txt b/python/test/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..1af9808f3b20e1f48d6652276c40a283c798ae3d
--- /dev/null
+++ b/python/test/CMakeLists.txt
@@ -0,0 +1,3 @@
+# Unit tests.
+# Builds: unit.
+add_subdirectory(cpp)
diff --git a/python/test/cpp/CMakeLists.txt b/python/test/cpp/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..291da50100446b139640fa0a2d6fe310878bba8e
--- /dev/null
+++ b/python/test/cpp/CMakeLists.txt
@@ -0,0 +1,27 @@
+set(py_unit_sources
+ s_expr.cpp
+
+ # unit test driver
+ test.cpp
+)
+
+add_executable(py_unit EXCLUDE_FROM_ALL ${py_unit_sources})
+add_dependencies(tests py_unit)
+
+add_library(py_unit_lib $<TARGET_OBJECTS:pyarb_obj>)
+target_link_libraries(py_unit_lib PRIVATE arbor pybind11::module)
+
+target_compile_options(py_unit PRIVATE ${ARB_CXXOPT_ARCH})
+target_include_directories(
+ py_unit PRIVATE
+ "${CMAKE_CURRENT_BINARY_DIR}"
+ "${PROJECT_SOURCE_DIR}/python"
+)
+
+target_link_libraries(
+ py_unit PRIVATE
+ gtest
+ arbor
+ py_unit_lib
+ pybind11::module
+)
diff --git a/python/test/cpp/s_expr.cpp b/python/test/cpp/s_expr.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..567ec754168b0c4eb07fc1816cbe35d998cf2aa4
--- /dev/null
+++ b/python/test/cpp/s_expr.cpp
@@ -0,0 +1,86 @@
+#include "gtest.h"
+
+#include <arbor/morph/region.hpp>
+#include <arbor/morph/locset.hpp>
+
+#include <morph_parse.hpp>
+#include <s_expr.hpp>
+
+using namespace pyarb;
+using namespace std::string_literals;
+
+TEST(s_expr, identifier) {
+ EXPECT_TRUE(test_identifier("foo"));
+ EXPECT_TRUE(test_identifier("f1"));
+ EXPECT_TRUE(test_identifier("f_"));
+ EXPECT_TRUE(test_identifier("f_1__"));
+ EXPECT_TRUE(test_identifier("A_1__"));
+
+ EXPECT_FALSE(test_identifier("_foobar"));
+ EXPECT_FALSE(test_identifier("2dogs"));
+ EXPECT_FALSE(test_identifier("1"));
+ EXPECT_FALSE(test_identifier("_"));
+ EXPECT_FALSE(test_identifier(""));
+ EXPECT_FALSE(test_identifier(" foo"));
+ EXPECT_FALSE(test_identifier("foo "));
+ EXPECT_FALSE(test_identifier("foo bar"));
+ EXPECT_FALSE(test_identifier("foo-bar"));
+ EXPECT_FALSE(test_identifier(""));
+}
+
+TEST(s_expr, atoms) {
+ auto get_atom = [](s_expr e) {
+ EXPECT_EQ(1u, length(e));
+ EXPECT_TRUE(e.head().is_atom());
+ return e.head().atom();
+ };
+
+ for (auto spelling: {"foo", "bar_", "car1", "car_1", "x_1__"}){
+ auto a = get_atom(parse("("s+spelling+")"));
+ EXPECT_EQ(a.kind, tok::name);
+ EXPECT_EQ(a.spelling, spelling);
+ }
+ // test parsing of integers
+ for (auto spelling: {"0", "-0", "1", "42", "-3287"}){
+ auto a = get_atom(parse("("s+spelling+")"));
+ EXPECT_EQ(a.kind, tok::integer);
+ EXPECT_EQ(a.spelling, spelling);
+ }
+ // test parsing of real numbers
+ for (auto spelling: {"0.", "-0.0", "1.21", "42.", "-3287.12"}){
+ auto a = get_atom(parse("("s+spelling+")"));
+ EXPECT_EQ(a.kind, tok::real);
+ EXPECT_EQ(a.spelling, spelling);
+ }
+ // test parsing of strings
+ for (auto spelling: {"foo", "dog cat", ""}) {
+ auto s = "(\""s+spelling+"\")";
+ auto a = get_atom(parse(s));
+ EXPECT_EQ(a.kind, tok::string);
+ }
+}
+
+TEST(s_expr, parse) {
+ auto round_trip_reg = [](const char* in) {
+ auto x = eval(parse(in));
+ return util::pprintf("{}", arb::util::any_cast<arb::region>(*x));
+ };
+ auto round_trip_loc = [](const char* in) {
+ auto x = eval(parse(in));
+ return util::pprintf("{}", arb::util::any_cast<arb::locset>(*x));
+ };
+
+ EXPECT_EQ("(cable 3 0 1)", round_trip_reg("(branch 3)"));
+ EXPECT_EQ("(cable 2 0.1 0.4)", round_trip_reg("(cable 2 0.1 0.4)"));
+ EXPECT_EQ("(all)", round_trip_reg("(all)"));
+ EXPECT_EQ("(region \"foo\")", round_trip_reg("(region \"foo\")"));
+
+ EXPECT_EQ("(terminal)", round_trip_loc("(terminal)"));
+ EXPECT_EQ("(root)", round_trip_loc("(root)"));
+ EXPECT_EQ("(locset \"cat_burgler\")", round_trip_loc("(locset \"cat_burgler\")"));
+
+ auto lhs = arb::util::any_cast<arb::region>(*eval(parse("(region \"dend\")")));
+ auto rhs = arb::util::any_cast<arb::region>(*eval(parse("(all)")));
+
+ EXPECT_EQ(util::pprintf("{}", join(lhs,rhs)), "(join (region \"dend\") (all))");
+}
diff --git a/python/test/cpp/test.cpp b/python/test/cpp/test.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..f2eab0931a0b4b2a5b17ed12ba6c2f3d1d09d27c
--- /dev/null
+++ b/python/test/cpp/test.cpp
@@ -0,0 +1,11 @@
+#include <iostream>
+#include <fstream>
+#include <numeric>
+#include <vector>
+
+#include "gtest.h"
+
+int main(int argc, char **argv) {
+ ::testing::InitGoogleTest(&argc, argv);
+ return RUN_ALL_TESTS();
+}
diff --git a/scripts/travis/build.sh b/scripts/travis/build.sh
index bacba3110db1945bbbb64ce1dfba67e13e4f926a..000e5e11662ddb9914b35408d43091f38f82e8b4 100755
--- a/scripts/travis/build.sh
+++ b/scripts/travis/build.sh
@@ -20,6 +20,8 @@ echo "compiler : ${compiler_version}"
echo "cmake : ${cmake_version}"
echo "build path : ${build_path}"
echo "base path : ${base_path}"
+echo "python3 : $(which python3)"
+echo "python3ver : $(python3 --version)"
if [[ "${WITH_DISTRIBUTED}" == "mpi" ]]; then
echo "mpi : on"
@@ -50,6 +52,13 @@ if [[ "${WITH_PYTHON}" == "true" ]]; then
python_path=$base_path/python
echo "python src : ${python_path}"
echo "PYTHONPATH : ${PYTHONPATH}"
+
+ if [[ "$TRAVIS_OS_NAME" = "linux" ]]; then
+ pypref=$(python-config --prefix)
+ PY_FLAGS="-DPYTHON_EXECUTABLE=$pypref/bin/python3.6 -DPYTHON_INCLUDE_DIR=$pypref/include/python3.6m -DPYTHON_LIBRARY=$pypref/lib/python3.6/config-3.6m-x86_64-linux-gnu/libpython3.6m.so"
+ else
+ PY_FLAGS=""
+ fi
else
echo "python : off"
ARB_WITH_PYTHON="OFF"
@@ -66,11 +75,11 @@ cd $build_path
#
progress "Configuring with cmake"
-cmake_flags="-DARB_WITH_ASSERTIONS=ON -DARB_WITH_MPI=${WITH_MPI} -DARB_WITH_PYTHON=${ARB_WITH_PYTHON} ${CXX_FLAGS}"
+cmake_flags="-DARB_WITH_ASSERTIONS=ON -DARB_WITH_MPI=${WITH_MPI} -DARB_WITH_PYTHON=${ARB_WITH_PYTHON} ${CXX_FLAGS} ${PY_FLAGS}"
echo "cmake flags: ${cmake_flags}"
cmake .. ${cmake_flags} || error "unable to configure cmake"
-export NMC_NUM_THREADS=2
+export ARB_NUM_THREADS=2
progress "C++ unit tests"
make unit -j4 || error "building unit tests"
@@ -96,8 +105,12 @@ if [[ "${WITH_PYTHON}" == "true" ]]; then
progress "Python unit tests"
python$PY $python_path/test/unit/runner.py -v2 || error "running python unit tests (serial)"
if [[ "${WITH_DISTRIBUTED}" = "mpi" ]]; then
- progress "Python distributed unit tests (MPI)"
- ${launch} python$PY $python_path/test/unit_distributed/runner.py -v2 || error "running python distributed unit tests (MPI)"
+ if [[ "$TRAVIS_OS_NAME" = "osx" ]]; then
+ progress "Python distributed unit tests (MPI)"
+ ${launch} python$PY $python_path/test/unit_distributed/runner.py -v2 || error "running python distributed unit tests (MPI)"
+ else
+ progress "Python distributed unit tests (MPI) -- skipping on Linux due to Travis config issues on Bionic."
+ fi
fi
fi