diff --git a/doc/cpp/cable_cell.rst b/doc/cpp/cable_cell.rst
index 6f42ccd2d70f704f345a4e4def4ff7f604b59beb..14a2b412056e4ef823e2162441dca78b2e51b1bc 100644
--- a/doc/cpp/cable_cell.rst
+++ b/doc/cpp/cable_cell.rst
@@ -41,198 +41,6 @@ on the cell via the ``paint`` method; while synapses, current clamps,
gap junction locations, and the site for testing the threshold potential
are specified via the ``place`` method. See :ref:`cppcablecell-dynamics`, below.
-.. _cppcablecell-morphology-construction:
-
-Constructing cell morphologies
-------------------------------
-
-Arbor provides an interface for constructing arbitrary cell morphologies
-composed of *segments*, see :ref:`the general documentation <morph>`.
-**TODO** C++ documentation for segment trees is not yet complete.
-
-The stitch-builder interface
-^^^^^^^^^^^^^^^^^^^^^^^^^^^^
-
-Like the segment tree, the :cpp:type:`stich_builder` class constructs morphologies
-through attaching simple components described by a pair of :cpp:type:`mpoint` values,
-proximal and distal. These components are :cpp:type:`mstitch` objects, and
-they differ from segments in two regards:
-
-1. Stitches are identified by a unique string identifier, in addition to an optional tag value.
-
-2. Stitches can be attached to a parent stitch at either end, or anywhere in the middle.
-
-The ability to attach a stitch some way along another stitch dictates that one
-stitch may correspond to more than one morphological segment once the morphology
-is fully specified. When these attachment points are internal to a stitch, the
-corresponding geometrical point is determined by linearly interpolating between
-the proximal and distal points.
-
-The required header file is ``arbor/morph/stitch.hpp``.
-
-:cpp:type:`mstitch` has two constructors:
-
-.. code::
-
- mstitch::mstitch(std::string id, mpoint prox, mpoint dist, int tag = 0)
- mstitch::mstitch(std::string id, mpoint dist, int tag = 0)
-
-If the proximal point is omitted, it will be inferred from the point at which
-the stitch is attached to its parent.
-
-The :cpp:type:`stitch_builder` class collects the stitches with the ``add`` method:
-
-.. code::
-
- stitch_builder::add(mstitch, const std::string& parent_id, double along = 1.)
- stitch_builder::add(mstitch, double along = 1.)
-
-The first stitch will have no parent. If no parent id is specified for a subsequent
-stitch, the last stitch added will be used as parent. The ``along`` parameter
-must lie between zero and one inclusive, and determines the point of attachment
-as a relative position between the parent's proximal and distal points.
-
-A :cpp:type:`stitched_morphology` is constructed from a :cpp:type:`stitch_builder`,
-and provides both the :cpp:type:`morphology` built from the stitches, and methods
-for querying the extent of individual stitches.
-
-.. cpp:class:: stitched_morphology
-
- .. cpp:function:: stitched_morphology(const stitch_builder&)
- .. cpp:function:: stitched_morphology(stitch_builder&&)
-
- Construct from a ``stitch_builder``. Note that constructing from an
- rvalue is more efficient, as it avoids making a copy of the underlying
- tree structure.
-
- .. cpp:function:: arb::morphology morphology() const
-
- Return the constructed morphology object.
-
- .. cpp:function:: region stitch(const std::string& id) const
-
- Return the region expression corresponding to the specified stitch.
-
- .. cpp:function:: std::vector<msize_t> segments(const std::string& id) const
-
- Return the collection of segments by index comprising the specified stitch.
-
- .. cpp:function:: label_dict labels(const std::string& prefix="") const
-
- Provide a :cpp:type:`label_dict` with a region entry for each stitch; if
- a prefix is provided, this prefix is applied to each segment id to determine
- the region labels.
-
-Example code, constructing a cable cell from a T-shaped morphology specified
-by two stitches:
-
-.. code::
-
- using namespace arb;
-
- mpoint soma0{0, 0, 0, 10};
- mpoint soma1{20, 0, 0, 10};
- mpoint dend_end{10, 100, 0, 1};
-
- stitch_builder builder;
- builder.add({"soma", soma0, soma1, 1});
- builder.add({"dend", dend_end, 4}, "soma", 0.5);
-
- stitched_morphology stitched(std::move(builder));
- cable_cell cell(stitched.morphology(), stitched.labels());
-
- cell.paint("\"soma\"", "hh");
-
-
-Supported morphology formats
-^^^^^^^^^^^^^^^^^^^^^^^^^^^^
-
-Arbor supports morphologies described using the SWC file format and the NeuroML file format.
-
-SWC
-"""
-
-Arbor supports reading morphologies described using the
-`SWC <http://www.neuronland.org/NLMorphologyConverter/MorphologyFormats/SWC/Spec.html>`_ file format. And
-has three different interpretation of that format.
-
-A :cpp:func:`parse_swc()` function is used to parse the SWC file and generate a :cpp:type:`swc_data` object.
-This object contains a vector of :cpp:type:`swc_record` objects that represent the SWC samples, with a number of
-basic checks performed on them. The :cpp:type:`swc_data` object can then be used to generate a
-:cpp:type:`morphology` object using one of the following functions: (See the morphology concepts
-:ref:`page <morph-formats>` for more details).
-
- * :cpp:func:`load_swc_arbor`
- * :cpp:func:`load_swc_allen`
- * :cpp:func:`load_swc_neuron`
-
-.. cpp:class:: swc_record
-
- .. cpp:member:: int id
-
- ID of the record
-
- .. cpp:member:: int tag
-
- Structure identifier (tag).
-
- .. cpp:member:: double x
-
- x coordinate in space.
-
- .. cpp:member:: double y
-
- y coordinate in space.
-
- .. cpp:member:: double z
-
- z coordinate in space.
-
- .. cpp:member:: double r
-
- Sample radius.
-
- .. cpp:member:: int parent_id
-
- Record parent's sample ID.
-
-.. cpp:class:: swc_data
-
- .. cpp:member:: std::string metadata
-
- Contains the comments of an SWC file.
-
- .. cpp:member:: std::vector<swc_record> records
-
- Stored the list of samples from an SWC file, after performing some checks.
-
-.. cpp:function:: swc_data parse_swc(std::istream&)
-
- Returns an :cpp:type:`swc_data` object given an std::istream object.
-
-.. cpp:function:: morphology load_swc_arbor(const swc_data& data)
-
- Returns a :cpp:type:`morphology` constructed according to Arbor's SWC specifications.
-
-.. cpp:function:: morphology load_swc_allen(const swc_data& data, bool no_gaps=false)
-
- Returns a :cpp:type:`morphology` constructed according to the Allen Institute's SWC
- specifications. By default, gaps in the morphology are allowed, this can be toggled
- using the ``no_gaps`` argument.
-
-.. cpp:function:: morphology load_swc_neuron(const swc_data& data)
-
- Returns a :cpp:type:`morphology` constructed according to NEURON's SWC specifications.
-
-.. _cppcablecell-locsets-and-regions:
-
-Identifying sites and subsets of the morphology
------------------------------------------------
-
-.. todo::
-
- TODO: Region and locset documentation is under development.
-
.. _cppcablecell-dynamics:
Cell dynamics
diff --git a/doc/cpp/index.rst b/doc/cpp/index.rst
index b7a7ee0c2a6d3dba4669bcc5b9f0ab0ecebb4e7a..ac5799ff74d60c57825df78d19f10c92e52158a2 100644
--- a/doc/cpp/index.rst
+++ b/doc/cpp/index.rst
@@ -17,6 +17,8 @@ A :cpp:type:`arb::recipe` describes a model, and a :cpp:type:`arb::simulation` i
recipe
cell
+ cable_cell
+ morphology
interconnectivity
hardware
domdec
diff --git a/doc/cpp/morphology.rst b/doc/cpp/morphology.rst
new file mode 100644
index 0000000000000000000000000000000000000000..46c939f72a6ec2bf5be115b27190a1fcd6970e88
--- /dev/null
+++ b/doc/cpp/morphology.rst
@@ -0,0 +1,316 @@
+.. _cppmorphology:
+
+Cable cell morphologies
+=======================
+
+Cell morphologies are required to describe a :ref:`cppcablecell`.
+Morphologies can be constructed directly, or read from a number of
+file formats; see :ref:`cppcablecell-morphology-construction` for details.
+
+Morphology API
+--------------
+
+.. todo::
+
+ TODO: Describe morphology methods.
+
+.. _cppcablecell-morphology-construction:
+
+Constructing cell morphologies
+------------------------------
+
+.. todo::
+
+ TODO: Description of segment trees.
+
+
+The stitch-builder interface
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+Like the segment tree, the :cpp:type:`stich_builder` class constructs morphologies
+through attaching simple components described by a pair of :cpp:type:`mpoint` values,
+proximal and distal. These components are :cpp:type:`mstitch` objects, and
+they differ from segments in two regards:
+
+1. Stitches are identified by a unique string identifier, in addition to an optional tag value.
+
+2. Stitches can be attached to a parent stitch at either end, or anywhere in the middle.
+
+The ability to attach a stitch some way along another stitch dictates that one
+stitch may correspond to more than one morphological segment once the morphology
+is fully specified. When these attachment points are internal to a stitch, the
+corresponding geometrical point is determined by linearly interpolating between
+the proximal and distal points.
+
+The required header file is ``arbor/morph/stitch.hpp``.
+
+:cpp:type:`mstitch` has two constructors:
+
+.. code::
+
+ mstitch::mstitch(std::string id, mpoint prox, mpoint dist, int tag = 0)
+ mstitch::mstitch(std::string id, mpoint dist, int tag = 0)
+
+If the proximal point is omitted, it will be inferred from the point at which
+the stitch is attached to its parent.
+
+The :cpp:type:`stitch_builder` class collects the stitches with the ``add`` method:
+
+.. code::
+
+ stitch_builder::add(mstitch, const std::string& parent_id, double along = 1.)
+ stitch_builder::add(mstitch, double along = 1.)
+
+The first stitch will have no parent. If no parent id is specified for a subsequent
+stitch, the last stitch added will be used as parent. The ``along`` parameter
+must lie between zero and one inclusive, and determines the point of attachment
+as a relative position between the parent's proximal and distal points.
+
+A :cpp:type:`stitched_morphology` is constructed from a :cpp:type:`stitch_builder`,
+and provides both the :cpp:type:`morphology` built from the stitches, and methods
+for querying the extent of individual stitches.
+
+.. cpp:class:: stitched_morphology
+
+ .. cpp:function:: stitched_morphology(const stitch_builder&)
+ .. cpp:function:: stitched_morphology(stitch_builder&&)
+
+ Construct from a ``stitch_builder``. Note that constructing from an
+ rvalue is more efficient, as it avoids making a copy of the underlying
+ tree structure.
+
+ .. cpp:function:: arb::morphology morphology() const
+
+ Return the constructed morphology object.
+
+ .. cpp:function:: region stitch(const std::string& id) const
+
+ Return the region expression corresponding to the specified stitch.
+
+ .. cpp:function:: std::vector<msize_t> segments(const std::string& id) const
+
+ Return the collection of segments by index comprising the specified stitch.
+
+ .. cpp:function:: label_dict labels(const std::string& prefix="") const
+
+ Provide a :cpp:type:`label_dict` with a region entry for each stitch; if
+ a prefix is provided, this prefix is applied to each segment id to determine
+ the region labels.
+
+Example code, constructing a cable cell from a T-shaped morphology specified
+by two stitches:
+
+.. code::
+
+ using namespace arb;
+
+ mpoint soma0{0, 0, 0, 10};
+ mpoint soma1{20, 0, 0, 10};
+ mpoint dend_end{10, 100, 0, 1};
+
+ stitch_builder builder;
+ builder.add({"soma", soma0, soma1, 1});
+ builder.add({"dend", dend_end, 4}, "soma", 0.5);
+
+ stitched_morphology stitched(std::move(builder));
+ cable_cell cell(stitched.morphology(), stitched.labels());
+
+ cell.paint("\"soma\"", "hh");
+
+
+Supported morphology formats
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+Arbor supports morphologies described using the SWC file format and the NeuroML file format.
+
+SWC
+"""
+
+Arbor supports reading morphologies described using the
+`SWC <http://www.neuronland.org/NLMorphologyConverter/MorphologyFormats/SWC/Spec.html>`_ file format. And
+has three different interpretation of that format.
+
+A :cpp:func:`parse_swc()` function is used to parse the SWC file and generate a :cpp:type:`swc_data` object.
+This object contains a vector of :cpp:type:`swc_record` objects that represent the SWC samples, with a number of
+basic checks performed on them. The :cpp:type:`swc_data` object can then be used to generate a
+:cpp:type:`morphology` object using one of the following functions: (See the morphology concepts
+:ref:`page <morph-formats>` for more details).
+
+ * :cpp:func:`load_swc_arbor`
+ * :cpp:func:`load_swc_allen`
+ * :cpp:func:`load_swc_neuron`
+
+.. cpp:class:: swc_record
+
+ .. cpp:member:: int id
+
+ ID of the record
+
+ .. cpp:member:: int tag
+
+ Structure identifier (tag).
+
+ .. cpp:member:: double x
+
+ x coordinate in space.
+
+ .. cpp:member:: double y
+
+ y coordinate in space.
+
+ .. cpp:member:: double z
+
+ z coordinate in space.
+
+ .. cpp:member:: double r
+
+ Sample radius.
+
+ .. cpp:member:: int parent_id
+
+ Record parent's sample ID.
+
+.. cpp:class:: swc_data
+
+ .. cpp:member:: std::string metadata
+
+ Contains the comments of an SWC file.
+
+ .. cpp:member:: std::vector<swc_record> records
+
+ Stored the list of samples from an SWC file, after performing some checks.
+
+.. cpp:function:: swc_data parse_swc(std::istream&)
+
+ Returns an :cpp:type:`swc_data` object given an std::istream object.
+
+.. cpp:function:: morphology load_swc_arbor(const swc_data& data)
+
+ Returns a :cpp:type:`morphology` constructed according to Arbor's SWC specifications.
+
+.. cpp:function:: morphology load_swc_allen(const swc_data& data, bool no_gaps=false)
+
+ Returns a :cpp:type:`morphology` constructed according to the Allen Institute's SWC
+ specifications. By default, gaps in the morphology are allowed, this can be toggled
+ using the ``no_gaps`` argument.
+
+.. cpp:function:: morphology load_swc_neuron(const swc_data& data)
+
+ Returns a :cpp:type:`morphology` constructed according to NEURON's SWC specifications.
+
+.. _locsets-and-regions:
+
+Identifying sites and subsets of the morphology
+-----------------------------------------------
+
+.. todo::
+
+ TODO: Region and locset documentation.
+
+
+Translating regions and locsets to cables and locations
+-------------------------------------------------------
+
+.. todo::
+
+ TODO: ``mprovider``, ``mextent`` and ``thingify``.
+
+
+From morphologies to points and segments
+----------------------------------------
+
+The :cpp:type:`morphology` class has the ``branch_segments`` method for
+returning a vector of :cpp:type:`msegment` objects that describe the geometry
+of that branch. However, determining the position in space of an
+:cpp:type:`mlocation`, for example, requires some assumptions about how to
+position points which fall inside a morphological segment.
+
+The :cpp:type:`place_pwlin` class takes a :cpp:type:`morphology` (and
+optionally an :cpp:type:`isometry`) and interprets it as describing a
+piecewise-linear object in space. It can then be queried to find the 3-d
+positions in space of points on the morphology and the extents in space of
+morphological sub-regions.
+
+Because the morphology need not be contiguous in space, a position query can
+potentially give more than one possible answer. Similarly, a description of a
+cable in terms of segments or partial segments in space may include multiple
+zero-length components as a result of such discontinuities.
+
+.. cpp:class:: place_pwlin
+
+ .. cpp:function:: place_pwlin(const morphology&, const isometry& = isometry())
+
+ Construct a piecewise linear placement of the morphology in space,
+ optionally applying the given isometry.
+
+ .. cpp:function:: mpoint at(mlocation) const
+
+ Return any single point corresponding to the given :cpp:class:`mlocation`
+ in the placement.
+
+ .. cpp:function:: std::vector<mpoint> all_at(mlocation) const
+
+ Return all points corresponding to the given :cpp:class:`mlocation` in
+ the placement.
+
+ .. cpp:function:: std::vector<msegment> segments(const mextent&) const
+
+ Return any minimal collection of segments and partial segments whose
+ union is coterminous with the given :cpp:class:`mextent` in the placement.
+
+ .. cpp:function:: std::vector<msegment> all_segments(const mextent&) const
+
+ Return the maximal set of segments and partial segments whose
+ union is coterminous with the given :cpp:class:`mextent` in the placement.
+
+Isometries
+^^^^^^^^^^
+
+The one cellular morphology may be used to represent multiple cable cells
+which are in principle sited in different locations and orientations.
+An explicit isometry allows the one morphology to be repositioned so as
+to answer location queries on such cells.
+
+An isometry consists of a rotation and a translation. Isometries can be
+composed; as interpreted by Arbor, translations are always regarded as
+being relative to the absolute, extrinsic co-ordinate system, while
+rotations are interpreted as *intrinsic rotations*: rotations are always
+applied with respect to the coordinate system carried with the object,
+not the absolute co-ordinate axes.
+
+.. cpp:class:: isometry
+
+ .. cpp:function:: isometry()
+
+ Construct an identity isometry.
+
+ .. cpp:function:: static isometry translate(double x, double y, double z)
+
+ Construct a translation (x, y, z) with respect to the extrinsic coordinate system.
+
+ .. cpp:function:: template <typename Point> static isometry translate(const Point& p)
+
+ Construct a translation (p.x, p.y, p.z) from an arbitrary object with the corresponding
+ public member variables.
+
+ .. cpp:function:: static isometry rotate(double theta, double x, double y, double z)
+
+ Construct a rotation of theta radians about the axis (x, y, z) with respect to the intrinsic coordinate system.
+
+ .. cpp:function:: template <typename Point> static isometry translate(double theta, const Point& p)
+
+ Construct a rotation of theta radians about the (p.x, p.y, p.z) from an arbitrary object with the corresponding
+ public member variables.
+
+ .. cpp:function:: template <typename Point> Point apply(Point p) const
+
+ The Point object is interpreted as a point in space given by public member variables x, y, and z.
+ The isometry is applied to the point (x, y, z), and a copy of ``p`` is returned with the new
+ coordinate values.
+
+.. cpp:function:: isometry operator*(const isometry& a, const isometry& b)
+
+ Compose two isometries to form a new isometry which applies the intrinsic rotation of *b*, and
+ then the intrinsic rotation of *a*, together with the translations of both *a* and *b*.
+
+
diff --git a/doc/python/morphology.rst b/doc/python/morphology.rst
index 65b83032238f5221996de501b8ed279600fca061..c4e20bc1d381443f25ac6b203f25bf1e54a9e9b5 100644
--- a/doc/python/morphology.rst
+++ b/doc/python/morphology.rst
@@ -394,3 +394,91 @@ Cell morphology
:param str filename: the name of the SWC file.
:param bool no_gaps: enforce that distance between soma center and branches attached to soma is the soma radius.
:rtype: morphology
+
+.. py:class:: place_pwlin
+
+ A :class:`place_pwlin` object allows the querying of the 3-d location of locations and cables
+ in a morphology. Refer to the C++ documentation for :cpp:type:`place_pwlin` for more details.
+
+ .. py:function:: place_pwlin(morphology, isometry)
+ .. py:function:: place_pwlin(morphology)
+ :noindex:
+
+ Construct a piecewise linear placement of the morphology in space,
+ optionally applying the given isometry.
+
+ .. py:method:: at(loc: location) -> location
+
+ Return any single point corresponding to the :class:`location` ``loc``
+ in the placement.
+
+ .. py:method:: all_at(loc: location) -> list[location]
+
+ Return all points corresponding to the given :class:`location` ``loc``
+ the placement.
+
+ .. py:method:: segments(cables: list[cable]) -> list[segment]
+
+ Return any minimal collection of segments and partial segments whose
+ union is coterminous with the sub-region of the morphology covered by
+ the given cables in the placement.
+
+ .. py:method:: all_segments(cables: list[cable]) -> list[segment]
+
+ Return the maximal set of segments and partial segments whose
+ union is coterminous with the sub-region of the morphology covered by
+ the given cables in the placement.
+
+.. py:class:: isometry
+
+ Isometries represent rotations and translations in space, and can be used with
+ :class:`place_pwlin` to position a morphology in an arbitrary spatial location
+ and orientation. Refer to the C++ documentation for :cpp:type:`isometry` for
+ more details.
+
+ .. py::function:: isometry()
+
+ Construct an identity isometry.
+
+ .. py:method:: translate(x: float, y: float, z: float) -> isometry
+ :staticmethod:
+
+ Construct a translation (x, y, z) with respect to the extrinsic coordinate system.
+
+ .. py:method:: translate(displacement: Tuple[float, float, float]) -> isometry
+ :staticmethod:
+ :noindex:
+
+ Construct a translation from the elements of the given tuple.
+
+ .. py:method:: translate(displacement: mpoint) -> isometry
+ :staticmethod:
+ :noindex:
+
+ Construct a translation from the (x, y, z) components of the given :py:class:`mpoint`.
+
+ .. py:method:: rotate(theta: float, x: float, y: float, z: float) -> isometry
+ :staticmethod:
+
+ Construct a rotation of ``theta`` radians about the axis (x, y, z) with respect to the intrinsic coordinate system.
+
+ .. py:method:: rotate(theta: float, axiss: Tuple[float, float, float]) -> isometry
+ :staticmethod:
+ :noindex:
+
+ Construct a rotation of ``theta`` radians about the axis given by the ``axis`` tuple.
+
+ .. py:method:: __call__(point: mpoint) -> mpoint
+
+ Apply the isometry to given point.
+
+ .. py:method:: __call__(point: Tuple[float, float, float, ...]) -> Tuple[float, float, float, ...]
+ :noindex:
+
+ Apply the isometry to the first three components of the given tuple, interpreted as a point.
+
+ .. py:function:: __mul__(a: isometry, b: isometry) -> isometry
+
+ Compose the two isometries to form a new isometry that applies *b* and then applies *a*.
+ Note that rotations are composed as being with respect to the *intrinsic* coordinate system,
+ while translations are always taken to be with respect to the *extrinsic* absolute coordinate system.
diff --git a/doc/requirements.txt b/doc/requirements.txt
index c4e363475fc7c5aa3fadbe0113a1f914ccc12ce0..7eb5a6f21c5e3231c0a2204def1956eeac892c2c 100644
--- a/doc/requirements.txt
+++ b/doc/requirements.txt
@@ -1,2 +1,2 @@
-sphinx>=2.3.0
+sphinx>=2.4.0
svgwrite
diff --git a/python/morphology.cpp b/python/morphology.cpp
index 9b89c0fca69021e8ab69e5899b34bdcbad6d6758..0412b3d0ec74cb7d76c775bf24878dbfbbbb9663 100644
--- a/python/morphology.cpp
+++ b/python/morphology.cpp
@@ -1,11 +1,13 @@
+#include <fstream>
+#include <tuple>
+
#include <pybind11/operators.h>
#include <pybind11/pybind11.h>
#include <pybind11/pytypes.h>
#include <pybind11/stl.h>
-#include <fstream>
-
#include <arbor/morph/morphology.hpp>
+#include <arbor/morph/place_pwlin.hpp>
#include <arbor/morph/primitives.hpp>
#include <arbor/morph/segment_tree.hpp>
@@ -18,6 +20,10 @@ namespace py = pybind11;
namespace pyarb {
+static inline bool cable_lt(const arb::mcable& a, const arb::mcable& b) {
+ return std::tuple(a.branch, a.prox_pos, a.dist_pos)<std::tuple(b.branch, b.prox_pos, b.dist_pos);
+}
+
void check_trailing(std::istream& in, std::string fname) {
if (!(in >> std::ws).eof()) {
throw pyarb_error(util::pprintf("Trailing data found at end of file '{}'", fname));
@@ -60,26 +66,25 @@ void register_morphology(py::module& m) {
// arb::mpoint
py::class_<arb::mpoint> mpoint(m, "mpoint");
mpoint
- .def(py::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(py::init<double, double, double, double>(),
+ "x"_a, "y"_a, "z"_a, "radius"_a,
+ "Create an mpoint object from parameters x, y, z, and radius, specified in µm.")
.def(py::init([](py::tuple t) {
if (py::len(t)!=4) throw std::runtime_error("tuple length != 4");
- return arb::mpoint{t[0].cast<double>(), t[1].cast<double>(), t[2].cast<double>(), t[3].cast<double>()};
- }))
+ return arb::mpoint{t[0].cast<double>(), t[1].cast<double>(), t[2].cast<double>(), t[3].cast<double>()}; }),
+ "Create an mpoint object from a tuple (x, y, z, radius), specified 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(py::self==py::self)
.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);});
+ [](const arb::mpoint& p) {return util::pprintf("{}", p);});
py::implicitly_convertible<py::tuple, arb::mpoint>();
@@ -112,6 +117,82 @@ void register_morphology(py::module& m) {
.def("__str__", [](const arb::mcable& c) { return util::pprintf("{}", c); })
.def("__repr__", [](const arb::mcable& c) { return util::pprintf("{}", c); });
+ // arb::isometry
+ py::class_<arb::isometry> isometry(m, "isometry");
+ isometry
+ .def(py::init<>(), "Construct a trivial isometry.")
+ .def("__call__", [](arb::isometry& iso, arb::mpoint& p) {
+ return iso.apply(p);
+ },
+ "Apply isometry to mpoint argument.")
+ .def("__call__", [](arb::isometry& iso, py::tuple t) {
+ int len = py::len(t);
+ if (len<3) throw std::runtime_error("tuple length < 3");
+
+ arb::mpoint p{t[0].cast<double>(), t[1].cast<double>(), t[2].cast<double>(), 0.};
+ p = iso.apply(p);
+
+ py::tuple result(len);
+ result[0] = p.x;
+ result[1] = p.y;
+ result[2] = p.z;
+ for (int i = 3; i<len; ++i) {
+ result[i] = t[i];
+ }
+ return result;
+ },
+ "Apply isometry to first three components of tuple argument.")
+ .def(py::self*py::self)
+ .def_static("translate",
+ [](double x, double y, double z) { return arb::isometry::translate(x, y, z); },
+ "x"_a, "y"_a, "z"_a,
+ "Construct a translation isometry from displacements x, y, and z.")
+ .def_static("translate",
+ [](py::tuple t) {
+ if (py::len(t)!=3) throw std::runtime_error("tuple length != 3");
+ return arb::isometry::translate(t[0].cast<double>(), t[1].cast<double>(), t[2].cast<double>());
+ },
+ "Construct a translation isometry from the first three components of a tuple.")
+ .def_static("translate",
+ [](arb::mpoint p) { return arb::isometry::translate(p.x, p.y, p.z); },
+ "Construct a translation isometry from the x, y, and z components of an mpoint.")
+ .def_static("rotate",
+ [](double theta, double x, double y, double z) { return arb::isometry::rotate(theta, x, y, z); },
+ "theta"_a, "x"_a, "y"_a, "z"_a,
+ "Construct a rotation isometry of angle theta about the axis in direction (x, y, z).")
+ .def_static("rotate",
+ [](double theta, py::tuple t) {
+ if (py::len(t)!=3) throw std::runtime_error("tuple length != 3");
+ return arb::isometry::rotate(theta, t[0].cast<double>(), t[1].cast<double>(), t[2].cast<double>());
+ },
+ "theta"_a, "axis"_a,
+ "Construct a rotation isometry of angle theta about the given axis in the direction described by a tuple.");
+
+ // arb::place_pwlin
+ py::class_<arb::place_pwlin> place(m, "place_pwlin");
+ place
+ .def(py::init<const arb::morphology&, const arb::isometry&>(),
+ "morphology"_a, "isometry"_a=arb::isometry{},
+ "Construct a piecewise-linear placement object from the given morphology and optional isometry.")
+ .def("at", &arb::place_pwlin::at, "location"_a,
+ "Return an interpolated mpoint corresponding to the location argument.")
+ .def("all_at", &arb::place_pwlin::all_at, "location"_a,
+ "Return list of all possible interpolated mpoints corresponding to the location argument.")
+ .def("segments",
+ [](const arb::place_pwlin& self, std::vector<arb::mcable> cables) {
+ std::sort(cables.begin(), cables.end(), cable_lt);
+ return self.segments(cables);
+ },
+ "Return minimal list of full or partial msegments whose union is coterminous "
+ "with the extent of the given list of cables.")
+ .def("all_segments",
+ [](const arb::place_pwlin& self, std::vector<arb::mcable> cables) {
+ std::sort(cables.begin(), cables.end(), cable_lt);
+ return self.all_segments(cables);
+ },
+ "Return maximal list of non-overlapping full or partial msegments whose union is coterminous "
+ "with the extent of the given list of cables.");
+
//
// Higher-level data structures (segment_tree, morphology)
//
diff --git a/python/test/unit/runner.py b/python/test/unit/runner.py
index 231b88b16ca7c6d041a00ac5252911960907bf35..b727c46647d251114cfef81d85c36d2703cf9cad 100644
--- a/python/test/unit/runner.py
+++ b/python/test/unit/runner.py
@@ -17,6 +17,7 @@ try:
import test_tests
import test_schedules
import test_cable_probes
+ import test_morphology
# add more if needed
except ModuleNotFoundError:
from test import options
@@ -26,6 +27,7 @@ except ModuleNotFoundError:
from test.unit import test_identifiers
from test.unit import test_schedules
from test.unit import test_cable_probes
+ from test.unit import test_morphology
# add more if needed
test_modules = [\
@@ -34,7 +36,8 @@ test_modules = [\
test_event_generators,\
test_identifiers,\
test_schedules,\
- test_cable_probes\
+ test_cable_probes,\
+ test_morphology\
] # add more if needed
def suite():
diff --git a/python/test/unit/test_morphology.py b/python/test/unit/test_morphology.py
new file mode 100644
index 0000000000000000000000000000000000000000..bbc37da98a002a91149138890d6ce332fc4e85ea
--- /dev/null
+++ b/python/test/unit/test_morphology.py
@@ -0,0 +1,145 @@
+# -*- coding: utf-8 -*-
+#
+# test_morphology.py
+
+import unittest
+import arbor as A
+import numpy as N
+import math
+
+# to be able to run .py file from child directory
+import sys, os
+sys.path.insert(0, os.path.abspath(os.path.join(os.path.dirname(__file__), '../../')))
+
+try:
+ import options
+except ModuleNotFoundError:
+ from test import options
+
+"""
+tests for morphology-related classes
+"""
+
+def as_matrix(iso):
+ trans = N.array(iso((0, 0, 0)))
+ return N.c_[N.array([iso(v) for v in [(1,0,0),(0,1,0),(0,0,1)]]).transpose()-N.c_[trans, trans, trans], trans]
+
+class PlacePwlin(unittest.TestCase):
+ def test_identity(self):
+ self.assertTrue(N.isclose(as_matrix(A.isometry()), N.eye(3, 4)).all())
+
+ def test_translation(self):
+ displacement=(4, 5, 6)
+ iso = A.isometry.translate(displacement)
+ expected = N.c_[N.eye(3), displacement]
+ self.assertTrue(N.isclose(as_matrix(iso), expected).all())
+
+ def test_rotation(self):
+ # 90 degrees about y axis.
+ iso = A.isometry.rotate(theta=math.pi/2, axis=(0, 1, 0))
+ expected = N.array([[0, 0, 1, 0],[0, 1, 0, 0],[-1, 0, 0, 0]])
+ self.assertTrue(N.isclose(as_matrix(iso), expected).all())
+
+ def test_compose(self):
+ # Translations are always extrinsic, rotations are intrinsic.
+ y90 = A.isometry.rotate(theta=math.pi/2, axis=(0, 1, 0))
+ z90 = A.isometry.rotate(theta=math.pi/2, axis=(0, 0, 1))
+ t123 = A.isometry.translate(1, 2, 3)
+ t456 = A.isometry.translate(4, 5, 6)
+ iso = t123*z90*t456*y90 # rot about y, then translate, then rot about z, then translate.
+ expected = N.array([[0, 0, 1, 5],[1, 0, 0, 7],[0, 1, 0, 9]])
+ self.assertTrue(N.isclose(as_matrix(iso), expected).all())
+
+ def test_mpoint(self):
+ # Translations can be built from mpoints, and isometry can act upon mpoints. Radius is ignored.
+ y90 = A.isometry.rotate(theta=math.pi/2, axis=(0, 1, 0))
+ z90 = A.isometry.rotate(theta=math.pi/2, axis=(0, 0, 1))
+ t123 = A.isometry.translate(A.mpoint(1, 2, 3, 20))
+ t456 = A.isometry.translate(A.mpoint(4, 5, 6, 30))
+ iso = t123*z90*t456*y90
+ expected = N.array([[0, 0, 1, 5],[1, 0, 0, 7],[0, 1, 0, 9]])
+ self.assertTrue(N.isclose(as_matrix(iso), expected).all())
+
+ q = iso(A.mpoint(2, 3, 4, 10))
+ q_arr = N.array((q.x, q.y, q.z, q.radius))
+ q_expected = N.array([4+5, 2+7, 3+9, 10])
+ self.assertTrue(N.isclose(q_arr, q_expected).all())
+
+ def test_place_pwlin_id(self):
+ # Single branch, discontiguous segments.
+ s0p = A.mpoint(0, 0, 0, 10)
+ s0d = A.mpoint(1, 0, 0, 10)
+ s1p = A.mpoint(3, 0, 0, 10)
+ s1d = A.mpoint(4, 0, 0, 10)
+
+ tree = A.segment_tree()
+ i = tree.append(A.mnpos, s0p, s0d, 1)
+ tree.append(i, s1p, s1d, 2)
+
+ m = A.morphology(tree)
+ place = A.place_pwlin(m)
+
+ L0 = place.at(A.location(0, 0))
+ L0s = place.all_at(A.location(0, 0))
+ self.assertEqual(s0p, L0)
+ self.assertEqual([s0p], L0s)
+
+ Lhalf = place.at(A.location(0, 0.5))
+ Lhalfs = place.all_at(A.location(0, 0.5))
+ self.assertTrue(s0d==Lhalf or s1p==Lhalf)
+ self.assertTrue([s0d, s1p]==Lhalfs)
+
+ Chalf = [(s.prox, s.dist) for s in place.segments([A.cable(0, 0., 0.5)])]
+ self.assertEqual([(s0p, s0d)], Chalf)
+
+ Chalf_all = [(s.prox, s.dist) for s in place.all_segments([A.cable(0, 0., 0.5)])]
+ self.assertEqual([(s0p, s0d), (s1p, s1p)], Chalf_all)
+
+ def test_place_pwlin_isometry(self):
+ # Single branch, discontiguous segments.
+ s0p = A.mpoint(0, 0, 0, 10)
+ s0d = A.mpoint(1, 0, 0, 10)
+ s1p = A.mpoint(3, 0, 0, 10)
+ s1d = A.mpoint(4, 0, 0, 10)
+
+ tree = A.segment_tree()
+ i = tree.append(A.mnpos, s0p, s0d, 1)
+ tree.append(i, s1p, s1d, 2)
+
+ m = A.morphology(tree)
+ iso = A.isometry.translate(2, 3, 4)
+ place = A.place_pwlin(m, iso)
+
+ x0p = iso(s0p)
+ x0d = iso(s0d)
+ x1p = iso(s1p)
+ x1d = iso(s1d)
+
+ L0 = place.at(A.location(0, 0))
+ L0s = place.all_at(A.location(0, 0))
+ self.assertEqual(x0p, L0)
+ self.assertEqual([x0p], L0s)
+
+ Lhalf = place.at(A.location(0, 0.5))
+ Lhalfs = place.all_at(A.location(0, 0.5))
+ self.assertTrue(x0d==Lhalf or x1p==Lhalf)
+ self.assertTrue([x0d, x1p]==Lhalfs)
+
+ Chalf = [(s.prox, s.dist) for s in place.segments([A.cable(0, 0., 0.5)])]
+ self.assertEqual([(x0p, x0d)], Chalf)
+
+ Chalf_all = [(s.prox, s.dist) for s in place.all_segments([A.cable(0, 0., 0.5)])]
+ self.assertEqual([(x0p, x0d), (x1p, x1p)], Chalf_all)
+
+def suite():
+ # specify class and test functions in tuple (here: all tests starting with 'test' from class Contexts
+ suite = unittest.makeSuite(PlacePwlin, ('test'))
+ return suite
+
+def run():
+ v = options.parse_arguments().verbosity
+ runner = unittest.TextTestRunner(verbosity = v)
+ runner.run(suite())
+
+if __name__ == "__main__":
+ run()