Skip to content
Snippets Groups Projects
Unverified Commit 8c85f006 authored by Thorsten Hater's avatar Thorsten Hater Committed by GitHub
Browse files

Add Developer Documentation (#1639)

Add a first round of developer documentation, evolving from the older `internal` part of the docs.
parent a9e4224b
No related branches found
No related tags found
No related merge requests found
Showing
with 689 additions and 47 deletions
......@@ -12,8 +12,12 @@ sphinx:
formats: []
build:
os: ubuntu-20.04
tools:
python: "3.7"
python:
version: 3.7
system_packages: true
install:
- requirements: doc/requirements.txt
\ No newline at end of file
set(GRAPHVIZ_IGNORE_TARGETS "bench;brunel;drybench;dryrun;gap_junctions;generators;lfp;probe-demo;ring;single;unit")
<div class="diag-cont">
<div id="recipe_src">Recipe
<div>Cell descriptions
<div>gid 1
<div>type <code>lif_cell_0</code>
<div>...</div>
</div>
</div>
<div>...</div>
<div> gid 37
<div>type <code>cable_cell_A</code>
<div>Morphology</div>
<div>Decor
<div>Mechanisms</div>
<div>Connection sites</div>
</div>
<div>Label dictionary</div>
</div>
</div>
</div>
<div>Cell kinds
<div>gid 1
<div>type <code>arbor.lif_cell</code>
</div>
</div>
<div>...</div>
<div>gid 37
<div>type <code>arbor.cable_cell</code></div>
</div>
</div>
<div>...
</div>
</div>
<div class="diag-col">
<div class="diag-right">
Simulator: requests cell description of <code>gid=37</code>
</div>
<div class="diag-left">
Recipe: <code>type=cable_cell_A</code>.
</div>
<div class="diag-right">
Simulator: requests cell kind of <code>gid=37</code>
</div>
<div class="diag-left">
Recipe: <code>kind=arbor.cable_cell</code>.
</div>
<div class="diag-right">
Simulator: look up cell group implementation for kind <code>arbor.cable_cell</code>.
</div>
<div class="diag-right2">
Domain decomposition: <code>cable_cell_group_gpu</code>
</div>
<div class="diag-right">
Simulator: <code>cable_cell_group_gpu</code> construct <code>cable_cell_A</code> object.
</div>
<div class="diag-right3">
<code>cable_cell_group_gpu</code>: Construction complete.
</div>
<div class="diag-right">
Simulator: <code>cable_cell_group_gpu</code> simulate for <code>t .. t + dt</code>
</div>
</div>
<div>Simulation
<div id="recipe_dst">Recipe</div>
<div id="conc_src">Context
<div>
12 threads
</div>
<div>
1 GPU
</div>
</div>
<div id="domdec_src">Domain decomposition
<div>
cable_cell_group_gpu
<div>1 GPU</div>
</div>
<div>
...
</div>
</div>
</div>
</div>
<div class="diag-caption">
An illustration of the cell-specific portion of the recipe, and how it is
used during the lifetime of the simulation: the simulation object will,
depending on its configuration, query the recipe for the neuroscientific
components it describes. This demonstration also show why the recipe separates
cell descriptions from cell types. The latter is, as you might expect,
shorthand, and is used in the allocation of the cell to a particular cell group.
A cell group implementation is a handler for a certain kind of cell, and Arbor
comes with these for all it's included cell kinds. However, users can develop
their own specialized cell group implementations. More on that in the internal
developer documentation.
</div>
......@@ -3,6 +3,11 @@
Domain decomposition
====================
An Arbor simulation requires a :ref:`modelrecipe`, a :ref:`(hardware) context <modelhardware>`, and a domain decomposition. The Recipe contains the neuroscientific model, the hardware context describes the computational resources you are going to execute the simulation on, and the domain decomposition describes how Arbor will use the hardware. Since the context and domain decomposition may seem closely related at first, it might be instructive to see how recipes are used by Arbor:
.. raw:: html
:file: domdec-diag-1.html
A *domain decomposition* describes the distribution of the model over the available computational resources.
The description partitions the cells in the model as follows:
......
.. _modelhardware:
Hardware management
===================
Hardware context
================
An Arbor simulation requires a :ref:`modelrecipe`, a (hardware) context, and a :ref:`modeldomdec`. The Recipe contains the neuroscientific model, the hardware context describes the computational resources you are going to execute the simulation on, and the domain decomposition describes how Arbor will use the hardware. Since the context and domain decomposition may seem closely related at first, it might be instructive to see how recipes are used by Arbor:
.. raw:: html
:file: domdec-diag-1.html
*Local resources* are locally available computational resources, specifically the number of hardware threads and the number of GPUs.
An *allocation* enumerates the computational resources to be used for a simulation, typically a subset of the resources available on a physical hardware node.
.. Note::
New users can find using contexts a little verbose.
The design is very deliberate, to allow fine-grained control over which
computational resources an Arbor simulation should use.
As a result Arbor is much easier to integrate into workflows that
run multiple applications or libraries on the same node, because
Arbor has a direct API for using on node resources (threads and GPU)
and distributed resources (MPI) that have been partitioned between
applications/libraries.
New users can find using contexts a little verbose.
The design is very deliberate, to allow fine-grained control over which
computational resources an Arbor simulation should use.
As a result Arbor is much easier to integrate into workflows that
run multiple applications or libraries on the same node, because
Arbor has a direct API for using on node resources (threads and GPU)
and distributed resources (MPI) that have been partitioned between
applications/libraries.
.. _modelcontext:
......
<div class="diag-spacer"></div>
<div class="diag-spacer"></div>
<div class="diag-cont">
<div id="recipe_src">Recipe
<div id="recipe_dsc" class="diag-note">describe the neuroscience</div>
<div>Cells
</div>
<div>Network
</div>
<div>In- & outputs
</div>
<div>...
</div>
</div>
<div id="sim_src">Simulation
<div id="sim_dsc" class="diag-note">describe the execution</div>
<div id="recipe_dst">Recipe</div>
<div class="diag-spacer"></div>
<div id="conc_src">Context
<div id="conc_dsc" class="diag-note">describe hardware</div>
</div>
<div class="diag-spacer"></div>
<div id="domdec_src">Domain decomposition
<div id="domdec_dsc" class="diag-note">describe how to use hardware</div>
</div>
</div>
</div>
<div class="diag-caption">An overview of Arbor. The two columns reflect a division that is core to Arbor: the neuroscience is described entirely separate from the execution of the simulation. The Recipe ties the neuroscience together: the user can provide any number of cells, each with a particular type, morphology and set of mechanisms; a network configuration; and in- & outputs like event generators and probes. This description, the recipe, is self-contained and can be executed by any configuration of execution. Execution is determined by the available hardware and instructions for how certain parts of a recipe (e.g. cell groups) are executed on that hardware.</div>
<connection width="1" from="#recipe_src" to="#recipe_dst" fromX="1" fromY="0" toX="0" toY="0" onlyVisible></connection>
<connection width="1" from="#recipe_src" to="#recipe_dst" fromX="1" fromY="1" toX="0" toY="1" onlyVisible></connection>
<connection from="#recipe_src" to="#recipe_dsc" color="rgb(221, 191, 146)" fromX="0.1" fromY="0" toX="0" tail onlyVisible></connection>
<connection from="#sim_src" to="#sim_dsc" color="rgb(221, 191, 146)" fromX="0.1" fromY="0" toX="0" tail onlyVisible></connection>
<connection from="#conc_src" to="#conc_dsc" color="rgb(221, 191, 146)" fromX="0.1" fromY="0" toX="0" tail onlyVisible></connection>
<connection from="#domdec_src" to="#domdec_dsc" color="rgb(221, 191, 146)" fromX="0.1" fromY="0" toX="0" tail onlyVisible></connection>
......@@ -3,6 +3,14 @@
Concepts overview
=================
Arbor is a library that lets you to model neural networks with morphologically
detailed cells; which it then executes the resulting simulation on a variety of
hardware. The execution can optionally be configured in high detail but comes
with sensible defaults.
.. raw:: html
:file: index-diag-1.html
To learn how to use Arbor, it is helpful to understand some of its concepts.
Arbor's design aims to enable scalability through abstraction.
To achieve this, Arbor makes a distinction between the **description** of a model, and the
......
<div class="diag-cont">
<div id="recipe_src">Recipe
<div>
<input id="collapsible1" type="checkbox" checked><label for="collapsible1">Cell descriptions</label>
<div>
<input id="collapsible10" type="checkbox"><label for="collapsible10">Cell 0</label>
<div>Morphology
<div>Segment Tree
<div>Segment 0 (soma)</div>
</div>
</div>
<div>Decor
<div>Paintables
<div>hh
<code>(all)</code>
</div>
</div>
<div>Placables
<div>iclamp <code>"midsoma"</code></div>
<div>detector_0 <code>"midsoma"</code></div>
</div>
</div>
<div>Label dictionary
<div>midsoma
<code>(location 0 0.5)</code>
</div>
</div>
</div>
<div>
<input id="collapsible11" type="checkbox"><label for="collapsible11">Cell 1</label>
<div>Morphology
<div>Segment Tree
<div>Segment 0 (soma)</div>
<div>Segment 1 (dendrite)</div>
<div>Segment 2 (dendrite)</div>
<div>Segment ...</div>
</div>
</div>
<div>Decor
<div>Paintables
<div>passive
<code>(all)</code>
</div>
</div>
<div>Placables
<div>synapse_1 <code>"endpoint"</code></div>
<div>gap_junction_1 <code>"midsoma"</code></div>
</div>
</div>
<div>Label dictionary
<div>midsoma
<code>(location 0 0.5)</code>
</div>
<div>endpoint
<code>(terminal)</code>
</div>
</div>
</div>
<div>
<input id="collapsible12" type="checkbox"><label for="collapsible12">Cell 2</label>
<div>Morphology
<div>Segment Tree
<div>Segment 0 (soma)</div>
<div>Segment 1 (dendrite)</div>
<div>Segment 2 (dendrite)</div>
<div>Segment ...</div>
</div>
</div>
<div>Decor
<div>Paintables
<div>passive
<code>(all)</code>
</div>
</div>
<div>Placables
<div>gap_junction_2 <code>"midsoma"</code></div>
</div>
</div>
<div>Label dictionary
<div>midsoma
<code>(location 0 0.5)</code>
</div>
</div>
</div>
</div>
<div>
<input id="collapsible2" type="checkbox"><label for="collapsible2">Cell kinds</label>
<div>Cell 0
<div>type <code>arbor.cable_cell</code></div>
</div>
<div>Cell 1
<div>type <code>arbor.cable_cell</code></div>
</div>
<div>Cell 2
<div>type <code>arbor.cable_cell</code></div>
</div>
</div>
<div>
<input id="collapsible3" type="checkbox"><label for="collapsible3">Connections</label>
<div>Connection 1
<div>from <code>detector_0</code></div>
<div>to <code>synapse_1</code></div>
</div>
<div>Gap Junction 1
<div>from <code>gap_junction_1</code></div>
<div>to <code>gap_junction_2</code></div>
</div>
</div>
<div>
<input id="collapsible4" type="checkbox"><label for="collapsible4">Probes</label>
<div>Probe 1
<div>gid
<code>2</code>
</div>
</div>
</div>
</div>
</div>
<div class="diag-caption">
An example recipe. Click "+" and "-" to (de)collapse contents.
</div>
......@@ -14,11 +14,11 @@ building phase to provide information about individual cells in the model, such
* **Gap junction connections** on each cell.
* **Probes** on each cell.
Recipes are structured to provide a consistent interface for describing each cell in the
network using their global identifier (`gid`).
This allows the simulator to be able to quickly look-up properties related to the connections
going in and out of a cell (think of synapses, gap junctions, but also probes and spike inputs);
which helps make Arbor fast and easily distributable over many nodes.
An example model
----------------
.. raw:: html
:file: recipe-diag-1.html
To better illustrate the content of a recipe, let's consider the following network of
three cells:
......@@ -33,7 +33,7 @@ three cells:
(because cable cells allow complex dynamics such as ``hh``). This is referred to as
the **kind** of the cell.
- ``Cell 1``: Is a soma and a single dendrite, with ``passive`` dynamics everywhere.
It has a single synapse at the end of the dendrite labeled "syanpse_1" and a gap
It has a single synapse at the end of the dendrite labeled "synapse_1" and a gap
junction mechanism in the middle of the soma labeled "gap_junction_1".
This is the **description** of the cell. It's also a cable cell, which is its **cell kind**.
- ``Cell 2``: Is a soma and a single dendrite, with ``passive`` dynamics everywhere.
......@@ -44,7 +44,7 @@ The total **number of cells** in the model is 3. The **kind**, and **description
is known and can be registered in the recipe. Next is the cell interaction.
The model is designed such that each cell has labeled source, target and gap junction sites.
A **network connection** can be formed from ``detector_0`` to ``synpase_1``; and a
A **network connection** can be formed from ``detector_0`` to ``synapse_1``; and a
**gap junction connection** between ``gap_junction_1`` and ``gap_junction_2``.
If ``detector_0`` spikes, a spike should be observed on ``gap_junction_2`` after some delay.
To monitor the voltage on ``gap_junction_2`` and record the spike, a **probe** can be set up
......@@ -59,12 +59,17 @@ Technical details of the recipe class are presented in the :ref:`Python <pyreci
:ref:`C++ <cpprecipe>` APIs.
Are recipes always necessary?
------------------------------
-----------------------------
Yes. They are a fundamental building block in Arbor simulations. Recipes are structured to provide
a consistent interface for describing each cell in the network using their global identifier (`gid`).
This allows the simulator to be able to quickly look-up properties related to the connections
going in and out of a cell (think of synapses, gap junctions, but also probes and spike inputs);
which helps make Arbor fast and easily distributable over many nodes.
Yes. However, we provide a python :class:`single_cell_model <py_single_cell_model>`
that abstracts away the details of a recipe for simulations of single, stand-alone
:ref:`cable cells<modelcablecell>`, which absolves the users from having to create the
recipe themselves. This is possible because the number of cells is fixed and known,
For single, stand-alone :ref:`cable cells<modelcablecell>` simulations, the Python API provides
a :class:`single_cell_model <py_single_cell_model>` that abstracts away the details of a recipe.
This is possible because the number of cells is fixed and known,
and it is guaranteed that there can be no connections or gap junctions in a model of a
single cell. The single cell model is able to fill out the details of the recipe under
the hood, and the user need only provide the cell description, and any probes they wish
......
......@@ -10,6 +10,8 @@ script_path=this_path+'/scripts'
sys.path.append(script_path)
html_static_path = ['static']
html_css_files = ['htmldiag.css']
html_js_files = ['domarrow.js']
def setup(app):
app.add_object_type('generic', 'gen', 'pair: %s; generic')
......
......@@ -130,25 +130,28 @@ Release
Post Release
------------
0. Update and submit Zenodo release if necessary.
1. Announce on our website
2. Announce on HBP newsletter newsletter@humanbrainproject.eu, HBP Twitter/socials evan.hancock@ebrains.eu
3. [AUTOMATED] Add tagged version of docs on ReadTheDocs
4. HBP internal admin
- Plus: https://plus.humanbrainproject.eu/components/2691/
- TC Wiki: https://wiki.ebrains.eu/bin/view/Collabs/technical-coordination/EBRAINS%20components/Arbor/
- KG: https://search.kg.ebrains.eu/instances/5cf4e24b-b0eb-4d05-96e5-a7751134a061
- Update howto: https://github.com/bweyers/HBPVisCatalogue/wiki/How-to-start-software-meta-data-curation%3F#update-curated-software
- Previous update as template: https://github.com/bweyers/HBPVisCatalogue/issues/480
- Supported file formats
- ContentTypes: https://humanbrainproject.github.io/openMINDS/v3/core/v4/data/contentType.html
- details: https://github.com/HumanBrainProject/openMINDS_core/tree/v3/instances/data/contentTypes
- Send an update to the folk in charge of HBP Twitter if we want to shout about it
5. FZJ admin
- https://juser.fz-juelich.de/submit
#. Update and submit Zenodo release if necessary.
#. Announce on our website
#. Announce on HBP newsletter newsletter@humanbrainproject.eu, HBP Twitter/socials evan.hancock@ebrains.eu
#. [AUTOMATED] Add tagged version of docs on ReadTheDocs
#. HBP internal admin
- Plus: https://plus.humanbrainproject.eu/components/2691/
- TC Wiki: https://wiki.ebrains.eu/bin/view/Collabs/technical-coordination/EBRAINS%20components/Arbor/
- KG: https://search.kg.ebrains.eu/instances/5cf4e24b-b0eb-4d05-96e5-a7751134a061
- Update howto: https://github.com/bweyers/HBPVisCatalogue/wiki/How-to-start-software-meta-data-curation%3F#update-curated-software
- Previous update as template: https://github.com/bweyers/HBPVisCatalogue/issues/480
- Supported file formats
- ContentTypes: https://humanbrainproject.github.io/openMINDS/v3/core/v4/data/contentType.html
- details: https://github.com/HumanBrainProject/openMINDS_core/tree/v3/instances/data/contentTypes
- Send an update to the folk in charge of HBP Twitter if we want to shout about it
#. FZJ admin
- https://juser.fz-juelich.de/submit
.. _GH tags: https://github.com/arbor-sim/arbor/tags
.. _AUTOMATED: https://github.com/arbor-sim/arbor/blob/master/.github/workflows/ebrains.yml
.. _cpphardware:
Hardware management
===================
Hardware context
================
Arbor provides two library APIs for working with hardware resources:
......
doc/dev/arbor-dep-graph.png

258 KiB

.. _cable_cell:
How the Cable Cell is made
==========================
This is a short tour through the inner workings of a cable cell intended for new
Arbor developers. Cable cells are not the sole cell type supported in Arbor, but
both the most common and most complex kind.
This is the introduction from which more detailed descriptions will branch out.
As such, we will start with a simple cable cell simulation and how the user input
is turned into an executable object.
Terminology
-----------
In Arbor's codebase some prefixes are used as a low-key namespacing
- ``arb_``:: Mech ABI types, in general use through out Arbor, eg
``arb_mechanism_type``.
- ``fvm_``:: Concerning use by the Finite Volume Method (FVM), eg
``fvm_lowered_cell``.
- ``mc_``:: Related to Multi-Compartment (Cells), identical to cable cells the
difference is purely historical, eg ``mc_cell_group``.
Setting up a Cable Cell simulation
----------------------------------
Arbor constructs a runnable simulation from three ingredients:
- ``recipe``:: a queryable structure that collects cells, connections, gap
junctions, etc.
- ``context``:: a hardware configuration, summarising threads, GPU, and MPI
resources to be used.
- ``domain_decomposition``:: Distribution of cells over the ``context``, made
from ``context`` and ``recipe``.
The interesting part here is the ``recipe``, which is used to lazily produce the
data required by the ``simulation`` and ``domain_decomposition``. A simple example
might be this model of a single cell
.. code:: c++
struct recipe: public arb::recipe {
recipe(const arb::morphology& m, const arb::decor& d): cell{m, d} {}
arb::cell_size_type num_cells() const override { return 1; }
std::vector<arb::probe_info> get_probes(arb::cell_gid_type) const override { return {}; }
arb::cell_kind get_cell_kind(arb::cell_gid_type) const override { return arb::cell_kind::cable; }
std::any get_global_properties(arb::cell_kind) const override { return gprop; }
arb::util::unique_any get_cell_description(arb::cell_gid_type) const override { return cell; }
arb::cable_cell cell
arb::cable_cell_global_properties gprop;
};
As you can see, data is produced on request by feeding the recipe a
``cell_gid_type``. Finally, we need to have a ``morphology`` and a ``decor`` to
generate a ``cable_cell``. Please refer to the documentation on how to construct
these objects. For now, it is sufficient to note that a ``cable_cell`` is a
description of a cell, consisting of a cell geometry and a mapping of
sub-geometries to properties, rather an object that can be simulated. At this point
ion channels are described by a structure ``(name, [(parameter, value)])``.
Lowered Cells, Shared State, and the Discretisation
---------------------------------------------------
To obtain a simulation we need to turn the ``cable_cell`` description object
into a ``fvm_lowered_cell``. However, multiple cells are collected into a
``mc_cell_group`` and ``fvm_lowered_cell`` is the lowered representation of a
full cell group. The ``fvm_lowered_cell`` is responsible for holding the
backend-specific data of a cell group, managing sampling and stimuli, facilitate
event processing, and driving time integration.
Discretisation splits the segments described by the morphology into *control
volumes* (CV; sometimes called *compartments*) according to a ``cv_policy``.
This allows us to construct a system of linear equations, the Hines matrix, to
describe the evolution of the CV voltages according to the cable equation. Refer
to :ref:`Discretisation <discretisation>` and :ref:`Cable equation
<cable_equation>`.
Backend-dependent data is stored in ``shared_state`` as per-compartment data and
indices into these arrays. Upon ``fvm_lowered_cell::initialize`` these are
populated using the concrete discretisation and the ``cable_cell`` description.
Also, mechanisms are concretised using the provided ``mechanism_catalogue`` and
their internal data is set up in ``shared_state``. See :ref:`Shared state <shared_state>`
for more details
Main integration loop
---------------------
Now we are in a state to simulate a cell group by calling
``simulation::run(t_end, dt)``.
The integration in Arbor proceeds in *epochs* with a length less than half a
time constant :math:`T_{min}`, which is the maximum time over which cell groups
can evolve independently. The length :math:`T_{min}` is derived as the minimum over all
axonal/synaptic delays. This works since we know that an event at time :math:`t`
can have an effect at time :math:`t + T_{min}` at the soonest. The factor of one
half stems from double-buffering to overlap communication and computation. So,
Arbor collects all events in an epoch and transmits them in bulk, see
:ref:`Communication <mpi_comm>` for details.
Integration in Arbor is then split into three parts:
1. apply effect of events to mechanisms :ref:`Event Handling <events>`
2. evolve mechanisms and apply currents :ref:`Mechanisms <mechanisms>`
3. solve voltage equations, see :ref:`Solver <matrix_solver>`
Integration proceeds as far as possible without needing to process an event, but
at most with the given time step `dt`.
.. _cell_groups:
Cell groups
===========
Cell groups represent a union of cells of a single *kind* simulated in lockstep.
In a sense, their existence is an optimisation, since parts of the internal
state and computations can be shared between cell in single group. The currently
most complicated cell group is the one for cable cells, called ``mc_cell_group``
(``mc`` stands for multi-compartment, used in older parts of Arbor), so we will
focus on this type here.
Cell groups are created by domain decomposition methods on consideration of soft
(like performance optimisation) and hard (cells connected by gap junctions must
be in the same group) constraints.
Cable Cell group ``mc_cell_group``
----------------------------------
Cable cell groups have backing store in ``shared_state`` (given the
introduction, we now understand that the ``shared`` stands for 'shared' between
cell in a group). In this data set, we also collect the private data of
mechanisms. One thing to watch out for here is that instances of the same
mechanism on a cell group will be collated.
Let us examine an example of such a cell group; assume the following
1. the group comprises two cells with a total of 9 CVs
1. cell 0 has 4 CVs
2. cell 1 has 5 CVs
2. ``pas`` has been painted on three regions and collated
1. region 0 covers CVs [0 1 2]
2. region 1 covers CV 5
3. region 2 covers CV 7
3. ``pas`` has two parameters ``g`` (conductance density) and ``E`` (resting potential)
.. code::
- shared_state
- mechanisms
- id 0
- name "pas"
- width 5
- parameters
- id 0
- name "g"
- values [* * * * *]
- id 1
- name "E"
- values [* * * * *]
- index [0 1 2 5 7]
/ | | | \
/ / / | |
/ / / \ \
- voltage [* * * * * * * * *]
Mechanisms access their view of the cell group data via the
``arb_mechanism_ppack`` structure. To continue with our example, the ``pas``
mechanism would iterate through its view on the cell group voltage to
compute the current density ``i`` like this
.. code::
for ix in 0..width
# Obtain parameters
g = ppack.parameters["g"][ix]
E = ppack.parameters["E"][ix]
# Fetch voltage, note the indirection
cv = ppack.index[ix]
u = ppack.voltage[cv]
# Write outgoing current
ppack.i[ix] = g*(u - E)
In general, cell group wide quantities (like ``voltage`` here) need to be
indexed via ``ppack.index``. Note, that the layout of parameters in ``ppack`` is
this in reality:
.. code::
- ppack
- parameters [g g g g g E E E E E]
When using NMODL, we translate names like ``g`` to offsets into the parameter array
at compile time. Handwritten mechanisms need to do this manually.
......@@ -23,7 +23,7 @@ However, we currently do not handle this case in our build scripts as it is not
.. Note::
When linking an application with **static** Arbor libraries the linker may issue warnings (particularly on macos). Thus, if you encounter problems, try building shared Arbor libraries (cmake option ``-DBUILD_SHARED_LIBS=ON``) instead.
Macro Descripiton
Macro Description
-----------------
.. c:macro:: ARB_LIBNAME_API
......
.. _internals-overview:
.. _dev-overview:
Internals
=========
Developers Guide
================
Here we document internal components of Arbor. These pages can be useful if you're interested in developing on Arbor itself.
.. figure:: arbor-dep-graph.png
:align: center
Arbor dependency graph, generated by ``cmake --graphviz=graphviz/arbor-dep-graph.dot .``
.. toctree::
:caption: Arbor Internals:
:caption: Arbor Developers Guide:
:maxdepth: 2
export
util
cable_cell
cell_groups
matrix_solver
simd_api
shared_state
export
extending_catalogues
mechanism_abi
util
.. numerics
.. _matrix_solver:
Matrix Solvers
==============
Cable Equation
--------------
At the heart of the time evolution step in Arbor we find a linear system that must
be solved once per time step. This system arises from the cable equation
.. math::
C \partial_t V = \frac{\sigma}{2\pi a}\partial_x(a^2\partial_x V) + I
I: \mbox{External currents}
V: \mbox{Membrane potential}
a: \mbox{cable radius}
\sigma: \mbox{conductivity}
after discretisation into CVs, application of the FVM, and choosing an implicit
Euler time step as
.. math::
\left(\frac{\sigma_i C_i}{\Delta\,t} + \sum_{\delta(i, j)} a_{ij}\right)V_i^{k+1} - \sum_{\delta(i, j)} a_ij V_i^{k+1}
= \frac{\sigma_i C_i}{\Delta\,t}V_i^k + \sigma_i I_i
where :math:`\delta(i, j)` indicates whether two CVs are adjacent. It is written
in form of a sparse matrix, symmetric by construction.
The currents :math:`I` originate from the ion channels on the CVs in question,
see the discussion on mechanisms for further details. As :math:`I` potentially
depends on :math:`V`, the cable equation is non-linear. We model these
dependencies up to first order as :math:`I = gV + J` and collect all higher
orders into :math:`J`. This is done to improve accuracy and stability of the
solver. Finding :math:`I` requires the computation of the symbolic derivative
:math:`g = \partial_V I` during compilation of the mechanisms. At runtime
:math:`g` is updated alongside with the currents :math:`I` using that symbolic
expression.
Each *branch* in the morphology leads to a tri-diagonal block in the matrix
describing the system, since *branches* do not contain interior branching
points. Thus, an interior CV couples to only its neighbours (and itself).
However, at branch points, we need to factor in the branch's parents, which
couple blocks via entries outside the tri-diagonal structure. To ensure
un-problematic data dependencies for use of a substitution algorithm, ie each
row depends only on those of larger indices, we enumerate CVs in breadth-first
ordering. This particular form of matrix is called a *Hines matrix*.
CPU
---
.. note:: See ``arbor/backends/multicore/matrix_state.hpp``:
* ``struct matrix_state``
* the ``matrix_state`` constructor sets up the static parts
* the dynamic part is found in ``assemble``
* the solver lives in ``solve``.
The matrix solver proceeds in two phases: assembly and the actual solving. Since
we are working on cell groups, not individual cells, this is executed for each
cell's matrix.
Assembly
^^^^^^^^
We store the matrix in compressed form, as its upper and main diagonals. The
static parts -- foremost the main diagonal -- are computed once at construction
time and stored. The dynamic parts of the matrix and the right-hand side of the
equation are initialised by calling ``assemble``.
Solving
^^^^^^^
The CPU implementation is a straight-forward implemenation of a modified
Thomas-algorithm, using an extra input for the parent relationship. If each
parent is simply the previous CV, we recover the Thomas algorithm.
.. code:: c++
void hines(const arb_value_type* diagonal, // main diagonal
const arb_value_type* upper, // upper diagonal
arb_value_type* rhs, // rhs / solution
const arb_index_type* parents, // CV's parent
int N) {
// backward substitution
for (int i = N-1; i>0; --i) {
const auto parent = parents[i];
const auto factor = upper[parent] / diagonal[i];
diagonal[parent] -= factor * upper[parent];
rhs[parent] -= factor * rhs[i];
}
// solve root
b[0] = b[0] / d[0];
// forward substitution
for(int i=1; i<N; ++i) {
const auto parent = parents[i];
rhs[i] -= upper[i] * rhs[parent];
rhs[i] /= diagonal[i];
}
}
GPU
---
.. note:: See ``arbor/backends/gpu/matrix_fine.hpp``:
* ``struct matrix_state``
* the ``matrix_state`` constructor sets up the static parts
* the dynamic part is found in ``assemble``
* the solver lives in ``solve``.
There is a simple solver in ``arbor/backends/gpu/matrix_flat.hpp``,
which is only used to test/verify the optimised solver described
below.
The GPU implementation of the matrix solver is more complex to improve
performance and make reasonable use of the hardware's capabilities.
In particular it trades a more complex assembly (and structure) for improved
performance.
Looking back at the structure of the Hines matrix, we find that we can solve
blocks in parallel, as long as their parents have been processed. Therefore,
starting at the root, we parallelise over the children of each branching point
and synchronise execution at each such branching point. Each such step is called
a *level*. Execution time is further optimised by packing blocks into threads by
size and splitting overly large blocks to minimise divergence.
A detailled description can be found `here
<https://arxiv.org/ftp/arxiv/papers/1810/1810.12742.pdf>`_ and the references
therein are worthwhile further reading.
File moved
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment