core.py

#!/usr/bin/env
# encoding: utf-8
"""
Author:     Daniel Boeckenhoff
Mail:       daniel.boeckenhoff@ipp.mpg.de

core of tfields library
contains numpy ndarray derived bases of the tfields package

Notes:
    It could be worthwhile concidering `np.li.mixins.NDArrayOperatorsMixin ...
    <https://docs.scipy.org/doc/numpy-1.15.1/reference/generated/...
    ... numpy.lib.mixins.NDArrayOperatorsMixin.html>`_
"""
# builtin
import warnings
import pathlib
from six import string_types
from contextlib import contextmanager
from collections import Counter
from copy import deepcopy

# 3rd party
import numpy as np
import sympy
import scipy as sp
import sortedcontainers
import rna

import tfields.bases

np.seterr(all="warn", over="raise")


def rank(tensor):
    """
    Tensor rank
    """
    tensor = np.asarray(tensor)
    return len(tensor.shape) - 1


def dim(tensor):
    """
    Manifold dimension
    """
    tensor = np.asarray(tensor)
    if rank(tensor) == 0:
        return 1
    return tensor.shape[1]


class AbstractObject(object):
    def save(self, path, *args, **kwargs):
        """
        Saving by redirecting to the correct save method depending on path

        Args:
            path (str or buffer)
            *args:
                forwarded to extension specific method
            **kwargs:
                extension (str): only needed if path is buffer
                ... remaining:forwarded to extension specific method
        """
        # get the extension
        if isinstance(path, (string_types, pathlib.Path)):
            extension = pathlib.Path(path).suffix.lstrip(".")
        else:
            raise ValueError("Wrong path type {0}".format(type(path)))
        path = str(path)

        # get the save method
        try:
            save_method = getattr(self, "_save_{extension}".format(**locals()))
        except:
            raise NotImplementedError(
                "Can not find save method for extension: "
                "{extension}.".format(**locals())
            )

        path = rna.path.resolve(path)
        return save_method(path, **kwargs)

    @classmethod
    def load(cls, path, *args, **kwargs):
        """
        load a file as a tensors object.

        Args:
            path (str or buffer)
            *args:
                forwarded to extension specific method
            **kwargs:
                extension (str): only needed if path is buffer
                ... remaining:forwarded to extension specific method
        """
        if isinstance(path, (string_types, pathlib.Path)):
            extension = pathlib.Path(path).suffix.lstrip(".")
            path = str(path)
            path = rna.path.resolve(path)
        else:
            extension = kwargs.pop("extension", "npz")

        try:
            load_method = getattr(cls, "_load_{e}".format(e=extension))
        except:
            raise NotImplementedError(
                "Can not find load method for extension: "
                "{extension}.".format(**locals())
            )
        return load_method(path, *args, **kwargs)

    def _save_npz(self, path, **kwargs):
        """
        Args:
            path (open file or str/unicode): destination to save file to.

        Examples:
            Build some dummies:
            >>> import tfields
            >>> from tempfile import NamedTemporaryFile
            >>> out_file = NamedTemporaryFile(suffix='.npz')
            >>> p = tfields.Points3D([[1., 2., 3.], [4., 5., 6.], [1, 2, -6]],
            ...                      name='my_points')
            >>> scalars = tfields.Tensors([0, 1, 2], name=42)
            >>> vectors = tfields.Tensors([[0, 0, 0], [0, 0, 1], [0, -1, 0]])
            >>> maps = [tfields.TensorFields([[0, 1, 2], [0, 1, 2]], [42, 21]),
            ...         tfields.TensorFields([[1], [2]], [-42, -21])]
            >>> m = tfields.TensorMaps(vectors, scalars,
            ...                        maps=maps)

            Simply give the file name to save
            >>> p.save(out_file.name)
            >>> _ = out_file.seek(0)  # this is only necessary in the test
            >>> p1 = tfields.Points3D.load(out_file.name)
            >>> assert p.equal(p1)
            >>> assert p.coord_sys == p1.coord_sys

            The fully nested structure of a TensorMaps object is reconstructed
            >>> out_file_maps = NamedTemporaryFile(suffix='.npz')
            >>> m.save(out_file_maps.name)
            >>> _ = out_file_maps.seek(0)
            >>> m1 = tfields.TensorMaps.load(out_file_maps.name,
            ...                              allow_pickle=True)
            >>> assert m.equal(m1)
            >>> assert m.maps[3].dtype == m1.maps[3].dtype

            Names are preserved
            >>> assert p.name == 'my_points'
            >>> m.names
            [42]

        """
        content_dict = self._as_dict()
        content_dict['tfields_version'] = tfields.__version__
        np.savez(path, **content_dict)

    @classmethod
    def _load_npz(cls, path, **load_kwargs):
        """
        Factory method
        Given a path to a npz file, construct the object
        """
        # TODO: think about allow_pickle, wheter it really should be True or
        # wheter we could avoid pickling (potential security issue)
        load_kwargs.setdefault('allow_pickle', True)
        np_file = np.load(path, **load_kwargs)
        d = dict(np_file)
        d.pop('tfields_version', None)
        return cls._from_dict(d)

    def _args(self) -> tuple:
        return tuple()

    def _kwargs(self) -> dict:
        return dict()

    _HIERARCHY_SEPARATOR = '::'

    def _as_dict(self):
        d = {}

        # type
        d["type"] = type(self).__name__

        # args and kwargs
        for base_attr, iterable in [
                ('args', ((str(i), arg)
                          for i, arg in enumerate(self._args()))),
                ('kwargs', self._kwargs().items())]:
            for attr, value in iterable:
                attr = base_attr + self._HIERARCHY_SEPARATOR + attr
                if hasattr(value, '_as_dict'):
                    part_dict = value._as_dict()
                    for part_attr, part_value in part_dict.items():
                        d[
                            attr + self._HIERARCHY_SEPARATOR + part_attr
                        ] = part_value
                else:
                    d[attr] = value
        return d

    @classmethod
    def _from_dict(cls, d: dict):
        try:
            d.pop('type')
        except KeyError:
            # legacy
            return cls._from_dict_legacy(**d)

        here = {}
        for string in d:  # TOO no sortelist
            value = d[string]

            attr, _, end = string.partition(cls._HIERARCHY_SEPARATOR)
            key, _, end = end.partition(cls._HIERARCHY_SEPARATOR)
            if attr not in here:
                here[attr] = {}
            if key not in here[attr]:
                here[attr][key] = {}
            here[attr][key][end] = value

        """
        Do the recursion
        """
        for attr in here:
            for key in here[attr]:
                if 'type' in here[attr][key]:
                    obj_type = here[attr][key].get("type")
                    if isinstance(obj_type, np.ndarray):  # happens on np.load
                        obj_type = obj_type.tolist()
                    if isinstance(obj_type, bytes):
                        # asthonishingly, this is not necessary under linux.
                        # Found under nt. ???
                        obj_type = obj_type.decode("UTF-8")
                    obj_type = getattr(tfields, obj_type)
                    attr_value = obj_type._from_dict(here[attr][key])
                else:  # if len(here[attr][key]) == 1:
                    attr_value = here[attr][key].pop('')
                here[attr][key] = attr_value

        '''
        Build the generic way
        '''
        args = here.pop('args', tuple())
        args = tuple(args[key] for key in sorted(args))
        kwargs = here.pop('kwargs', {})
        assert len(here) == 0
        obj = cls(*args, **kwargs)
        return obj

    @classmethod
    def _from_dict_legacy(cls, **d):
        """
        legacy method of _from_dict - Opposite of old _as_dict method
        which is overridden in this version
        """
        list_dict = {}
        kwargs = {}
        """
        De-Flatten the first layer of lists
        """
        for key in sorted(list(d)):
            if "::" in key:
                attr, _, end = key.partition("::")
                if attr not in list_dict:
                    list_dict[attr] = {}

                index, _, end = end.partition("::")
                if not index.isdigit():
                    raise ValueError("None digit index given")
                index = int(index)
                if index not in list_dict[attr]:
                    list_dict[attr][index] = {}
                list_dict[attr][index][end] = d[key]
            else:
                kwargs[key] = d[key]

        """
        Build the lists (recursively)
        """
        for key in list(list_dict):
            sub_dict = list_dict[key]
            list_dict[key] = []
            for index in sorted(list(sub_dict)):
                bulk_type = sub_dict[index].get("bulk_type")
                bulk_type = bulk_type.tolist()
                if isinstance(bulk_type, bytes):
                    # asthonishingly, this is not necessary under linux.
                    # Found under nt. ???
                    bulk_type = bulk_type.decode("UTF-8")
                bulk_type = getattr(tfields, bulk_type)
                list_dict[key].append(bulk_type._from_dict_legacy(**sub_dict[index]))

        with cls._bypass_setters('fields', demand_existence=False):
            '''
            Build the normal way
            '''
            bulk = kwargs.pop('bulk')
            bulk_type = kwargs.pop('bulk_type')
            obj = cls.__new__(cls, bulk, **kwargs)

            '''
            Set list attributes
            '''
            for attr, list_value in list_dict.items():
                setattr(obj, attr, list_value)
        return obj


class AbstractNdarray(np.ndarray, AbstractObject):
    """
    All tensors and subclasses should derive from AbstractNdarray.
    AbstractNdarray implements all the inheritance specifics for np.ndarray
    Whene inheriting, three attributes are of interest:

    Attributes:
        __slots__ (List(str)): If you want to add attributes to
            your AbstractNdarray subclass, add the attribute name to __slots__
        __slot_defaults__ (list): if __slot_defaults__ is None, the
            defaults for the attributes in __slots__ will be None
            other values will be treaded as defaults to the corresponding
            arg at the same position in the __slots__ list.
        __slot_dtype__ (List(dtypes)): for the conversion of the
            args in __slots__ to numpy arrays. None values mean no
            conversion.
        __slot_setters__ (List(callable)): Because __slots__ and properties are
            mutually exclusive this is a possibility to take care of proper
            attribute handling. None will be passed for 'not set'.

    Args:
        array (array-like): input array
        **kwargs: arguments corresponding to __slots__

    TODO:
        equality check

    """

    __slots__ = []
    __slot_defaults__ = []
    __slot_dtypes__ = []
    __slot_setters__ = []

    def __new__(cls, array, **kwargs):  # pragma: no cover
        raise NotImplementedError(
            "{clsType} type must implement '__new__'".format(clsType=type(cls))
        )

    def __array_finalize__(self, obj):
        if obj is None:
            return
        for attr in self._iter_slots():
            setattr(self, attr, getattr(obj, attr, None))

    def __array_wrap__(self, out_arr, context=None):
        return np.ndarray.__array_wrap__(self, out_arr, context)

    @classmethod
    def _iter_slots(cls):
        return [att for att in cls.__slots__ if att != "_cache"]

    @classmethod
    def _update_slot_kwargs(cls, kwargs):
        """
        set the defaults in kwargs according to __slot_defaults__
        and convert the kwargs according to __slot_dtypes__
        """
        slot_defaults = cls.__slot_defaults__ + [None] * (
            len(cls.__slots__) - len(cls.__slot_defaults__)
        )
        slot_dtypes = cls.__slot_dtypes__ + [None] * (
            len(cls.__slots__) - len(cls.__slot_dtypes__)
        )
        for attr, default, dtype in zip(
            cls.__slots__, slot_defaults, slot_dtypes
        ):
            if attr == "_cache":
                continue
            if attr not in kwargs:
                kwargs[attr] = default
            if dtype is not None:
                try:
                    kwargs[attr] = np.array(kwargs[attr], dtype=dtype)
                except Exception as err:
                    raise ValueError(
                        str(attr) + str(dtype) + str(kwargs[attr]) + str(err)
                    )

    def __setattr__(self, name, value):
        if name in self.__slots__:
            index = self.__slots__.index(name)
            try:
                setter = self.__slot_setters__[index]
            except IndexError:
                setter = None
            if isinstance(setter, str):
                setter = getattr(self, setter)
            if setter is not None:
                value = setter(value)
        super(AbstractNdarray, self).__setattr__(name, value)

    def _args(self):
        return (np.array(self),)

    def _kwargs(self):
        return dict((attr, getattr(self, attr)) for attr in self._iter_slots())

    def __reduce__(self):
        """
        important for pickling (see `here <https://stackoverflow.com/questions/26598109/preserve-custom-attributes-when-pickling-subclass-of-numpy-array>`_)

        Examples:
            >>> from tempfile import NamedTemporaryFile
            >>> import pickle
            >>> import tfields

            Build a dummy scalar field

            >>> from tfields import Tensors, TensorFields
            >>> scalars = Tensors([0, 1, 2])
            >>> vectors = Tensors([[0, 0, 0], [0, 0, 1], [0, -1, 0]])
            >>> scalar_field = TensorFields(vectors,
            ...                             scalars,
            ...                             coord_sys='cylinder')

            Save it and restore it

            >>> out_file = NamedTemporaryFile(suffix='.pickle')

            >>> pickle.dump(scalar_field,
            ...             out_file)
            >>> _ = out_file.seek(0)

            >>> sf = pickle.load(out_file)
            >>> sf.coord_sys == 'cylinder'
            True
            >>> sf.fields[0][2] == 2.
            True

        """
        # Get the parent's __reduce__ tuple
        pickled_state = super(AbstractNdarray, self).__reduce__()

        # Create our own tuple to pass to __setstate__
        new_state = pickled_state[2] + tuple(
            [getattr(self, slot) for slot in self._iter_slots()]
        )

        # Return a tuple that replaces the parent's __setstate__
        # tuple with our own
        return (pickled_state[0], pickled_state[1], new_state)

    def __setstate__(self, state):
        """
        Counterpart to __reduce__. Important for unpickling.
        """
        # Call the parent's __setstate__ with the other tuple elements.
        super(AbstractNdarray, self).__setstate__(
            state[0:-len(self._iter_slots())]
        )

        # set the __slot__ attributes
        valid_slot_attrs = list(self._iter_slots())
        added_slot_attrs = ['name']  # attributes that have been added later
                                     # have not been pickled with the full
                                     # information and thus need to be
                                     # excluded from the __setstate__
                                     # need to be in the same order as they
                                     # have been added to __slots__
        n_np = 5  # number of numpy array states
        n_old = len(valid_slot_attrs) - len(state[n_np:])
        if n_old > 0:
            for latest_index in range(n_old):
                new_slot = added_slot_attrs[-latest_index]
                warnings.warn("Slots with names '{new_slot}' appears to have "
                              "been added after the creation of the reduced "
                              "state. No corresponding state found in "
                              "__setstate__."
                              .format(**locals()))
                valid_slot_attrs.pop(valid_slot_attrs.index(new_slot))
                setattr(self, new_slot, None)

        for slot_index, slot in enumerate(valid_slot_attrs):
            state_index = n_np + slot_index
            setattr(self, slot, state[state_index])

    @property
    def bulk(self):
        """
        The pure ndarray version of the actual state
            -> nothing attached
        """
        return np.array(self)

    @classmethod
    @contextmanager
    def _bypass_setters(cls, *slots,
                        empty_means_all=True,
                        demand_existence=False):
        """
        Temporarily remove the setter in __slot_setters__ corresponding to slot
        position in __slot__. You should know what you do, when using this.

        Args:
            *slots (str): attribute names in __slots__
            empty_means_all (bool): defines behaviour when slots is empty.
                When True: if slots is empty mute all slots in __slots__
            demand_existence (bool): if false do not check the existence of the
                slot in __slots__ - do nothing for that slot. Handle with care!
        """
        if not slots and empty_means_all:
            slots = cls.__slots__
        slot_indices = []
        setters = []
        for slot in slots:
            slot_index = cls.__slots__.index(slot)\
                if slot in cls.__slots__ else None
            if slot_index is None:
                # slot not in cls.__slots__.
                if demand_existence:
                    raise ValueError(
                        "Slot {slot} not existing".format(**locals()))
                continue
            if len(cls.__slot_setters__) < slot_index + 1:
                # no setter to be found
                continue
            slot_indices.append(slot_index)
            setter = cls.__slot_setters__[slot_index]
            setters.append(setter)
            cls.__slot_setters__[slot_index] = None
        yield
        for slot_index, setter in zip(slot_indices, setters):
            cls.__slot_setters__[slot_index] = setter

    def copy(self, *args, **kwargs):
        """
        The standard ndarray copy does not copy slots. Correct for this.

        Examples:
            >>> import tfields
            >>> m = tfields.TensorMaps(
            ...     [[1,2,3], [3,3,3], [0,0,0], [5,6,7]],
            ...     [[1], [3], [0], [5]],
            ...     maps=[
            ...         ([[0, 1, 2], [1, 2, 3]], [21, 42]),
            ...         [[1]],
            ...         [[0, 1, 2, 3]]
            ...     ])
            >>> mc = m.copy()
            >>> mc.equal(m)
            True
            >>> mc is m
            False
            >>> mc.fields is m.fields
            False
            >>> mc.fields[0] is m.fields[0]
            False
            >>> mc.maps[3].fields[0] is m.maps[3].fields[0]
            False

        """
        # works with __reduce__ / __setstate__
        return deepcopy(self)


class Tensors(AbstractNdarray):
    """
    Set of tensors with the same basis.

    TODO:
        all slot args should be protected -> _base

    Args:
        tensors: np.ndarray or AbstractNdarray subclass
        **kwargs:
            name: optional - custom name, can be anything

    Examples:
        >>> import numpy as np
        >>> import tfields

        Initialize a scalar range

        >>> scalars = tfields.Tensors([0, 1, 2])
        >>> scalars.rank == 0
        True

        Initialize vectors

        >>> vectors = tfields.Tensors([[0, 0, 0], [0, 0, 1], [0, -1, 0]])
        >>> vectors.rank == 1
        True
        >>> vectors.dim == 3
        True
        >>> assert vectors.coord_sys == 'cartesian'

        Initialize the Levi-Zivita Tensor

        >>> matrices = tfields.Tensors(
        ...                     [[[0, 0, 0], [0, 0, 1], [0, -1, 0]],
        ...                      [[0, 0, -1], [0, 0, 0], [1, 0, 0]],
        ...                      [[0, 1, 0], [-1, 0, 0], [0, 0, 0]]])
        >>> matrices.shape == (3, 3, 3)
        True
        >>> matrices.rank == 2
        True
        >>> matrices.dim == 3
        True

        Initializing in different start coordinate system

        >>> cyl = tfields.Tensors([[5, np.arctan(4. / 3.), 42]],
        ...                       coord_sys='cylinder')
        >>> assert cyl.coord_sys == 'cylinder'
        >>> cyl.transform('cartesian')
        >>> assert cyl.coord_sys == 'cartesian'
        >>> cart = cyl
        >>> assert round(cart[0, 0], 10) == 3.
        >>> assert round(cart[0, 1], 10) == 4.
        >>> assert cart[0, 2] == 42

        Initialize with copy constructor keeps the coordinate system

        >>> with vectors.tmp_transform('cylinder'):
        ...     vect_cyl = tfields.Tensors(vectors)
        ...     assert vect_cyl.coord_sys == vectors.coord_sys
        >>> assert vect_cyl.coord_sys == 'cylinder'

        You can demand a special dimension.

        >>> _ = tfields.Tensors([[1, 2, 3]], dim=3)
        >>> _ = tfields.Tensors([[1, 2, 3]], dim=2)  # doctest: +ELLIPSIS
        Traceback (most recent call last):
            ...
        ValueError: Incorrect dimension: 3 given, 2 demanded.

        The dimension argument (dim) becomes necessary if you want to initialize
        an empty array

        >>> _ = tfields.Tensors([])  # doctest: +ELLIPSIS
        Traceback (most recent call last):
            ...
        ValueError: Empty tensors need dimension parameter 'dim'.
        >>> tfields.Tensors([], dim=7)
        Tensors([], shape=(0, 7), dtype=float64)

    """
    __slots__ = ['coord_sys', 'name']
    __slot_defaults__ = ['cartesian']
    __slot_setters__ = [tfields.bases.get_coord_system_name]

    def __new__(cls, tensors, **kwargs):
        dtype = kwargs.pop("dtype", None)
        order = kwargs.pop("order", None)
        dim = kwargs.pop("dim", None)

        """ copy constructor extracts the kwargs from tensors"""
        if issubclass(type(tensors), Tensors):
            if dim is not None:
                dim = tensors.dim
            coord_sys = kwargs.pop("coord_sys", tensors.coord_sys)
            tensors = tensors.copy()
            tensors.transform(coord_sys)
            kwargs['coord_sys'] = coord_sys
            kwargs['name'] = kwargs.pop('name', tensors.name)
            if dtype is None:
                dtype = tensors.dtype
        else:
            if dtype is None:
                if hasattr(tensors, "dtype"):
                    dtype = tensors.dtype
                else:
                    dtype = np.float64

        """ demand iterable structure """
        try:
            len(tensors)
        except TypeError as err:
            raise TypeError(
                "Iterable structure necessary."
                " Got {tensors}".format(**locals())
            )

        """ process empty inputs """
        if len(tensors) == 0:
            if issubclass(type(tensors), tfields.Tensors):
                tensors = np.empty(tensors.shape, dtype=tensors.dtype)
            elif dim is not None:
                tensors = np.empty((0, dim))
            if issubclass(type(tensors), np.ndarray):
                # np.empty
                pass
            elif hasattr(tensors, 'shape'):
                dim = dim(tensors)
            else:
                raise ValueError(
                    "Empty tensors need dimension parameter 'dim'."
                )

        tensors = np.asarray(tensors, dtype=dtype, order=order)
        obj = tensors.view(cls)

        """ check dimension(s) """
        for d in obj.shape[1:]:
            if not d == obj.dim:
                raise ValueError(
                    "Dimensions are inconstistent. "
                    "Manifold dimension is {obj.dim}. "
                    "Found dimensions {found} in {obj}.".format(
                        found=obj.shape[1:], **locals()
                    )
                )
        if dim is not None:
            if dim != obj.dim:
                raise ValueError(
                    "Incorrect dimension: {obj.dim} given,"
                    " {dim} demanded.".format(**locals())
                )

        """ update kwargs with defaults from slots """
        cls._update_slot_kwargs(kwargs)

        """ set kwargs to slots attributes """
        for attr in kwargs:
            if attr not in cls._iter_slots():
                raise AttributeError(
                    "Keyword argument {attr} not accepted "
                    "for class {cls}".format(**locals())
                )
            setattr(obj, attr, kwargs[attr])

        return obj

    def __iter__(self):
        """
        Forwarding iterations to the bulk array. Otherwise __getitem__ would
        kick in and slow down imensely.

        Examples:
            >>> import tfields
            >>> vectors = tfields.Tensors([[0, 0, 0], [0, 0, 1], [0, -1, 0]])
            >>> scalar_field = tfields.TensorFields(
            ...     vectors, [42, 21, 10.5], [1, 2, 3])
            >>> [(point.rank, point.dim) for point in scalar_field]
            [(0, 1), (0, 1), (0, 1)]

        """
        for index in range(len(self)):
            yield super(Tensors, self).__getitem__(index).view(Tensors)

    @classmethod
    def merged(cls, *objects, **kwargs):
        """
        Factory method
        Merges all input arguments to one object

        Args:
            return_templates (bool): return the templates which can be used
                together with cut to retrieve the original objects
            dim (int):
            **kwargs: passed to cls

        Examples:
            >>> import numpy as np
            >>> import tfields
            >>> import tfields.bases

            The new object with turn out in the most frequent coordinate
            system if not specified explicitly

            >>> vec_a = tfields.Tensors([[0, 0, 0], [0, 0, 1], [0, -1, 0]])
            >>> vec_b = tfields.Tensors([[5, 4, 1]],
            ...     coord_sys=tfields.bases.cylinder)
            >>> vec_c = tfields.Tensors([[4, 2, 3]],
            ...     coord_sys=tfields.bases.cylinder)
            >>> merge = tfields.Tensors.merged(
            ...     vec_a, vec_b, vec_c, [[2, 0, 1]])
            >>> assert merge.coord_sys == 'cylinder'
            >>> assert merge.equal([[0, 0, 0],
            ...                     [0, 0, 1],
            ...                     [1, -np.pi / 2, 0],
            ...                     [5, 4, 1],
            ...                     [4, 2, 3],
            ...                     [2, 0, 1]])

            Merge also shifts the maps to still refer to the same tensors

            >>> tm_a = tfields.TensorMaps(merge, maps=[[[0, 1, 2]]])
            >>> tm_b = tm_a.copy()
            >>> assert tm_a.coord_sys == 'cylinder'
            >>> tm_merge = tfields.TensorMaps.merged(tm_a, tm_b)
            >>> assert tm_merge.coord_sys == 'cylinder'
            >>> assert tm_merge.maps[3].equal([[0, 1, 2],
            ...                               list(range(len(merge),
            ...                                          len(merge) + 3,
            ...                                          1))])

            >>> obj_list = [tfields.Tensors([[1, 2, 3]],
            ...             coord_sys=tfields.bases.CYLINDER),
            ...             tfields.Tensors([[3] * 3]),
            ...             tfields.Tensors([[5, 1, 3]])]
            >>> merge2 = tfields.Tensors.merged(
            ...     *obj_list, coord_sys=tfields.bases.CARTESIAN)
            >>> assert merge2.equal([[-0.41614684, 0.90929743, 3.],
            ...                      [3, 3, 3], [5, 1, 3]], atol=1e-8)

            The return_templates argument allows to retrieve a template which
            can be used with the cut method.

            >>> merge, templates = tfields.Tensors.merged(
            ...     vec_a, vec_b, vec_c, return_templates=True)
            >>> assert merge.cut(templates[0]).equal(vec_a)
            >>> assert merge.cut(templates[1]).equal(vec_b)
            >>> assert merge.cut(templates[2]).equal(vec_c)

        """

        """ get most frequent coord_sys or predefined coord_sys """
        coord_sys = kwargs.get("coord_sys", None)
        return_templates = kwargs.pop("return_templates", False)
        if coord_sys is None:
            bases = []
            for t in objects:
                try:
                    bases.append(t.coord_sys)
                except AttributeError:
                    pass
            if bases:
                # get most frequent coord_sys
                coord_sys = sorted(
                    bases, key=Counter(bases).get, reverse=True
                )[0]
                kwargs["coord_sys"] = coord_sys
            else:
                default = cls.__slot_defaults__[
                    cls.__slots__.index("coord_sys")
                ]
                kwargs["coord_sys"] = default

        """ transform all raw inputs to cls type with correct coord_sys. Also
        automatically make a copy of those instances that are of the correct
        type already."""
        objects = [cls.__new__(cls, t, **kwargs) for t in objects]

        """ check rank and dimension equality """
        if not len(set([t.rank for t in objects])) == 1:
            raise TypeError("Tensors must have the same rank for merging.")
        if not len(set([t.dim for t in objects])) == 1:
            raise TypeError("Tensors must have the same dimension for merging.")

        """ merge all objects """
        remainingObjects = objects[1:] or []
        tensors = objects[0]

        for i, obj in enumerate(remainingObjects):
            tensors = np.append(tensors, obj, axis=0)

        if len(tensors) == 0 and not kwargs.get('dim', None):
            # if you can not determine the tensor dimension, search for the
            # first object with some entries
            kwargs['dim'] = dim(objects[0])

        inst = cls.__new__(cls, tensors, **kwargs)
        if not return_templates:
            return inst
        else:
            tensor_lengths = [len(o) for o in objects]
            cum_tensor_lengths = [sum(tensor_lengths[:i])
                                  for i in range(len(objects))]
            templates = [
                tfields.TensorFields(
                    np.empty((len(obj), 0)),
                    np.arange(tensor_lengths[i]) + cum_tensor_lengths[i])
                for i, obj in enumerate(objects)]
            return inst, templates

    @classmethod
    def grid(cls, *base_vectors, **kwargs):
        """
        Args:
            *base_vectors (Iterable): base coordinates. The amount of base
                vectors defines the dimension

            **kwargs:
                iter_order (list): order in which the iteration will be done.
                    Frequency rises with position in list. default is [0, 1, 2]
                    iteration will be done like::

                    for v0 in base_vectors[iter_order[0]]:
                        for v1 in base_vectors[iter_order[1]]:
                            for v2 in base_vectors[iter_order[2]]:
                                coords0.append(locals()['v%i' % iter_order[0]])
                                coords1.append(locals()['v%i' % iter_order[1]])
                                coords2.append(locals()['v%i' % iter_order[2]])

        Examples:
            Initilaize using the mgrid notation

            >>> import tfields
            >>> mgrid = tfields.Tensors.grid((0, 1, 2j), (3, 4, 2j), (6, 7, 2j))
            >>> mgrid.equal([[0, 3, 6],
            ...              [0, 3, 7],
            ...              [0, 4, 6],
            ...              [0, 4, 7],
            ...              [1, 3, 6],
            ...              [1, 3, 7],
            ...              [1, 4, 6],
            ...              [1, 4, 7]])
            True

            Lists or arrays are accepted also.
            Furthermore, the iteration order can be changed

            >>> lins = tfields.Tensors.grid(
            ...     np.linspace(3, 4, 2), np.linspace(0, 1, 2),
            ...     np.linspace(6, 7, 2), iter_order=[1, 0, 2])
            >>> lins.equal([[3, 0, 6],
            ...             [3, 0, 7],
            ...             [4, 0, 6],
            ...             [4, 0, 7],
            ...             [3, 1, 6],
            ...             [3, 1, 7],
            ...             [4, 1, 6],
            ...             [4, 1, 7]])
            True
            >>> lins2 = tfields.Tensors.grid(np.linspace(0, 1, 2),
            ...                              np.linspace(3, 4, 2),
            ...                              np.linspace(6, 7, 2),
            ...                              iter_order=[2, 0, 1])
            >>> lins2.equal([[0, 3, 6],
            ...              [0, 4, 6],
            ...              [1, 3, 6],
            ...              [1, 4, 6],
            ...              [0, 3, 7],
            ...              [0, 4, 7],
            ...              [1, 3, 7],
            ...              [1, 4, 7]])
            True

            When given the coord_sys argument, the grid is performed in the
            given coorinate system:

            >>> lins3 = tfields.Tensors.grid(np.linspace(4, 9, 2),
            ...                              np.linspace(np.pi/2, np.pi/2, 1),
            ...                              np.linspace(4, 4, 1),
            ...                              iter_order=[2, 0, 1],
            ...                              coord_sys=tfields.bases.CYLINDER)
            >>> assert lins3.coord_sys == 'cylinder'
            >>> lins3.transform('cartesian')
            >>> assert np.array_equal(lins3[:, 1], [4, 9])

        """
        cls_kwargs = {
            attr: kwargs.pop(attr)
            for attr in list(kwargs)
            if attr in cls.__slots__
        }
        inst = cls.__new__(
            cls, tfields.lib.grid.igrid(*base_vectors, **kwargs), **cls_kwargs
        )
        return inst

    @property
    def rank(self):
        """
        Tensor rank
        """
        return rank(self)

    @property
    def dim(self):
        """
        Manifold dimension
        """
        return dim(self)

    def transform(self, coord_sys):
        """
        Args:
            coord_sys (str)

        Examples:
            >>> import numpy as np
            >>> import tfields

            CARTESIAN to SPHERICAL
            >>> t = tfields.Tensors([[1, 2, 2], [1, 0, 0], [0, 0, -1],
            ...                      [0, 0, 1], [0, 0, 0]])
            >>> t.transform('spherical')

            r

            >>> assert t[0, 0] == 3

            phi

            >>> assert t[1, 1] == 0.
            >>> assert t[2, 1] == 0.

            theta is 0 at (0, 0, 1) and pi / 2 at (0, 0, -1)

            >>> assert round(t[1, 2], 10) == round(0, 10)
            >>> assert t[2, 2] == -np.pi / 2
            >>> assert t[3, 2] == np.pi / 2

            theta is defined 0 for R == 0

            >>> assert t[4, 0] == 0.
            >>> assert t[4, 2] == 0.


            CARTESIAN to CYLINDER

            >>> tCart = tfields.Tensors([[3, 4, 42], [1, 0, 0], [0, 1, -1],
            ...                          [-1, 0, 1], [0, 0, 0]])
            >>> t_cyl = tCart.copy()
            >>> t_cyl.transform('cylinder')
            >>> assert t_cyl.coord_sys == 'cylinder'

            R

            >>> assert t_cyl[0, 0] == 5
            >>> assert t_cyl[1, 0] == 1
            >>> assert t_cyl[2, 0] == 1
            >>> assert t_cyl[4, 0] == 0

            Phi

            >>> assert round(t_cyl[0, 1], 10) == round(np.arctan(4. / 3), 10)
            >>> assert t_cyl[1, 1] == 0
            >>> assert round(t_cyl[2, 1], 10) == round(np.pi / 2, 10)
            >>> assert t_cyl[1, 1] == 0

            Z

            >>> assert t_cyl[0, 2] == 42
            >>> assert t_cyl[2, 2] == -1

            >>> t_cyl.transform('cartesian')
            >>> assert t_cyl.coord_sys == 'cartesian'
            >>> assert t_cyl[0, 0] == 3

        """
        #           scalars                 empty             already there
        if self.rank == 0 or self.shape[0] == 0 or self.coord_sys == coord_sys:
            self.coord_sys = coord_sys
            return

        tfields.bases.transform(self, self.coord_sys, coord_sys)
        # self[:] = tfields.bases.transform(self, self.coord_sys, coord_sys)
        self.coord_sys = coord_sys

    @contextmanager
    def tmp_transform(self, coord_sys):
        """
        Temporarily change the coord_sys to another coord_sys and change it back at exit
        This method is for cleaner code only.
        No speed improvements go with this.

        Args:
            see transform

        Examples:
            >>> import tfields
            >>> p = tfields.Tensors([[1,2,3]], coord_sys=tfields.bases.SPHERICAL)
            >>> with p.tmp_transform(tfields.bases.CYLINDER):
            ...     assert p.coord_sys == tfields.bases.CYLINDER
            >>> assert p.coord_sys == tfields.bases.SPHERICAL

        """
        baseBefore = self.coord_sys
        if baseBefore == coord_sys:
            yield
        else:
            self.transform(coord_sys)

            yield

            self.transform(baseBefore)

    def mirror(self, coordinate, condition=None):
        """
        Reflect/Mirror the entries meeting <condition> at <coordinate> = 0

        Args:
            coordinate (int): coordinate index

        Examples:
            >>> import tfields
            >>> p = tfields.Tensors([[1., 2., 3.], [4., 5., 6.], [1, 2, -6]])
            >>> p.mirror(1)
            >>> assert p.equal([[1, -2, 3], [4, -5,  6], [1, -2, -6]])

            multiple coordinates can be mirrored at the same time
            i.e. a point mirrorion would be

            >>> p = tfields.Tensors([[1., 2., 3.], [4., 5., 6.], [1, 2, -6]])
            >>> p.mirror([0,2])
            >>> assert p.equal([[-1, 2, -3], [-4, 5, -6], [-1, 2., 6.]])

            You can give a condition as mask or as str.
            The mirroring will only be applied to the points meeting the
            condition.

            >>> import sympy
            >>> x, y, z = sympy.symbols('x y z')
            >>> p.mirror([0, 2], y > 3)
            >>> p.equal([[-1, 2, -3], [4, 5, 6], [-1, 2, 6]])
            True

        """
        if condition is None:
            condition = np.array([True for i in range(len(self))])
        elif isinstance(condition, sympy.Basic):
            condition = self.evalf(condition)
        if isinstance(coordinate, list) or isinstance(coordinate, tuple):
            for c in coordinate:
                self.mirror(c, condition=condition)
        elif isinstance(coordinate, int):
            self[:, coordinate][condition] *= -1
        else:
            raise TypeError()

    def to_segment(
        self,
        segment,
        num_segments,
        coordinate,
        periodicity=2 * np.pi,
        offset=0.0,
        coord_sys=None,
    ):
        """
        For circular (close into themself after
        <periodicity>) coordinates at index <coordinate> assume
        <num_segments> segments and transform all values to
        segment number <segment>

        Args:
            segment (int): segment index (starting at 0)
            num_segments (int): number of segments
            coordinate (int): coordinate index
            periodicity (float): after what lenght, the coordiante repeats
            offset (float): offset in the mapping
            coord_sys (str or sympy.CoordinateSystem): in which coord sys the
                transformation should be done

        Examples:
            >>> import tfields
            >>> import numpy as np
            >>> pStart = tfields.Points3D([[6, 2 * np.pi, 1],
            ...                            [6, 2 * np.pi / 5 * 3, 1]],
            ...                           coord_sys='cylinder')
            >>> p = tfields.Points3D(pStart)
            >>> p.to_segment(0, 5, 1, offset=-2 * np.pi / 10)
            >>> assert np.array_equal(p[:, 1], [0, 0])

            >>> p2 = tfields.Points3D(pStart)
            >>> p2.to_segment(1, 5, 1, offset=-2 * np.pi / 10)
            >>> assert np.array_equal(np.round(p2[:, 1], 4), [1.2566] * 2)

        """
        if segment > num_segments - 1:
            raise ValueError("Segment {0} not existent.".format(segment))

        if coord_sys is None:
            coord_sys = self.coord_sys
        with self.tmp_transform(coord_sys):
            # map all values to first segment
            self[:, coordinate] = (
                (self[:, coordinate] - offset) % (periodicity / num_segments)
                + offset
                + segment * periodicity / num_segments
            )

    def equal(
        self, other, rtol=None, atol=None, equal_nan=False, return_bool=True
    ):
        """
        Evaluate, whether the instance has the same content as other.

        Args:
            optional:
                rtol (float)
                atol (float)
                equal_nan (bool)
            see numpy.isclose
        """
        if (
            issubclass(type(other), Tensors)
            and self.coord_sys != other.coord_sys
        ):
            other = other.copy()
            other.transform(self.coord_sys)
        x, y = np.asarray(self), np.asarray(other)
        if rtol is None and atol is None:
            mask = x == y
            if equal_nan:
                both_nan = np.isnan(x) & np.isnan(y)
                mask[both_nan] = both_nan[both_nan]
        else:
            if rtol is None:
                rtol = 0.0
            if atol is None:
                atol = 0.0
            mask = np.isclose(x, y, rtol=rtol, atol=atol, equal_nan=equal_nan)
        if return_bool:
            return bool(np.all(mask))
        return mask

    def contains(self, other, **kwargs):
        """
        Inspired by a speed argument @
        stackoverflow.com/questions/14766194/testing-whether-a-numpy-array-contains-a-given-row

        Examples:
            >>> import tfields
            >>> p = tfields.Tensors([[1,2,3], [4,5,6], [6,7,8]])
            >>> p.contains([4,5,6])
            True

        """
        return any(self.equal(other, return_bool=False).all(1))

    def indices(self, tensor, rtol=None, atol=None):
        """
        Returns:
            list of int: indices of tensor occuring

        Examples:
            Rank 1 Tensors

            >>> import tfields
            >>> p = tfields.Tensors([[1,2,3], [4,5,6], [6,7,8], [4,5,6],
            ...                      [4.1, 5, 6]])
            >>> p.indices([4,5,6])
            array([1, 3])
            >>> p.indices([4,5,6.1], rtol=1e-5, atol=1e-1)
            array([1, 3, 4])

            Rank 0 Tensors

            >>> p = tfields.Tensors([2, 3, 6, 3.01])
            >>> p.indices(3)
            array([1])
            >>> p.indices(3, rtol=1e-5, atol=1e-1)
            array([1, 3])

        """
        x, y = np.asarray(self), np.asarray(tensor)
        if rtol is None and atol is None:
            equal_method = np.equal
        else:
            equal_method = lambda a, b: np.isclose(a, b, rtol=rtol, atol=atol)

        # inspired by https://stackoverflow.com/questions/19228295/find-ordered-vector-in-numpy-array
        if self.rank == 0:
            indices = np.where(equal_method((x - y), 0))[0]