# Source code for brian2.spatialneuron.morphology

    """
    Neuronal morphology module.
    This module defines classes to load and build neuronal morphologies.
    """

    import abc
    import numbers
    import os
    from abc import abstractmethod
    from collections import OrderedDict, defaultdict, namedtuple

    from brian2 import numpy as np
    from brian2.units.allunits import meter
    from brian2.units.fundamentalunits import Quantity, check_units, have_same_dimensions
    from brian2.units.stdunits import um
    from brian2.utils.logger import get_logger

    logger = get_logger(__name__)

    __all__ = ["Morphology", "Section", "Cylinder", "Soma"]

    Node = namedtuple("Node", field_names="index,comp_name,x,y,z,diameter,parent,children")

    def _to_meters(value):
        """
        Helper function to convert a floating point value (or array) to a `Quantity`
        in units of "meter", but also allow for ``None`` and return it as it is.
        """
        if value is None:
            return None
        else:
            return Quantity(value, dim=meter.dim)

    def _from_morphology(variable, i, j):
        """
        Helper function to return coordinates from a main morphology (used by
        `SubMorphology`), dealing with ``None``.
        """
        if variable is None:
            return None
        return variable[i:j]

    class MorphologyIndexWrapper:
        """
        A simpler version of `~brian2.groups.group.IndexWrapper`, not allowing for
        string indexing (`Morphology` is not a `Group`). It allows to use
        ``morphology.indices[...]`` instead of ``morphology[...]._indices()``.
        """

        def __init__(self, morphology):
            self.morphology = morphology

        def __getitem__(self, item):
            if isinstance(item, str):
                raise NotImplementedError("Morphologies do not support string indexing")
            assert isinstance(self.morphology, (SubMorphology, Morphology))
            return self.morphology._indices(item)

    def _calc_start_idx(section):
        """
        Calculate the absolute start index that will be used by a flattened
        representation.
        """
        # calculate the absolute start index of this section
        # 1. find the root of the tree
        root = section
        while root._parent is not None:
            root = root._parent
        # 2. go down from the root and advance the indices until we find
        # the current section
        start_idx, found = _find_start_index(root, section)
        assert found
        return start_idx

    def _find_start_index(current, target_section, index=0):
        if current == target_section:
            return index, True
        index += current.n
        for child in current.children:
            if child == target_section:
                return index, True
            else:
                index, found = _find_start_index(child, target_section, index)
                if found:
                    return index, True
        return index, False

    class Topology:
        """
        A representation of the topology of a `Morphology`. Has a useful string
        representation, inspired by NEURON's ``topology`` function.
        """

        def __init__(self, morphology):
            self.morphology = morphology

        def __str__(self):
            # TODO: Make sure that the shown compartments do not get out of hand
            divisor = 1
            return Topology._str_topology(self.morphology, compartments_divisor=divisor)

        @staticmethod
        def _str_topology(
            morphology, indent=0, named_path="", compartments_divisor=1, parent=None
        ):
            """
            A simple string-based representation of a morphology. Inspired by
            NEURON's ``topology`` function.
            """
            description = " " * indent
            length = max([1, morphology.n // compartments_divisor])
            if parent is not None:
                description += "`"
            if isinstance(morphology, Soma):
                description += "( )"
            else:
                description += "-" * length
                description += "|"
            if len(named_path) == 0:
                description += "  [root] \n"
            else:
                description += f"  {named_path}\n"
            for child in morphology.children:
                name = morphology.children.name(child)
                description += Topology._str_topology(
                    child,
                    indent=indent + 2 + length,
                    named_path=f"{named_path}.{name}",
                    compartments_divisor=compartments_divisor,
                    parent=morphology,
                )
            return description

        __repr__ = __str__

    def _rotate(vec, axis, angle):
        """
        Rotate a vector around an arbitrary axis.

        Parameters
        ----------
        vec : `ndarray`
            The vector to rotate.
        axis : `ndarray`
            The axis around which the vector should be rotated.
        angle : float
            The rotation angle (in radians).

        Returns
        -------
        rotated : `ndarray`
            The rotated vector.
        """
        return (
            vec * np.cos(angle)
            - np.cross(axis, vec) * np.sin(angle)
            + axis * np.dot(axis, vec) * (1 - np.cos(angle))
        )

    def _perturb(vec, sigma):
        if sigma == 0:
            return vec
        # Get an arbitrary orthogonal vector
        if vec[1] != 0 or vec[0] != 0:
            orthogonal = np.hstack([vec[1], vec[0], 0])
        else:  # special case for the [0, 0, 1] vector
            orthogonal = np.array([1, 0, 0])

        # Rotate the orthogonal vector
        orthogonal = _rotate(orthogonal, vec, np.random.rand() * np.pi * 2)

        # Use an exponentially distributed angle for the perturbation
        perturbation = np.random.exponential(sigma, 1)
        return _rotate(vec, orthogonal, perturbation)

    def _add_coordinates(
        orig_morphology,
        root=None,
        parent=None,
        name=None,
        section_randomness=0.0,
        compartment_randomness=0.0,
        n_th_child=0,
        total_children=0,
        overwrite_existing=False,
    ):
        # Note that in the following, all values are without physical units

        # The new direction is based on the direction of the parent section
        if parent is None:
            section_dir = np.array([0, 0, 0])
        else:
            section_dir = np.hstack(
                [
                    np.asarray(parent.end_x[-1] - parent.start_x[0]),
                    np.asarray(parent.end_y[-1] - parent.start_y[0]),
                    np.asarray(parent.end_z[-1] - parent.start_z[0]),
                ]
            )
            parent_dir_norm = np.sqrt(np.sum(section_dir**2))
            if parent_dir_norm != 0:
                section_dir /= parent_dir_norm
            else:
                section_dir = np.array([0, 0, 0])
        if not overwrite_existing and orig_morphology.x is not None:
            section = orig_morphology.copy_section()
        elif isinstance(orig_morphology, Soma):
            # No perturbation for the soma
            section = Soma(
                diameter=orig_morphology.diameter,
                x=section_dir[0] * meter,
                y=section_dir[1] * meter,
                z=section_dir[2] * meter,
            )
        else:
            if np.sum(section_dir**2) == 0:
                # We don't have any direction to base this section on (most common
                # case is that the root section is a soma)
                # We stay in the x-y plane and distribute all children in a 360 degree
                # circle around (0, 0, 0)
                section_dir = np.array([1, 0, 0])
                rotation_axis = np.array([0, 0, 1])
                angle_increment = 2 * np.pi / total_children
                rotation_angle = np.pi / 2 + angle_increment * n_th_child
                section_dir = _rotate(section_dir, rotation_axis, rotation_angle)
            else:
                if (
                    section_randomness == 0 and section_dir[2] == 0
                ):  # If we are in the x-y plane, stay there
                    rotation_axis = np.array([0, 0, 1])
                else:
                    rotation_axis = np.array([-section_dir[1], section_dir[2], 0])
                if section_randomness == 0:
                    angle_increment = np.pi / (total_children + 1)
                    rotation_angle = -np.pi / 2 + angle_increment * (n_th_child + 1)
                    section_dir = _rotate(section_dir, rotation_axis, rotation_angle)
            if section_randomness > 0:
                # Rotate randomly
                section_dir = _perturb(section_dir, section_randomness)

            section_dir_norm = np.sqrt(np.sum(section_dir**2))
            section_dir /= section_dir_norm

            # For a soma, we let child sections begin at the surface of the sphere
            if isinstance(parent, Soma):
                origin = parent.diameter / 2 * section_dir
            else:
                origin = (0, 0, 0) * um
            coordinates = np.zeros((orig_morphology.n + 1, 3)) * meter
            start_coords = origin
            coordinates[0, :] = origin
            # Perturb individual compartments as well
            for idx, length in enumerate(orig_morphology.length):
                compartment_dir = _perturb(section_dir, compartment_randomness)
                compartment_dir_norm = np.sqrt(np.sum(compartment_dir**2))
                compartment_dir /= compartment_dir_norm
                current_coords = start_coords + length * compartment_dir
                coordinates[idx + 1, :] = current_coords
                start_coords = current_coords

            if isinstance(orig_morphology, Cylinder) and compartment_randomness == 0:
                section = Cylinder(
                    n=orig_morphology.n,
                    diameter=orig_morphology.diameter[0],
                    x=coordinates[[0, -1], 0],
                    y=coordinates[[0, -1], 1],
                    z=coordinates[[0, -1], 2],
                    type=orig_morphology.type,
                )
            elif isinstance(orig_morphology, Section):
                section = Section(
                    n=orig_morphology.n,
                    diameter=np.hstack(
                        [orig_morphology.start_diameter[0], orig_morphology.end_diameter]
                    )
                    * meter,
                    x=coordinates[:, 0],
                    y=coordinates[:, 1],
                    z=coordinates[:, 2],
                    type=orig_morphology.type,
                )
            else:
                raise NotImplementedError(
                    f"Do not know how to deal with section of type {type(orig_morphology)}."
                )
        if parent is None:
            root = section
        else:
            parent.children.add(name, section)

        for idx, child in enumerate(orig_morphology.children):
            _add_coordinates(
                child,
                root=root,
                parent=section,
                name=orig_morphology.children.name(child),
                n_th_child=idx,
                total_children=len(orig_morphology.children),
                section_randomness=section_randomness,
                compartment_randomness=compartment_randomness,
                overwrite_existing=overwrite_existing,
            )
        return section

    class Children:
        """
        Helper class to represent the children (sub trees) of a section. Can be
        used like a dictionary (mapping names to `Morphology` objects), but iterates
        over the values (sub trees) instead of over the keys (names).
        """

        def __init__(self, owner):
            self._owner = owner
            self._counter = 0
            self._children = []
            self._named_children = {}
            self._given_name = defaultdict(lambda: None)

        def __iter__(self):
            return iter(self._children)

        def __len__(self):
            return len(self._children)

        def __contains__(self, item):
            return item in self._named_children

        def name(self, child):
            """
            Return the given name (i.e. not the automatic name such as ``1``) for a
            child subtree.

            Parameters
            ----------
            child : `Morphology`

            Returns
            -------
            name : str
                The given name for the ``child``.
            """
            return self._given_name[child]

        def __getitem__(self, item):
            if isinstance(item, str):
                return self._named_children[item]
            else:
                raise TypeError("Index has to be an integer or a string.")

        def add(self, name, subtree, automatic_name=False):
            """
            Add a new child to the morphology.

            Parameters
            ----------
            name : str
                The name (e.g. ``"axon"``, ``"soma"``) to use for this sub tree.
            subtree : `Morphology`
                The subtree to link as a child.
            automatic_name : bool, optional
                Whether to chose a new name automatically, if a subtree of the same
                name already exists (uses e.g. ``"dend2"`` instead ``"dend"``).
                Defaults to ``False`` and will raise an error instead.
            """
            if name in self._named_children and self._named_children[name] is not subtree:
                if automatic_name:
                    basename = name
                    counter = 1
                    while name in self._named_children:
                        counter += 1
                        name = basename + str(counter)
                else:
                    raise AttributeError(f"The name {name} is already used for a subtree.")

            if subtree not in self._children:
                self._counter += 1
                self._children.append(subtree)
                self._named_children[str(self._counter)] = subtree
                self._given_name[subtree] = name

            if name is not None:
                self._named_children[name] = subtree

            subtree._parent = self._owner

        def remove(self, name):
            """
            Remove a subtree from this morphology.

            Parameters
            ----------
            name : str
                The name of the sub tree to remove.
            """
            if name not in self:
                raise AttributeError(f"The subtree {name} does not exist")
            subtree = self._named_children[name]
            del self._named_children[name]
            self._children.remove(subtree)
            subtree._parent = None

        def __repr__(self):
            n = len(self._children)
            s = f"<{int(n)} children"
            if n > 0:
                name_dict = {self.name(sec): sec for sec in self._children}
                s += f": {name_dict!r}"
            return f"{s}>"

    class Morphology(metaclass=abc.ABCMeta):
        """
        Neuronal morphology (tree structure).

        The data structure is a tree where each node is an un-branched section
        consisting of a number of connected compartments, each one defined by its
        geometrical properties (length, area, diameter, position).

        Notes
        -----
        You cannot create objects of this class, create a `Soma`, a `Section`, or
        a `Cylinder` instead.
        """

        @check_units(n=1)
        def __init__(self, n, type=None):
            if isinstance(n, str):
                raise TypeError(
                    "Need the number of compartments, not a string. "
                    "If you want to load a morphology from a file, "
                    "use 'Morphology.from_file' instead."
                )
            self._n = int(n)
            if self._n != n:
                raise TypeError("The number of compartments n has to be an integer value.")
            if n <= 0:
                raise ValueError("The number of compartments n has to be at least 1.")
            self.type = type
            self._children = Children(self)
            self._parent = None
            self.indices = MorphologyIndexWrapper(self)

        def __getitem__(self, item):
            """
            Return the subtree with the given name/index.

            Ex.: ```neuron['axon']``` or ```neuron['11213']```
            ```neuron[10*um:20*um]``` returns the subbranch from 10 um to 20 um.
            ```neuron[10*um]``` returns one compartment.
            ```neuron[5]``` returns compartment number 5.
            """
            if isinstance(item, slice):  # neuron[10*um:20*um] or neuron[1:3]
                using_lengths = all(
                    [
                        arg is None or have_same_dimensions(arg, meter)
                        for arg in [item.start, item.stop]
                    ]
                )
                using_ints = all(
                    [
                        arg is None or int(arg) == float(arg)
                        for arg in [item.start, item.stop]
                    ]
                )
                if not (using_lengths or using_ints):
                    raise TypeError("Index slice has to use lengths or integers")

                if using_lengths:
                    if item.step is not None:
                        raise TypeError(
                            "Cannot provide a step argument when slicing with lengths"
                        )
                    pos = np.cumsum(np.asarray(self.length))  # coordinate on the section
                    # We use a special handling for values very close to the points
                    # between the compartments to avoid non-intuitive rounding
                    # effects: a point closer than 1e-12*length of section will be
                    # considered to be within the following section (for a start
                    # index), respectively within the previous section (for an end
                    # index)
                    if item.start is None:
                        i = 0
                    else:
                        diff = np.abs(float(item.start) - pos)
                        if min(diff) < 1e-12 * pos[-1]:
                            i = np.argmin(diff) + 1
                        else:
                            i = np.searchsorted(pos, item.start)
                    if item.stop is None:
                        j = len(pos)
                    else:
                        diff = np.abs(float(item.stop) - pos)
                        if min(diff) < 1e-12 * pos[-1]:
                            j = np.argmin(diff) + 1
                        else:
                            j = np.searchsorted(pos, item.stop) + 1
                else:  # integers
                    i, j, step = item.indices(self.n)
                    if step != 1:
                        raise TypeError("Can only slice a contiguous segment")
            elif isinstance(item, Quantity) and have_same_dimensions(item, meter):
                pos = np.hstack(
                    [0, np.cumsum(np.asarray(self.length))]
                )  # coordinate on the section
                if float(item) < 0 or float(item) > (1 + 1e-12) * pos[-1]:
                    raise IndexError(
                        f"Invalid index {item}, has to be in the interval "
                        f"[{0 * meter!s}, {pos[-1] * meter!s}]."
                    )
                diff = np.abs(float(item) - pos)
                if min(diff) < 1e-12 * pos[-1]:
                    i = np.argmin(diff)
                else:
                    i = np.searchsorted(pos, item) - 1
                j = i + 1
            elif isinstance(item, numbers.Integral):  # int: returns one compartment
                if item < 0:  # allows e.g. to use -1 to get the last compartment
                    item += self.n
                if item >= self.n:
                    raise IndexError(f"Invalid index {item} for {self.n} compartments")
                i = item
                j = i + 1
            elif isinstance(item, str):
                item = str(item)  # convert int to string
                if (len(item) > 1) and all(
                    [c in "LR123456789" for c in item]
                ):  # binary string of the form LLLRLR or 1213 (or mixed)
                    return self._children[item[0]][item[1:]]
                elif item in self._children:
                    return self._children[item]
                else:
                    raise AttributeError(f"The subtree {item} does not exist")
            else:
                raise TypeError(f"Index of type {type(item)} not understood")

            return SubMorphology(self, i, j)

        def __setitem__(self, item, child):
            """
            Inserts the subtree and name it ``item``.
            Ex.: ``neuron['axon']`` or ``neuron['11213']``
            """
            item = str(item)  # convert int to string
            if (len(item) > 1) and all([c in "LR123456789" for c in item]):
                # binary string of the form LLLRLR or 1213 (or mixed)
                self.children[item[0]][item[1:]] = child
            else:
                self.children.add(item, child)

        def __delitem__(self, item):
            """
            Remove the subtree ``item``.
            """
            item = str(item)  # convert int to string
            if (len(item) > 1) and all([c in "LR123456789" for c in item]):
                # binary string of the form LLLRLR or 1213 (or mixed)
                del self._children[item[0]][item[1:]]
            else:
                self._children.remove(item)

        def __getattr__(self, item):
            """
            Return the subtree named ``item``.

            Ex.: ``axon = neuron.axon``
            """
            if item.startswith("_"):
                return super(object, self).__getattr__(item)
            else:
                return self[item]

        def __setattr__(self, item, child):
            """
            Attach a subtree and name it ``item``.

            Ex.: ``neuron.axon = Soma(diameter=10*um)``
            """
            if isinstance(child, Morphology) and not item.startswith("_"):
                self[item] = child
            else:  # If it is not a subtree, then it's a normal class attribute
                object.__setattr__(self, item, child)

        def __delattr__(self, item):
            """
            Remove the subtree ``item``.
            """
            del self[item]

        def _indices(self, item=None, index_var="_idx"):
            """
            Return compartment indices for the main section, relative to the
            original morphology.
            """
            if index_var != "_idx":
                raise AssertionError(f"Unexpected index {index_var}")
            if not (item is None or item == slice(None)):
                if isinstance(item, slice):
                    # So that this always returns an array of values, even if it is
                    # just a single value
                    return self[item]._indices(slice(None))
                else:
                    return self[item]._indices(None)
            else:
                start_idx = _calc_start_idx(self)
                if self.n == 1 and item is None:
                    return start_idx
                else:
                    return np.arange(start_idx, start_idx + self.n)

        def topology(self):
            """
            Return a representation of the topology

            Returns
            -------
            topology : `Topology`
                An object representing the topology (can be converted to a string
                by using ``str(...)`` or simply by printing it with `print`.)
            """
            return Topology(self)

        def generate_coordinates(
            self,
            section_randomness=0.0,
            compartment_randomness=0.0,
            overwrite_existing=False,
        ):
            r"""
            Create a new `Morphology`, with coordinates filled in place where the
            previous morphology did not have any. This is mostly useful for
            plotting a morphology, it does not affect its electrical properties.

            Parameters
            ----------
            section_randomness : float, optional
                The randomness when deciding the direction vector for each new
                section. The given number is the :math:`\beta` parameter of an
                exponential distribution (in degrees) which will be used to
                determine the deviation from the direction of the parent section.
                If the given value equals 0 (the default), then a deterministic
                algorithm will be used instead.
            compartment_randomness : float, optional
                The randomness when deciding the direction vector for each
                compartment within a section. The given number is the :math:`\beta`
                parameter of an exponential distribution (in degrees) which will be
                used to determine the deviation from the main direction of the
                current section. If the given value equals 0 (the default), then all
                compartments will be along a straight line.
            overwrite_existing : bool, optional
                Whether to overwrite existing coordinates in the morphology. This
                is by default set to ``False``, meaning that only sections that do
                not currently have any coordinates set will get new coordinates.
                This allows to conveniently generate a morphology that can be
                plotted for a morphology that is based on points but also has
                artificially added sections (the most common case: an axon added
                to a reconstructed morphology). If set to ``True``, all sections
                will get new coordinates. This can be useful to either get a
                schematic representation of the morphology (with
                ``section_randomness`` and ``compartment_randomness`` both 0) or to
                simply generate a new random variation of a morphology (which will
                still be electrically equivalent, of course).

            Returns
            -------
            morpho_with_coordinates : `Morphology`
                The same morphology, but with coordinates
            """
            # Convert to radians
            section_randomness *= np.pi / 180
            compartment_randomness *= np.pi / 180
            return _add_coordinates(
                self,
                section_randomness=section_randomness,
                compartment_randomness=compartment_randomness,
                overwrite_existing=overwrite_existing,
            )

        @abstractmethod
        def copy_section(self):
            """
            Create a copy of the current section (attributes of this section only,
            not re-creating the parent/children relation)

            Returns
            -------
            copy : `Morphology`
                A copy of this section (without the links to the parent/children)
            """
            raise NotImplementedError()

        @property
        def n(self):
            """
            The number of compartments in this section.
            """
            return self._n

        def __len__(self):
            """
            This is not well-defined, use `Morphology.n` or
            `Morphology.total_compartments` instead.
            """
            raise TypeError(
                "The 'length' of a Morphology is ambiguous, use its "
                "'n' attribute for the number of compartments in this "
                "section or the 'total_compartments' attribute for the "
                "total number of compartments in the whole sub-tree."
            )

        @property
        def total_compartments(self):
            """
            The total number of compartments in this subtree (i.e. the number of
            compartments in this section plus all the compartments in the sections
            deeper in the tree).
            """
            return self.n + sum(c.total_compartments for c in self.children)

        @property
        def total_sections(self):
            """
            The total number of sections in this subtree.
            """
            return 1 + sum(c.total_sections for c in self.children)

        @property
        def parent(self):
            """
            The parent section of this section.
            """
            return self._parent

        @property
        def children(self):
            """
            The children (as a `Children` object) of this section.
            """
            return self._children

        @property
        @abc.abstractmethod
        def end_distance(self):
            """
            The distance to the root of the morphology at the end of this section.
            """
            raise NotImplementedError()

        # Per-compartment attributes
        @property
        @abc.abstractmethod
        def area(self):
            """
            The membrane surface area of each compartment in this section.
            """
            raise NotImplementedError()

        @property
        @abc.abstractmethod
        def volume(self):
            """
            The volume of each compartment in this section.
            """
            raise NotImplementedError()

        @property
        @abc.abstractmethod
        def length(self):
            """
            The length of each compartment in this section.
            """
            raise NotImplementedError()

        @property
        @abc.abstractmethod
        def r_length_1(self):
            """
            The geometry-dependent term to calculate the conductance between the
            start and the midpoint of each compartment. Dividing this value by the
            Intracellular resistivity gives the conductance.
            """
            raise NotImplementedError()

        @property
        @abc.abstractmethod
        def r_length_2(self):
            """
            The geometry-dependent term to calculate the conductance between the
            midpoint and the end of each compartment. Dividing this value by the
            Intracellular resistivity gives the conductance.
            """
            raise NotImplementedError()

        # At-midpoint attributes
        @property
        @abc.abstractmethod
        def diameter(self):
            """
            The diameter at the middle of each compartment in this section.
            """
            raise NotImplementedError()

        @property
        @abc.abstractmethod
        def distance(self):
            """
            The total distance between the midpoint of each compartment and the root
            of the morphology.
            """
            raise NotImplementedError()

        @property
        def start_x(self):
            """
            The x coordinate at the beginning of each compartment. Returns ``None``
            for morphologies without coordinates.
            """
            return _to_meters(self.start_x_)

        @property
        def start_y(self):
            """
            The y coordinate at the beginning of each compartment. Returns ``None``
            for morphologies without coordinates.
            """
            return _to_meters(self.start_y_)

        @property
        def start_z(self):
            """
            The z coordinate at the beginning of each compartment. Returns ``None``
            for morphologies without coordinates.
            """
            return _to_meters(self.start_z_)

        @property
        @abc.abstractmethod
        def start_x_(self):
            """
            The x coordinate (as a unitless floating point number) at the beginning
            of each compartment. Returns ``None`` for morphologies without
            coordinates.
            """
            raise NotImplementedError()

        @property
        @abc.abstractmethod
        def start_y_(self):
            """
            The y coordinate (as a unitless floating point number) at the beginning
            of each compartment. Returns ``None`` for morphologies without
            coordinates.
            """
            raise NotImplementedError()

        @property
        @abc.abstractmethod
        def start_z_(self):
            """
            The z coordinate (as a unitless floating point number) at the beginning
            of each compartment. Returns ``None`` for morphologies without
            coordinates.
            """
            raise NotImplementedError()

        @property
        def x(self):
            """
            The x coordinate at the midpoint of each compartment. Returns ``None``
            for morphologies without coordinates.
            """
            return _to_meters(self.x_)

        @property
        def y(self):
            """
            The y coordinate at the midpoint of each compartment. Returns ``None``
            for morphologies without coordinates.
            """
            return _to_meters(self.y_)

        @property
        def z(self):
            """
            The y coordinate at the midpoint of each compartment. Returns ``None``
            for morphologies without coordinates.
            """
            return _to_meters(self.z_)

        @property
        @abc.abstractmethod
        def x_(self):
            """
            The x coordinate (as a unitless floating point number) at the midpoint
            of each compartment. Returns ``None`` for morphologies without
            coordinates.
            """
            raise NotImplementedError()

        @property
        @abc.abstractmethod
        def y_(self):
            """
            The y coordinate (as a unitless floating point number) at the midpoint
            of each compartment. Returns ``None`` for morphologies without
            coordinates.
            """
            raise NotImplementedError()

        @property
        @abc.abstractmethod
        def z_(self):
            """
            The z coordinate (as a unitless floating point number) at the midpoint
            of each compartment. Returns ``None`` for morphologies without
            coordinates.
            """
            raise NotImplementedError()

        @property
        def end_x(self):
            """
            The x coordinate at the end of each compartment. Returns ``None``
            for morphologies without coordinates.
            """
            return _to_meters(self.end_x_)

        @property
        def end_y(self):
            """
            The y coordinate at the end of each compartment. Returns ``None``
            for morphologies without coordinates.
            """
            return _to_meters(self.end_y_)

        @property
        def end_z(self):
            """
            The z coordinate at the end of each compartment. Returns ``None``
            for morphologies without coordinates.
            """
            return _to_meters(self.end_z_)

        @property
        @abc.abstractmethod
        def end_x_(self):
            """
            The x coordinate (as a unitless floating point number) at the end of
            each compartment. Returns ``None`` for morphologies without coordinates.
            """
            raise NotImplementedError()

        @property
        @abc.abstractmethod
        def end_y_(self):
            """
            The y coordinate (as a unitless floating point number) at the end of
            each compartment. Returns ``None`` for morphologies without coordinates.
            """
            raise NotImplementedError()

        @property
        @abc.abstractmethod
        def end_z_(self):
            """
            The z coordinate (as a unitless floating point number) at the end of
            each compartment. Returns ``None`` for morphologies without coordinates.
            """
            raise NotImplementedError()

        @property
        def coordinates(self):
            r"""
            Array with all coordinates at the start- and end-points of each
            compartment in this section. The array has size :math:`(n+1) \times 3`,
            where :math:`n` is the number of compartments in this section. Each
            row is one point (start point of first compartment, end point of first
            compartment, end point of second compartment, ...), with the columns
            being the x, y, and z coordinates. Returns ``None`` for morphologies
            without coordinates.
            """
            if self.x_ is None:
                return None
            else:
                return Quantity(self.coordinates_, dim=meter.dim)

        @property
        def coordinates_(self):
            r"""
            Array with all coordinates (as unitless floating point numbers) at the
            start- and end-points of each compartment in this section. The array has
            size :math:`(n+1) \times 3`, where :math:`n` is the number of
            compartments in this section. Each row is one point (start point of
            first compartment, end point of first compartment, end point of second
            compartment, ...), with the columns being the x, y, and z coordinates.
            Returns ``None`` for morphologies without coordinates.
            """
            if self.x_ is None:
                return None
            else:
                return np.vstack(
                    [
                        np.hstack([self.start_x_[0], self.end_x_[:]]),
                        np.hstack([self.start_y_[0], self.end_y_[:]]),
                        np.hstack([self.start_z_[0], self.end_z_[:]]),
                    ]
                ).T

        @staticmethod
        def _create_section(compartments, name, parent, sections, spherical_soma):
            if (
                spherical_soma
                and len(compartments) == 1
                and compartments[0].comp_name == "soma"
            ):
                soma = compartments[0]
                section = Soma(
                    diameter=soma.diameter * um, x=soma.x * um, y=soma.y * um, z=soma.z * um
                )
            else:
                sec_x, sec_y, sec_z, sec_diameter = zip(
                    *[(c.x, c.y, c.z, c.diameter) for c in compartments]
                )
                # Add a point for the end of the parent_idx compartment
                if parent is not None:
                    n = len(compartments)
                    if parent.comp_name is not None and parent.comp_name.lower() == "soma":
                        # For a Soma, we don't use its diameter
                        start_diameter = sec_diameter[0]
                    else:
                        start_diameter = parent.diameter
                    # Use relative coordinates
                    sec_x = np.array(sec_x) - parent.x
                    sec_y = np.array(sec_y) - parent.y
                    sec_z = np.array(sec_z) - parent.z
                    start_x = start_y = start_z = 0.0
                else:
                    n = len(compartments) - 1
                    start_diameter = sec_diameter[0]
                    sec_diameter = sec_diameter[1:]
                    start_x = sec_x[0]
                    start_y = sec_y[0]
                    start_z = sec_z[0]
                    sec_x = sec_x[1:]
                    sec_y = sec_y[1:]
                    sec_z = sec_z[1:]

                diameter = np.hstack([start_diameter, sec_diameter]) * um
                x = np.hstack([start_x, sec_x]) * um
                y = np.hstack([start_y, sec_y]) * um
                z = np.hstack([start_z, sec_z]) * um
                section = Section(n=n, diameter=diameter, x=x, y=y, z=z, type=name)

            # Add the section as a child to its parent
            if parent is not None:
                parent_sec = sections[parent.index]
                parent_sec.children.add(name, section, automatic_name=True)

            return section

        @staticmethod
        def _compartments_to_sections(
            compartment, spherical_soma, current_compartments=None, sections=None
        ):
            # Merge all unbranched compartments of the same type into a single
            # section
            if sections is None:
                sections = OrderedDict()
            if current_compartments is None:
                current_compartments = []

            # We have to create a new section, if we are either
            # 1. at a leaf of the tree or at a branching point, or
            # 2. if the compartment type changes
            while True:
                current_compartments.append(compartment)
                if (
                    len(compartment.children) != 1
                    or compartment.children[0].comp_name != compartment.comp_name
                ):
                    break
                compartment = compartment.children[0]

            parent = current_compartments[0].parent
            section = Morphology._create_section(
                current_compartments,
                compartment.comp_name,
                parent=parent,
                sections=sections,
                spherical_soma=spherical_soma,
            )
            sections[current_compartments[-1].index] = section
            # If we are at a branching point, recurse into all subtrees
            for child in compartment.children:
                Morphology._compartments_to_sections(
                    child,
                    spherical_soma=spherical_soma,
                    current_compartments=None,
                    sections=sections,
                )

            return sections

        @staticmethod
        def _replace_three_point_soma(compartment, all_compartments):
            # Replace a three-point/two-cylinder soma by a single spherical soma
            # if possible (see http://neuromorpho.org/SomaFormat.html for some
            # details)

            # We are looking for a node with two children of the soma type (and
            # other childen of other types), where the two children don't have any
            # children of their own
            # To avoid recursion errors we are not searching any further
            soma_children = [c for c in compartment.children if c.comp_name == "soma"]
            if (
                compartment.comp_name == "soma"
                and len(soma_children) == 2
                and all(len(c.children) == 0 for c in soma_children)
            ):
                # This is a 3-point soma to replace
                soma_c = [compartment] + soma_children
                if not all(abs(c.diameter - soma_c[0].diameter) < 1e-15 for c in soma_c):
                    indices = ", ".join(str(c.index) for c in soma_c)
                    raise ValueError(
                        f"Found a '3-point-soma' (lines: {indices}), but not "
                        "all the diameters are "
                        "identical."
                    )
                diameter = soma_c[0].diameter
                point_0 = np.array([soma_c[0].x, soma_c[0].y, soma_c[0].z])
                point_1 = np.array([soma_c[1].x, soma_c[1].y, soma_c[1].z])
                point_2 = np.array([soma_c[2].x, soma_c[2].y, soma_c[2].z])
                length_1 = np.sqrt(np.sum((point_1 - point_0) ** 2))
                length_2 = np.sqrt(np.sum((point_2 - point_0) ** 2))
                if (
                    np.abs(length_1 - diameter / 2) > 0.1
                    or np.abs(length_2 - diameter / 2) > 0.1
                ):
                    raise ValueError(
                        "Cannot replace '3-point-soma' by a single "
                        "point, the second and third points should "
                        "be positioned one radius away from the "
                        f"first point. Distances are {length_1:.3f}um and "
                        f"{length_2:.3f}um, respectively, while the "
                        f"radius is {diameter / 2:.3f}um."
                    )
                children = [c for c in compartment.children if c not in soma_c]
                compartment = Node(
                    index=compartment.index,
                    comp_name="soma",
                    x=point_0[0],
                    y=point_0[1],
                    z=point_0[2],
                    diameter=diameter,
                    parent=compartment.parent,
                    children=children,
                )
                all_compartments[compartment.index] = compartment
                del all_compartments[soma_children[0].index]
                del all_compartments[soma_children[1].index]
            else:
                # Raise a warning if there is more than one soma compartment
                somata = [c for c in all_compartments.values() if c.comp_name == "soma"]
                if len(somata) > 1:
                    logger.warning(
                        f"Found {len(somata)} soma compartments. If you have a 3-point soma, "
                        "make sure it is at the beginning.",
                        name_suffix="soma_compartments",
                    )

        @staticmethod
        def from_points(points, spherical_soma=True):
            """
            Create a morphology from a sequence of points (similar to the ``SWC``
            format, see `Morphology.from_swc_file`). Each point has to be
            a 7-tuple: ``(index, name, x, y, z, diameter, parent)``

            Note that the values should not use units, but are instead all taken
            to be in micrometers.

            Parameters
            ----------
            points : sequence of 7-tuples
                The points of the morphology.
            spherical_soma : bool, optional
                Whether to model a soma as a sphere.
            Returns
            -------
            morphology : `Morphology`

            Notes
            -----
            This format closely follows the SWC format (see
            `Morphology.from_swc_file`) with two differences: the ``type`` should
            be a string (e.g. ``'soma'``) instead of an integer and the 6-th element
            should be the diameter and not the radius.
            """
            # First pass through all points to get the dependency structure
            compartments = OrderedDict()
            for counter, point in enumerate(points):
                if len(point) != 7:
                    raise ValueError(
                        "Each point needs to be described by 7 values, got"
                        f" {len(point)} instead."
                    )
                index, name, x, y, z, diameter, parent_idx = point

                if index in compartments:
                    raise ValueError(f"Two compartments with index {int(index)}")
                if parent_idx == index:
                    raise ValueError(
                        f"Compartment {int(index)} lists itself as the parent compartment."
                    )

                if counter == 0 and parent_idx == -1:
                    parent = None  # The first compartment does not have a parent
                elif parent_idx not in compartments:
                    raise ValueError(
                        f"Did not find the compartment {parent_idx} (parent "
                        f"compartment of compartment {index}). Make sure "
                        "that parent compartments are listed before "
                        "their children."
                    )
                else:
                    parent = compartments[parent_idx]
                children = []
                node = Node(index, name, x, y, z, diameter, parent, children)
                compartments[index] = node
                if parent is not None:
                    parent.children.append(node)

            if spherical_soma:
                Morphology._replace_three_point_soma(
                    list(compartments.values())[0], compartments
                )

            sections = Morphology._compartments_to_sections(
                list(compartments.values())[0], spherical_soma
            )

            # Go through all the sections again and add standard names for all
            # sections (potentially in addition to the name they already have):
            # "L" + "R" for one or two children, "1", "2", "3", etc. otherwise
            children_counter = defaultdict(int)
            for section in sections.values():
                parent = section.parent
                if parent is not None:
                    children_counter[parent] += 1
                    children = parent.children
                    nth_child = children_counter[parent]
                    if len(children) <= 2:
                        name = "L" if nth_child == 1 else "R"
                    else:
                        name = f"{int(nth_child)}"
                    children.add(name, section)

            # There should only be one section without parents
            root = [sec for sec in sections.values() if sec.parent is None]
            assert len(root) == 1
            return root[0]

        @staticmethod
        def from_swc_file(filename, spherical_soma=True):
            """
            Load a morphology from a ``SWC`` file. A large database of morphologies
            in this format can be found at http://neuromorpho.org

            The format consists of an optional header of lines starting with ``#``
            (ignored), followed by a sequence of points, each described in a line
            following the format::

                index type x y z radius parent

            ``index`` is an integer label (starting at 1) that identifies the
            current point and increases by one each line. ``type`` is an integer
            representing the type of the neural segment. The only type that changes
            the interpretation by Brian is the type ``1`` which signals a soma.
            Types ``2`` (axon), ``3`` (dendrite), and ``4`` (apical dendrite) are
            used to give corresponding names to the respective sections. All other
            types are ignored. ``x``, ``y``, and ``z`` are the cartesian coordinates
            at each point and ``r`` is its radius. ``parent`` refers to the index
            of the parent point or is ``-1`` for the root point.

            Parameters
            ----------
            filename : str
                The name of the ``SWC`` file.
            spherical_soma : bool, optional
                Whether to model the soma as a sphere.

            Returns
            -------
            morpho : `Morphology`
                The morphology stored in the given file.
            """
            swc_types = defaultdict(lambda: None)
            # The following names will be translated into names, all other will be
            # ignored
            swc_types.update({"1": "soma", "2": "axon", "3": "dend", "4": "apic"})

            with open(filename) as f:
                points = []
                for line_no, line in enumerate(f):
                    line = line.strip()
                    if line.startswith("#") or len(line) == 0:
                        # Ignore comments or empty lines
                        continue
                    splitted = line.split()
                    if len(splitted) != 7:
                        raise ValueError(
                            "Each line of an SWC file has to contain "
                            "7 space-separated entries, but line "
                            f"{line_no + 1} contains {len(splitted)}."
                        )
                    index, comp_type, x, y, z, radius, parent = splitted
                    points.append(
                        (
                            int(index),
                            swc_types[comp_type],
                            float(x),
                            float(y),
                            float(z),
                            2 * float(radius),
                            int(parent),
                        )
                    )

            return Morphology.from_points(points, spherical_soma=spherical_soma)

        @staticmethod
        def from_file(filename, spherical_soma=True):
            """
            Convencience method to load a morphology from a given file. At the
            moment, only ``SWC`` files are supported, calling this function is
            therefore equivalent to calling `Morphology.from_swc_file` directly.

            Parameters
            ----------
            filename : str
                The name of a file storing a morphology.
            spherical_soma : bool, optional
                Whether to model the soma as a sphere.
            Returns
            -------
            morphology : `Morphology`
                The morphology stored in the given file.
            """
            _, ext = os.path.splitext(filename)
            if ext.lower() == ".swc":
                return Morphology.from_swc_file(filename, spherical_soma=spherical_soma)
            else:
                raise NotImplementedError(
                    "Currently, SWC is the only supported file format."
                )

    class SubMorphology:
        """
        A view on a subset of a section in a morphology.
        """

        def __init__(self, morphology, i, j):
            self._morphology = morphology
            self.indices = MorphologyIndexWrapper(self)
            self._i = i
            self._j = j

        def _indices(self, item=None):
            if not (item is None or item == slice(None)):
                raise IndexError("Cannot index a view on a subset of a section further")
            # Start index of the main section
            start_idx = _calc_start_idx(self._morphology)
            if item is None and self.n == 1:
                return start_idx + self._i
            else:
                return np.arange(start_idx + self._i, start_idx + self._j)

        @property
        def n(self):
            """
            The number of compartments in this sub-section.
            """
            return self._j - self._i

        def __len__(self):
            return self.n

        @property
        def n_sections(self):
            """
            The number of sections in this sub-section (always 1).
            """
            return 1

        # Per-compartment attributes
        @property
        def area(self):
            """
            The membrane surface area of each compartment in this sub-section.
            """
            return self._morphology.area[self._i : self._j]

        @property
        def volume(self):
            """
            The volume of each compartment in this sub-section.
            """
            return self._morphology.volume[self._i : self._j]

        @property
        def length(self):
            """
            The length of each compartment in this sub-section.
            """
            return self._morphology.length[self._i : self._j]

        @property
        def r_length_1(self):
            """
            The geometry-dependent term to calculate the conductance between the
            start and the midpoint of each compartment in this sub-section.
            Dividing this value by the Intracellular resistivity gives the
            conductance.
            """
            return self._morphology.r_length_1[self._i : self._j]

        @property
        def r_length_2(self):
            """
            The geometry-dependent term to calculate the conductance between the
            midpoint and the end of each compartment in this sub-section. Dividing
            this value by the Intracellular resistivity gives the conductance.
            """
            return self._morphology.r_length_2[self._i : self._j]

        # At-midpoint attributes
        @property
        def diameter(self):
            """
            The diameter at the middle of each compartment in this sub-section.
            """
            return self._morphology.diameter[self._i : self._j]

        @property
        def distance(self):
            """
            The total distance between the midpoint of each compartment in this
            sub-section and the root of the morphology.
            """
            return self._morphology.distance[self._i : self._j]

        @property
        def start_x(self):
            """
            The x coordinate at the beginning of each compartment in this
            sub-section. Returns ``None`` for morphologies without coordinates.
            """
            return _to_meters(self.start_x_)

        @property
        def start_y(self):
            """
            The y coordinate at the beginning of each compartment in this
            sub-section. Returns ``None`` for morphologies without coordinates.
            """
            return _to_meters(self.start_y_)

        @property
        def start_z(self):
            """
            The x coordinate at the beginning of each compartment in this
            sub-section. Returns ``None`` for morphologies without coordinates.
            """
            return _to_meters(self.start_z_)

        @property
        def start_x_(self):
            """
            The x coordinate (as a unitless floating point number) at the beginning
            of each compartment in this sub-section. Returns ``None`` for
            morphologies without coordinates.
            """
            return _from_morphology(self._morphology.start_x_, self._i, self._j)

        @property
        def start_y_(self):
            """
            The y coordinate (as a unitless floating point number) at the beginning
            of each compartment in this sub-section. Returns ``None`` for
            morphologies without coordinates.
            """
            return _from_morphology(self._morphology.start_y_, self._i, self._j)

        @property
        def start_z_(self):
            """
            The z coordinate (as a unitless floating point number) at the beginning
            of each compartment in this sub-section. Returns ``None`` for
            morphologies without coordinates.
            """
            return _from_morphology(self._morphology.start_z_, self._i, self._j)

        @property
        def x(self):
            """
            The x coordinate at the midpoint of each compartment in this
            sub-section. Returns ``None`` for morphologies without coordinates.
            """
            return _to_meters(self.x_)

        @property
        def y(self):
            """
            The y coordinate at the midpoint of each compartment in this
            sub-section. Returns ``None`` for morphologies without coordinates.
            """
            return _to_meters(self.y_)

        @property
        def z(self):
            """
            The z coordinate at the midpoint of each compartment in this
            sub-section. Returns ``None`` for morphologies without coordinates.
            """
            return _to_meters(self.z_)

        @property
        def x_(self):
            """
            The x coordinate (as a unitless floating point number) at the midpoint
            of each compartment in this sub-section. Returns ``None`` for
            morphologies without coordinates.
            """
            return _from_morphology(self._morphology.x_, self._i, self._j)

        @property
        def y_(self):
            """
            The y coordinate (as a unitless floating point number) at the midpoint
            of each compartment in this sub-section. Returns ``None`` for
            morphologies without coordinates.
            """
            return _from_morphology(self._morphology.y_, self._i, self._j)

        @property
        def z_(self):
            """
            The z coordinate (as a unitless floating point number) at the midpoint
            of each compartment in this sub-section. Returns ``None`` for
            morphologies without coordinates.
            """
            return _from_morphology(self._morphology.z_, self._i, self._j)

        @property
        def end_x(self):
            """
            The x coordinate at the end of each compartment in this sub-section.
            Returns ``None`` for morphologies without coordinates.
            """
            return _to_meters(self.end_x_)

        @property
        def end_y(self):
            """
            The y coordinate at the end of each compartment in this sub-section.
            Returns ``None`` for morphologies without coordinates.
            """
            return _to_meters(self.end_y_)

        @property
        def end_z(self):
            """
            The z coordinate at the end of each compartment in this sub-section.
            Returns ``None`` for morphologies without coordinates.
            """
            return _to_meters(self.end_z_)

        @property
        def end_x_(self):
            """
            The x coordinate (as a unitless floating point number) at the end of
            each compartment in this sub-section. Returns ``None`` for morphologies
            without coordinates.
            """
            return _from_morphology(self._morphology.end_x_, self._i, self._j)

        @property
        def end_y_(self):
            """
            The y coordinate (as a unitless floating point number) at the end of
            each compartment in this sub-section. Returns ``None`` for morphologies
            without coordinates.
            """
            return _from_morphology(self._morphology.end_y_, self._i, self._j)

        @property
        def end_z_(self):
            """
            The z coordinate (as a unitless floating point number) at the end of
            each compartment in this sub-section. Returns ``None`` for morphologies
            without coordinates.
            """
            return _from_morphology(self._morphology.end_z_, self._i, self._j)

    class Soma(Morphology):
        """
        A spherical, iso-potential soma.

        Parameters
        ----------
        diameter : `Quantity`
            Diameter of the sphere.
        x : `Quantity`, optional
            The x coordinate of the position of the soma.
        y : `Quantity`, optional
            The y coordinate of the position of the soma.
        z : `Quantity`, optional
            The z coordinate of the position of the soma.
        type : str, optional
            The ``type`` of this section, defaults to ``'soma'``.
        """

        @check_units(diameter=meter, x=meter, y=meter, z=meter)
        def __init__(self, diameter, x=None, y=None, z=None, type="soma"):
            Morphology.__init__(self, n=1, type=type)
            if diameter.shape != () and len(diameter) != 1:
                raise TypeError("Diameter has to be a scalar value.")
            for coord in [x, y, z]:
                if coord is not None and coord.shape != () and len(coord) != 1:
                    raise TypeError("Coordinates have to be scalar values.")
            self._diameter = np.ones(1) * diameter
            if any(coord is not None for coord in (x, y, z)):
                default_value = np.array([0.0])
            else:
                default_value = None
            self._x = np.atleast_1d(np.asarray(x)) if x is not None else default_value
            self._y = np.atleast_1d(np.asarray(y)) if y is not None else default_value
            self._z = np.atleast_1d(np.asarray(z)) if z is not None else default_value

        def __repr__(self):
            s = f"{self.__class__.__name__}(diameter={self.diameter[0]!r}"
            if self._x is not None:
                s += f", x={self.x[0]!r}, y={self.y[0]!r}, z={self.z[0]!r}"
            if self.type != "soma":
                s += f", type={self.type!r}"
            return f"{s})"

        def copy_section(self):
            return Soma(self.diameter, x=self.x, y=self.y, z=self.z, type=self.type)

        # Note that the per-compartment properties should always return 1D arrays,
        # i.e. for the soma arrays of length 1 instead of scalar values
        @property
        def area(self):
            """
            The membrane surface area of this section (as an array of length 1).
            """
            return np.pi * self.diameter**2

        @property
        def volume(self):
            """
            The volume of this section (as an array of length 1).
            """
            return (np.pi * self.diameter**3) / 6

        @property
        def length(self):
            """
            The "length" (equal to `diameter`) of this section (as an array of
            length 1).
            """
            return self.diameter

        @property
        def r_length_1(self):
            """
            The geometry-dependent term to calculate the conductance between the
            start and the midpoint of each compartment. Returns a fixed (high)
            value for a `Soma`, corresponding to a section with very low
            intracellular resistance.
            """
            return [1] * meter

        @property
        def r_length_2(self):
            """
            The geometry-dependent term to calculate the conductance between the
            midpoint and the end of each compartment. Returns a fixed (high)
            value for a `Soma`, corresponding to a section with very low
            intracellular resistance.
            """
            return [1] * meter

        @property
        def diameter(self):
            """
            The diameter of this section (as an array of length 1).
            """
            return self._diameter

        @property
        def distance(self):
            """
            The total distance between the midpoint of this section and the root
            of the morphology. The `Soma` is most likely the root of the
            morphology, and therefore the `distance` is 0.
            """
            dist = self._parent.distance[-1:] if self._parent is not None else [0] * um
            return dist

        @property
        def start_x_(self):
            """
            The x-coordinate of the current section (as an array of length 1). Note
            that a `Soma` is modelled as a "point" with finite surface/volume,
            equivalent to that of a sphere with the given `diameter`. It's start-,
            midpoint-, and end-coordinates are therefore identical.
            """
            return self._x

        @property
        def start_y_(self):
            """
            The y-coordinate of the current section (as an array of length 1). Note
            that a `Soma` is modelled as a "point" with finite surface/volume,
            equivalent to that of a sphere with the given `diameter`. It's start-,
            midpoint-, and end-coordinates are therefore identical.
            """
            return self._y

        @property
        def start_z_(self):
            """
            The z-coordinate of the current section (as an array of length 1). Note
            that a `Soma` is modelled as a "point" with finite surface/volume,
            equivalent to that of a sphere with the given `diameter`. It's start-,
            midpoint-, and end-coordinates are therefore identical.
            """
            return self._z

        @property
        def x_(self):
            """
            The x-coordinate of the current section (as an array of length 1). Note
            that a `Soma` is modelled as a "point" with finite surface/volume,
            equivalent to that of a sphere with the given `diameter`. It's start-,
            midpoint-, and end-coordinates are therefore identical.
            """
            return self._x

        @property
        def y_(self):
            """
            The y-coordinate of the current section (as an array of length 1). Note
            that a `Soma` is modelled as a "point" with finite surface/volume,
            equivalent to that of a sphere with the given `diameter`. It's start-,
            midpoint-, and end-coordinates are therefore identical.
            """
            return self._y

        @property
        def z_(self):
            """
            The z-coordinate of the current section (as an array of length 1). Note
            that a `Soma` is modelled as a "point" with finite surface/volume,
            equivalent to that of a sphere with the given `diameter`. It's start-,
            midpoint-, and end-coordinates are therefore identical.
            """
            return self._z

        @property
        def end_x_(self):
            """
            The x-coordinate of the current section (as an array of length 1). Note
            that a `Soma` is modelled as a "point" with finite surface/volume,
            equivalent to that of a sphere with the given `diameter`. It's start-,
            midpoint-, and end-coordinates are therefore identical.
            """
            return self._x

        @property
        def end_y_(self):
            """
            The y-coordinate of the current section (as an array of length 1). Note
            that a `Soma` is modelled as a "point" with finite surface/volume,
            equivalent to that of a sphere with the given `diameter`. It's start-,
            midpoint-, and end-coordinates are therefore identical.
            """
            return self._y

        @property
        def end_z_(self):
            """
            The z-coordinate of the current section (as an array of length 1). Note
            that a `Soma` is modelled as a "point" with finite surface/volume,
            equivalent to that of a sphere with the given `diameter`. It's start-,
            midpoint-, and end-coordinates are therefore identical.
            """
            return self._z

        @property
        def end_distance(self):
            """
            The distance to the root of the morphology at the end of this section.
            Note that since a `Soma` is modeled as a point (see docs of `x`, etc.),
            it does not add anything to the total distance, e.g. a section
            connecting to a `Soma` has a `distance` of 0 um at its start.
            """
            dist = self._parent.end_distance if self._parent is not None else 0 * um
            return dist

    class Section(Morphology):
        """
        A section (unbranched structure), described as a sequence of truncated
        cones with potentially varying diameters and lengths per compartment.

        Parameters
        ----------
        diameter : `Quantity`
            Either a single value (the constant diameter along the whole section),
            or a value of length ``n+1``. When ``n+1`` values are given, they
            will be interpreted as the diameters at the start of the first
            compartment and the diameters at the end of each compartment (which is
            equivalent to: the diameter at the start of each compartment and the
            diameter at the end of the last compartment.
        n : int, optional
            The number of compartments in this section. Defaults to 1.
        length : `Quantity`, optional
            Either a single value (the total length of the section), or a value of
            length ``n``, the length of each individual compartment. Cannot be
            combined with the specification of coordinates.
        x : `Quantity`, optional
            ``n+1`` values, specifying the x coordinates of the start point of the
            first compartment and the end-points of all compartments (which is
            equivalent to: the start point of all compartments and the end point of
            the last compartment). The coordinates are interpreted as relative to
            the end point of the parent compartment (if any), so in most cases the
            start point should be ``0*um``. The common exception is a cylinder
            connecting to a `Soma`, here the start point can be used to make the
            cylinder start at the surface of the sphere instead of at its center.
            You can specify all of ``x``, ``y``, or ``z`` to specify
            a morphology in 3D, or only one or two out of them to specify a
            morphology in 1D or 2D.
        y : `Quantity`, optional
            See ``x``
        z : `Quantity`, optional
            See ``x``
        type : str, optional
            The type (e.g. ``"axon"``) of this `Section`.
        """

        @check_units(
            n=1,
            length=meter,
            diameter=meter,
            start_diameter=meter,
            x=meter,
            y=meter,
            z=meter,
        )
        def __init__(
            self,
            diameter,
            n=1,
            length=None,
            x=None,
            y=None,
            z=None,
            start_diameter=None,
            origin=None,
            type=None,
        ):
            n = int(n)
            Morphology.__init__(self, n=n, type=type)

            if diameter.ndim != 1 or len(diameter) != n + 1:
                raise TypeError(
                    "The diameter argument has to be a one-dimensional array of length "
                    f"{int(n + 1)}"
                )
            self._diameter = Quantity(diameter, copy=True).reshape((n + 1,))

            if (x is not None or y is not None or z is not None) and length is not None:
                raise TypeError("Cannot specify coordinates and length at the same time.")

            if length is not None:
                # Length
                if length.ndim != 1 or len(length) != n:
                    raise TypeError(
                        "The length argument has to be a one-dimensional array of length "
                        f"{int(n)}"
                    )
                self._length = Quantity(length, copy=True).reshape((n,))
                self._x = self._y = self._z = None
            else:
                # Coordinates
                if x is None and y is None and z is None:
                    raise TypeError(
                        "No length specified, need to specify at least "
                        "one out of x, y, or z."
                    )
                for name, value in [("x", x), ("y", y), ("z", z)]:
                    if value is not None and (value.ndim != 1 or len(value) != n + 1):
                        raise TypeError(
                            f"'{name}' needs to be a 1-dimensional array of length {n + 1}."
                        )
                self._x = (
                    np.asarray(x).reshape((n + 1,)) if x is not None else np.zeros(n + 1)
                )
                self._y = (
                    np.asarray(y).reshape((n + 1,)) if y is not None else np.zeros(n + 1)
                )
                self._z = (
                    np.asarray(z).reshape((n + 1,)) if z is not None else np.zeros(n + 1)
                )

                length = np.sqrt(
                    (self.end_x - self.start_x) ** 2
                    + (self.end_y - self.start_y) ** 2
                    + (self.end_z - self.start_z) ** 2
                )
                self._length = length

        def __repr__(self):
            if all(
                np.abs(self.end_diameter - self.end_diameter[0])
                < self.end_diameter[0] * 1e-12
            ):
                # Constant diameter
                diam = self.end_diameter[0]
            else:
                diam = (
                    np.hstack(
                        [np.asarray(self.start_diameter[0]), np.asarray(self.end_diameter)]
                    )
                    * meter
                )
            s = f"{self.__class__.__name__}(diameter={diam!r}"
            if self.n != 1:
                s += f", n={self.n}"
            if self._x is not None:
                s += f", x={self._x!r}, y={self._y!r}, z={self._z!r}"
            else:
                s += f", length={sum(self._length)!r}"
            if self.type is not None:
                s += f", type={self.type!r}"
            return f"{s})"

        def copy_section(self):
            if self.x is None:
                x, y, z = None, None, None
                length = self.length
            else:
                x, y, z = self._x * meter, self._y * meter, self._z * meter
                length = None
            return Section(
                diameter=self._diameter,
                n=self.n,
                x=x,
                y=y,
                z=z,
                length=length,
                type=self.type,
            )

        @property
        def area(self):
            r"""
            The membrane surface area of each compartment in this section. The
            surface area of each compartment is calculated as
            :math:`\frac{\pi}{2}(d_1 + d_2)\sqrt{\frac{(d_1 - d_2)^2}{4} + l^2)}`,
            where :math:`l` is the length of the compartment, and :math:`d_1` and
            :math:`d_2` are the diameter at the start and end of the compartment,
            respectively. Note that this surface area does not contain the area of
            the two disks at the two sides of the truncated cone.
            """
            d_1 = self.start_diameter
            d_2 = self.end_diameter
            return (
                np.pi / 2 * (d_1 + d_2) * np.sqrt(((d_1 - d_2) ** 2) / 4 + self._length**2)
            )

        @property
        def volume(self):
            r"""
            The volume of each compartment in this section. The volume of each
            compartment is calculated as
            :math:`\frac{\pi}{12} l (d_1^2 + d_1 d_2 + d_2^2)`,
            where :math:`l` is the length of the compartment, and :math:`d_1` and
            :math:`d_2` are the diameter at the start and end of the compartment,
            respectively.
            """
            d_1 = self.start_diameter
            d_2 = self.end_diameter
            return np.pi * self._length * (d_1**2 + d_1 * d_2 + d_2**2) / 12

        @property
        def length(self):
            """
            The length of each compartment in this section.
            """
            return self._length

        @property
        def start_diameter(self):
            """
            The diameter at the start of each compartment in this section.
            """
            return Quantity(self._diameter[:-1], copy=True)

        @property
        def end_diameter(self):
            """
            The diameter at the end of each compartment in this section.
            """
            return Quantity(self._diameter[1:], copy=True)

        @property
        def diameter(self):
            """
            The diameter at the middle of each compartment in this section.
            """
            d_1 = self.start_diameter
            d_2 = self.end_diameter
            # Diameter at the center
            return 0.5 * (d_1 + d_2)

        @property
        def distance(self):
            """
            The total distance between the midpoint of each compartment and the root
            of the morphology.
            """
            dist = self._parent.end_distance if self._parent is not None else 0 * um
            return dist + np.cumsum(self.length) - 0.5 * self.length

        @property
        def end_distance(self):
            """
            The distance to the root of the morphology at the end of this section.
            """
            return self.distance[-1] + 0.5 * self.length[-1]

        @property
        def r_length_1(self):
            """
            The geometry-dependent term to calculate the conductance between the
            start and the midpoint of each compartment. Dividing this value by the
            Intracellular resistivity gives the conductance.
            """
            d_1 = self.start_diameter
            d_2 = (self.start_diameter + self.end_diameter) * 0.5
            return np.pi / 2 * (d_1 * d_2) / self._length

        @property
        def r_length_2(self):
            """
            The geometry-dependent term to calculate the conductance between the
            midpoint and the end of each compartment. Dividing this value by the
            Intracellular resistivity gives the conductance.
            """
            d_1 = (self.start_diameter + self.end_diameter) * 0.5
            d_2 = self.end_diameter
            return np.pi / 2 * (d_1 * d_2) / self._length

        @property
        def start_x_(self):
            """
            The x coordinate (as a unitless floating point number) at the beginning
            of each compartment. Returns ``None`` for morphologies without
            coordinates.
            """
            if self._x is None:
                return None
            parent_end_x = self.parent.end_x_ if self.parent is not None else None
            if parent_end_x is not None:
                return parent_end_x[-1] + self._x[:-1]
            else:
                return self._x[:-1]

        @property
        def start_y_(self):
            """
            The y coordinate (as a unitless floating point number) at the beginning
            of each compartment. Returns ``None`` for morphologies without
            coordinates.
            """
            if self._y is None:
                return None
            parent_end_y = self.parent.end_y_ if self.parent is not None else None
            if parent_end_y is not None:
                return parent_end_y[-1] + self._y[:-1]
            else:
                return self._y[:-1]

        @property
        def start_z_(self):
            """
            The z coordinate (as a unitless floating point number) at the beginning
            of each compartment. Returns ``None`` for morphologies without
            coordinates.
            """
            if self._z is None:
                return None
            parent_end_z = self.parent.end_z_ if self.parent is not None else None
            if parent_end_z is not None:
                return parent_end_z[-1] + self._z[:-1]
            else:
                return self._z[:-1]

        @property
        def x_(self):
            """
            The x coordinate (as a unitless floating point number) at the midpoint
            of each compartment. Returns ``None`` for morphologies without
            coordinates.
            """
            if self._x is None:
                return None
            start_x = self.start_x_
            diff_x = self.end_x_ - start_x
            return start_x + 0.5 * diff_x

        @property
        def y_(self):
            """
            The y coordinate (as a unitless floating point number) at the midpoint
            of each compartment. Returns ``None`` for morphologies without
            coordinates.
            """
            if self._y is None:
                return None
            start_y = self.start_y_
            diff_y = self.end_y_ - start_y
            return start_y + 0.5 * diff_y

        @property
        def z_(self):
            """
            The z coordinate (as a unitless floating point number) at the midpoint
            of each compartment. Returns ``None`` for morphologies without
            coordinates.
            """
            if self._z is None:
                return None
            start_z = self.start_z_
            diff_z = self.end_z_ - start_z
            return start_z + 0.5 * diff_z

        @property
        def end_x_(self):
            """
            The x coordinate (as a unitless floating point number) at the end of
            each compartment. Returns ``None`` for morphologies without coordinates.
            """
            if self._x is None:
                return None
            parent_end_x = self.parent.end_x_ if self.parent is not None else None
            if parent_end_x is not None:
                return parent_end_x[-1] + self._x[1:]
            else:
                return self._x[1:]

        @property
        def end_y_(self):
            """
            The y coordinate (as a unitless floating point number) at the end of
            each compartment. Returns ``None`` for morphologies without coordinates.
            """
            if self._y is None:
                return None
            parent_end_y = self.parent.end_y_ if self.parent is not None else None
            if parent_end_y is not None:
                return parent_end_y[-1] + self._y[1:]
            else:
                return self._y[1:]

        @property
        def end_z_(self):
            """
            The z coordinate (as a unitless floating point number) at the end of
            each compartment. Returns ``None`` for morphologies without coordinates.
            """
            if self._z is None:
                return None
            parent_end_z = self.parent.end_z_ if self.parent is not None else None
            if parent_end_z is not None:
                return parent_end_z[-1] + self._z[1:]
            else:
                return self._z[1:]

    class Cylinder(Section):
        """
        A cylindrical section. For sections with more complex geometry (varying
        length and/or diameter of each compartment), use the `Section` class.

        Parameters
        ----------
        diameter : `Quantity`
            The diameter of the cylinder.
        n : int, optional
            The number of compartments in this section. Defaults to 1.
        length : `Quantity`, optional
            The length of the cylinder. Cannot be combined with the specification
            of coordinates.
        x : `Quantity`, optional
            A sequence of two values, the start and the end point of the cylinder.
            The coordinates are interpreted as relative to the end point of the
            parent compartment (if any), so in most cases the start point should
            be ``0*um``. The common exception is a cylinder connecting to a `Soma`,
            here the start point can be used to make the cylinder start at the
            surface of the sphere instead of at its center.
            You can specify all of ``x``, ``y``, or ``z`` to specify
            a morphology in 3D, or only one or two out of them to specify a
            morphology in 1D or 2D.
        y : `Quantity`, optional
            See ``x``
        z : `Quantity`, optional
            See ``x``
        type : str, optional
            The type (e.g. ``"axon"``) of this `Cylinder`.
        """

        @check_units(n=1, length=meter, diameter=meter, x=meter, y=meter, z=meter)
        def __init__(self, diameter, n=1, length=None, x=None, y=None, z=None, type=None):
            n = int(n)
            Morphology.__init__(self, n=n, type=type)

            # Diameter
            if diameter.shape != () and (diameter.ndim > 1 or len(diameter) != 1):
                raise TypeError("The diameter argument has to be a single value.")
            diameter = np.ones(n) * diameter
            self._diameter = diameter

            if (x is not None or y is not None or z is not None) and length is not None:
                raise TypeError("Cannot specify coordinates and length at the same time.")

            if length is not None:
                # Length
                if length.shape != () and (length.ndim > 1 or len(length) != 1):
                    raise TypeError("The length argument has to be a single value.")
                self._length = np.ones(n) * (length / n)  # length was total length
                self._x = self._y = self._z = None
            else:
                # Coordinates
                if x is None and y is None and z is None:
                    raise TypeError(
                        "No length specified, need to specify at least "
                        "one out of x, y, or z."
                    )
                for name, value in [("x", x), ("y", y), ("z", z)]:
                    if value is not None and (value.ndim != 1 or len(value) != 2):
                        raise TypeError(
                            f"{name} needs to be a 1-dimensional array of length 2 (start "
                            "and end point)"
                        )
                self._x = (
                    np.asarray(np.linspace(x[0], x[1], n + 1))
                    if x is not None
                    else np.zeros(n + 1)
                )
                self._y = (
                    np.asarray(np.linspace(y[0], y[1], n + 1))
                    if y is not None
                    else np.zeros(n + 1)
                )
                self._z = (
                    np.asarray(np.linspace(z[0], z[1], n + 1))
                    if z is not None
                    else np.zeros(n + 1)
                )
                length = np.sqrt(
                    (self.end_x - self.start_x) ** 2
                    + (self.end_y - self.start_y) ** 2
                    + (self.end_z - self.start_z) ** 2
                )
                self._length = length

        def __repr__(self):
            s = f"{self.__class__.__name__}(diameter={self.diameter[0]!r}"
            if self.n != 1:
                s += f", n={self.n}"
            if self._x is not None:
                s += (
                    f", x={self._x[[0, -1]]!r}, y={self._y[[0, -1]]!r},"
                    f" z={self._z[[0, -1]]!r}"
                )
            else:
                s += f", length={sum(self._length)!r}"
            if self.type is not None:
                s += f", type={self.type!r}"
            return f"{s})"

        def copy_section(self):
            if self.x is None:
                return Cylinder(
                    self.diameter[0], n=self.n, length=self.length, type=self.type
                )
            else:
                return Cylinder(
                    self.diameter[0],
                    n=self.n,
                    x=self._x[[0, -1]],
                    y=self._y[[0, -1]],
                    z=self._z[[0, -1]],
                    type=self.type,
                )

        # Overwrite the properties that differ from `Section`
        @property
        def area(self):
            r"""
            The membrane surface area of each compartment in this section. The
            surface area of each compartment is calculated as
            :math:`\pi d l`,
            where :math:`l` is the length of the compartment, and :math:`d` is its
            diameter. Note that this surface area does not contain the area of
            the two disks at the two sides of the cylinder.
            """
            return np.pi * self._diameter * self.length

        @property
        def start_diameter(self):
            """
            The diameter at the start of each compartment in this section.
            """
            return self._diameter

        @property
        def diameter(self):
            """
            The diameter at the middle of each compartment in this section.
            """
            return self._diameter

        @property
        def end_diameter(self):
            """
            The diameter at the end of each compartment in this section.
            """
            return self._diameter

        @property
        def volume(self):
            r"""
            The volume of each compartment in this section. The volume of each
            compartment is calculated as
            :math:`\pi \frac{d}{2}^2 l` ,
            where :math:`l` is the length of the compartment, and :math:`d` is its
            diameter.
            """
            return np.pi * (self._diameter / 2) ** 2 * self.length

        @property
        def r_length_1(self):
            """
            The geometry-dependent term to calculate the conductance between the
            start and the midpoint of each compartment. Dividing this value by the
            Intracellular resistivity gives the conductance.
            """
            return np.pi / 2 * (self._diameter**2) / self.length

        @property
        def r_length_2(self):
            """
            The geometry-dependent term to calculate the conductance between the
            midpoint and the end of each compartment. Dividing this value by the
            Intracellular resistivity gives the conductance.
            """
            return np.pi / 2 * (self._diameter**2) / self.length