diff --git a/geometer/__init__.py b/geometer/__init__.py index 2a349e8..36cc4de 100644 --- a/geometer/__init__.py +++ b/geometer/__init__.py @@ -1,5 +1,5 @@ from geometer.__version__ import __version__ -from geometer.curve import Circle, Cone, Conic, Cylinder, Ellipse, Quadric, Sphere +from geometer.curve import Circle, Cone, Conic, Cylinder, Ellipse, Quadric, QuadricCollection, Sphere from geometer.operators import ( angle, angle_bisectors, @@ -12,10 +12,35 @@ is_coplanar, is_perpendicular, ) -from geometer.point import I, J, Line, Plane, Point, infty, infty_plane, join, meet -from geometer.shapes import Cuboid, Polygon, Polyhedron, Polytope, Rectangle, RegularPolygon, Segment, Simplex, Triangle +from geometer.point import ( + I, + J, + Line, + LineCollection, + Plane, + PlaneCollection, + Point, + PointCollection, + infty, + infty_plane, + join, + meet, +) +from geometer.shapes import ( + Cuboid, + Polygon, + PolygonCollection, + Polyhedron, + Rectangle, + RegularPolygon, + Segment, + SegmentCollection, + Simplex, + Triangle, +) from geometer.transformation import ( Transformation, + TransformationCollection, affine_transform, identity, reflection, @@ -23,11 +48,3 @@ scaling, translation, ) - -PointCollection = Point -LineCollection = Line -PlaneCollection = Plane -QuadricCollection = Quadric -TransformationCollection = Transformation -SegmentCollection = Segment -PolygonCollection = Polygon diff --git a/geometer/base.py b/geometer/base.py index 2c6b7e3..e973d1b 100644 --- a/geometer/base.py +++ b/geometer/base.py @@ -1,14 +1,14 @@ from __future__ import annotations from abc import ABC -from collections.abc import Iterable, Iterator, Sequence, Sized +from collections.abc import Iterable, Iterator, Sized from itertools import permutations -from typing import TYPE_CHECKING, Any, ClassVar, Literal, TypedDict, Union +from typing import TYPE_CHECKING, Any, Callable, ClassVar, Generic, Literal, TypeVar import numpy as np import numpy.typing as npt -from geometer.exceptions import TensorComputationError +from geometer.exceptions import IncompatibleShapeError, TensorComputationError from geometer.utils import ( is_multiple, is_numerical_dtype, @@ -17,32 +17,19 @@ posify_index, sanitize_index, ) -from geometer.utils.math import NumericalDType +from geometer.utils.ops_dispatch import maybe_dispatch_ufunc_to_dunder_op +from geometer.utils.typing import NDArrayParameters, NumericalDType, Shape, TensorIndex if TYPE_CHECKING: - from typing_extensions import Self, TypeAlias, Unpack + from typing_extensions import Self, Unpack - from geometer.transformation import Transformation + from geometer.transformation import TransformationTensor EQ_TOL_REL = 1e-15 EQ_TOL_ABS = 1e-8 -IntegerIndex1D: TypeAlias = Union[int, np.int_, slice, Sequence[int], Sequence[np.int_], npt.NDArray[np.int_]] -BooleanIndex1D: TypeAlias = Union[bool, np.bool_, slice, Sequence[bool], Sequence[np.bool_], npt.NDArray[np.bool_]] -TensorIndex: TypeAlias = Union[IntegerIndex1D, BooleanIndex1D, tuple[IntegerIndex1D, ...], tuple[BooleanIndex1D, ...]] -Shape: TypeAlias = tuple[int, ...] - -class NDArrayParameters(TypedDict, total=False): - dtype: npt.DTypeLike - copy: bool - order: Literal["K", "A", "C", "F"] - subok: bool - ndim: int - like: npt.ArrayLike - - -class Tensor(Sized, Iterable["Tensor"]): +class Tensor: """Wrapper class around a numpy array that keeps track of covariant and contravariant indices. Covariant indices are the lower indices (subscripts) and contravariant indices are the upper indices (superscripts) @@ -54,7 +41,7 @@ class Tensor(Sized, Iterable["Tensor"]): indices of the array will be covariant indices and all others contravariant indices. By default, all indices are covariant. If the first argument is a tensor then its indices are copied. tensor_rank: If the Tensor object contains multiple tensors, this parameter specifies the rank of the tensors - contained in the collection. By default, only a single tensor is contained in a Tensor object. + contained in it. By default, only a single tensor is contained in a Tensor object. **kwargs: Additional keyword arguments for the constructor of the numpy array as defined in `numpy.array`. Attributes: @@ -67,7 +54,6 @@ class Tensor(Sized, Iterable["Tensor"]): """ array: npt.NDArray[NumericalDType] - _collection_indices: set[int] _covariant_indices: set[int] _contravariant_indices: set[int] @@ -84,7 +70,6 @@ def __init__( if len(args) == 1: if isinstance(args[0], Tensor): self.array = np.array(args[0].array, **kwargs) - self._collection_indices = args[0]._collection_indices self._covariant_indices = args[0]._covariant_indices self._contravariant_indices = args[0]._contravariant_indices return @@ -93,9 +78,6 @@ def __init__( else: self.array = np.array(args, **kwargs) - if not is_numerical_dtype(self.dtype): - raise TypeError(f"The dtype of a Tensor must be a numeric dtype not {self.dtype.name}") - if tensor_rank is None: tensor_rank = self.rank elif tensor_rank < 0: @@ -104,11 +86,9 @@ def __init__( if self.rank < tensor_rank: raise ValueError("tensor_rank must be smaller than or equal to the number of array dimensions.") - n_col = self.rank - tensor_rank - self._collection_indices = set(range(n_col)) - + n_free_indices = self.rank - tensor_rank if covariant is True: - self._covariant_indices = set(range(n_col, self.rank)) + self._covariant_indices = set(range(n_free_indices, self.rank)) elif covariant is False: self._covariant_indices = set() else: @@ -118,11 +98,18 @@ def __init__( raise IndexError(f"Index {idx} out of range [{-tensor_rank}, {tensor_rank})") idx = sanitize_index(idx) # type: ignore[no-untyped-call] idx = posify_index(tensor_rank, idx) # type: ignore[no-untyped-call] - self._covariant_indices.add(n_col + idx) + self._covariant_indices.add(n_free_indices + idx) - self._contravariant_indices = set(range(self.rank)) - self._covariant_indices - self._collection_indices + free_indices = set(range(n_free_indices)) + self._contravariant_indices = set(range(self.rank)) - self._covariant_indices - free_indices - def __apply__(self, transformation: Transformation) -> Self: + self._validate_tensor() + + def _validate_tensor(self) -> None: + if not is_numerical_dtype(self.dtype): + raise TypeError(f"The dtype of a Tensor must be a numeric dtype not {self.dtype.name}") + + def __apply__(self, transformation: TransformationTensor) -> Self: ts = self.tensor_shape edges: list[tuple[Tensor, Tensor]] = [(self, transformation.copy()) for _ in range(ts[0])] if ts[1] > 0: @@ -144,9 +131,10 @@ def dtype(self) -> np.dtype[NumericalDType]: return self.array.dtype @property - def cdim(self) -> int: - """Number of collection indices. For single tensors this is always zero.""" - return len(self._collection_indices) + def free_indices(self) -> int: + """Number of free indices, i.e. indices that are not covariant and not contravariant.""" + cov, con = self.tensor_shape + return self.rank - (cov + con) @property def tensor_shape(self) -> tuple[int, int]: @@ -172,16 +160,18 @@ def tensor_product(self, other: Tensor) -> Tensor: The tensor product. """ + if self.free_indices > 0 or other.free_indices > 0: + raise NotImplementedError("tensor_product is only implemented for tensors without free indices") offset = self.rank covariant = list(self._covariant_indices) + [offset + i for i in other._covariant_indices] contravariant = list(self._contravariant_indices) + [offset + i for i in other._contravariant_indices] - result = np.tensordot(self.array, other.array, 0) + result = np.tensordot(self.array, other.array, 0) # type: ignore[arg-type] result = np.transpose(result, axes=covariant + contravariant) return Tensor(result, covariant=range(len(covariant)), copy=False) def transpose(self, perm: Iterable[int] | None = None) -> Tensor: - """Permute the indices of the tensor. + """Permute the indices of the tensor. Free indices are not permuted. Args: perm: A list of permuted indices or a shorter list representing a permutation in cycle notation. @@ -192,10 +182,11 @@ def transpose(self, perm: Iterable[int] | None = None) -> Tensor: """ if perm is None: - # TODO: don't permute collection indices - perm = reversed(range(self.rank)) + perm = list(range(self.free_indices)) + list(reversed(range(self.free_indices, self.rank))) + else: + perm = list(perm) - perm = list(perm) + # TODO: check that free indices are not permuted if len(perm) < self.rank: a = list(range(self.rank)) @@ -206,19 +197,15 @@ def transpose(self, perm: Iterable[int] | None = None) -> Tensor: covariant_indices = [] contravariant_indices = [] - collection_indices = [] for i, j in enumerate(perm): if j in self._covariant_indices: covariant_indices.append(i) elif j in self._contravariant_indices: contravariant_indices.append(i) - elif j in self._collection_indices: - collection_indices.append(i) result = Tensor(self.array.transpose(perm), copy=False) result._covariant_indices = set(covariant_indices) result._contravariant_indices = set(contravariant_indices) - result._collection_indices = set(collection_indices) return result @@ -233,49 +220,6 @@ def copy(self) -> Self: result.__dict__.update(self.__dict__) return result - def expand_dims(self, axis: int) -> Self: - """Add a new index as collection index. - - Args: - axis: Position in the new shape where the new axis is placed. - - Returns: - The tensor with an additional index. - - """ - # TODO: return tensor - result = self.copy() - result.array = np.expand_dims(self.array, axis) - - axis = sanitize_index(axis) # type: ignore[no-untyped-call] - axis = posify_index(self.rank, axis) # type: ignore[no-untyped-call] - result._collection_indices = {i + 1 if i >= axis else i for i in self._collection_indices} - result._covariant_indices = {i + 1 if i >= axis else i for i in self._covariant_indices} - result._contravariant_indices = {i + 1 if i >= axis else i for i in self._contravariant_indices} - result._collection_indices.add(axis) - - return result - - def squeeze(self, axis: int) -> Tensor: - """Remove an axis of length one. - - Args: - axis: Axis to remove. - - Returns: - The tensor without the removed axis. - - """ - result = self.copy() - result.array = np.squeeze(self.array, axis) - - axis = sanitize_index(axis) # type: ignore[no-untyped-call] - axis = posify_index(self.rank, axis) # type: ignore[no-untyped-call] - result._collection_indices = {i - 1 if i > axis else i for i in self._collection_indices if i != axis} - result._covariant_indices = {i - 1 if i > axis else i for i in self._covariant_indices if i != axis} - result._contravariant_indices = {i - 1 if i > axis else i for i in self._contravariant_indices if i != axis} - return result - def is_zero(self, tol: float = EQ_TOL_ABS) -> npt.NDArray[np.bool_]: """Test whether the tensor is zero with respect to covariant and contravariant indices. @@ -288,13 +232,10 @@ def is_zero(self, tol: float = EQ_TOL_ABS) -> npt.NDArray[np.bool_]: """ axes = tuple(self._covariant_indices) + tuple(self._contravariant_indices) - return np.all(np.isclose(self.array, 0, atol=tol), axis=axes) # type: ignore[no-any-return] + return np.all(np.isclose(self.array, 0, atol=tol), axis=axes) def __repr__(self) -> str: - class_name = self.__class__.__name__ - if self.cdim > 0: - class_name += "Collection" - return f"{class_name}({self.array.tolist()})" + return f"{self.__class__.__name__}({self.array.tolist()})" def _get_index_mapping(self, index: TensorIndex) -> list[int | None]: normalized_index = normalize_index(index, self.shape) # type: ignore[no-untyped-call] @@ -340,14 +281,11 @@ def __getitem__(self, index: TensorIndex) -> Tensor | np.generic: covariant_indices = [] contravariant_indices = [] - collection_indices = [] index_mapping = self._get_index_mapping(index) for new_axis, old_axis in enumerate(index_mapping): - if old_axis is None or old_axis in self._collection_indices: - collection_indices.append(new_axis) - elif old_axis in self._covariant_indices: + if old_axis in self._covariant_indices: covariant_indices.append(new_axis) elif old_axis in self._contravariant_indices: contravariant_indices.append(new_axis) @@ -355,7 +293,6 @@ def __getitem__(self, index: TensorIndex) -> Tensor | np.generic: result_tensor = Tensor(result, copy=False) result_tensor._covariant_indices = set(covariant_indices) result_tensor._contravariant_indices = set(contravariant_indices) - result_tensor._collection_indices = set(collection_indices) return result_tensor @@ -364,23 +301,19 @@ def __setitem__(self, key: TensorIndex, value: Tensor | npt.ArrayLike) -> None: value = value.array self.array[key] = value - def __iter__(self) -> Iterator[Tensor]: - for i in range(len(self)): - yield self[i] - - def __copy__(self) -> Tensor: + def __copy__(self) -> Self: return self.copy() def __mul__(self, other: Tensor | npt.ArrayLike) -> Tensor: if is_numerical_scalar(other): - return Tensor(self.array * other, covariant=self._covariant_indices, copy=False) + return Tensor(self.array * other, covariant=self._covariant_indices, copy=False) # type: ignore[operator] if not isinstance(other, Tensor): other = Tensor(other, copy=False) return TensorDiagram((other, self)).calculate() def __rmul__(self, other: Tensor | npt.ArrayLike) -> Tensor: if is_numerical_scalar(other): - return self * other + return self * other # type: ignore[operator] if not isinstance(other, Tensor): other = Tensor(other, copy=False) return TensorDiagram((self, other)).calculate() @@ -403,13 +336,13 @@ def __pow__(self, power: int, modulo: int | None = None) -> Tensor: def __truediv__(self, other: Tensor | npt.ArrayLike) -> Tensor: if is_numerical_scalar(other): - return Tensor(self.array / other, covariant=self._covariant_indices, copy=False) + return Tensor(self.array / other, covariant=self._covariant_indices, copy=False) # type: ignore[operator] return NotImplemented def __add__(self, other: Tensor | npt.ArrayLike) -> Tensor: if isinstance(other, Tensor): other = other.array - return Tensor(self.array + other, covariant=self._covariant_indices, copy=False) + return Tensor(self.array + other, covariant=self._covariant_indices, copy=False) # type: ignore[operator] def __radd__(self, other: Tensor | npt.ArrayLike) -> Tensor: return self + other @@ -417,7 +350,7 @@ def __radd__(self, other: Tensor | npt.ArrayLike) -> Tensor: def __sub__(self, other: Tensor | npt.ArrayLike) -> Tensor: if isinstance(other, Tensor): other = other.array - return Tensor(self.array - other, covariant=self._covariant_indices, copy=False) + return Tensor(self.array - other, covariant=self._covariant_indices, copy=False) # type: ignore[operator] def __rsub__(self, other: Tensor | npt.ArrayLike) -> Tensor: return -self + other @@ -433,27 +366,147 @@ def __eq__(self, other: object) -> bool: except TypeError: return NotImplemented - def __len__(self) -> int: - return len(self.array) - def __array__(self, dtype: npt.DTypeLike | None = None) -> np.ndarray: if dtype and dtype != self.dtype: return self.array.astype(dtype) return self.array - def __array_ufunc__(self, ufunc: np.ufunc, method: str, *inputs: Any, **kwargs: Any): + def __array_ufunc__( + self, + ufunc: np.ufunc, + method: Literal["__call__", "reduce", "reduceat", "accumulate", "outer", "inner"], + *inputs: Any, + **kwargs: Any, + ) -> Any: + return maybe_dispatch_ufunc_to_dunder_op(self, ufunc, method, *inputs, **kwargs) + + def __array_function__( + self, func: Callable[..., Any], types: Iterable[type], args: tuple[Any, ...], kwargs: dict[str, Any] + ) -> Any: return NotImplemented - def __array_function__(self, func, types, args, kwargs): - return NotImplemented + +class BoundTensor(Tensor): + """A tensor without free indices.""" + + def _validate_tensor(self) -> None: + super()._validate_tensor() + if self.free_indices != 0: + raise IncompatibleShapeError( + f"Only tensors without free indices can be used in a {self.__class__.__name__}" + ) + + +T = TypeVar("T", bound=Tensor) + + +class TensorCollection(Tensor, Generic[T], Sized, Iterable[T]): + """A collection of tensors.""" + + _element_class: ClassVar[type[Tensor]] = Tensor + + def __init__( + self, + *args: Tensor | npt.ArrayLike, + covariant: bool | Iterable[int] = True, + tensor_rank: int | None = 1, + **kwargs: Unpack[NDArrayParameters], + ) -> None: + super().__init__(*args, covariant=covariant, tensor_rank=tensor_rank, **kwargs) + + def _validate_tensor(self) -> None: + super()._validate_tensor() + if self.free_indices == 0: + raise IncompatibleShapeError( + f"Tensor has no free indices and cannot be used in a {self.__class__.__name__}." + ) + + @classmethod + def from_tensor(cls, tensor: Tensor, **kwargs: Unpack[NDArrayParameters]) -> Self | T: + """Construct an object from another tensor. If the tensor has no free indices, an object of type T is returned. + + By default the array in the tensor is not copied. + + Args: + tensor: A tensor to use for the new tensor. + **kwargs: Additional keyword arguments for the constructor of the numpy array as defined in `numpy.array`. + + Returns: + A new tensor. + + """ + kwargs.setdefault("copy", False) + if tensor.free_indices > 0: + return cls(tensor, **kwargs) + else: + return cls._element_class(tensor, **kwargs) + + @classmethod + def from_array(cls, array: npt.ArrayLike, **kwargs: Unpack[NDArrayParameters]) -> Self | T: + """Try to construct a new collection from an array. If the rank is too low, an object of type T is returned. + + By default the array is not copied. + + Args: + array: A numpy array to use for the new tensor. + **kwargs: Additional keyword arguments for the constructor of the numpy array as defined in `numpy.array`. + + Returns: + A new tensor. + + """ + kwargs.setdefault("copy", False) + try: + return cls(array, **kwargs) + except IncompatibleShapeError: + return cls._element_class(array, **kwargs) + + def expand_dims(self, axis: int) -> Self: + """Add a new index as a free index. + + Args: + axis: Position in the new shape where the new axis is placed. + + Returns: + The tensor collection with an additional free index. + + """ + axis = sanitize_index(axis) # type: ignore[no-untyped-call] + axis = posify_index(self.rank + 1, axis) # type: ignore[no-untyped-call] + if axis > self.free_indices: + raise ValueError("Free indices can only be inserted at the beginning.") + + result = self.copy() + result.array = np.expand_dims(self.array, axis) + result._covariant_indices = {i + 1 if i >= axis else i for i in self._covariant_indices} + result._contravariant_indices = {i + 1 if i >= axis else i for i in self._contravariant_indices} + return result @property def size(self) -> npt.NDArray[np.int_]: - """The number of tensors in the tensor, i.e. the product of the size of all collection axes.""" - return np.prod([self.shape[i] for i in self._collection_indices], dtype=int) + """The number of tensors in the tensor collection, i.e. the product of the size of all collection axes.""" + return np.prod([self.shape[i] for i in range(self.free_indices)], dtype=int) + + def __getitem__(self, index: TensorIndex) -> Tensor | np.generic: + result = super().__getitem__(index) + + if not isinstance(result, Tensor): + return result + + if result.free_indices > 0: + return TensorCollection(result, copy=False) + + return self._element_class(result, copy=False) + + def __len__(self) -> int: + return len(self.array) + + def __iter__(self) -> Iterator[T]: + for i in range(len(self)): + yield self[i] -class LeviCivitaTensor(Tensor): +class LeviCivitaTensor(BoundTensor): r"""This class can be used to construct a tensor representing the Levi-Civita symbol. The Levi-Civita symbol is also called :math:`\varepsilon`-Tensor and is defined as follows: @@ -494,7 +547,7 @@ def __init__(self, size: int, covariant: bool = True) -> None: super().__init__(array, covariant=bool(covariant), copy=False) -class KroneckerDelta(Tensor): +class KroneckerDelta(BoundTensor): r"""This class can be used to construct a (p, p)-tensor representing the Kronecker delta tensor. The following generalized definition of the Kronecker delta is used: @@ -526,7 +579,7 @@ def __init__(self, n: int, p: int = 1) -> None: array = self._cache[(p, n)] elif p == n: e = LeviCivitaTensor(n) - array = np.tensordot(e.array, e.array, 0) + array = np.tensordot(e.array, e.array, 0) # type: ignore[arg-type] self._cache[(p, n)] = array else: @@ -567,7 +620,7 @@ class TensorDiagram: def __init__(self, *edges: tuple[Tensor, Tensor]) -> None: self._nodes: list[Tensor] = [] - self._free_indices: list[tuple[list[int], list[int]]] = [] + self._unused_indices: list[tuple[list[int], list[int]]] = [] self._node_positions: list[int] = [] self._contraction_list: list[tuple[int, int, int, int]] = [] self._index_count: int = 0 @@ -589,7 +642,7 @@ def add_node(self, node: Tensor) -> tuple[list[int], list[int]]: self._node_positions.append(self._index_count) self._index_count += node.rank ind = (list(node._covariant_indices), list(node._contravariant_indices)) - self._free_indices.append(ind) + self._unused_indices.append(ind) return ind def add_edge(self, source: Tensor, target: Tensor) -> None: @@ -600,7 +653,7 @@ def add_edge(self, source: Tensor, target: Tensor) -> None: target: The target tensor of the edge in the diagram. Raises: - TensorComputationError: If the tensors have no free indices or the dimensions do not match. + TensorComputationError: If the tensors have no unused indices or the dimensions do not match. """ @@ -609,10 +662,10 @@ def add_edge(self, source: Tensor, target: Tensor) -> None: for index, node in enumerate(self._nodes): if node is source: source_index = index - free_source = self._free_indices[index][0] + free_source = self._unused_indices[index][0] if node is target: target_index = index - free_target = self._free_indices[index][1] + free_target = self._unused_indices[index][1] if source_index is not None and target_index is not None: break @@ -629,7 +682,7 @@ def add_edge(self, source: Tensor, target: Tensor) -> None: free_target = self.add_node(target)[1] if len(free_source) == 0 or len(free_target) == 0: - raise TensorComputationError("Could not add the edge because no free indices are left.") + raise TensorComputationError("Could not add the edge because no indices are left.") # Third step: Pick some free indices i = free_source.pop(0) @@ -659,13 +712,13 @@ def calculate(self) -> Tensor: # Split the indices and insert the arrays args: list[npt.ArrayLike] = [] result_indices: tuple[list[int], list[int], list[int]] = ([], [], []) - for i, (node, ind, offset) in enumerate(zip(self._nodes, self._free_indices, self._node_positions)): - if len(node._collection_indices) > 0: - col_ind = sorted(node._collection_indices, reverse=True) - for j, k in enumerate(col_ind[: len(result_indices[0])]): + for i, (node, ind, offset) in enumerate(zip(self._nodes, self._unused_indices, self._node_positions)): + if node.free_indices > 0: + free_ind = list(reversed(range(node.free_indices))) + for j, k in enumerate(free_ind[: len(result_indices[0])]): indices[offset + k] = result_indices[0][-j - 1] - if len(node._collection_indices) > len(result_indices[0]): - for k in col_ind[len(result_indices[0]) :]: + if node.free_indices > len(result_indices[0]): + for k in free_ind[len(result_indices[0]) :]: result_indices[0].insert(0, offset + k) result_indices[1].extend(offset + x for x in ind[0]) result_indices[2].extend(offset + x for x in ind[1]) @@ -673,17 +726,17 @@ def calculate(self) -> Tensor: s = slice(offset, self._node_positions[i + 1] if i + 1 < len(self._node_positions) else None) args.append(indices[s]) - result = np.einsum(*args, result_indices[0] + result_indices[1] + result_indices[2]) + result = np.einsum(*args, result_indices[0] + result_indices[1] + result_indices[2]) # type: ignore[arg-type] - n_col = len(result_indices[0]) + n_free = len(result_indices[0]) n_cov = len(result_indices[1]) - return Tensor(result, covariant=range(n_cov), tensor_rank=result.ndim - n_col, copy=False) + return Tensor(result, covariant=range(n_cov), tensor_rank=result.ndim - n_free, copy=False) def copy(self) -> TensorDiagram: result = TensorDiagram() result._nodes = self._nodes.copy() - result._free_indices = self._free_indices.copy() + result._unused_indices = self._unused_indices.copy() result._node_positions = self._node_positions.copy() result._contraction_list = self._contraction_list.copy() result._index_count = self._index_count @@ -715,7 +768,7 @@ def __eq__(self, other: object) -> bool: if not isinstance(other, Tensor): try: - other = Tensor(other) + other = Tensor(other) # type: ignore[arg-type] except TypeError: return NotImplemented diff --git a/geometer/curve.py b/geometer/curve.py index 03df24f..8c65d32 100644 --- a/geometer/curve.py +++ b/geometer/curve.py @@ -1,6 +1,7 @@ from __future__ import annotations import math +from abc import ABC from itertools import combinations from typing import TYPE_CHECKING, Union, cast @@ -8,17 +9,42 @@ import numpy.typing as npt from numpy.lib.scimath import sqrt as csqrt -from geometer.base import EQ_TOL_ABS, EQ_TOL_REL, NDArrayParameters, ProjectiveTensor, Tensor, TensorDiagram +from geometer.base import ( + EQ_TOL_ABS, + EQ_TOL_REL, + BoundTensor, + ProjectiveTensor, + Tensor, + TensorCollection, + TensorDiagram, +) from geometer.exceptions import IncidenceError, NotReducible -from geometer.point import I, J, Line, Plane, Point, Subspace, infty_plane, join +from geometer.point import ( + I, + J, + Line, + LineCollection, + LineTensor, + Plane, + PlaneCollection, + PlaneTensor, + Point, + PointCollection, + PointTensor, + SubspaceTensor, + infty_plane, + join, +) from geometer.transformation import rotation, translation from geometer.utils import adjugate, det, hat_matrix, inv, is_multiple, matmul, matvec, outer, roots if TYPE_CHECKING: from typing_extensions import Unpack + from geometer.utils.typing import NDArrayParameters -class Quadric(ProjectiveTensor): + +class QuadricTensor(ProjectiveTensor, ABC): r"""Represents a quadric, i.e. the zero set of a polynomial of degree 2, in any dimension. The quadric is defined by a symmetric matrix of size :math:`n+1` where :math:`n` is the dimension of the projective @@ -59,19 +85,19 @@ def __init__( super().__init__(matrix, tensor_rank=2, **kwargs) def __add__(self, other: Tensor | npt.ArrayLike) -> Tensor: - if not isinstance(other, Point): + if not isinstance(other, PointTensor): return super().__add__(other) return translation(other).apply(self) def __sub__(self, other: Tensor | npt.ArrayLike) -> Tensor: - if not isinstance(other, Point): + if not isinstance(other, PointTensor): return super().__add__(other) return translation(-other).apply(self) @classmethod - def from_planes(cls, e: Plane, f: Plane) -> Quadric: + def from_planes(cls, e: PlaneTensor, f: PlaneTensor) -> QuadricTensor: """Construct a degenerate quadric from two hyperplanes. Args: @@ -83,9 +109,9 @@ def from_planes(cls, e: Plane, f: Plane) -> Quadric: """ m = outer(e.array, f.array) m += m.T - return Quadric(m, normalize_matrix=True) + return cls(m, normalize_matrix=True) - def tangent(self, at: Point) -> Plane: + def tangent(self, at: PointTensor) -> PlaneTensor: """Returns the hyperplane defining the tangent space at a given point. Args: @@ -95,9 +121,9 @@ def tangent(self, at: Point) -> Plane: The tangent plane at the given point. """ - return Plane(matvec(self.array, at.array), copy=False) + return PlaneCollection.from_array(matvec(self.array, at.array)) - def is_tangent(self, plane: Subspace) -> npt.NDArray[np.bool_]: + def is_tangent(self, plane: SubspaceTensor) -> npt.NDArray[np.bool_]: """Tests if a given hyperplane is tangent to the quadric. Args: @@ -109,7 +135,7 @@ def is_tangent(self, plane: Subspace) -> npt.NDArray[np.bool_]: """ return self.dual.contains(plane) - def contains(self, other: Point | Subspace, tol: float = EQ_TOL_ABS) -> npt.NDArray[np.bool_]: + def contains(self, other: PointTensor | SubspaceTensor, tol: float = EQ_TOL_ABS) -> npt.NDArray[np.bool_]: """Tests if a given point lies on the quadric. Args: @@ -132,7 +158,7 @@ def is_degenerate(self) -> npt.NDArray[np.bool_]: return np.isclose(det(self.array), 0, atol=EQ_TOL_ABS) @property - def components(self) -> list[Point | Line | Plane]: + def components(self) -> list[PointTensor | LineTensor | PlaneTensor]: """The components of a degenerate quadric.""" # Algorithm adapted from Perspectives on Projective Geometry, Section 11.1 n = self.shape[-1] @@ -150,7 +176,7 @@ def components(self) -> list[Point | Line | Plane]: [np.delete(np.delete(ind, i, axis=1), i, axis=2) for i in combinations(range(n), n - 2)], axis=1 ) minors = det(self.array[..., ind[0], ind[1]]) - p = csqrt(-minors) + p = csqrt(-minors) # type: ignore[arg-type] # use the skew symmetric matrix m to get a matrix of rank 1 defining the same quadric m = hat_matrix(p) @@ -168,12 +194,12 @@ def components(self) -> list[Point | Line | Plane]: p, q = np.real_if_close(p), np.real_if_close(q) if self.is_dual: - return [Point(p, copy=False), Point(q, copy=False)] + return [PointCollection.from_array(p), PointCollection.from_array(q)] elif n == 3: - return [Line(p, copy=False), Line(q, copy=False)] - return [Plane(p, copy=False), Plane(q, copy=False)] + return [LineCollection.from_array(p), LineCollection.from_array(q)] + return [PlaneCollection.from_array(p), PlaneCollection.from_array(q)] - def intersect(self, other: Line) -> list[Point] | list[Line]: + def intersect(self, other: LineTensor) -> list[PointTensor] | list[LineTensor]: """Calculates points of intersection of a line with the quadric. This method also returns complex points of intersection, even if the quadric and the line do not intersect in @@ -189,9 +215,6 @@ def intersect(self, other: Line) -> list[Point] | list[Line]: - J. Richter-Gebert: Perspectives on Projective Geometry, Section 11.3 """ - if not isinstance(other, Line): - raise TypeError(f"Expected Line but got {type(other)}") - reducible: bool | np.bool_ = np.all(self.is_degenerate) if reducible: try: @@ -203,39 +226,47 @@ def intersect(self, other: Line) -> list[Point] | list[Line]: if self.dim > 2: arr = other.array.reshape(other.shape[: -other.tensor_shape[1]] + (-1, self.dim + 1)) - if other.cdim == 0: + if isinstance(other, Line): i = arr.nonzero()[0][0] m = Plane(arr[i], copy=False).basis_matrix else: i = np.any(arr, axis=-1).argmax(-1) - m = Plane(arr[(*tuple(np.indices(i.shape)), i)], copy=False).basis_matrix + m = PlaneCollection(arr[(*tuple(np.indices(i.shape)), i)], copy=False).basis_matrix line = other._matrix_transform(m) - projected_quadric = Quadric(matmul(matmul(m, self.array), m, transpose_b=True), copy=False) + projected_quadric = QuadricCollection.from_array(matmul(matmul(m, self.array), m, transpose_b=True)) return [ - Point(np.squeeze(matmul(np.expand_dims(point.array, -2), m), -2), copy=False) + PointCollection.from_array(np.squeeze(matmul(np.expand_dims(point.array, -2), m), -2)) for point in projected_quadric.intersect(line) ] else: m = hat_matrix(other.array) b = matmul(matmul(m, self.array, transpose_a=True), m) - p, q = Quadric(b, is_dual=not self.is_dual, copy=False).components + p, q = QuadricCollection.from_array(b, is_dual=not self.is_dual).components else: if self.is_dual: - e, f = cast(Point, e), cast(Point, f) + e, f = cast(PointTensor, e), cast(PointTensor, f) p, q = e.join(other), f.join(other) else: - e, f = cast(Union[Line, Plane], e), cast(Union[Line, Plane], f) + e, f = cast(Union[LineTensor, PlaneTensor], e), cast(Union[LineTensor, PlaneTensor], f) p, q = e.meet(other), f.meet(other) - if p.cdim == 0 and p == q: + if p.free_indices == 0 and p == q: return [p] return [p, q] @property - def dual(self) -> Quadric: + def dual(self) -> QuadricTensor: """The dual quadric.""" - return Quadric(inv(self.array), is_dual=not self.is_dual, copy=False) + return type(self)(inv(self.array), is_dual=not self.is_dual, copy=False) + + +class Quadric(QuadricTensor, BoundTensor): + pass + + +class QuadricCollection(QuadricTensor, TensorCollection[Quadric]): + _element_class = Quadric class Conic(Quadric): @@ -339,7 +370,7 @@ def from_foci(cls, f1: Point, f2: Point, bound: Point) -> Conic: p1, p2 = Point(t1.array, copy=False), Point(t2.array, copy=False) p3, p4 = Point(t3.array, copy=False), Point(t4.array, copy=False) c = cls.from_tangent(Line(bound.array, copy=False), p1, p2, p3, p4) - return Conic(np.linalg.inv(c.array), copy=False) + return cls(np.linalg.inv(c.array), copy=False) @classmethod def from_crossratio(cls, cr: float, a: Point, b: Point, c: Point, d: Point) -> Conic: @@ -467,7 +498,13 @@ class Ellipse(Conic): """ - def __init__(self, center: Point = Point(0, 0), hradius: float = 1, vradius: float = 1, **kwargs) -> None: + def __init__( + self, + center: Point = Point(0, 0), + hradius: float = 1, + vradius: float = 1, + **kwargs: Unpack[NDArrayParameters], + ) -> None: if hradius == vradius == 0: raise ValueError("hradius and vradius can not both be zero.") @@ -497,7 +534,7 @@ class Circle(Ellipse): """ - def __init__(self, center: Point = Point(0, 0), radius: float = 1, **kwargs) -> None: + def __init__(self, center: Point = Point(0, 0), radius: float = 1, **kwargs: Unpack[NDArrayParameters]) -> None: if radius <= 0: raise ValueError(f"radius must be greater than 0, but is {radius}") super().__init__(center, radius, radius, **kwargs) diff --git a/geometer/exceptions.py b/geometer/exceptions.py index b447da7..af0a238 100644 --- a/geometer/exceptions.py +++ b/geometer/exceptions.py @@ -45,3 +45,7 @@ def __init__(self, message: str, dependent_values: npt.ArrayLike = True) -> None class NotReducible(GeometryException, ValueError): """The given geometric object is not reducible.""" + + +class IncompatibleShapeError(ValueError): + """The given tensor has a shape that is not compatible.""" diff --git a/geometer/operators.py b/geometer/operators.py index 575075c..59b597c 100644 --- a/geometer/operators.py +++ b/geometer/operators.py @@ -1,6 +1,6 @@ from __future__ import annotations -from typing import cast +from typing import TYPE_CHECKING, cast import numpy as np import numpy.typing as npt @@ -8,16 +8,34 @@ from geometer.base import EQ_TOL_ABS, EQ_TOL_REL, LeviCivitaTensor, TensorDiagram from geometer.curve import absolute_conic from geometer.exceptions import NotCollinear, NotConcurrent -from geometer.point import I, J, Line, Plane, Point, infty, infty_plane, join, meet +from geometer.point import ( + I, + J, + LineCollection, + LineTensor, + PlaneCollection, + PlaneTensor, + Point, + PointCollection, + PointTensor, + SubspaceTensor, + infty, + infty_plane, + join, + meet, +) from geometer.utils import det, matvec, orth +if TYPE_CHECKING: + from geometer.shapes import PolytopeTensor + def crossratio( - a: Point | Line | Plane, - b: Point | Line | Plane, - c: Point | Line | Plane, - d: Point | Line | Plane, - from_point: Point | None = None, + a: PointTensor | SubspaceTensor, + b: PointTensor | SubspaceTensor, + c: PointTensor | SubspaceTensor, + d: PointTensor | SubspaceTensor, + from_point: PointTensor | None = None, ) -> np.ndarray: """Calculates the cross ratio of points or lines. @@ -35,19 +53,29 @@ def crossratio( """ if a == b: - return np.ones(a.shape[: a.cdim]) - - if isinstance(a, Line): + return np.ones(a.shape[: a.free_indices]) + + if ( + isinstance(a, LineTensor) + and isinstance(b, LineTensor) + and isinstance(c, LineTensor) + and isinstance(d, LineTensor) + ): if not np.all(is_concurrent(a, b, c, d)): raise NotConcurrent("The lines are not concurrent: " + str([a, b, c, d])) from_point = a.meet(b) a, b, c, d = a.base_point, b.base_point, c.base_point, d.base_point - elif isinstance(a, Plane): + elif ( + isinstance(a, PlaneTensor) + and isinstance(b, PlaneTensor) + and isinstance(c, PlaneTensor) + and isinstance(d, PlaneTensor) + ): l = a.meet(b) - l = cast(Line, l) - e = Plane(l.direction.array, copy=False) + l = cast(LineTensor, l) + e = PlaneCollection.from_array(l.direction.array) a, b, c, d = e.meet(a), e.meet(b), e.meet(c), e.meet(d) m = e.basis_matrix p = e.meet(l) @@ -57,6 +85,14 @@ def crossratio( c = (p + c.direction)._matrix_transform(m) d = (p + d.direction)._matrix_transform(m) + elif not ( + isinstance(a, PointTensor) + and isinstance(b, PointTensor) + and isinstance(c, PointTensor) + and isinstance(d, PointTensor) + ): + raise TypeError(f"Unsupported combination of types: a: {type(a)}, b: {type(b)}, c: {type(c)}, d: {type(d)}") + if a.dim > 2 or (from_point is None and a.dim == 2): if not np.all(is_collinear(a, b, c, d)): raise NotCollinear("The points are not collinear: " + str([a, b, c, d])) @@ -84,7 +120,7 @@ def crossratio( return ac * bd / (ad * bc) -def harmonic_set(a: Point, b: Point, c: Point) -> Point: +def harmonic_set(a: PointTensor, b: PointTensor, c: PointTensor) -> PointTensor: """Constructs a fourth point that forms a harmonic set with the given points. The three given points must be collinear. @@ -122,7 +158,7 @@ def harmonic_set(a: Point, b: Point, c: Point) -> Point: return result -def angle(*args: Point | Line | Plane) -> npt.NDArray[np.float_]: +def angle(*args: PointTensor | LineTensor | PlaneTensor) -> npt.NDArray[np.float_]: r"""Calculates the (oriented) angle between given points, lines or planes. The function uses the Laguerre formula to calculate angles in two or three-dimensional projective space @@ -159,25 +195,26 @@ def angle(*args: Point | Line | Plane) -> npt.NDArray[np.float_]: elif len(args) == 2: x, y = args - if isinstance(x, Plane) and isinstance(y, Plane): + if isinstance(x, PlaneTensor) and isinstance(y, PlaneTensor): l = x.meet(y) p = l.meet(infty_plane) - polar = Line(p.array[..., :-1], copy=False) + polar = LineCollection.from_array(p.array[..., :-1]) tangent_points = absolute_conic.intersect(polar) tangent_points = [ - Point(np.append(p.array, np.zeros(p.shape[:-1] + (1,)), axis=-1), copy=False) for p in tangent_points + PointCollection.from_array(np.append(p.array, np.zeros(p.shape[:-1] + (1,)), axis=-1)) + for p in tangent_points ] i = l.join(p.join(tangent_points[0])) j = l.join(p.join(tangent_points[1])) return 1 / 2j * np.log(crossratio(x, y, i, j)) - if isinstance(x, Line) and isinstance(y, Line): + if isinstance(x, LineTensor) and isinstance(y, LineTensor): a = x.meet(y) else: a = Point(*(x.dim * [0])) - if isinstance(x, Point): + if isinstance(x, PointTensor): x = a.join(x) - if isinstance(y, Point): + if isinstance(y, PointTensor): y = a.join(y) if a.dim > 2: @@ -195,7 +232,7 @@ def angle(*args: Point | Line | Plane) -> npt.NDArray[np.float_]: return np.real(1 / 2j * np.log(crossratio(b, c, I, J, a))) -def angle_bisectors(l: Line, m: Line) -> tuple[Line, Line]: +def angle_bisectors(l: LineTensor, m: LineTensor) -> tuple[LineTensor, LineTensor]: """Constructs the angle bisectors of two given lines. Args: @@ -225,8 +262,8 @@ def angle_bisectors(l: Line, m: Line) -> tuple[Line, Line]: a, b = np.expand_dims(a, -1), np.expand_dims(b, -1) r, s = a * i + b * j, a * i - b * j - r = Point(r, copy=False) - s = Point(s, copy=False) + r = PointCollection.from_array(r) + s = PointCollection.from_array(s) if o.dim > 2: r = r._matrix_transform(np.swapaxes(basis, -1, -2)) @@ -235,7 +272,30 @@ def angle_bisectors(l: Line, m: Line) -> tuple[Line, Line]: return join(o, r), join(o, s) -def dist(p: Point | Line | Plane, q: Point | Line | Plane) -> npt.NDArray[np.float_]: +def _point_dist(p: PointTensor, q: PointTensor) -> npt.NDArray[np.float_]: + if p.dim > 2: + x = np.stack(np.broadcast_arrays(p.normalized_array, q.normalized_array), axis=-2) + z = x[..., -1] + + m = orth(np.swapaxes(x, -1, -2), 2) + x = np.matmul(x, m) + x = np.concatenate([x, np.expand_dims(z, -1)], axis=-1) + p, q, i, j = np.broadcast_arrays(x[..., 0, :], x[..., 1, :], I.array, J.array) + else: + p, q, i, j = np.broadcast_arrays(p.array, q.array, I.array, J.array) + + pqi = det(np.stack([p, q, i], axis=-2)) + pqj = det(np.stack([p, q, j], axis=-2)) + pij = det(np.stack([p, i, j], axis=-2)) + qij = det(np.stack([q, i, j], axis=-2)) + + with np.errstate(divide="ignore", invalid="ignore"): + return 4 * np.abs(np.sqrt(pqi * pqj) / (pij * qij)) + + +def dist( + p: PointTensor | SubspaceTensor | PolytopeTensor, q: PointTensor | SubspaceTensor | PolytopeTensor +) -> npt.NDArray[np.float_]: r"""Calculates the (euclidean) distance between two objects. Instead of the usual formula for the euclidean distance this function uses @@ -255,65 +315,47 @@ def dist(p: Point | Line | Plane, q: Point | Line | Plane) -> npt.NDArray[np.flo """ if p == q: - return np.zeros(p.shape[: p.cdim]) + return np.zeros(p.shape[: p.free_indices]) - subspace_types = (Line, Plane) - - if isinstance(p, subspace_types) and isinstance(q, Point): + if isinstance(p, PointTensor) and isinstance(q, PointTensor): + return _point_dist(p, q) + if isinstance(p, SubspaceTensor) and isinstance(q, PointTensor): return dist(p.project(q), q) - if isinstance(p, Point) and isinstance(q, subspace_types): + if isinstance(p, PointTensor) and isinstance(q, SubspaceTensor): return dist(q.project(p), p) - if isinstance(p, Plane) and isinstance(q, subspace_types): - return dist(p, Point(q.basis_matrix[0, :], copy=False)) - if isinstance(q, Plane) and isinstance(p, Line): - return dist(q, p.base_point) + if isinstance(p, SubspaceTensor) and isinstance(q, PlaneTensor): + return dist(q, p) + if isinstance(p, PlaneTensor) and isinstance(q, LineTensor): + return dist(p, q.base_point) + if isinstance(p, PlaneTensor) and isinstance(q, SubspaceTensor): + return dist(p, PointCollection.from_array(q.basis_matrix[0, :])) - from geometer.shapes import Polygon, Polyhedron, Segment + from geometer.shapes import PolygonTensor, Polyhedron, SegmentTensor - if isinstance(p, Point) and isinstance(q, Polygon): + if isinstance(p, PointTensor) and isinstance(q, SegmentTensor): + return dist(q, p) + if isinstance(p, SegmentTensor) and isinstance(q, PointTensor): + result = np.min([dist(v, q) for v in p.vertices], axis=0) + r = p._line.project(q) + return np.where(p.contains(r), dist(r, q), result) + if isinstance(p, PointTensor) and isinstance(q, PolygonTensor): return dist(q, p) - if isinstance(p, Polygon) and isinstance(q, Point): + if isinstance(p, PolygonTensor) and isinstance(q, PointTensor): result = np.min([dist(e, q) for e in p.edges], axis=0) if p.dim > 2: r = p._plane.project(q) return np.where(p.contains(r), dist(r, q), result) return result - if isinstance(p, Point) and isinstance(q, Polyhedron): + if isinstance(p, PointTensor) and isinstance(q, Polyhedron): return dist(q, p) - if isinstance(p, Polyhedron) and isinstance(q, Point): + if isinstance(p, Polyhedron) and isinstance(q, PointTensor): return np.min([dist(f, q) for f in p.faces], axis=0) - if isinstance(p, Point) and isinstance(q, Segment): - return dist(q, p) - if isinstance(p, Segment) and isinstance(q, Point): - result = np.min([dist(v, q) for v in p.vertices], axis=0) - r = p._line.project(q) - return np.where(p.contains(r), dist(r, q), result) - - if not isinstance(p, Point) or not isinstance(q, Point): - raise TypeError(f"Unsupported types {type(p)} and {type(q)}.") - - if p.dim > 2: - x = np.stack(np.broadcast_arrays(p.normalized_array, q.normalized_array), axis=-2) - z = x[..., -1] - - m = orth(np.swapaxes(x, -1, -2), 2) - x = np.matmul(x, m) - x = np.concatenate([x, np.expand_dims(z, -1)], axis=-1) - p, q, i, j = np.broadcast_arrays(x[..., 0, :], x[..., 1, :], I.array, J.array) - else: - p, q, i, j = np.broadcast_arrays(p.array, q.array, I.array, J.array) - pqi = det(np.stack([p, q, i], axis=-2)) - pqj = det(np.stack([p, q, j], axis=-2)) - pij = det(np.stack([p, i, j], axis=-2)) - qij = det(np.stack([q, i, j], axis=-2)) - - with np.errstate(divide="ignore", invalid="ignore"): - return 4 * np.abs(np.sqrt(pqi * pqj) / (pij * qij)) + raise TypeError(f"Unsupported types {type(p)} and {type(q)}.") def is_cocircular( - a: Point, b: Point, c: Point, d: Point, rtol: float = EQ_TOL_REL, atol: float = EQ_TOL_ABS + a: PointTensor, b: PointTensor, c: PointTensor, d: PointTensor, rtol: float = EQ_TOL_REL, atol: float = EQ_TOL_ABS ) -> npt.NDArray[np.bool_]: """Tests whether four points lie on a circle. @@ -343,7 +385,7 @@ def is_cocircular( def is_perpendicular( - l: Line | Plane, m: Line | Plane, rtol: float = EQ_TOL_REL, atol: float = EQ_TOL_ABS + l: LineTensor | PlaneTensor, m: LineTensor | PlaneTensor, rtol: float = EQ_TOL_REL, atol: float = EQ_TOL_ABS ) -> npt.NDArray[np.bool_]: """Tests whether two lines/planes are perpendicular. @@ -360,19 +402,20 @@ def is_perpendicular( L = l.meet(infty) M = m.meet(infty) - elif isinstance(l, Line) and isinstance(m, Line): + elif isinstance(l, LineTensor) and isinstance(m, LineTensor): e = join(l, m) basis = e.basis_matrix L = l.meet(infty_plane)._matrix_transform(basis) M = m.meet(infty_plane)._matrix_transform(basis) - elif isinstance(l, Plane) and isinstance(m, Plane): + elif isinstance(l, PlaneTensor) and isinstance(m, PlaneTensor): x = l.meet(m) p = x.meet(infty_plane) - polar = Line(p.array[..., :-1], copy=False) + polar = LineCollection.from_array(p.array[..., :-1]) tangent_points = absolute_conic.intersect(polar) tangent_points = [ - Point(np.append(p.array, np.zeros(p.shape[:-1] + (1,)), axis=-1), copy=False) for p in tangent_points + PointCollection.from_array(np.append(p.array, np.zeros(p.shape[:-1] + (1,)), axis=-1)) + for p in tangent_points ] i = x.join(p.join(tangent_points[0])) j = x.join(p.join(tangent_points[1])) @@ -384,7 +427,7 @@ def is_perpendicular( return np.isclose(crossratio(L, M, I, J, Point(1, 1)), -1, rtol, atol) -def is_coplanar(*args: Point | Line, tol: float = EQ_TOL_ABS) -> npt.NDArray[np.bool_]: +def is_coplanar(*args: PointTensor | LineTensor, tol: float = EQ_TOL_ABS) -> npt.NDArray[np.bool_]: """Tests whether the given points or lines are collinear, coplanar or concurrent. Works in any dimension. Due to line point duality this function has dual versions :obj:`is_collinear` and :obj:`is_concurrent`. diff --git a/geometer/point.py b/geometer/point.py index 85e8bbf..71fac34 100644 --- a/geometer/point.py +++ b/geometer/point.py @@ -1,6 +1,7 @@ from __future__ import annotations -from typing import TYPE_CHECKING, cast, overload +from abc import ABC, abstractmethod +from typing import TYPE_CHECKING, TypeVar, cast, overload import numpy as np import numpy.typing as npt @@ -8,69 +9,68 @@ if TYPE_CHECKING: from typing_extensions import Literal, Unpack + from geometer.utils.typing import NDArrayParameters, TensorIndex + from geometer.base import ( EQ_TOL_ABS, + EQ_TOL_REL, + BoundTensor, LeviCivitaTensor, - NDArrayParameters, ProjectiveTensor, Tensor, + TensorCollection, TensorDiagram, - TensorIndex, ) from geometer.exceptions import GeometryException, LinearDependenceError, NotCoplanar -from geometer.utils import is_numerical_scalar, matmul, matvec, null_space - - -@overload -def _join_meet_duality( - *args: Line, intersect_lines: Literal[True], check_dependence: bool = True, normalize_result: bool = True -) -> Point: - ... +from geometer.utils import is_multiple, is_numerical_scalar, matmul, matvec, null_space @overload def _join_meet_duality( - *args: Unpack[tuple[Subspace, Line]], + *args: Unpack[tuple[SubspaceTensor, LineTensor]], intersect_lines: Literal[True], check_dependence: bool = True, normalize_result: bool = True, -) -> Point: +) -> PointTensor: ... @overload def _join_meet_duality( - *args: Unpack[tuple[Line, Subspace]], + *args: Unpack[tuple[LineTensor, SubspaceTensor]], intersect_lines: Literal[True], check_dependence: bool = True, normalize_result: bool = True, -) -> Point: +) -> PointTensor: ... @overload def _join_meet_duality( - *args: Unpack[tuple[Point, Point]], + *args: Unpack[tuple[PointTensor, PointTensor]], intersect_lines: Literal[True], check_dependence: bool = True, normalize_result: bool = True, -) -> Line: +) -> LineTensor: ... @overload def _join_meet_duality( - *args: Point | Subspace, + *args: PointTensor | SubspaceTensor, intersect_lines: Literal[False], check_dependence: bool = True, normalize_result: bool = True, -) -> Subspace: +) -> SubspaceTensor: ... def _join_meet_duality( - *args: Point | Subspace, intersect_lines: bool = True, check_dependence: bool = True, normalize_result: bool = True -) -> Point | Subspace: + *args: PointTensor | SubspaceTensor, + intersect_lines: bool = True, + check_dependence: bool = True, + normalize_result: bool = True, +) -> PointTensor | SubspaceTensor: if len(args) < 2: raise ValueError(f"Expected at least 2 arguments, got {len(args)}.") @@ -87,14 +87,19 @@ def _join_meet_duality( # two lines/planes elif len(args) == 2: a, b = args - if isinstance(a, Line) and isinstance(b, Plane) or isinstance(b, Line) and isinstance(a, Plane): + if ( + isinstance(a, LineTensor) + and isinstance(b, PlaneTensor) + or isinstance(b, LineTensor) + and isinstance(a, PlaneTensor) + ): e = LeviCivitaTensor(n) result = TensorDiagram(*[(e, a)] * a.tensor_shape[1], *[(e, b)] * b.tensor_shape[1]).calculate() - elif isinstance(a, Subspace) and isinstance(b, Point): + elif isinstance(a, SubspaceTensor) and isinstance(b, PointTensor): result = a * b - elif isinstance(a, Point) and isinstance(b, Subspace): + elif isinstance(a, PointTensor) and isinstance(b, SubspaceTensor): result = b * a - elif isinstance(a, Line) and isinstance(b, Line): + elif isinstance(a, LineTensor) and isinstance(b, LineTensor): # can assume that n >= 4, because for n = 3 lines are 1-tensors e = LeviCivitaTensor(n) @@ -130,7 +135,7 @@ def _join_meet_duality( elif intersect_lines or n == 4: # can't intersect lines that are not coplanar and can't join skew lines in 3D raise NotCoplanar("The given lines are not all coplanar.") - elif np.any(coplanar) and (a.cdim > 0 or b.cdim > 0): + elif np.any(coplanar) and (a.free_indices > 0 or b.free_indices > 0): raise GeometryException("Can only join tensors that are either all coplanar or all not coplanar.") else: @@ -143,7 +148,7 @@ def _join_meet_duality( if check_dependence: is_zero = result.is_zero() - if result.cdim == 0 and is_zero: + if result.free_indices == 0 and is_zero: raise LinearDependenceError("Arguments are not linearly independent.") elif np.any(is_zero): raise LinearDependenceError("Some arguments are not linearly independent.", is_zero) @@ -155,15 +160,15 @@ def _join_meet_duality( result.array = _divide_by_power_of_two(result.array, max_exponent) if result.tensor_shape == (0, 1): - return Line(result, copy=False) if n == 3 else Plane(result, copy=False) + return LineCollection.from_tensor(result) if n == 3 else PlaneCollection.from_tensor(result) if result.tensor_shape == (1, 0): - return Point(result, copy=False) + return PointCollection.from_tensor(result) if result.tensor_shape == (2, 0): - return Line(result, copy=False).contravariant_tensor + return LineCollection.from_tensor(result).contravariant_tensor if result.tensor_shape == (0, n - 2): - return Line(result, copy=False) + return LineCollection.from_tensor(result) - return Subspace(result, copy=False) + raise RuntimeError(f"Unexpected tensor of type {result.tensor_shape}") def _divide_by_power_of_two(array: np.ndarray, power: int) -> np.ndarray: @@ -183,16 +188,22 @@ def _divide_by_power_of_two(array: np.ndarray, power: int) -> np.ndarray: @overload -def join(*args: Unpack[tuple[Point, Point]], _check_dependence: bool = True, _normalize_result: bool = True) -> Line: +def join( + *args: Unpack[tuple[PointTensor, PointTensor]], _check_dependence: bool = True, _normalize_result: bool = True +) -> LineTensor: ... @overload -def join(*args: Point | Subspace, _check_dependence: bool = True, _normalize_result: bool = True) -> Subspace: +def join( + *args: PointTensor | SubspaceTensor, _check_dependence: bool = True, _normalize_result: bool = True +) -> SubspaceTensor: ... -def join(*args: Point | Subspace, _check_dependence: bool = True, _normalize_result: bool = True) -> Subspace: +def join( + *args: PointTensor | SubspaceTensor, _check_dependence: bool = True, _normalize_result: bool = True +) -> SubspaceTensor: """Joins a number of objects to form a line, plane or subspace. Args: @@ -212,21 +223,27 @@ def join(*args: Point | Subspace, _check_dependence: bool = True, _normalize_res @overload -def meet(*args: Line, _check_dependence: bool = True, _normalize_result: bool = True) -> Point: +def meet(*args: LineTensor, _check_dependence: bool = True, _normalize_result: bool = True) -> PointTensor: ... @overload -def meet(*args: Unpack[tuple[Subspace, Line]], _check_dependence: bool = True, _normalize_result: bool = True) -> Point: +def meet( + *args: Unpack[tuple[SubspaceTensor, LineTensor]], _check_dependence: bool = True, _normalize_result: bool = True +) -> PointTensor: ... @overload -def meet(*args: Unpack[tuple[Line, Subspace]], _check_dependence: bool = True, _normalize_result: bool = True) -> Point: +def meet( + *args: Unpack[tuple[LineTensor, SubspaceTensor]], _check_dependence: bool = True, _normalize_result: bool = True +) -> PointTensor: ... -def meet(*args: Subspace, _check_dependence: bool = True, _normalize_result: bool = True) -> Point | Subspace: +def meet( + *args: SubspaceTensor, _check_dependence: bool = True, _normalize_result: bool = True +) -> PointTensor | SubspaceTensor: """Intersects a number of given objects. Args: @@ -245,68 +262,102 @@ def meet(*args: Subspace, _check_dependence: bool = True, _normalize_result: boo ) -class Point(ProjectiveTensor): - """Represents points in a projective space of arbitrary dimension. - - The number of supplied coordinates determines the dimension of the space that the point lives in. - If the coordinates are given as scalar arguments (not in a single iterable), the coordinates will automatically be - transformed into homogeneous coordinates, i.e. a one added as an additional coordinate. - - Addition and subtraction of finite and infinite points will always give a finite result if one of the points - was finite beforehand. - - Args: - *args: A single iterable object or tensor or multiple (affine) coordinates. - homogenize: If True, and the first argument is an array of points, all points in the array will be converted to - homogeneous coordinates, i.e. 1 will be added to the coordinates of each point. - By default, homogenize is False. - **kwargs: Additional keyword arguments for the constructor of the numpy array as defined in `numpy.array`. - - """ +class PointLikeTensor(ProjectiveTensor, ABC): + """Base class for point tensors and polytopes implementing arithmetic operations.""" def __init__( - self, *args: Tensor | npt.ArrayLike, homogenize=False, tensor_rank=1, **kwargs: Unpack[NDArrayParameters] + self, + *args: Tensor | npt.ArrayLike, + homogenize: bool = False, + tensor_rank: int = 1, + **kwargs: Unpack[NDArrayParameters], ) -> None: - if np.isscalar(args[0]): - super().__init__(*args, 1, tensor_rank=tensor_rank, **kwargs) - homogenize = False - else: - super().__init__(*args, tensor_rank=tensor_rank, **kwargs) + super().__init__(*args, tensor_rank=tensor_rank, **kwargs) if homogenize is True: self.array = np.append(self.array, np.ones(self.shape[:-1] + (1,), self.dtype), axis=-1) + if self.tensor_shape != (1, 0): + raise ValueError(f"Expected tensor of type (1, 0), but is {self.tensor_shape}") + + @staticmethod + def _normalize_array(array: np.ndarray) -> np.ndarray: + z = array[..., -1, None] + isinf = np.isclose(z, 0, atol=EQ_TOL_ABS) + if np.all(isinf | (z == 1)): + return array + dtype = np.promote_types(np.float64, array.dtype) + result = array.astype(dtype) + np.divide(result, z, where=~isinf, out=result) + return np.real_if_close(result) + + @property + def normalized_array(self) -> np.ndarray: + """The coordinate array of the points with the last coordinates normalized to 1.""" + return self._normalize_array(self.array) def __add__(self, other: Tensor | npt.ArrayLike) -> Tensor: - if not isinstance(other, Point): + if not isinstance(other, PointLikeTensor): return super().__add__(other) a, b = self.normalized_array, other.normalized_array result = a[..., :-1] + b[..., :-1] result = np.append(result, np.maximum(a[..., -1:], b[..., -1:]), axis=-1) - return Point(result, copy=False) + return PointCollection.from_array(result) def __sub__(self, other: Tensor | npt.ArrayLike) -> Tensor: - if not isinstance(other, Point): + if not isinstance(other, PointLikeTensor): return super().__add__(other) a, b = self.normalized_array, other.normalized_array result = a[..., :-1] - b[..., :-1] result = np.append(result, np.maximum(a[..., -1:], b[..., -1:]), axis=-1) - return Point(result, copy=False) + return PointCollection.from_array(result) def __mul__(self, other: Tensor | npt.ArrayLike) -> Tensor: if not is_numerical_scalar(other): return super().__mul__(other) result = self.normalized_array[..., :-1] * other result = np.append(result, self.array[..., -1:] != 0, axis=-1) - return Point(result, copy=False) + return PointCollection.from_array(result) def __truediv__(self, other: Tensor | npt.ArrayLike) -> Tensor: if not is_numerical_scalar(other): return super().__truediv__(other) result = self.normalized_array[..., :-1] / other result = np.append(result, self.array[..., -1:] != 0, axis=-1) - return Point(result, copy=False) + return PointCollection.from_array(result) + + +class PointTensor(PointLikeTensor, ABC): + """Represents points in a projective space of arbitrary dimension. + + The number of supplied coordinates determines the dimension of the space that the point lives in. + If the coordinates are given as scalar arguments (not in a single iterable), the coordinates will automatically be + transformed into homogeneous coordinates, i.e. a one added as an additional coordinate. + + Addition and subtraction of finite and infinite points will always give a finite result if one of the points + was finite beforehand. + + Args: + *args: A single iterable object or tensor or multiple (affine) coordinates. + homogenize: If True, and the first argument is an array of points, all points in the array will be converted to + homogeneous coordinates, i.e. 1 will be added to the coordinates of each point. + By default, homogenize is False. + **kwargs: Additional keyword arguments for the constructor of the numpy array as defined in `numpy.array`. + + """ + + def __init__( + self, + *args: Tensor | npt.ArrayLike, + homogenize: bool = False, + tensor_rank: int = 1, + **kwargs: Unpack[NDArrayParameters], + ) -> None: + if np.isscalar(args[0]): + super().__init__(*args, 1, homogenize=False, tensor_rank=tensor_rank, **kwargs) + else: + super().__init__(*args, homogenize=homogenize, tensor_rank=tensor_rank, **kwargs) def __repr__(self) -> str: - if self.cdim == 0: + if self.free_indices == 0: result = f"Point({', '.join(self.normalized_array[:-1].astype(str))})" if self.isinf: result += " at Infinity" @@ -314,14 +365,16 @@ def __repr__(self) -> str: return f"PointCollection({self.normalized_array.tolist()})" @overload - def join(self, *others: Unpack[tuple[Point, Point]]) -> Line: + def join(self, *others: Unpack[tuple[PointTensor, PointTensor]]) -> LineTensor: ... @overload - def join(self, *others: Unpack[tuple[Point, Subspace]] | Unpack[tuple[Subspace, Point]]) -> Subspace: + def join( + self, *others: Unpack[tuple[PointTensor, SubspaceTensor]] | Unpack[tuple[SubspaceTensor, PointTensor]] + ) -> SubspaceTensor: ... - def join(self, *others: Point | Subspace) -> Subspace: + def join(self, *others: PointTensor | LineTensor) -> SubspaceTensor: return join(self, *others) def __getitem__(self, index: TensorIndex) -> Tensor | np.generic: @@ -330,28 +383,12 @@ def __getitem__(self, index: TensorIndex) -> Tensor | np.generic: if not isinstance(result, Tensor) or result.tensor_shape != (1, 0): return result - return Point(result, copy=False) - - @staticmethod - def _normalize_array(array: np.ndarray) -> np.ndarray: - z = array[..., -1, None] - isinf = np.isclose(z, 0, atol=EQ_TOL_ABS) - if np.all(isinf | (z == 1)): - return array - dtype = np.promote_types(np.float64, array.dtype) - result = array.astype(dtype) - np.divide(result, z, where=~isinf, out=result) - return np.real_if_close(result) + return PointCollection.from_array(result) - def _matrix_transform(self, m: npt.ArrayLike) -> Point: - if self.cdim == 0: - return Point(np.dot(m, self.array), copy=False) - return Point(matvec(m, self.array), copy=False) - - @property - def normalized_array(self) -> np.ndarray: - """The coordinate array of the points with the last coordinates normalized to 1.""" - return self._normalize_array(self.array) + def _matrix_transform(self, m: npt.ArrayLike) -> PointTensor: + if self.free_indices == 0: + return PointCollection.from_array(np.dot(m, self.array)) + return PointCollection.from_array(matvec(m, self.array)) @property def isinf(self) -> npt.NDArray[np.bool_]: @@ -364,8 +401,16 @@ def isreal(self) -> npt.NDArray[np.bool_]: return np.all(np.isreal(self.normalized_array), axis=-1) -class Subspace(ProjectiveTensor): - """Represents a general subspace of a projective space. Line and Plane are subclasses. +class Point(PointTensor, BoundTensor): + pass + + +class PointCollection(PointTensor, TensorCollection[Point]): + _element_class = Point + + +class SubspaceTensor(ProjectiveTensor, ABC): + """Abstract base class for subspaces of a projective space. Line and Plane are subclasses. Args: *args: The coordinates of the subspace. Instead of separate coordinates, a single iterable can be supplied. @@ -381,29 +426,21 @@ def __init__( super().__init__(*args, tensor_rank=tensor_rank, **kwargs) @overload - def meet(self, other: Line) -> Point: + def meet(self, other: LineTensor) -> PointTensor: ... @overload - def meet(self, other: Subspace) -> Point | Subspace: + def meet(self, other: SubspaceTensor) -> PointTensor | SubspaceTensor: ... - def meet(self, other: Subspace) -> Point | Subspace: + def meet(self, other: SubspaceTensor) -> PointTensor | SubspaceTensor: return meet(self, other) - def join(self, *others: Point | Subspace) -> Subspace: + def join(self, *others: PointTensor | SubspaceTensor) -> SubspaceTensor: return join(self, *others) - def __getitem__(self, index: TensorIndex) -> Tensor | np.generic: - result = super().__getitem__(index) - - if not isinstance(result, Tensor) or result.tensor_shape != self.tensor_shape: - return result - - return type(self)(result, copy=False) - def __add__(self, other: Tensor | npt.ArrayLike) -> Tensor: - if not isinstance(other, Point): + if not isinstance(other, PointTensor): return super().__add__(other) from geometer.transformation import translation @@ -418,20 +455,20 @@ def __sub__(self, other: Tensor | npt.ArrayLike) -> Tensor: @property def basis_matrix(self) -> np.ndarray: x = self.array - x = x.reshape(x.shape[: self.cdim] + (-1, x.shape[-1])) - return np.swapaxes(null_space(x, self.shape[-1] - self.rank + self.cdim), -1, -2) + x = x.reshape(x.shape[: self.free_indices] + (-1, x.shape[-1])) + return np.swapaxes(null_space(x, self.shape[-1] - self.rank + self.free_indices), -1, -2) - def _matrix_transform(self, m: npt.ArrayLike) -> Subspace: + def _matrix_transform(self, m: npt.ArrayLike) -> SubspaceTensor: transformed_basis = matmul(self.basis_matrix, m, transpose_b=True) - transformed_basis_points = [Point(p, copy=False) for p in np.swapaxes(transformed_basis, 0, -2)] + transformed_basis_points = [PointCollection.from_array(p) for p in np.swapaxes(transformed_basis, 0, -2)] return join(*transformed_basis_points) @property - def general_point(self) -> Point: - """Points in general position i.e. not in the subspaces.""" + def general_point(self) -> PointTensor: + """Point in general position i.e. not in the subspaces.""" n = self.dim + 1 - s = [self.shape[i] for i in self._collection_indices] - p = Point(np.zeros([*s, n], dtype=int), copy=False) + s = [self.shape[i] for i in range(self.free_indices)] + p = PointCollection.from_array(np.zeros([*s, n], dtype=int)) ind = np.ones(s, dtype=bool) for i in range(n): p[ind, -i - 1] = 1 @@ -440,20 +477,20 @@ def general_point(self) -> Point: break return p - def parallel(self, through: Point) -> Subspace: - """Returns the subspaces through given points that are parallel to this collection of subspaces. + def parallel(self, through: PointTensor) -> SubspaceTensor: + """Returns the subspace through given point that is parallel to this subspace. Args: - through: The point through which the parallel subspaces are to be constructed. + through: The point through which the parallel subspace is to be constructed. Returns: - The parallel subspaces. + The parallel subspace. """ x = self.meet(infty_hyperplane(self.dim)) return join(x, through) - def contains(self, other: Point | Subspace, tol: float = EQ_TOL_ABS) -> npt.NDArray[np.bool_]: + def contains(self, other: PointTensor | LineTensor, tol: float = EQ_TOL_ABS) -> npt.NDArray[np.bool_]: """Tests whether given points or lines lie in the subspaces. Args: @@ -464,20 +501,17 @@ def contains(self, other: Point | Subspace, tol: float = EQ_TOL_ABS) -> npt.NDAr Boolean array that indicates which of given points/lines lies in the subspaces. """ - if isinstance(other, Point): + if isinstance(other, PointTensor): result = self * other - - elif isinstance(other, Line): + elif isinstance(other, LineTensor): result = self * other.covariant_tensor - else: - # TODO: test subspace - raise ValueError(f"argument of type {type(other)} not supported") + raise TypeError(f"argument of type {type(other)} not supported") axes = tuple(result._covariant_indices) + tuple(result._contravariant_indices) return np.all(np.isclose(result.array, 0, atol=tol), axis=axes) - def is_parallel(self, other: Subspace) -> npt.NDArray[np.bool_]: + def is_parallel(self, other: SubspaceTensor) -> npt.NDArray[np.bool_]: """Tests whether the given subspace is parallel to this subspace. Args: @@ -490,8 +524,56 @@ def is_parallel(self, other: Subspace) -> npt.NDArray[np.bool_]: x = self.meet(other) return infty_hyperplane(self.dim).contains(x) + @abstractmethod + def mirror(self, pt: PointTensor) -> PointTensor: + """Construct the reflection of a point at the subspace. + + Args: + pt: The point to reflect. + + Returns: + The mirror point. + + """ + + @abstractmethod + def perpendicular(self, through: PointTensor | LineTensor) -> SubspaceTensor: + """Construct the perpendicular subspace though the given point or line. + + Args: + through: The point or line through which the perpendicular is constructed. + + Returns: + The perpendicular subspace. + + """ + + def project(self, pt: PointTensor) -> PointTensor: + """The orthogonal projection of a point onto the subspace. + + Args: + pt: The point to project. + + Returns: + The projected point. + + """ + l = self.perpendicular(pt) + return self.meet(l) + -class Line(Subspace): +class Subspace(SubspaceTensor, BoundTensor, ABC): + pass + + +SubspaceT = TypeVar("SubspaceT", bound=Subspace) + + +class SubspaceCollection(SubspaceTensor, TensorCollection[SubspaceT], ABC): + pass + + +class LineTensor(SubspaceTensor, ABC): """Represents a line in a projective space of arbitrary dimension. Args: @@ -505,27 +587,41 @@ def __init__(self, *args: Tensor | npt.ArrayLike, **kwargs: Unpack[NDArrayParame if len(args) == 2: kwargs["copy"] = False p, q = args - super().__init__(join(Point(p, copy=False), Point(q, copy=False)), tensor_rank=-2, **kwargs) + super().__init__( + join(PointCollection.from_array(p), PointCollection.from_array(q)), tensor_rank=-2, **kwargs + ) else: super().__init__(*args, tensor_rank=-2, **kwargs) + if self.tensor_shape not in {(0, self.dim - 1), (self.dim - 1, 0)}: + raise ValueError(f"Unexpected tensor of type {self.tensor_shape}") + if self.dim == 3 and self.shape[-1] != self.shape[-2]: + raise ValueError(f"Expected quadratic matrix, but last two dimensions are {self.shape[-2:]}") + + def __getitem__(self, index: TensorIndex) -> Tensor | np.generic: + result = super().__getitem__(index) + + if not isinstance(result, Tensor) or result.tensor_shape != self.tensor_shape: + return result + + return LineCollection.from_tensor(result) - def _matrix_transform(self, m: npt.ArrayLike) -> Line: + def _matrix_transform(self, m: npt.ArrayLike) -> LineTensor: result = super()._matrix_transform(m) - return cast(Line, result) + return cast(LineTensor, result) - def meet(self, other: Subspace) -> Point: + def meet(self, other: SubspaceTensor) -> PointTensor: result = super().meet(other) - return cast(Point, result) + return cast(PointTensor, result) @property - def base_point(self) -> Point: + def base_point(self) -> PointTensor: """Base points for the lines, arbitrarily chosen.""" if self.dim > 2: base = self.basis_matrix p, q = base[..., 0, :], base[..., 1, :] isinf = np.isclose(p[..., -1, None], 0, atol=EQ_TOL_ABS) result = np.where(isinf, q, p) - return Point(result, copy=False) + return PointCollection.from_array(result) y_zero = np.isclose(self.array[..., 1], 0, atol=EQ_TOL_ABS) z_zero = np.isclose(self.array[..., 2], 0, atol=EQ_TOL_ABS) @@ -539,10 +635,10 @@ def base_point(self) -> Point: result[y_zero & ~z_zero, 0] = self.array[y_zero & ~z_zero, 2] result[y_zero & ~z_zero, 2] = -self.array[y_zero & ~z_zero, 0] - return Point(result, copy=False) + return PointCollection.from_array(result) @property - def direction(self) -> Point: + def direction(self) -> PointTensor: """The direction of the lines (not normalized).""" if self.dim > 2: return meet(self, infty_hyperplane(self.dim), _normalize_result=False) @@ -556,20 +652,20 @@ def direction(self) -> Point: result[~(x_zero & y_zero), 0] = self.array[~(x_zero & y_zero), 1] result[~(x_zero & y_zero), 1] = -self.array[~(x_zero & y_zero), 0] - return Point(result, copy=False) + return PointCollection.from_array(result) @property def basis_matrix(self) -> np.ndarray: """A matrix with orthonormal basis vectors as rows.""" if self.dim == 2: a = self.base_point.array - b = np.cross(self.array, a) + b = np.cross(self.array, a) # type: ignore[arg-type] m = [a / np.linalg.norm(a, axis=-1, keepdims=True), b / np.linalg.norm(b, axis=-1, keepdims=True)] return np.stack(m, axis=-2) return super().basis_matrix @property - def covariant_tensor(self) -> Line: + def covariant_tensor(self) -> LineTensor: """The covariant tensors of lines in 3D.""" if self.dim != 3: raise NotImplementedError(f"Expected dimension 3 but found dimension {self.dim}.") @@ -577,10 +673,10 @@ def covariant_tensor(self) -> Line: return self e = LeviCivitaTensor(4) diagram = TensorDiagram((e, self), (e, self)) - return Line(diagram.calculate(), copy=False) + return LineCollection.from_tensor(diagram.calculate()) @property - def contravariant_tensor(self) -> Line: + def contravariant_tensor(self) -> LineTensor: """The contravariant tensors of lines in 3D.""" if self.dim != 3: raise NotImplementedError(f"Expected dimension 3 but found dimension {self.dim}.") @@ -588,9 +684,9 @@ def contravariant_tensor(self) -> Line: return self e = LeviCivitaTensor(4, False) diagram = TensorDiagram((self, e), (self, e)) - return Line(diagram.calculate(), copy=False) + return LineCollection.from_tensor(diagram.calculate()) - def is_coplanar(self, other: Line) -> npt.NDArray[np.bool_]: + def is_coplanar(self, other: LineTensor) -> npt.NDArray[np.bool_]: """Tests whether another line lies in the same plane as this line, i.e. whether two lines intersect. Args: @@ -604,13 +700,13 @@ def is_coplanar(self, other: Line) -> npt.NDArray[np.bool_]: """ if self.dim == 2: - return np.ones(self.shape[: self.cdim], dtype=bool) + return np.ones(self.shape[: self.free_indices], dtype=bool) e = LeviCivitaTensor(self.dim + 1) d = TensorDiagram(*[(e, self)] * (self.dim - 1), *[(e, other)] * (self.dim - 1)) return d.calculate().is_zero() - def mirror(self, pt: Point) -> Point: + def mirror(self, pt: PointTensor) -> PointTensor: """Construct the reflection of points at the lines. Args: @@ -629,7 +725,7 @@ def mirror(self, pt: Point) -> Point: e = join(self, pt) if self.dim == 3: - e = cast(Plane, e) + e = cast(PlaneTensor, e) f = e.perpendicular(self) return f.mirror(pt) @@ -638,7 +734,7 @@ def mirror(self, pt: Point) -> Point: arr_sort = np.argsort(np.abs(arr), axis=-1) arr_ind = tuple(np.indices(arr.shape)[:-1]) m = m[(*arr_ind, arr_sort, slice(None))] - pt = Point(arr[(*arr_ind, arr_sort)], copy=False) + pt = PointCollection(arr[(*arr_ind, arr_sort)], copy=False) l = self._matrix_transform(m) l1 = join(I, pt, _normalize_result=False) l2 = join(J, pt, _normalize_result=False) @@ -651,7 +747,7 @@ def mirror(self, pt: Point) -> Point: return result._matrix_transform(np.swapaxes(m, -1, -2)) return result - def perpendicular(self, through: Point, plane: Subspace | None = None) -> Line: + def perpendicular(self, through: PointTensor, plane: PlaneTensor | None = None) -> LineTensor: """Construct the perpendicular line though a point. Args: @@ -665,55 +761,52 @@ def perpendicular(self, through: Point, plane: Subspace | None = None) -> Line: """ n = self.dim + 1 contains = self.contains(through) - result = Line(np.empty(contains.shape + (n,) * (n - 2), np.complex128)) + result = LineCollection.from_array(np.empty(contains.shape + (n,) * (n - 2), np.complex128)) if np.any(contains): l = self - if self.cdim > 0: + if self.free_indices > 0: l = self[contains] if n > 3: if plane is None: # additional point is required to determine the exact line plane = join(l, l.general_point) - elif plane.cdim > 0: + elif isinstance(plane, PlaneCollection): plane = plane[contains] basis = plane.basis_matrix line_pts = matmul(l.basis_matrix, basis, transpose_b=True) - l = Line(np.cross(line_pts[..., 0, :], line_pts[..., 1, :]), copy=False) + l = LineCollection.from_array(np.cross(line_pts[..., 0, :], line_pts[..., 1, :])) - p = Point(np.append(l.array[..., :-1], np.zeros(l.shape[:-1] + (1,), dtype=l.dtype), axis=-1), copy=False) + p = PointCollection.from_array( + np.append(l.array[..., :-1], np.zeros(l.shape[:-1] + (1,), dtype=l.dtype), axis=-1) + ) if n > 3: p = p._matrix_transform(np.swapaxes(basis, -1, -2)) - result[contains] = join(through if through.cdim == 0 else through[contains], p) + result[contains] = join(through if through.free_indices == 0 else through[contains], p) if np.any(~contains): - if through.cdim > 0: + if through.free_indices > 0: through = through[~contains] - if self.cdim > 0: - result[~contains] = cast(Line, self[~contains]).mirror(through).join(through) + if self.free_indices > 0: + result[~contains] = cast(LineTensor, self[~contains]).mirror(through).join(through) else: result[~contains] = self.mirror(through).join(through) - return Line(np.real_if_close(result.array), copy=False) + return LineCollection.from_array(np.real_if_close(result.array)) - def project(self, pt: Point) -> Point: - """The orthogonal projection of points onto the lines. - Args: - pt: The points to project. +class Line(LineTensor, Subspace): + pass - Returns: - The projected points. - """ - l = self.perpendicular(pt) - return self.meet(l) +class LineCollection(LineTensor, SubspaceCollection[Line]): + _element_class = Line -class Plane(Subspace): +class PlaneTensor(SubspaceTensor): """Represents a hyperplane in a projective space of arbitrary dimension. Args: @@ -724,11 +817,21 @@ class Plane(Subspace): """ def __init__(self, *args: Tensor | npt.ArrayLike, **kwargs: Unpack[NDArrayParameters]) -> None: - if all(isinstance(o, (Line, Point)) for o in args): + if all(isinstance(o, (LineTensor, PointTensor)) for o in args): kwargs["copy"] = False super().__init__(join(*args), **kwargs) else: super().__init__(*args, **kwargs) + if self.tensor_shape != (0, 1): + raise ValueError(f"Expected tensor of type (0, 1), but is {self.tensor_shape}") + + def __getitem__(self, index: TensorIndex) -> Tensor | np.generic: + result = super().__getitem__(index) + + if not isinstance(result, Tensor) or result.tensor_shape != (0, 1): + return result + + return PlaneCollection.from_tensor(result) @property def basis_matrix(self) -> np.ndarray: @@ -745,7 +848,7 @@ def basis_matrix(self) -> np.ndarray: q, r = np.linalg.qr(result) return np.swapaxes(q, -1, -2) - def mirror(self, pt: Point) -> Point: + def mirror(self, pt: PointTensor) -> PointTensor: """Construct the reflections of a point at the plane. Only works in 3D. @@ -761,15 +864,15 @@ def mirror(self, pt: Point) -> Point: raise NotImplementedError(f"Expected dimension 3 but found dimension {self.dim}.") l = self.meet(infty_plane) basis = l.basis_matrix - l = Line(np.cross(basis[..., 0, :-1], basis[..., 1, :-1]), copy=False) + l = LineCollection.from_array(np.cross(basis[..., 0, :-1], basis[..., 1, :-1])) p = l.base_point - polars = Line(p.array, copy=False) + polars = LineCollection.from_array(p.array) from geometer.curve import absolute_conic p1, p2 = absolute_conic.intersect(polars) - p1 = Point(np.append(p1.array, np.zeros(p1.shape[:-1] + (1,)), axis=-1), copy=False) - p2 = Point(np.append(p2.array, np.zeros(p2.shape[:-1] + (1,)), axis=-1), copy=False) + p1 = PointCollection.from_array(np.append(p1.array, np.zeros(p1.shape[:-1] + (1,)), axis=-1)) + p2 = PointCollection.from_array(np.append(p2.array, np.zeros(p2.shape[:-1] + (1,)), axis=-1)) l1 = p1.join(pt) l2 = p2.join(pt) @@ -779,28 +882,15 @@ def mirror(self, pt: Point) -> Point: m2 = q2.join(p1) return m1.meet(m2) - def project(self, pt: Point) -> Point: - """The orthogonal projection of points onto the planes. - - Args: - pt: The points to project. - - Returns: - The projected points. - - """ - l = self.perpendicular(pt) - return self.meet(l) - @overload - def perpendicular(self, through: Point) -> Line: + def perpendicular(self, through: PointTensor) -> LineTensor: ... @overload - def perpendicular(self, through: Line) -> Plane: + def perpendicular(self, through: LineTensor) -> PlaneTensor: ... - def perpendicular(self, through: Point | Line) -> Line | Plane: + def perpendicular(self, through: PointTensor | LineTensor) -> LineTensor | PlaneTensor: """Construct the perpendicular lines though the given points or the perpendicular planes through the given lines. @@ -813,13 +903,26 @@ def perpendicular(self, through: Point | Line) -> Line | Plane: The perpendicular lines or planes. """ - if self.dim != 3 and isinstance(through, Line): + if self.dim != 3 and isinstance(through, LineTensor): raise NotImplementedError(f"Expected dimension 3 but found dimension {self.dim}.") p = self.array[..., :-1] - p = Point(np.append(p, np.zeros(p.shape[:-1] + (1,), dtype=p.dtype), axis=-1), copy=False) + p = PointCollection.from_array(np.append(p, np.zeros(p.shape[:-1] + (1,), dtype=p.dtype), axis=-1)) return through.join(p) + @property + def isinf(self) -> npt.NDArray[np.bool_]: + """Boolean array that indicates whether the plane is the hyperplane at infinity.""" + return is_multiple(self.array, infty_hyperplane(self.dim).array, axis=-1, rtol=EQ_TOL_REL, atol=EQ_TOL_ABS) + + +class Plane(PlaneTensor, Subspace): + pass + + +class PlaneCollection(PlaneTensor, SubspaceCollection[Plane]): + _element_class = Plane + I = Point([-1j, 1, 0]) J = Point([1j, 1, 0]) diff --git a/geometer/shapes.py b/geometer/shapes.py index 5646254..2046089 100644 --- a/geometer/shapes.py +++ b/geometer/shapes.py @@ -1,24 +1,38 @@ from __future__ import annotations import math +from abc import ABC from itertools import combinations -from typing import TYPE_CHECKING, cast +from typing import TYPE_CHECKING, TypeVar, cast import numpy as np import numpy.typing as npt -from geometer.base import EQ_TOL_ABS, EQ_TOL_REL, NDArrayParameters, Tensor, TensorIndex -from geometer.exceptions import LinearDependenceError, NotCoplanar +from geometer.base import EQ_TOL_ABS, EQ_TOL_REL, Tensor, TensorCollection +from geometer.exceptions import IncompatibleShapeError, LinearDependenceError, NotCoplanar from geometer.operators import angle, dist, harmonic_set -from geometer.point import Line, Plane, Point, infty_hyperplane, join, meet -from geometer.transformation import Transformation, rotation, translation +from geometer.point import ( + LineTensor, + Plane, + PlaneTensor, + Point, + PointCollection, + PointLikeTensor, + PointTensor, + infty_hyperplane, + join, + meet, +) +from geometer.transformation import TransformationTensor, rotation, translation from geometer.utils import det, distinct, is_multiple, matmul, matvec if TYPE_CHECKING: from typing_extensions import Unpack + from geometer.utils.typing import NDArrayParameters, TensorIndex -class Polytope(Point): + +class PolytopeTensor(PointLikeTensor, ABC): """A class representing polytopes in arbitrary dimension. A (n+1)-polytope is a collection of n-polytopes that have some (n-1)-polytopes in common, where 3-polytopes are polyhedra, 2-polytopes are polygons, 1-polytopes are line segments and 0-polytopes are points. @@ -44,21 +58,18 @@ def __init__(self, *args: Tensor | npt.ArrayLike, pdim: int = 0, **kwargs: Unpac super().__init__(*args, **kwargs) def __repr__(self) -> str: - class_name = self.__class__.__name__ - if self.cdim - max(self.pdim - 1, 1) > 0: - class_name += "Collection" - return f"{class_name}({', '.join(str(v) for v in self.vertices)})" + return f"{self.__class__.__name__}({', '.join(str(v) for v in self.vertices)})" @property - def vertices(self) -> list[Point]: + def vertices(self) -> list[PointTensor]: """The vertices of the polytope.""" first_polygon_index = self.rank - max(self.pdim - 1, 1) - 1 new_shape = self.shape[:first_polygon_index] + (-1, self.shape[-1]) array = self.array.reshape(new_shape) - return list(distinct(Point(x, copy=False) for x in np.moveaxis(array, -2, 0))) + return list(distinct(PointCollection.from_array(x) for x in np.moveaxis(array, -2, 0))) @property - def facets(self) -> list[Polytope]: + def facets(self) -> list[PolytopeTensor]: """The facets of the polytope.""" first_polygon_index = self.rank - max(self.pdim - 1, 1) - 1 slices = (slice(None),) * first_polygon_index @@ -70,8 +81,8 @@ def _edges(self) -> np.ndarray: v2 = np.roll(v1, -1, axis=-2) return np.stack([v1, v2], axis=-2) - def __eq__(self, other: Tensor | npt.ArrayLike) -> bool: - if not isinstance(other, Polytope): + def __eq__(self, other: object) -> bool: + if not isinstance(other, PolytopeTensor): return super().__eq__(other) if self.shape != other.shape: @@ -98,7 +109,7 @@ def __eq__(self, other: Tensor | npt.ArrayLike) -> bool: return False def __add__(self, other: Tensor | npt.ArrayLike) -> Tensor: - if not isinstance(other, Point): + if not isinstance(other, PointTensor): return super().__add__(other) return translation(other).apply(self) @@ -106,35 +117,36 @@ def __sub__(self, other: Tensor | npt.ArrayLike) -> Tensor: return self + (-other) @staticmethod - def _cast_polytope(tensor: Tensor, pdim: int) -> Polytope: + def _cast_polytope(tensor: Tensor, pdim: int) -> PolytopeTensor: if pdim == 1: - return Segment(tensor, copy=False) + return SegmentCollection.from_tensor(tensor) if pdim == 2: if tensor.shape[-2] == 3: + # TODO: check if a collection can be returned here return Triangle(tensor, copy=False) if tensor.shape[-2] == 4: try: return Rectangle(tensor, copy=False) except NotCoplanar: - return Polytope(tensor, copy=False) + return PolytopeCollection.from_tensor(tensor) try: - return Polygon(tensor, copy=False) + return PolygonCollection.from_tensor(tensor) except NotCoplanar: - return Polytope(tensor, copy=False) + return PolytopeCollection.from_tensor(tensor) if pdim == 3: return Polyhedron(tensor, copy=False) - return Polytope(tensor, copy=False) + return PolytopeCollection.from_tensor(tensor) def __getitem__(self, index: TensorIndex) -> Tensor | np.generic: result = super().__getitem__(index) - if not isinstance(result, Tensor) or result.cdim == 0 or result.tensor_shape != (1, 0): + if not isinstance(result, Tensor) or result.free_indices == 0 or result.tensor_shape != (1, 0): return result index_mapping = self._get_index_mapping(index) - removed_cdims = [i for i in self._collection_indices if i not in index_mapping] + removed_cdims = [i for i in range(self.free_indices) if i not in index_mapping] first_polygon_index = self.rank - max(self.pdim - 1, 1) - 1 result_pdim = self.pdim - len([i for i in removed_cdims if i >= first_polygon_index]) @@ -144,7 +156,23 @@ def __getitem__(self, index: TensorIndex) -> Tensor | np.generic: return self._cast_polytope(result, result_pdim) -class Segment(Polytope): +class Polytope(PolytopeTensor, ABC): + pass + + +SubspaceT = TypeVar("SubspaceT", bound=Polytope) + + +class PolytopeCollection(PolytopeTensor, TensorCollection[SubspaceT], ABC): + def _validate_tensor(self) -> None: + super()._validate_tensor() + if self.free_indices <= max(self.pdim - 1, 1): + raise IncompatibleShapeError( + f"Tensor has not enough free indices and cannot be used in a {self.__class__.__name__}." + ) + + +class SegmentTensor(PolytopeTensor): """Represents a line segment in an arbitrary projective space. Segments with one point at infinity represent rays/half-lines in a traditional sense. @@ -155,7 +183,7 @@ class Segment(Polytope): """ - _line: Line + _line: LineTensor def __init__(self, *args: Tensor | npt.ArrayLike, **kwargs: Unpack[NDArrayParameters]) -> None: if len(args) == 2: @@ -168,7 +196,7 @@ def __init__(self, *args: Tensor | npt.ArrayLike, **kwargs: Unpack[NDArrayParame self._line = join(*self.vertices) - def __apply__(self, transformation: Transformation) -> Segment: + def __apply__(self, transformation: TransformationTensor) -> SegmentTensor: result = super().__apply__(transformation) result._line = transformation.apply(result._line) return result @@ -179,29 +207,24 @@ def __getitem__(self, index: TensorIndex) -> Tensor | np.generic: if not isinstance(result, Tensor) or result.rank < 2 or result.shape[-2] != 2: return result - return Segment(result, copy=False) - - def expand_dims(self, axis: int) -> Segment: - result = super().expand_dims(axis) - result._line = result._line.expand_dims(axis - self.dim + 3 if axis < -1 else axis) - return result + return SegmentCollection.from_tensor(result) @property - def vertices(self) -> list[Point]: + def vertices(self) -> list[PointTensor]: """The start and endpoint of the line segment.""" - a = Point(self.array[..., 0, :], copy=False) - b = Point(self.array[..., 1, :], copy=False) + a = PointCollection.from_array(self.array[..., 0, :]) + b = PointCollection.from_array(self.array[..., 1, :]) return [a, b] @property - def facets(self) -> list[Point]: + def facets(self) -> list[PointTensor]: return self.vertices @property def _edges(self) -> np.ndarray: return self.array - def contains(self, other: Point, tol: float = EQ_TOL_ABS) -> npt.NDArray[np.bool_]: + def contains(self, other: PointTensor, tol: float = EQ_TOL_ABS) -> npt.NDArray[np.bool_]: """Tests whether a point is contained in the segment. Args: @@ -240,7 +263,9 @@ def contains(self, other: Point, tol: float = EQ_TOL_ABS) -> npt.NDArray[np.bool y_zero = np.isclose(y, 0, atol=EQ_TOL_ABS) return result & (~x_zero | ~y_zero) & (0 <= x + tol) & (x <= y + tol) - def intersect(self, other: Line | Plane | Segment | Polygon | Polyhedron) -> list[Point]: + def intersect( + self, other: LineTensor | PlaneTensor | SegmentTensor | PolygonTensor | Polyhedron + ) -> list[PointTensor]: """Intersect the line segment with another object. Args: @@ -250,24 +275,24 @@ def intersect(self, other: Line | Plane | Segment | Polygon | Polyhedron) -> lis The points of intersection. """ - if isinstance(other, (Polygon, Polyhedron)): + if isinstance(other, (PolygonTensor, Polyhedron)): return other.intersect(self) - if isinstance(other, Segment): + if isinstance(other, SegmentTensor): result = meet(self._line, other._line, _check_dependence=False) ind = ~result.is_zero() & self.contains(result) & other.contains(result) else: result = meet(self._line, other, _check_dependence=False) ind = ~result.is_zero() & self.contains(result) - if result.cdim > 0: + if result.free_indices > 0: return list(result[ind]) if ind: return [result] return [] @property - def midpoint(self) -> Point: + def midpoint(self) -> PointTensor: """The midpoint of the segment.""" l = self._line.meet(infty_hyperplane(self.dim)) return harmonic_set(*self.vertices, l) @@ -278,6 +303,19 @@ def length(self) -> npt.NDArray[np.float_]: return dist(*self.vertices) +class Segment(SegmentTensor, Polytope): + pass + + +class SegmentCollection(SegmentTensor, PolytopeCollection[Segment]): + _element_class = Segment + + def expand_dims(self, axis: int) -> SegmentCollection: + result = super().expand_dims(axis) + result._line = result._line.expand_dims(axis - self.dim + 3 if axis < -1 else axis) + return result + + class Simplex(Polytope): """Represents a simplex in any dimension, i.e. a k-polytope with k+1 vertices where k is the dimension. @@ -290,11 +328,11 @@ class Simplex(Polytope): """ # TODO: return only instances of type Simplex - def __new__(cls, *args: Point, **kwargs: Unpack[NDArrayParameters]) -> Polytope: + def __new__(cls, *args: Point, **kwargs: Unpack[NDArrayParameters]) -> PolytopeTensor: if len(args) == 2: return Segment(*args, **kwargs) - return super(Polytope, cls).__new__(cls) + return super(PolytopeTensor, cls).__new__(cls) def __init__(self, *args: Point, **kwargs: Unpack[NDArrayParameters]) -> None: kwargs.setdefault("pdim", len(args) - 1) @@ -324,7 +362,7 @@ def volume(self) -> float: return np.sqrt((-1) ** n / (math.factorial(n - 1) ** 2 * 2 ** (n - 1)) * det(m)) -class Polygon(Polytope): +class PolygonTensor(PolytopeTensor): """A flat polygon with vertices in any dimension. The vertices of the polygon must be given either in clockwise or counterclockwise order. @@ -336,10 +374,10 @@ class Polygon(Polytope): """ - _plane: Plane + _plane: PlaneTensor def __init__(self, *args: Tensor | npt.ArrayLike, **kwargs: Unpack[NDArrayParameters]) -> None: - if all(isinstance(x, Segment) for x in args): + if all(isinstance(x, SegmentTensor) for x in args): args = tuple(s.array[..., 0, :] for s in args) if len(args) > 1: args = tuple(a.array for a in args) @@ -347,28 +385,28 @@ def __init__(self, *args: Tensor | npt.ArrayLike, **kwargs: Unpack[NDArrayParame kwargs["copy"] = False args = (np.stack(args, axis=-2),) super().__init__(*args, pdim=2, **kwargs) - self._plane = Plane(*self.vertices[: self.dim]) if self.dim > 2 else None + self._plane = join(*self.vertices[: self.dim]) if self.dim > 2 else None - def __apply__(self, transformation: Transformation) -> Polygon: + def __apply__(self, transformation: TransformationTensor) -> PolygonTensor: result = super().__apply__(transformation) if result.dim > 2: - result._plane = Plane(*result.vertices[: result.dim]) + result._plane = join(*result.vertices[: result.dim]) return result @property - def vertices(self) -> list[Point]: - return [Point(self.array[..., i, :], copy=False) for i in range(self.shape[-2])] + def vertices(self) -> list[PointTensor]: + return [PointCollection.from_array(self.array[..., i, :]) for i in range(self.shape[-2])] @property - def facets(self) -> list[Segment]: + def facets(self) -> list[SegmentTensor]: return list(self.edges) @property - def edges(self) -> Segment: + def edges(self) -> SegmentCollection: """The edges of the polygon.""" - return Segment(self._edges, copy=False) + return SegmentCollection(self._edges, copy=False) - def contains(self, other: Point) -> npt.NDArray[np.bool_]: + def contains(self, other: PointTensor) -> npt.NDArray[np.bool_]: """Tests whether a point is contained in the polygon. Points on an edge of the polygon are considered True. @@ -392,40 +430,37 @@ def contains(self, other: Point) -> npt.NDArray[np.bool_]: if not np.any(coplanar): return coplanar - isinf = is_multiple(self._plane.array, infty_hyperplane(self.dim).array, axis=-1) - # remove coordinates with the largest absolute value in the normal vector i = np.argmax(np.abs(self._plane.array[..., :-1]), axis=-1) - i = np.where(isinf, self.dim, i) + i = np.where(self._plane.isinf, self.dim, i) i = np.expand_dims(i, -1) s = self.array.shape - arr = np.delete(self.array, np.ravel_multi_index((*tuple(np.indices(s[:-1])), i), s)).reshape( - s[:-1] + (-1,) - ) + arr = np.delete(self.array, np.ravel_multi_index((*tuple(np.indices(s[:-1])), i), s)) + arr = arr.reshape(s[:-1] + (-1,)) - if other.cdim == 0 and i.ndim == 1: + if isinstance(other, Point) and i.ndim == 1: other = Point(np.delete(other.array, i), copy=False) else: s = other.shape[:-1] + (1, other.shape[-1]) other = np.delete(other.array, np.ravel_multi_index((*tuple(np.indices(s[:-1])), i), s)) - other = Point(other.reshape(s[:-2] + (-1,)), copy=False) + other = PointCollection(other.reshape(s[:-2] + (-1,)), copy=False) # TODO: only test coplanar points - return coplanar & Polygon(arr, copy=False).contains(other) + return coplanar & PolygonCollection.from_array(arr).contains(other) edges = self.edges - edge_points = edges.contains(Point(np.expand_dims(other.array, -2), copy=False)) + edge_points = edges.contains(PointCollection(np.expand_dims(other.array, -2), copy=False)) - if other.cdim == 0: - direction = [other.array[1], -other.array[0], 0] if other.isinf else [1, 0, 0] - rays = Segment(np.stack([other.array, direction], axis=-2), copy=False) - else: + if isinstance(other, PointCollection): direction = np.zeros_like(other.array) ind = other.isinf direction[ind, 0] = other.array[ind, 1] direction[ind, 1] = -other.array[ind, 0] direction[~ind, 0] = 1 - rays = Segment(np.stack([other.array, direction], axis=-2), copy=False).expand_dims(-3) + rays = SegmentCollection(np.stack([other.array, direction], axis=-2), copy=False).expand_dims(-3) + else: + direction = [other.array[1], -other.array[0], 0] if other.isinf else [1, 0, 0] + rays = Segment(np.stack([other.array, direction], axis=-2), copy=False) intersections = meet(edges._line, rays._line, _check_dependence=False) @@ -450,7 +485,7 @@ def contains(self, other: Point) -> npt.NDArray[np.bool_]: return result - def intersect(self, other: Line | Segment) -> list[Point]: + def intersect(self, other: LineTensor | SegmentTensor) -> list[PointTensor]: """Intersect the polygon with another object. Args: @@ -463,15 +498,18 @@ def intersect(self, other: Line | Segment) -> list[Point]: if self.dim == 2: return list(distinct(self.edges.intersect(other))) - if isinstance(other, Segment): + if isinstance(other, SegmentTensor): try: result = self._plane.meet(other._line) except LinearDependenceError as e: - if isinstance(other, Segment): - other = cast(Segment, other[~e.dependent_values]) - result = cast(Plane, self._plane[~e.dependent_values]).meet(other._line) + if isinstance(other, SegmentTensor): + other = cast(SegmentTensor, other[~e.dependent_values]) + result = cast(PlaneTensor, self._plane[~e.dependent_values]).meet(other._line) return list( - result[Polygon(self[~e.dependent_values], copy=False).contains(result) & other.contains(result)] + result[ + PolygonCollection.from_tensor(self[~e.dependent_values]).contains(result) + & other.contains(result) + ] ) else: return list(result[self.contains(result) & other.contains(result)]) @@ -479,10 +517,10 @@ def intersect(self, other: Line | Segment) -> list[Point]: try: result = self._plane.meet(other) except LinearDependenceError as e: - if other.cdim > 0: + if other.free_indices > 0: other = other[~e.dependent_values] - result = cast(Plane, self._plane[~e.dependent_values]).meet(other) - return list(result[Polygon(self[~e.dependent_values], copy=False).contains(result)]) + result = cast(PlaneTensor, self._plane[~e.dependent_values]).meet(other) + return list(result[PolygonCollection.from_tensor(self[~e.dependent_values]).contains(result)]) else: return list(result[self.contains(result)]) @@ -492,9 +530,9 @@ def _normalized_projection(self) -> np.ndarray: if self.dim > 2: e = self._plane o = Point(*[0] * self.dim) - if e.cdim > 0: + if e.free_indices > 0: ind = ~e.contains(o) - e[ind] = cast(Plane, e[ind]).parallel(o) + e[ind] = cast(PlaneTensor, e[ind]).parallel(o) elif not e.contains(o): # use parallel hyperplane for projection to avoid rescaling e = e.parallel(o) @@ -510,27 +548,33 @@ def area(self) -> npt.NDArray[np.float_]: a = sum(det(points[..., [0, i, i + 1], :]) for i in range(1, points.shape[-2] - 1)) return 1 / 2 * np.abs(a) - @property - def centroid(self) -> Point: - """The centroid (center of mass) of the polygon.""" - points = self.normalized_array - centroids = [np.average(points[[0, i, i + 1], :-1], axis=0) for i in range(1, points.shape[0] - 1)] - weights = [det(self._normalized_projection()[[0, i, i + 1]]) / 2 for i in range(1, points.shape[0] - 1)] - return Point(*np.average(centroids, weights=weights, axis=0)) - @property def angles(self) -> list[npt.NDArray[np.float_]]: """The interior angles of the polygon.""" result = [] - a = cast(Segment, self.edges[-1]) + a = cast(SegmentTensor, self.edges[-1]) for b in self.edges: - b = cast(Segment, b) + b = cast(SegmentTensor, b) result.append(angle(a.vertices[1], a.vertices[0], b.vertices[1])) a = b return result +class Polygon(PolygonTensor, Polytope): + @property + def centroid(self) -> Point: + """The centroid (center of mass) of the polygon.""" + points = self.normalized_array + centroids = [np.average(points[[0, i, i + 1], :-1], axis=0) for i in range(1, points.shape[0] - 1)] + weights = [det(self._normalized_projection()[[0, i, i + 1]]) / 2 for i in range(1, points.shape[0] - 1)] + return Point(*np.average(centroids, weights=weights, axis=0)) + + +class PolygonCollection(PolygonTensor, PolytopeCollection[Polygon]): + _element_class = Polygon + + class RegularPolygon(Polygon): """A class that can be used to construct regular polygon from a radius and a center point. @@ -543,7 +587,9 @@ class RegularPolygon(Polygon): """ - def __init__(self, center: Point, radius: float, n: int, axis: Point | None = None, **kwargs) -> None: + def __init__( + self, center: Point, radius: float, n: int, axis: Point | None = None, **kwargs: Unpack[NDArrayParameters] + ) -> None: if axis is None: p = Point(1, 0) else: @@ -573,7 +619,7 @@ def center(self) -> Point: @property def inradius(self) -> npt.NDArray[np.float_]: """The inradius of the regular polygon.""" - return dist(self.center, cast(Segment, self.edges[0]).midpoint) + return dist(self.center, cast(SegmentTensor, self.edges[0]).midpoint) class Triangle(Polygon, Simplex): @@ -585,7 +631,7 @@ class Triangle(Polygon, Simplex): """ - def __init__(self, *args, **kwargs): + def __init__(self, *args: Tensor | npt.ArrayLike, **kwargs: Unpack[NDArrayParameters]): super().__init__(*args, **kwargs) assert self.shape[-2] == 3, "Unexpected number of vertices." @@ -597,7 +643,7 @@ def circumcenter(self) -> Point: bisector2 = e2._line.perpendicular(e2.midpoint, plane=self._plane) return bisector1.meet(bisector2) - def contains(self, other: Point) -> npt.NDArray[np.bool_]: + def contains(self, other: PointTensor) -> npt.NDArray[np.bool_]: # faster algorithm using barycentric coordinates # TODO: vectorize @@ -649,12 +695,12 @@ def __init__(self, *args: Tensor | npt.ArrayLike, **kwargs: Unpack[NDArrayParame super().__init__(*args, pdim=3, **kwargs) @property - def faces(self) -> Polygon: + def faces(self) -> PolygonCollection: """The faces of the polyhedron.""" - return Polygon(self.array, copy=False) + return PolygonCollection(self.array, copy=False) @property - def edges(self) -> list[Segment]: + def edges(self) -> list[SegmentTensor]: """The edges of the polyhedron.""" result = self._edges return list(distinct(Segment(result[idx], copy=False) for idx in np.ndindex(*self.shape[:2]))) @@ -664,7 +710,7 @@ def area(self) -> npt.NDArray[np.float_] | np.float_: """The surface area of the polyhedron.""" return np.sum(self.faces.area) - def intersect(self, other: Line | Segment) -> list[Point]: + def intersect(self, other: LineTensor | SegmentTensor) -> list[PointTensor]: """Intersect the polyhedron with another object. Args: diff --git a/geometer/transformation.py b/geometer/transformation.py index 8564334..058dfda 100644 --- a/geometer/transformation.py +++ b/geometer/transformation.py @@ -1,21 +1,42 @@ from __future__ import annotations +from abc import ABC from collections.abc import Sequence -from typing import TYPE_CHECKING, TypeVar, overload +from typing import TYPE_CHECKING, Literal, TypeVar, overload import numpy as np import numpy.typing as npt -from geometer.base import LeviCivitaTensor, ProjectiveTensor, Tensor, TensorDiagram, TensorIndex +from geometer.base import ( + BoundTensor, + LeviCivitaTensor, + ProjectiveTensor, + Tensor, + TensorCollection, + TensorDiagram, +) from geometer.exceptions import NoIncidence -from geometer.point import Line, Point, Subspace, infty_hyperplane +from geometer.point import LineTensor, Point, PointTensor, Subspace, infty_hyperplane from geometer.utils import inv, matmul, outer if TYPE_CHECKING: + from typing_extensions import Unpack + from geometer.curve import Conic + from geometer.utils.typing import NDArrayParameters, TensorIndex + + +@overload +def identity(dim: int, collection_dims: Literal[None] = None) -> Transformation: + ... + + +@overload +def identity(dim: int, collection_dims: tuple[int, ...] | None = None) -> TransformationCollection: + ... -def identity(dim: int, collection_dims: tuple[int, ...] | None = None) -> Transformation: +def identity(dim: int, collection_dims: tuple[int, ...] | None = None) -> TransformationTensor: """Returns the identity transformation. Args: @@ -31,7 +52,7 @@ def identity(dim: int, collection_dims: tuple[int, ...] | None = None) -> Transf e = np.eye(dim + 1) e = e.reshape((1,) * len(collection_dims) + e.shape) e = np.tile(e, (*collection_dims, 1, 1)) - return Transformation(e, copy=False) + return TransformationCollection(e, copy=False) return Transformation(np.eye(dim + 1), copy=False) @@ -68,7 +89,7 @@ def affine_transform(matrix: npt.ArrayLike | None = None, offset: npt.ArrayLike return Transformation(result, copy=False) -def rotation(angle: float, axis: Point | None = None) -> Transformation: +def rotation(angle: float | np.float_, axis: Point | None = None) -> Transformation: """Returns a projective transformation that represents a rotation by the specified angle (and axis). Args: @@ -122,7 +143,7 @@ def scaling(*factors: npt.ArrayLike) -> Transformation: The scaling transformation. """ - return affine_transform(np.diag(factors)) + return affine_transform(np.diag(factors)) # type: ignore[arg-type] def reflection(axis: Subspace) -> Transformation: @@ -142,7 +163,7 @@ def reflection(axis: Subspace) -> Transformation: return identity(axis.dim) v = axis.array[:-1] - v = v / np.linalg.norm(v) + v = v / np.linalg.norm(v) # type: ignore[operator] p = affine_transform(np.eye(axis.dim) - 2 * outer(v, v.conj())) @@ -154,7 +175,7 @@ def reflection(axis: Subspace) -> Transformation: return translation(x) * p * translation(-x) -class Transformation(ProjectiveTensor): +class TransformationTensor(ProjectiveTensor, ABC): """Represents a projective transformation in an arbitrary projective space. The underlying array is the matrix representation of the projective transformation. The matrix must be @@ -167,15 +188,19 @@ class Transformation(ProjectiveTensor): """ - def __init__(self, *args: Tensor | npt.ArrayLike, **kwargs) -> None: + def __init__(self, *args: Tensor | npt.ArrayLike, **kwargs: Unpack[NDArrayParameters]) -> None: kwargs.setdefault("covariant", [0]) super().__init__(*args, tensor_rank=2, **kwargs) + if self.tensor_shape != (1, 1): + raise ValueError(f"Expected tensor of type (1, 1), but is {self.tensor_shape}") + if self.shape[-1] != self.shape[-2]: + raise ValueError(f"Expected quadratic matrix, but last two dimensions are {self.shape[-2:]}") - def __apply__(self, transformation: Transformation) -> Transformation: - return Transformation(matmul(transformation.array, self.array), copy=False) + def __apply__(self, transformation: TransformationTensor) -> TransformationTensor: + return TransformationCollection.from_array(matmul(transformation.array, self.array)) @classmethod - def from_points(cls, *args: tuple[Point, Point]) -> Transformation: + def from_points(cls, *args: tuple[PointTensor, PointTensor]) -> TransformationTensor: """Constructs a projective transformation in n-dimensional projective space from the image of n + 2 points in general position. @@ -196,17 +221,17 @@ def from_points(cls, *args: tuple[Point, Point]) -> Transformation: b = [y.array for x, y in args] m1 = np.column_stack(a[:-1]) m2 = np.column_stack(b[:-1]) - d1 = np.linalg.solve(m1, a[-1]) - d2 = np.linalg.solve(m2, b[-1]) + d1 = np.linalg.solve(m1, a[-1]) # type: ignore[arg-type] + d2 = np.linalg.solve(m2, b[-1]) # type: ignore[arg-type] t1 = m1.dot(np.diag(d1)) t2 = m2.dot(np.diag(d2)) - return Transformation(t2.dot(np.linalg.inv(t1))) + return cls(t2.dot(np.linalg.inv(t1))) @classmethod def from_points_and_conics( - cls, points1: Sequence[Point], points2: Sequence[Point], conic1: Conic, conic2: Conic - ) -> Transformation: + cls, points1: Sequence[PointTensor], points2: Sequence[PointTensor], conic1: Conic, conic2: Conic + ) -> TransformationTensor: """Constructs a projective transformation from two conics and the image of pairs of 3 points on the conics. Args: @@ -225,9 +250,9 @@ def from_points_and_conics( a1, b1, c1 = points1 l1, l2 = conic1.tangent(a1), conic1.tangent(b1) - if not isinstance(l1, Line): + if not isinstance(l1, LineTensor): raise NoIncidence(f"Point {a1} does not lie on the conic {conic1}.") - if not isinstance(l2, Line): + if not isinstance(l2, LineTensor): raise NoIncidence(f"Point {b1} does not lie on the conic {conic1}.") m = l1.meet(l2).join(c1) @@ -237,16 +262,16 @@ def from_points_and_conics( a2, b2, c2 = points2 l1, l2 = conic2.tangent(a2), conic2.tangent(b2) - if not isinstance(l1, Line): + if not isinstance(l1, LineTensor): raise NoIncidence(f"Point {a2} does not lie on the conic {conic2}.") - if not isinstance(l2, Line): + if not isinstance(l2, LineTensor): raise NoIncidence(f"Point {b2} does not lie on the conic {conic2}.") m = l1.meet(l2).join(c2) p, q = conic2.intersect(m) d2 = p if q == c2 else q - return Transformation.from_points((a1, a2), (b1, b2), (c1, c2), (d1, d2)) + return cls.from_points((a1, a2), (b1, b2), (c1, c2), (d1, d2)) T = TypeVar("T", bound=Tensor) @@ -265,7 +290,7 @@ def apply(self, other: T) -> T: raise NotImplementedError(f"Object of type {type(other)} cannot be transformed.") @overload - def __mul__(self, other: Transformation) -> Transformation: + def __mul__(self, other: TransformationTensor) -> TransformationTensor: ... @overload @@ -284,14 +309,14 @@ def __mul__(self, other: Tensor | npt.ArrayLike) -> Tensor: def __pow__(self, power: int, modulo: int | None = None) -> Tensor: if power == 0: - if self.cdim == 0: + if self.free_indices == 0: return identity(self.dim) - return identity(self.dim, self.shape[: self.cdim]) + return identity(self.dim, self.shape[: self.free_indices]) if power < 0: return self.inverse().__pow__(-power, modulo) result = super().__pow__(power, modulo) - return Transformation(result, copy=False) + return type(self)(result, copy=False) def __getitem__(self, index: TensorIndex) -> Tensor | np.generic: result = super().__getitem__(index) @@ -299,13 +324,21 @@ def __getitem__(self, index: TensorIndex) -> Tensor | np.generic: if not isinstance(result, Tensor) or result.tensor_shape != (1, 1): return result - return Transformation(result, copy=False) + return TransformationCollection.from_tensor(result) - def inverse(self) -> Transformation: + def inverse(self) -> TransformationTensor: """Calculates the inverse projective transformation. Returns: The inverse transformation. """ - return Transformation(inv(self.array), copy=False) + return type(self)(inv(self.array), copy=False) + + +class Transformation(TransformationTensor, BoundTensor): + pass + + +class TransformationCollection(TransformationTensor, TensorCollection[Transformation]): + _element_class = Transformation diff --git a/geometer/utils/math.py b/geometer/utils/math.py index d6f89d9..cfdcfcb 100644 --- a/geometer/utils/math.py +++ b/geometer/utils/math.py @@ -2,18 +2,16 @@ import math from numbers import Number -from typing import TYPE_CHECKING, Literal, TypedDict, Union +from typing import TYPE_CHECKING, Literal, TypedDict import numpy as np import numpy.typing as npt from numpy.lib.scimath import sqrt as csqrt if TYPE_CHECKING: - from typing_extensions import TypeAlias, TypeGuard, Unpack + from typing_extensions import TypeGuard, Unpack -NumericalDType: TypeAlias = Union[np.number, np.bool_] -NumericalArray: TypeAlias = npt.NDArray[NumericalDType] -NumericalScalar: TypeAlias = Union[Number, np.number, np.bool_, NumericalArray] # TODO: restrict array shape + from geometer.utils.typing import NumericalDType, NumericalScalar def is_numerical_scalar(element: npt.ArrayLike) -> TypeGuard[NumericalScalar]: @@ -238,10 +236,10 @@ def det(A: npt.ArrayLike) -> npt.NDArray[np.number]: n = A.shape[-1] if n == 2: - return A[..., 0, 0] * A[..., 1, 1] - A[..., 1, 0] * A[..., 0, 1] # type: ignore[no-any-return] + return A[..., 0, 0] * A[..., 1, 1] - A[..., 1, 0] * A[..., 0, 1] if n == 3 and A.size >= 9 * 64: - return ( # type: ignore[no-any-return] + return ( A[..., 0, 0] * A[..., 1, 1] * A[..., 2, 2] + A[..., 0, 1] * A[..., 1, 2] * A[..., 2, 0] + A[..., 0, 2] * A[..., 1, 0] * A[..., 2, 1] @@ -250,7 +248,7 @@ def det(A: npt.ArrayLike) -> npt.NDArray[np.number]: - A[..., 2, 2] * A[..., 1, 0] * A[..., 0, 1] ) - return np.linalg.det(A) # type: ignore[no-any-return] + return np.linalg.det(A) def inv(A: npt.ArrayLike) -> npt.NDArray[np.number]: @@ -329,7 +327,7 @@ def orth(A: npt.ArrayLike, dim: int | None = None) -> npt.NDArray[np.number]: raise ValueError("Cannot calculate the image spaces of matrices when the spaces have different dimensions.") dim = dims.flat[0] - return u[..., :dim] # type: ignore[no-any-return] + return u[..., :dim] class UFuncParameters(TypedDict, total=False): diff --git a/geometer/utils/ops_dispatch.py b/geometer/utils/ops_dispatch.py new file mode 100644 index 0000000..30b5ced --- /dev/null +++ b/geometer/utils/ops_dispatch.py @@ -0,0 +1,123 @@ +# Adapted from https://github.com/pandas-dev/pandas/blob/main/pandas/_libs/ops_dispatch.pyx +# See license at https://github.com/pandas-dev/pandas/blob/main/LICENSE +from typing import Any, Literal + +import numpy as np + +DISPATCHED_UFUNCS = { + "add", + "sub", + "mul", + "pow", + "mod", + "floordiv", + "truediv", + "divmod", + "eq", + "ne", + "lt", + "gt", + "le", + "ge", + "remainder", + "matmul", + "or", + "xor", + "and", + "neg", + "pos", + "abs", +} +UNARY_UFUNCS = { + "neg", + "pos", + "abs", +} +UFUNC_ALIASES = { + "subtract": "sub", + "multiply": "mul", + "floor_divide": "floordiv", + "true_divide": "truediv", + "power": "pow", + "remainder": "mod", + "divide": "truediv", + "equal": "eq", + "not_equal": "ne", + "less": "lt", + "less_equal": "le", + "greater": "gt", + "greater_equal": "ge", + "bitwise_or": "or", + "bitwise_and": "and", + "bitwise_xor": "xor", + "negative": "neg", + "absolute": "abs", + "positive": "pos", +} + +# For op(., Array) -> Array.__r{op}__ +REVERSED_NAMES = { + "lt": "__gt__", + "le": "__ge__", + "gt": "__lt__", + "ge": "__le__", + "eq": "__eq__", + "ne": "__ne__", +} + + +def maybe_dispatch_ufunc_to_dunder_op( + obj: Any, + ufunc: np.ufunc, + method: Literal["__call__", "reduce", "reduceat", "accumulate", "outer", "inner"], + *inputs: Any, + **kwargs: Any, +) -> Any: + """ + Dispatch a ufunc to the equivalent dunder method. + + Args: + obj: The object whose dunder method we dispatch to. + ufunc: A NumPy ufunc. + method: How the ufunc was called. + inputs: The input arrays. + kwargs: The additional keyword arguments, e.g. ``out``. + + Returns: + The result of applying the ufunc + """ + # special has the ufuncs we dispatch to the dunder op on + + op_name = ufunc.__name__ + op_name = UFUNC_ALIASES.get(op_name, op_name) + + def not_implemented(*args, **kwargs): # type: ignore # noqa: ARG001 + return NotImplemented + + if kwargs or ufunc.nin > 2: + return NotImplemented + + if method == "__call__" and op_name in DISPATCHED_UFUNCS: + if inputs[0] is obj: + name = f"__{op_name}__" + meth = getattr(obj, name, not_implemented) + + if op_name in UNARY_UFUNCS: + assert len(inputs) == 1 + return meth() + + return meth(inputs[1]) + + elif inputs[1] is obj: + name = REVERSED_NAMES.get(op_name, f"__r{op_name}__") + + meth = getattr(obj, name, not_implemented) + result = meth(inputs[0]) + return result + + else: + # should not be reached, but covering our bases + return NotImplemented + + else: + return NotImplemented diff --git a/geometer/utils/typing.py b/geometer/utils/typing.py new file mode 100644 index 0000000..bc47ce2 --- /dev/null +++ b/geometer/utils/typing.py @@ -0,0 +1,30 @@ +from __future__ import annotations + +from collections.abc import Sequence +from numbers import Number +from typing import TYPE_CHECKING, Literal, TypedDict, Union + +import numpy as np +from numpy import typing as npt + +if TYPE_CHECKING: + from typing_extensions import TypeAlias + +IntegerIndex1D: TypeAlias = Union[int, np.int_, slice, Sequence[int], Sequence[np.int_], npt.NDArray[np.int_]] +BooleanIndex1D: TypeAlias = Union[bool, np.bool_, slice, Sequence[bool], Sequence[np.bool_], npt.NDArray[np.bool_]] +TensorIndex: TypeAlias = Union[IntegerIndex1D, BooleanIndex1D, tuple[IntegerIndex1D, ...], tuple[BooleanIndex1D, ...]] +Shape: TypeAlias = tuple[int, ...] + + +class NDArrayParameters(TypedDict, total=False): + dtype: npt.DTypeLike + copy: bool + order: Literal["K", "A", "C", "F"] + subok: bool + ndim: int + like: npt.ArrayLike + + +NumericalDType: TypeAlias = Union[np.number, np.bool_] +NumericalArray: TypeAlias = npt.NDArray[NumericalDType] +NumericalScalar: TypeAlias = Union[Number, np.number, np.bool_, NumericalArray] # TODO: restrict array shape diff --git a/pyproject.toml b/pyproject.toml index b654dfe..79b4498 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -113,6 +113,7 @@ plugins = ["numpy.typing.mypy_plugin"] strict = true implicit_reexport = true disallow_any_generics = false # TODO +warn_return_any = false # TODO [[tool.mypy.overrides]] module = "geometer.utils.indexing" diff --git a/tests/test_base.py b/tests/test_base.py index b9e9d44..c013038 100644 --- a/tests/test_base.py +++ b/tests/test_base.py @@ -1,6 +1,6 @@ import numpy as np -from geometer.base import KroneckerDelta, LeviCivitaTensor, Tensor, TensorDiagram +from geometer.base import KroneckerDelta, LeviCivitaTensor, Tensor, TensorCollection, TensorDiagram class TestTensor: @@ -42,39 +42,26 @@ def test_dtype(self): a = Tensor(2, 3, dtype=np.complex64) assert a.dtype == np.complex64 - def test_expand_dims(self): - a = Tensor([[1, 2], [3, 4]], covariant=[0]) - b = a.expand_dims(0) - - assert b.shape == (1, 2, 2) - assert b[0] == a - assert b.squeeze(0) == a - assert b.squeeze(0).tensor_shape == a.tensor_shape - class TestTensorCollection: def test_init(self): - # empty list - a = Tensor([], tensor_rank=1) - assert len(a) == 0 - # numpy array - a = Tensor(np.ones((1, 2, 3)), tensor_rank=1) + a = TensorCollection[Tensor](np.ones((1, 2, 3)), tensor_rank=1) assert len(a) == 1 assert a.size == 2 # nested list of numbers - a = Tensor([[1, 2], [3, 4]], tensor_rank=1) + a = TensorCollection[Tensor]([[1, 2], [3, 4]], tensor_rank=1) assert len(a) == 2 assert a.size == 2 # nested tuple of numbers - a = Tensor(((1, 2), (3, 4)), tensor_rank=1) + a = TensorCollection[Tensor](((1, 2), (3, 4)), tensor_rank=1) assert len(a) == 2 assert a.size == 2 # nested list of Tensor objects - a = Tensor([[Tensor(1, 2, 3), Tensor(3, 4, 5)]], tensor_rank=1) + a = TensorCollection[Tensor]([[Tensor(1, 2, 3), Tensor(3, 4, 5)]], tensor_rank=1) assert a.shape == (1, 2, 3) assert len(a) == 1 assert a.size == 2 @@ -84,14 +71,14 @@ class A: def __array__(self): return np.array([Tensor(1, 2), Tensor(3, 4)]) - a = Tensor(A(), tensor_rank=1) + a = TensorCollection[Tensor](A(), tensor_rank=1) assert len(a) == 2 assert a.size == 2 def test_getitem(self): - a = Tensor([[1, 2], [3, 4]]) - b = Tensor([[5, 6], [7, 8]]) - c = Tensor([a, b], tensor_rank=1) + a = TensorCollection[Tensor]([[1, 2], [3, 4]]) + b = TensorCollection[Tensor]([[5, 6], [7, 8]]) + c = TensorCollection[Tensor]([a, b], tensor_rank=1) assert c[0] == a assert c[1] == b @@ -99,6 +86,23 @@ def test_getitem(self): assert c[:, 1] == Tensor([Tensor([3, 4]), Tensor([7, 8])], tensor_rank=1) assert c[:, 0, 0] == [1, 5] + def test_expand_dims(self): + a = TensorCollection[Tensor]([[1, 2], [3, 4]]) + b = a.expand_dims(0) + + assert b.shape == (1, 2, 2) + assert b[0] == a + + c = a.expand_dims(-2) + + assert c.shape == (2, 1, 2) + assert c[:, 0, :] == a + + d = a.expand_dims(-3) + + assert d.shape == (1, 2, 2) + assert d[0] == a + class TestTensorDiagram: def test_add_edge(self):