diff options
Diffstat (limited to '')
| -rw-r--r-- | mdpoly/__init__.py | 14 | ||||
| -rw-r--r-- | mdpoly/abc.py | 21 | ||||
| -rw-r--r-- | mdpoly/algebra.py | 701 | ||||
| -rw-r--r-- | mdpoly/errors.py | 3 | ||||
| -rw-r--r-- | mdpoly/expressions/__init__.py | 66 | ||||
| -rw-r--r-- | mdpoly/expressions/matrix.py | 175 | ||||
| -rw-r--r-- | mdpoly/expressions/poly.py | 169 | ||||
| -rw-r--r-- | mdpoly/index.py | 13 | ||||
| -rw-r--r-- | mdpoly/operations/__init__.py | 0 | ||||
| -rw-r--r-- | mdpoly/operations/add.py | 81 | ||||
| -rw-r--r-- | mdpoly/operations/derivative.py | 62 | ||||
| -rw-r--r-- | mdpoly/operations/exp.py | 52 | ||||
| -rw-r--r-- | mdpoly/operations/mul.py | 186 | ||||
| -rw-r--r-- | mdpoly/operations/transpose.py | 20 | ||||
| -rw-r--r-- | mdpoly/representations.py | 12 | ||||
| -rw-r--r-- | mdpoly/state.py | 1 | ||||
| -rw-r--r-- | mdpoly/util.py | 9 |
17 files changed, 854 insertions, 731 deletions
diff --git a/mdpoly/__init__.py b/mdpoly/__init__.py index 15ab8c9..1ac17a5 100644 --- a/mdpoly/__init__.py +++ b/mdpoly/__init__.py @@ -61,9 +61,9 @@ Expression Within the code, we will often call something like :math:`x^2 + 1` and *expression* (for example :py:class:`mdpoly.abc.Expr` or -:py:class:`mdpoly.algebra.PolyRingExpr`) instead of polynomial. Expressions are -more general and can include any mathematical operation. Of course in the end -and expression should result in a polynomial. +:py:class:`mdpoly.expressions.poly.PolyRingExpr`) instead of polynomial. +Expressions are more general and can include any mathematical operation. Of +course in the end and expression should result in a polynomial. Decision Variables, Program @@ -107,12 +107,8 @@ TODO: data structure(s) that represent the polynomial from .abc import (Shape as _Shape) -from .algebra import (PolyConst as _PolyConst, - PolyVar as _PolyVar, - PolyParam as _PolyParam, - MatConst as _MatConst, - MatVar as _MatVar, - MatParam as _MatParam) +from .expressions.poly import (PolyConst as _PolyConst, PolyVar as _PolyVar, PolyParam as _PolyParam) +from .expressions.matrix import (MatConst as _MatConst, MatVar as _MatVar, MatParam as _MatParam) from .state import State as _State diff --git a/mdpoly/abc.py b/mdpoly/abc.py index 16b155b..cd90f85 100644 --- a/mdpoly/abc.py +++ b/mdpoly/abc.py @@ -2,23 +2,21 @@ from __future__ import annotations from typing import TYPE_CHECKING -from .index import Number, Shape, MatrixIndex, PolyIndex -from .constants import NUMERICS_EPS -from .errors import AlgebraicError -from .util import iszero - -if TYPE_CHECKING: - from .state import State - from typing import Self, Type, TypeVar, Generic, Any, Iterable, Sequence from enum import Enum, auto from copy import copy from abc import ABC, abstractmethod -from dataclassabc import dataclassabc from hashlib import sha256 +from dataclassabc import dataclassabc +from .constants import NUMERICS_EPS +from .errors import AlgebraicError +from .util import iszero + +if TYPE_CHECKING: + from .index import Number, Shape, MatrixIndex, PolyIndex + from .state import State -ReprT = TypeVar("ReprT", bound="Repr") # ┏━╸╻ ╻┏━┓┏━┓┏━╸┏━┓┏━┓╻┏━┓┏┓╻┏━┓ # ┣╸ ┏╋┛┣━┛┣┳┛┣╸ ┗━┓┗━┓┃┃ ┃┃┗┫┗━┓ @@ -380,6 +378,9 @@ class Rel(ABC): # ╹┗╸┗━╸╹ ╹┗╸┗━╸┗━┛┗━╸╹ ╹ ╹ ╹ ╹ ╹ ╹┗━┛╹ ╹┗━┛ +ReprT = TypeVar("ReprT", bound="Repr") + + class Repr(ABC): r""" Representation of a multivariate matrix polynomial expression. diff --git a/mdpoly/algebra.py b/mdpoly/algebra.py deleted file mode 100644 index 03a6fd5..0000000 --- a/mdpoly/algebra.py +++ /dev/null @@ -1,701 +0,0 @@ -""" Algebraic Structures for Expressions """ -from __future__ import annotations - -from .abc import Algebra, Expr, ReprT, Nothing, Const, Var, Param -from .errors import AlgebraicError, InvalidShape, MissingParameters -from .state import State -from .index import Shape, MatrixIndex, PolyIndex, PolyVarIndex, Number - -from typing import cast, Sequence, Iterable, Type, TypeVar, Self -from functools import reduce -from itertools import product, chain, combinations_with_replacement -from abc import abstractmethod -from dataclasses import field, dataclass -import operator - -from dataclassabc import dataclassabc - - -@dataclassabc(eq=False) -class BinaryOp(Expr): - left: Expr = field() - right: Expr = field() - - -@dataclassabc(eq=False) -class UnaryOp(Expr): - right: Expr = field() - - @property - def inner(self) -> Expr: - """ Inner expression on which the operator is acting, alias for right. """ - return self.right - - @property - def left(self) -> Expr: - return Nothing() - - @left.setter - def left(self, left) -> None: - if not isinstance(left, Nothing): - raise ValueError("Cannot set left of left-acting unary operator " - "to something that is not of type Nothing.") - - -class Reducible(Expr): - """ Reducible Expression. - - Algebraic expression that can be written in terms of other (existing) - expression objects, i.e. can be reduced to another expression made of - simpler operations. For example subtraction can be written in term of - addition and multiplication. - """ - - def to_repr(self, repr_type: Type[ReprT], state: State) -> tuple[ReprT, State]: - """ See :py:meth:`mdpoly.abc.Expr.to_repr`. """ - return self.reduce().to_repr(repr_type, state) - - @abstractmethod - def reduce(self) -> Expr: - """ Reduce the expression to its basic elements """ - - -# TODO: review this idea before implementing -# class Simplifiable(Expr): -# """ Simplifiable Expression. """ -# -# @abstractmethod -# def simplify(self) -> Expr: -# """ Simplify the expression """ - - -# ┏━┓┏━┓╻ ╻ ╻┏┓╻┏━┓┏┳┓╻┏━┓╻ ┏━┓╻ ┏━╸┏━╸┏┓ ┏━┓┏━┓ -# ┣━┛┃ ┃┃ ┗┳┛┃┗┫┃ ┃┃┃┃┃┣━┫┃ ┣━┫┃ ┃╺┓┣╸ ┣┻┓┣┳┛┣━┫ -# ╹ ┗━┛┗━╸ ╹ ╹ ╹┗━┛╹ ╹╹╹ ╹┗━╸ ╹ ╹┗━╸┗━┛┗━╸┗━┛╹┗╸╹ ╹ - - -class PolyRingExpr(Expr): - r""" Expression with the algebraic behaviour of a polynomial ring. - - This is the algebra of :math:`\mathbb{R}[x_1, \ldots, x_n]`. Note that the - polynomials are scalars. - """ - - # -- Properties --- - - @property - def algebra(self) -> Algebra: - return Algebra.poly_ring - - @property - def shape(self) -> Shape: - """ See :py:meth:`mdpoly.abc.Expr.shape`. """ - if self.left.shape != self.right.shape: - raise InvalidShape(f"Cannot perform operation {repr(self)} with " - f"shapes {self.left.shape} and {self.right.shape}.") - return self.left.shape - - # --- Helper methods for construction --- - - @classmethod - def from_powers(cls, variable: PolyRingExpr, exponents: Iterable[int]) -> PolyRingExpr: - """ Given a variable, say :math:`x`, and a list of exponents ``[0, 3, 5]`` - returns a polynomial :math:`x^5 + x^3 + 1`. - """ - nonzero_exponents = filter(lambda e: e != 0, exponents) - # Ignore typecheck because dataclassabc has no typing stub - constant_term = PolyConst(1 if 0 in exponents else 0) # type: ignore[call-arg] - # FIXME: remove annoying 0 - return sum((variable ** e for e in nonzero_exponents), constant_term) - - - @classmethod - def make_combinations(cls, variables: Iterable[PolyVar], max_degree: int) -> Iterable[PolyRingExpr]: - """ Make combinations of terms. - - For example given :math:`x, y` and *max_degree* of 2 generates - :math:`x^2, xy, x, y^2, y, 1`. - """ - # Ignore typecheck because dataclassabc has no typing stub - vars_and_const = chain(variables, (PolyConst(1),)) # type: ignore[call-arg] - for comb in combinations_with_replacement(vars_and_const, max_degree): - yield cast(PolyRingExpr, reduce(operator.mul, comb)) - - # --- Operator Overloading --- - - def __add__(self, other): - self._assert_same_algebra(self, other) - other = self._wrap(Number, PolyConst, other) - return PolyAdd(self, other) - - def __radd__(self, other): - self._assert_same_algebra(self, other) - other = self._wrap(Number, PolyConst, other) - return PolyAdd(other, self) - - def __sub__(self, other): - self._assert_same_algebra(self, other) - other = self._wrap(Number, PolyConst, other) - return PolySub(self, other) - - def __rsub__(self, other): - self._assert_same_algebra(self, other) - other = self._wrap(Number, PolyConst, other) - return PolyAdd(other, self) - - def __neg__(self): - # FIXME: Create PolyNeg? - return PolyMul(PolyConst(-1), self) - - def __mul__(self, other): - self._assert_same_algebra(self, other) - other = self._wrap(Number, PolyConst, other) - return PolyMul(self, other) - - def __rmul__(self, other): - self._assert_same_algebra(self, other) - other = self._wrap(Number, PolyConst, other) - return PolyMul(other, self) - - def __matmul__(self, other): - raise AlgebraicError("Cannot perform matrix multiplication in polynomial ring (they are scalars).") - - def __rmatmul__(self, other): - self.__rmatmul__(other) - - def __truediv__(self, other): - other = self._wrap(Number, PolyConst, other) - if not isinstance(other, PolyConst | PolyParam): - raise AlgebraicError("Cannot divide by variables in polynomial ring.") - - return PolyMul(self, other) - - def __rtruediv__(self, other): - raise AlgebraicError("Cannot perform right division in polynomial ring.") - - def __pow__(self, other): - other = self._wrap(Number, PolyConst, other) - if not isinstance(other, PolyConst | PolyParam): - raise AlgebraicError(f"Cannot raise to powers of type {type(other)} in " - "polynomial ring. Only constants and parameters are allowed.") - return PolyExp(left=self, right=other) - - # -- Other mathematical operations --- - - def diff(self, wrt: PolyVar) -> PolyPartialDiff: - return PolyPartialDiff(self, wrt) - - -@dataclassabc(frozen=True) -class PolyVar(Var, PolyRingExpr): - """ Variable TODO: desc """ - name: str # overloads Leaf.name - shape: Shape = Shape.scalar() # ovearloads PolyRingExpr.shape - - def to_repr(self, repr_type: Type[ReprT], state: State) -> tuple[ReprT, State]: - r = repr_type(self.shape) - idx = PolyVarIndex.from_var(self, state), # important comma! - r.set(MatrixIndex.scalar(), PolyIndex(idx), 1) - return r, state - - -@dataclassabc(frozen=True) -class PolyConst(Const, PolyRingExpr): - """ Constant TODO: desc """ - value: Number # overloads Const.value - name: str = "" # overloads Leaf.name - shape: Shape = Shape.scalar() # ovearloads PolyRingExpr.shape - - def to_repr(self, repr_type: Type[ReprT], state: State) -> tuple[ReprT, State]: - r = repr_type(self.shape) - r.set(MatrixIndex.scalar(), PolyIndex.constant(), self.value) - return r, state - - -@dataclassabc(frozen=True) -class PolyParam(Param, PolyRingExpr): - """ Polynomial parameter TODO: desc """ - name: str # overloads Leaf.name - shape: Shape = Shape.scalar() # overloads PolyRingExpr.shape - - def to_repr(self, repr_type: Type[ReprT], state: State) -> tuple[ReprT, State]: - if self not in state.parameters: - raise MissingParameters("Cannot construct representation because " - f"value for parameter {self} was not given.") - - return PolyConst(state.parameters[self]).to_repr(repr_type, state) # type: ignore[call-arg] - - -@dataclass(eq=False) -class PolyAdd(BinaryOp, PolyRingExpr): - """ Addition operator between scalar polynomials. """ - - def to_repr(self, repr_type: Type[ReprT], state: State) -> tuple[ReprT, State]: - """ See :py:meth:`mdpoly.abc.Expr.to_repr`. """ - # Make a new empty representation - r = repr_type(self.shape) - - # Make representations for existing stuff - for node in self.children(): - nrepr, state = node.to_repr(repr_type, state) - - # Process non-zero terms - for entry in nrepr.entries(): - for term in nrepr.terms(entry): - s = r.at(entry, term) + nrepr.at(entry, term) - r.set(entry, term, s) - - return r, state - - def __str__(self) -> str: - return f"({self.left} + {self.right})" - - -@dataclass(eq=False) -class PolySub(BinaryOp, PolyRingExpr, Reducible): - """ Subtraction operator between scalar polynomials. """ - - def reduce(self) -> Expr: - """ See :py:meth:`mdpoly.algebra.Reducible.reduce`. """ - return self.left + (-1 * self.right) - - def __str__(self) -> str: - return f"({self.left} - {self.right})" - - -@dataclass(eq=False) -class PolyMul(BinaryOp, PolyRingExpr): - """ Multiplication operator between scalar polynomials. """ - - def to_repr(self, repr_type: Type[ReprT], state: State) -> tuple[ReprT, State]: - """ See :py:meth:`mdpoly.abc.Expr.to_repr`. """ - r = repr_type(self.shape) - - lrepr, state = self.left.to_repr(repr_type, state) - rrepr, state = self.right.to_repr(repr_type, state) - - # Non zero entries are the intersection since if either is zero the - # result is zero - nonzero_entries = set(lrepr.entries()) & set(rrepr.entries()) - for entry in nonzero_entries: - # Compute polynomial product between non-zero entries - for lterm, rterm in product(lrepr.terms(entry), rrepr.terms(entry)): - # Compute where the results should go - term = PolyIndex.product(lterm, rterm) - - # Compute product - p = r.at(entry, term) + lrepr.at(entry, lterm) * rrepr.at(entry, rterm) - r.set(entry, term, p) - - return r, state - - def __str__(self) -> str: - return f"({self.left} * {self.right})" - - -@dataclass(eq=False) -class PolyExp(BinaryOp, PolyRingExpr, Reducible): - """ Exponentiation operator between scalar polynomials. """ - - @property # type: ignore[override] - def right(self) -> Const: # type: ignore[override] - if not isinstance(super().right, Const): - raise AlgebraicError(f"Cannot raise {self.left} to {self.right} because" - f"{self.right} is not a constant.") - - return cast(Const, super().right) - - @right.setter - def right(self, right: Expr) -> None: - # Workaround necessary because of a bug - # https://bugs.python.org/issue14965 - BinaryOp.right.fset(self, right) # type: ignore - - def reduce(self) -> Expr: - """ See :py:meth:`mdpoly.algebra.Reducible.reduce`. """ - var = self.left - if not isinstance(self.right.value, int): - raise NotImplementedError("Cannot raise to non-integer powers (yet).") - - ntimes = self.right.value - 1 - return reduce(operator.mul, (var for _ in range(ntimes)), var) - - def __str__(self) -> str: - return f"({self.left} ** {self.right})" - - -@dataclassabc(eq=False) -class PolyPartialDiff(UnaryOp, PolyRingExpr): - """ Partial differentiation of scalar polynomials. """ - wrt: PolyVar - - @property - def shape(self) -> Shape: - """ See :py:meth:`mdpoly.abc.Expr.shape`. """ - return self.inner.shape - - def to_repr(self, repr_type: Type[ReprT], state: State) -> tuple[ReprT, State]: - """ See :py:meth:`mdpoly.abc.Expr.to_repr`. """ - r = repr_type(self.shape) - lrepr, state = self.inner.to_repr(repr_type, state) - - entry = MatrixIndex.scalar() - wrt = state.index(self.wrt) - - for term in lrepr.terms(entry): - # sadly, walrus operator does not support unpacking - # https://bugs.python.org/issue43143 - # if ((newterm, fac) := PolyIndex.differentiate(term, wrt)): - if result := PolyIndex.differentiate(term, wrt): - newterm, fac = result - r.set(entry, newterm, fac * lrepr.at(entry, term)) - - return r, state - - def __str__(self) -> str: - return f"(∂_{self.wrt} {self.inner})" - - def replace(self, old: PolyRingExpr, new: PolyRingExpr) -> Self: - """ Overloads :py:meth:`mdpoly.abc.Expr.replace` """ - if self.wrt == old: - if not isinstance(new, PolyVar): - # FIXME: implement chain rule - raise AlgebraicError(f"Cannot take a derivative with respect to {new}.") - - self.wrt = cast(PolyVar, new) - - return cast(Self, Expr.replace(self, old, new)) - - -# ┏┳┓┏━┓╺┳╸┏━┓╻╻ ╻ ┏━┓╻ ┏━╸┏━╸┏┓ ┏━┓┏━┓ -# ┃┃┃┣━┫ ┃ ┣┳┛┃┏╋┛ ┣━┫┃ ┃╺┓┣╸ ┣┻┓┣┳┛┣━┫ -# ╹ ╹╹ ╹ ╹ ╹┗╸╹╹ ╹ ╹ ╹┗━╸┗━┛┗━╸┗━┛╹┗╸╹ ╹ - -class MatrixExpr(Expr): - r""" Expression with the algebraic properties of a matrix ring and / or - module (depending on the shape). - - We denote with :math:`R` a polynomial ring. - - - If the shape is square, like :math:`(n, n)` then this is :math:`M_n(R)` - the ring of matrices over :math:`R`. - - - If the shape is something else like a row or column (:math:`(m, p)`, - :math:`(1, n)` or :math:`(n, 1)`) this is a module, i.e. an algebra with - addition and scalar multiplication, where the "scalars" come from - :math:`R`. - - Furthermore some operators that are usually expected from matrices are - already included (eg. transposition). - """ - - @property - def algebra(self) -> Algebra: - return Algebra.matrix_ring - - def at(self, entry: MatrixIndex, scalar_type: type) -> Expr: - raise NotImplementedError - - def __add__(self, other): - self._assert_same_algebra(self, other) - other = self._wrap(Number, MatConst, other) - return MatAdd(self, other) - - def __sub__(self, other): - self._assert_same_algebra(self, other) - other = self._wrap(Number, MatConst, other) - return MatSub(self, other) - - def __rsub__(self, other): - self._assert_same_algebra(self, other) - other = self._wrap(Number, MatConst, other) - return MatSub(other, self) - - def __mul__(self, other): - self._assert_same_algebra(self, other) - other = self._wrap(Number, MatConst, other) - return MatScalarMul(other, self) - - def __rmul__(self, other): - self._assert_same_algebra(self, other) - other = self._wrap(Number, MatConst, other) - return MatScalarMul(other, self) - - def __matmul__(self, other): - self._assert_same_algebra(self, other) - other = self._wrap(Number, MatConst, other) - return MatMul(self, other) - - def __rmatmul(self, other): - self._assert_same_algebra(self, other) - other = self._wrap(Number, MatConst, other) - return MatMul(other, self) - - def __truediv__(self, scalar): - scalar = self._wrap_if_constant(scalar) - raise NotImplementedError - - def transpose(self) -> MatrixExpr: - """ Matrix transposition. """ - raise NotImplementedError - return MatTranspose(self) - - @property - def T(self) -> MatrixExpr: - """ Shorthand for :py:meth:`mdpoly.algebra.MatrixExpr.transpose`. """ - return self.transpose() - - def to_scalar(self, scalar_type: type): - """ Convert to a scalar expression. """ - raise NotImplementedError - - -@dataclassabc(frozen=True) -class MatConst(Const, MatrixExpr): - """ Matrix constant. TODO: desc. """ - value: Sequence[Sequence[Number]] # Row major, overloads Const.value - shape: Shape # overloads Expr.shape - name: str = "" # overloads Leaf.name - - def to_repr(self, repr_type: Type[ReprT], state: State) -> tuple[ReprT, State]: - r = repr_type(self.shape) - - for i, row in enumerate(self.value): - for j, val in enumerate(row): - r.set(MatrixIndex(row=i, col=j), PolyIndex.constant(), val) - - return r, state - - - def __str__(self) -> str: - if not self.name: - return repr(self.value) - - return self.name - - - -T = TypeVar("T", bound=Var) - -@dataclassabc(frozen=True) -class MatVar(Var, MatrixExpr): - """ Matrix of polynomial variables. TODO: desc """ - name: str # overloads Leaf.name - shape: Shape # overloads Expr.shape - - # TODO: review this API, can be moved elsewhere? - def to_scalars(self, scalar_var_type: Type[T]) -> Iterable[tuple[MatrixIndex, T]]: - for row in range(self.shape.rows): - for col in range(self.shape.cols): - var = scalar_var_type(name=f"{self.name}_[{row},{col}]") # type: ignore[call-arg] - entry = MatrixIndex(row, col) - - yield entry, var - - def to_repr(self, repr_type: Type[ReprT], state: State) -> tuple[ReprT, State]: - r = repr_type(self.shape) - - # FIXME: do not hardcode scalar type - for entry, var in self.to_scalars(PolyVar): - idx = PolyVarIndex.from_var(var, state), # important comma! - r.set(entry, PolyIndex(idx), 1) - - return r, state - - def __str__(self) -> str: - return self.name - - -@dataclassabc(frozen=True) -class MatParam(Param, MatrixExpr): - """ Matrix parameter. TODO: desc. """ - name: str - shape: Shape - - def to_repr(self, repr_type: Type[ReprT], state: State) -> tuple[ReprT, State]: - if self not in state.parameters: - raise MissingParameters("Cannot construct representation because " - f"value for parameter {self} was not given.") - - # FIXME: add conversion to scalar variables - # Ignore typecheck because dataclassabc has not type stub - return MatConst(state.parameters[self]).to_repr(repr_type, state) # type: ignore[call-arg] - - def __str__(self) -> str: - return self.name - - -@dataclassabc -class MatAdd(BinaryOp, MatrixExpr): - """ Addition operator between matrices. """ - - @property - def shape(self) -> Shape: - """ See :py:meth:`mdpoly.abc.Expr.shape`. """ - if not self.left.shape == self.right.shape: - raise AlgebraicError(f"Cannot add matrices {self.left} and {self.right} " - "with different shapes.") - return self.left.shape - - - def to_repr(self, repr_type: Type[ReprT], state: State) -> tuple[ReprT, State]: - """ See :py:meth:`mdpoly.abc.Expr.to_repr`. """ - # Make a new empty representation - r = repr_type(self.shape) - - # Make representations for existing stuff - for node in self.children(): - nrepr, state = node.to_repr(repr_type, state) - - # Process non-zero terms - for entry in nrepr.entries(): - for term in nrepr.terms(entry): - s = r.at(entry, term) + nrepr.at(entry, term) - r.set(entry, term, s) - - return r, state - - def __str__(self) -> str: - return f"({self.left} + {self.right})" - - -@dataclassabc -class MatSub(BinaryOp, MatrixExpr, Reducible): - """ Subtraction operator between matrices. """ - - @property - def shape(self) -> Shape: - """ See :py:meth:`mdpoly.abc.Expr.shape`. """ - if not self.left.shape == self.right.shape: - raise AlgebraicError(f"Cannot subtract matrices {self.left} and {self.right} " - f"with different shapes, {self.left.shape} and {self.right.shape}.") - return self.left.shape - - def reduce(self) -> Expr: - """ See :py:meth:`mdpoly.algebra.Reducible.reduce`. """ - return self.left + (-1 * self.right) - - def __str__(self) -> str: - return f"({self.left} - {self.right})" - - -@dataclassabc -class MatElemMul(BinaryOp, MatrixExpr): - """ Elementwise Matrix Multiplication. """ - - @property - def shape(self) -> Shape: - """ See :py:meth:`mdpoly.abc.Expr.shape`. """ - if not self.left.shape == self.right.shape: - raise AlgebraicError("Cannot perform element-wise multiplication of matrices " - f"{self.left} and {self.right} with different shapes, " - f"{self.left.shape} and {self.right.shape}") - return self.left.shape - - def to_repr(self, repr_type: Type[ReprT], state: State) -> tuple[ReprT, State]: - """ See :py:meth:`mdpoly.abc.Expr.to_repr`. """ - r = repr_type(self.shape) - - lrepr, state = self.left.to_repr(repr_type, state) - rrepr, state = self.right.to_repr(repr_type, state) - - # Non zero entries are the intersection since if either is zero the - # result is zero - nonzero_entries = set(lrepr.entries()) & set(rrepr.entries()) - for entry in nonzero_entries: - # Compute polynomial product between non-zero entries - for lterm, rterm in product(lrepr.terms(entry), rrepr.terms(entry)): - # Compute where the results should go - term = PolyIndex.product(lterm, rterm) - - # Compute product - p = r.at(entry, term) + lrepr.at(entry, lterm) * rrepr.at(entry, rterm) - r.set(entry, term, p) - - return r, state - - def __str__(self) -> str: - return f"({self.left} .* {self.right})" - - -@dataclassabc -class MatScalarMul(BinaryOp, MatrixExpr): - """ Matrix-Scalar Multiplication. Assumes scalar is on the left and matrix - on the right. """ - - @property - def shape(self) -> Shape: - """ See :py:meth:`mdpoly.abc.Expr.shape`. """ - if not self.left.shape == Shape.scalar(): - raise InvalidShape(f"Matrix-scalar product assumes that left argumet {self.left} " - f"but it has shape {self.left.shape}") - - return self.right.shape - - - def to_repr(self, repr_type: Type[ReprT], state: State) -> tuple[ReprT, State]: - """ See :py:meth:`mdpoly.abc.Expr.to_repr`. """ - r = repr_type(self.shape) - - scalar_repr, state = self.left.to_repr(repr_type, state) - mat_repr, state = self.right.to_repr(repr_type, state) - - for entry in mat_repr.entries(): - scalar_terms = scalar_repr.terms(MatrixIndex.scalar()) - mat_terms = mat_repr.terms(entry) - for scalar_term, mat_term in product(scalar_terms, mat_terms): - term = PolyIndex.product(scalar_term, mat_term) - - p = r.at(entry, term) + scalar_repr.at(entry, scalar_term) + mat_repr.at(entry, mat_term) - r.set(entry, term, p) - - return r, state - - def __str__(self) -> str: - return f"({self.left} * {self.right})" - - -@dataclass(eq=False) -class MatDotProd(BinaryOp, MatrixExpr): - """ Dot product. """ - - @property - def shape(self) -> Shape: - if not self.left.shape.is_row(): - raise AlgebraicError(f"Left operand {self.left} must be a row!") - - if not self.right.shape.is_col(): - raise AlgebraicError(f"Right operand {self.right} must be a column!") - - if self.left.shape.cols != self.right.shape.rows: - raise AlgebraicError(f"Rows of {self.right} and columns {self.left} do not match!") - - return Shape.scalar() - - -@dataclass(eq=False) -class MatMul(BinaryOp, MatrixExpr): - """ Matrix Multiplication. """ - - @property - def shape(self) -> Shape: - """ See :py:meth:`mdpoly.abc.Expr.shape`. """ - if not self.left.shape.rows == self.right.shape.cols: - raise AlgebraicError("Cannot perform matrix multiplication between " - f"{self.left} and {self.right} (shapes {self.left.shape} and {self.right.shape})") - - return Shape(self.left.shape.rows, self.right.shape.cols) - - - def __str__(self) -> str: - return f"({self.left} @ {self.right})" - - -@dataclass(eq=False) -class MatTranspose(UnaryOp, MatrixExpr): - """ Matrix transposition """ - - @property - def shape(self) -> Shape: - """ See :py:meth:`mdpoly.abc.Expr.shape`. """ - return Shape(self.inner.shape.cols, self.inner.shape.rows) diff --git a/mdpoly/errors.py b/mdpoly/errors.py index 85432b8..262f27e 100644 --- a/mdpoly/errors.py +++ b/mdpoly/errors.py @@ -1,7 +1,6 @@ - class AlgebraicError(Exception): """ This is raised whenever an invalid mathematical operation is performed. - See also :py:mod:`mdpoly.algebra`. + See also :py:mod:`mdpoly.operations`. """ class InvalidShape(Exception): diff --git a/mdpoly/expressions/__init__.py b/mdpoly/expressions/__init__.py new file mode 100644 index 0000000..2efe342 --- /dev/null +++ b/mdpoly/expressions/__init__.py @@ -0,0 +1,66 @@ +from __future__ import annotations +from typing import TYPE_CHECKING + +from abc import abstractmethod +from typing import Type +from dataclasses import field +from dataclassabc import dataclassabc + +from ..abc import Expr, Nothing + +if TYPE_CHECKING: + from ..abc import ReprT + from ..state import State + + +@dataclassabc(eq=False) +class BinaryOp(Expr): + left: Expr = field() + right: Expr = field() + + +@dataclassabc(eq=False) +class UnaryOp(Expr): + right: Expr = field() + + @property + def inner(self) -> Expr: + """ Inner expression on which the operator is acting, alias for right. """ + return self.right + + @property + def left(self) -> Expr: + return Nothing() + + @left.setter + def left(self, left) -> None: + if not isinstance(left, Nothing): + raise ValueError("Cannot set left of left-acting unary operator " + "to something that is not of type Nothing.") + + +class Reducible(Expr): + """ Reducible Expression. + + Algebraic expression that can be written in terms of other (existing) + expression objects, i.e. can be reduced to another expression made of + simpler operations. For example subtraction can be written in term of + addition and multiplication. + """ + + def to_repr(self, repr_type: Type[ReprT], state: State) -> tuple[ReprT, State]: + """ See :py:meth:`mdpoly.abc.Expr.to_repr`. """ + return self.reduce().to_repr(repr_type, state) + + @abstractmethod + def reduce(self) -> Expr: + """ Reduce the expression to its basic elements """ + + +# TODO: review this idea before implementing +# class Simplifiable(Expr): +# """ Simplifiable Expression. """ +# +# @abstractmethod +# def simplify(self) -> Expr: +# """ Simplify the expression """ diff --git a/mdpoly/expressions/matrix.py b/mdpoly/expressions/matrix.py new file mode 100644 index 0000000..65f86c2 --- /dev/null +++ b/mdpoly/expressions/matrix.py @@ -0,0 +1,175 @@ +from __future__ import annotations +from typing import TYPE_CHECKING + +from typing import Type, TypeVar, Iterable, Sequence +from dataclassabc import dataclassabc + +from .poly import PolyVar, PolyConst +from ..abc import Expr, Var, Const, Param, Algebra +from ..index import MatrixIndex, PolyIndex, PolyVarIndex +from ..errors import MissingParameters + +from .. import operations + +if TYPE_CHECKING: + from ..abc import ReprT + from ..index import Shape, Number + from ..state import State + + +class MatrixExpr(Expr): + r""" Expression with the algebraic properties of a matrix ring and / or + module (depending on the shape). + + We denote with :math:`R` a polynomial ring. + + - If the shape is square, like :math:`(n, n)` then this is :math:`M_n(R)` + the ring of matrices over :math:`R`. + + - If the shape is something else like a row or column (:math:`(m, p)`, + :math:`(1, n)` or :math:`(n, 1)`) this is a module, i.e. an algebra with + addition and scalar multiplication, where the "scalars" come from + :math:`R`. + + Furthermore some operators that are usually expected from matrices are + already included (eg. transposition). + """ + + @property + def algebra(self) -> Algebra: + return Algebra.matrix_ring + + def __add__(self, other): + self._assert_same_algebra(self, other) + other = self._wrap(Number, MatConst, other) + return operations.add.MatAdd(self, other) + + def __sub__(self, other): + self._assert_same_algebra(self, other) + other = self._wrap(Number, MatConst, other) + return operations.add.MatSub(self, other) + + def __rsub__(self, other): + self._assert_same_algebra(self, other) + other = self._wrap(Number, MatConst, other) + return operations.add.MatSub(other, self) + + def __neg__(self): + # FIXME: Create PolyNeg? + return operations.mul.MatScalarMul(PolyConst(-1), self) + + def __mul__(self, other): + self._assert_same_algebra(self, other) + other = self._wrap(Number, MatConst, other) + + # TODO: case distiction based on shapes + return operations.mul.MatScalarMul(other, self) + + + def __rmul__(self, other): + self._assert_same_algebra(self, other) + other = self._wrap(Number, MatConst, other) + return operations.mul.MatScalarMul(other, self) + + def __matmul__(self, other): + self._assert_same_algebra(self, other) + other = self._wrap(Number, MatConst, other) + return operations.MatMul(self, other) + + def __rmatmul(self, other): + self._assert_same_algebra(self, other) + other = self._wrap(Number, MatConst, other) + return operations.mul.MatMul(other, self) + + def __truediv__(self, scalar): + scalar = self._wrap_if_constant(scalar) + raise NotImplementedError + + def transpose(self) -> MatrixExpr: + """ Matrix transposition. """ + raise NotImplementedError + return operations.transpose.MatTranspose(self) + + @property + def T(self) -> MatrixExpr: + """ Shorthand for :py:meth:`mdpoly.expressions.matrix.MatrixExpr.transpose`. """ + return self.transpose() + + def to_scalar(self, scalar_type: type): + """ Convert to a scalar expression. """ + raise NotImplementedError + + +@dataclassabc(frozen=True) +class MatConst(Const, MatrixExpr): + """ Matrix constant. TODO: desc. """ + value: Sequence[Sequence[Number]] # Row major, overloads Const.value + shape: Shape # overloads Expr.shape + name: str = "" # overloads Leaf.name + + def to_repr(self, repr_type: Type[ReprT], state: State) -> tuple[ReprT, State]: + r = repr_type(self.shape) + + for i, row in enumerate(self.value): + for j, val in enumerate(row): + r.set(MatrixIndex(row=i, col=j), PolyIndex.constant(), val) + + return r, state + + + def __str__(self) -> str: + if not self.name: + return repr(self.value) + + return self.name + + + +T = TypeVar("T", bound=Var) + +@dataclassabc(frozen=True) +class MatVar(Var, MatrixExpr): + """ Matrix of polynomial variables. TODO: desc """ + name: str # overloads Leaf.name + shape: Shape # overloads Expr.shape + + # TODO: review this API, can be moved elsewhere? + def to_scalars(self, scalar_var_type: Type[T]) -> Iterable[tuple[MatrixIndex, T]]: + for row in range(self.shape.rows): + for col in range(self.shape.cols): + var = scalar_var_type(name=f"{self.name}_[{row},{col}]") # type: ignore[call-arg] + entry = MatrixIndex(row, col) + + yield entry, var + + def to_repr(self, repr_type: Type[ReprT], state: State) -> tuple[ReprT, State]: + r = repr_type(self.shape) + + # FIXME: do not hardcode scalar type + for entry, var in self.to_scalars(PolyVar): + idx = PolyVarIndex.from_var(var, state), # important comma! + r.set(entry, PolyIndex(idx), 1) + + return r, state + + def __str__(self) -> str: + return self.name + + +@dataclassabc(frozen=True) +class MatParam(Param, MatrixExpr): + """ Matrix parameter. TODO: desc. """ + name: str + shape: Shape + + def to_repr(self, repr_type: Type[ReprT], state: State) -> tuple[ReprT, State]: + if self not in state.parameters: + raise MissingParameters("Cannot construct representation because " + f"value for parameter {self} was not given.") + + # FIXME: add conversion to scalar variables + # Ignore typecheck because dataclassabc has not type stub + return MatConst(state.parameters[self]).to_repr(repr_type, state) # type: ignore[call-arg] + + def __str__(self) -> str: + return self.name diff --git a/mdpoly/expressions/poly.py b/mdpoly/expressions/poly.py new file mode 100644 index 0000000..3fe83ca --- /dev/null +++ b/mdpoly/expressions/poly.py @@ -0,0 +1,169 @@ +from __future__ import annotations +from typing import TYPE_CHECKING + +from typing import Type, Iterable +from functools import reduce +from itertools import chain, combinations_with_replacement +from dataclassabc import dataclassabc + +from ..abc import Expr, Var, Const, Param, Algebra, ReprT +from ..index import Shape, MatrixIndex, PolyIndex, PolyVarIndex, Number +from ..errors import AlgebraicError, MissingParameters, InvalidShape + +from .. import operations + +if TYPE_CHECKING: + from ..index import Number + from ..state import State + + +class PolyRingExpr(Expr): + r""" Expression with the algebraic behaviour of a polynomial ring. + + This is the algebra of :math:`\mathbb{R}[x_1, \ldots, x_n]`. Note that the + polynomials are scalars. + """ + + # -- Properties --- + + @property + def algebra(self) -> Algebra: + return Algebra.poly_ring + + @property + def shape(self) -> Shape: + """ See :py:meth:`mdpoly.abc.Expr.shape`. """ + if self.left.shape != self.right.shape: + raise InvalidShape(f"Cannot perform operation {repr(self)} with " + f"shapes {self.left.shape} and {self.right.shape}.") + return self.left.shape + + # --- Helper methods for construction --- + + @classmethod + def from_powers(cls, variable: PolyRingExpr, exponents: Iterable[int]) -> PolyRingExpr: + """ Given a variable, say :math:`x`, and a list of exponents ``[0, 3, 5]`` + returns a polynomial :math:`x^5 + x^3 + 1`. + """ + nonzero_exponents = filter(lambda e: e != 0, exponents) + # Ignore typecheck because dataclassabc has no typing stub + constant_term = PolyConst(1 if 0 in exponents else 0) # type: ignore[call-arg] + # FIXME: remove annoying 0 + return sum((variable ** e for e in nonzero_exponents), constant_term) + + + @classmethod + def make_combinations(cls, variables: Iterable[PolyVar], max_degree: int) -> Iterable[PolyRingExpr]: + """ Make combinations of terms. + + For example given :math:`x, y` and *max_degree* of 2 generates + :math:`x^2, xy, x, y^2, y, 1`. + """ + # Ignore typecheck because dataclassabc has no typing stub + vars_and_const = chain(variables, (PolyConst(1),)) # type: ignore[call-arg] + for comb in combinations_with_replacement(vars_and_const, max_degree): + yield cast(PolyRingExpr, reduce(operator.mul, comb)) + + # --- Operator Overloading --- + + def __add__(self, other): + self._assert_same_algebra(self, other) + other = self._wrap(Number, PolyConst, other) + return operations.add.PolyAdd(self, other) + + def __radd__(self, other): + self._assert_same_algebra(self, other) + other = self._wrap(Number, PolyConst, other) + return operations.add.PolyAdd(other, self) + + def __sub__(self, other): + self._assert_same_algebra(self, other) + other = self._wrap(Number, PolyConst, other) + return operations.add.PolySub(self, other) + + def __rsub__(self, other): + self._assert_same_algebra(self, other) + other = self._wrap(Number, PolyConst, other) + return operations.add.PolyAdd(other, self) + + def __neg__(self): + # FIXME: Create PolyNeg? + return operations.mul.PolyMul(PolyConst(-1), self) + + def __mul__(self, other): + self._assert_same_algebra(self, other) + other = self._wrap(Number, PolyConst, other) + return operations.mul.PolyMul(self, other) + + def __rmul__(self, other): + self._assert_same_algebra(self, other) + other = self._wrap(Number, PolyConst, other) + return operations.mul.PolyMul(other, self) + + def __matmul__(self, other): + raise AlgebraicError("Cannot perform matrix multiplication in polynomial ring (they are scalars).") + + def __rmatmul__(self, other): + self.__rmatmul__(other) + + def __truediv__(self, other): + other = self._wrap(Number, PolyConst, other) + if not isinstance(other, PolyConst | PolyParam): + raise AlgebraicError("Cannot divide by variables in polynomial ring.") + + return operations.mul.PolyMul(self, other) + + def __rtruediv__(self, other): + raise AlgebraicError("Cannot perform right division in polynomial ring.") + + def __pow__(self, other): + other = self._wrap(Number, PolyConst, other) + if not isinstance(other, PolyConst | PolyParam): + raise AlgebraicError(f"Cannot raise to powers of type {type(other)} in " + "polynomial ring. Only constants and parameters are allowed.") + return operations.exp.PolyExp(left=self, right=other) + + # -- Other mathematical operations --- + + def diff(self, wrt: PolyVar) -> operations.derivative.PolyPartialDiff: + return operations.derivative.PolyPartialDiff(right=self, wrt=wrt) + + +@dataclassabc(frozen=True) +class PolyVar(Var, PolyRingExpr): + """ Variable TODO: desc """ + name: str # overloads Leaf.name + shape: Shape = Shape.scalar() # ovearloads PolyRingExpr.shape + + def to_repr(self, repr_type: Type[ReprT], state: State) -> tuple[ReprT, State]: + r = repr_type(self.shape) + idx = PolyVarIndex.from_var(self, state), # important comma! + r.set(MatrixIndex.scalar(), PolyIndex(idx), 1) + return r, state + + +@dataclassabc(frozen=True) +class PolyConst(Const, PolyRingExpr): + """ Constant TODO: desc """ + value: Number # overloads Const.value + name: str = "" # overloads Leaf.name + shape: Shape = Shape.scalar() # ovearloads PolyRingExpr.shape + + def to_repr(self, repr_type: Type[ReprT], state: State) -> tuple[ReprT, State]: + r = repr_type(self.shape) + r.set(MatrixIndex.scalar(), PolyIndex.constant(), self.value) + return r, state + + +@dataclassabc(frozen=True) +class PolyParam(Param, PolyRingExpr): + """ Polynomial parameter TODO: desc """ + name: str # overloads Leaf.name + shape: Shape = Shape.scalar() # overloads PolyRingExpr.shape + + def to_repr(self, repr_type: Type[ReprT], state: State) -> tuple[ReprT, State]: + if self not in state.parameters: + raise MissingParameters("Cannot construct representation because " + f"value for parameter {self} was not given.") + + return PolyConst(state.parameters[self]).to_repr(repr_type, state) # type: ignore[call-arg] diff --git a/mdpoly/index.py b/mdpoly/index.py index 5211245..070412c 100644 --- a/mdpoly/index.py +++ b/mdpoly/index.py @@ -1,15 +1,17 @@ from __future__ import annotations +from typing import TYPE_CHECKING + +from typing import Self, NamedTuple, Optional, Sequence +from itertools import filterfalse from .errors import InvalidShape from .util import partition, isclose -from itertools import filterfalse -from typing import Self, NamedTuple, Optional, Sequence, TYPE_CHECKING - if TYPE_CHECKING: from .state import State from .abc import Var +# FIXME: this is not really an appropriate place for this -> constants? Number = int | float @@ -206,6 +208,11 @@ class PolyIndex(tuple[PolyVarIndex]): return cls(PolyVarIndex(k, v) for k, v in d.items()) @classmethod + def from_expr(cls, product_of_vars) -> Self: + # TODO: implement product of variables + raise NotImplementedError + + @classmethod def constant(cls) -> Self: """ Index of the constant term. """ return cls((PolyVarIndex.constant(),)) diff --git a/mdpoly/operations/__init__.py b/mdpoly/operations/__init__.py new file mode 100644 index 0000000..e69de29 --- /dev/null +++ b/mdpoly/operations/__init__.py diff --git a/mdpoly/operations/add.py b/mdpoly/operations/add.py new file mode 100644 index 0000000..bf0c5a0 --- /dev/null +++ b/mdpoly/operations/add.py @@ -0,0 +1,81 @@ +from __future__ import annotations +from typing import TYPE_CHECKING + +from typing import Type +from dataclassabc import dataclassabc + +from ..abc import Expr +from ..errors import AlgebraicError + +from ..expressions import BinaryOp, Reducible +from ..expressions.matrix import MatrixExpr +from ..expressions.poly import PolyRingExpr + +if TYPE_CHECKING: + from ..abc import ReprT + from ..index import Shape + from ..state import State + + +@dataclassabc(eq=False) +class MatAdd(BinaryOp, MatrixExpr): + """ Addition operator between matrices. """ + + @property + def shape(self) -> Shape: + """ See :py:meth:`mdpoly.abc.Expr.shape`. """ + if not self.left.shape == self.right.shape: + raise AlgebraicError(f"Cannot add matrices {self.left} and {self.right} " + "with different shapes.") + return self.left.shape + + + def to_repr(self, repr_type: Type[ReprT], state: State) -> tuple[ReprT, State]: + """ See :py:meth:`mdpoly.abc.Expr.to_repr`. """ + # Make a new empty representation + r = repr_type(self.shape) + + # Make representations for existing stuff + for node in self.children(): + nrepr, state = node.to_repr(repr_type, state) + + # Process non-zero terms + for entry in nrepr.entries(): + for term in nrepr.terms(entry): + s = r.at(entry, term) + nrepr.at(entry, term) + r.set(entry, term, s) + + return r, state + + def __str__(self) -> str: + return f"({self.left} + {self.right})" + + +@dataclassabc(eq=False) +class MatSub(BinaryOp, MatrixExpr, Reducible): + """ Subtraction operator between matrices. """ + + @property + def shape(self) -> Shape: + """ See :py:meth:`mdpoly.abc.Expr.shape`. """ + if not self.left.shape == self.right.shape: + raise AlgebraicError(f"Cannot subtract matrices {self.left} and {self.right} " + f"with different shapes, {self.left.shape} and {self.right.shape}.") + return self.left.shape + + def reduce(self) -> Expr: + """ See :py:meth:`mdpoly.expressions.Reducible.reduce`. """ + return self.left + (-1 * self.right) + + def __str__(self) -> str: + return f"({self.left} - {self.right})" + + +@dataclassabc(eq=False) +class PolyAdd(MatAdd, PolyRingExpr): + ... + + +@dataclassabc(eq=False) +class PolySub(MatSub, PolyRingExpr): + ... diff --git a/mdpoly/operations/derivative.py b/mdpoly/operations/derivative.py new file mode 100644 index 0000000..6e02b0e --- /dev/null +++ b/mdpoly/operations/derivative.py @@ -0,0 +1,62 @@ +from __future__ import annotations +from typing import TYPE_CHECKING + +from typing import Type, Self, cast +from dataclassabc import dataclassabc + +from ..abc import Expr +from ..errors import AlgebraicError +from ..index import Shape, MatrixIndex, PolyIndex + +from ..expressions import UnaryOp +from ..expressions.poly import PolyRingExpr, PolyVar + +if TYPE_CHECKING: + from ..abc import ReprT + from ..state import State + + +# TODO: implement vector derivatives + + +@dataclassabc(eq=False) +class PolyPartialDiff(UnaryOp, PolyRingExpr): + """ Partial differentiation of scalar polynomials. """ + wrt: PolyVar + + @property + def shape(self) -> Shape: + """ See :py:meth:`mdpoly.abc.Expr.shape`. """ + return self.inner.shape + + def to_repr(self, repr_type: Type[ReprT], state: State) -> tuple[ReprT, State]: + """ See :py:meth:`mdpoly.abc.Expr.to_repr`. """ + r = repr_type(self.shape) + lrepr, state = self.inner.to_repr(repr_type, state) + + entry = MatrixIndex.scalar() + wrt = state.index(self.wrt) + + for term in lrepr.terms(entry): + # sadly, walrus operator does not support unpacking + # https://bugs.python.org/issue43143 + # if ((newterm, fac) := PolyIndex.differentiate(term, wrt)): + if result := PolyIndex.differentiate(term, wrt): + newterm, fac = result + r.set(entry, newterm, fac * lrepr.at(entry, term)) + + return r, state + + def __str__(self) -> str: + return f"(∂_{self.wrt} {self.inner})" + + def replace(self, old: PolyRingExpr, new: PolyRingExpr) -> Self: + """ Overloads :py:meth:`mdpoly.abc.Expr.replace` """ + if self.wrt == old: + if not isinstance(new, PolyVar): + # FIXME: implement chain rule + raise AlgebraicError(f"Cannot take a derivative with respect to {new}.") + + self.wrt = cast(PolyVar, new) + + return cast(Self, Expr.replace(self, old, new)) diff --git a/mdpoly/operations/exp.py b/mdpoly/operations/exp.py new file mode 100644 index 0000000..fbc67f5 --- /dev/null +++ b/mdpoly/operations/exp.py @@ -0,0 +1,52 @@ +from __future__ import annotations +from typing import TYPE_CHECKING + +from operator import mul as opmul +from typing import cast +from functools import reduce +from dataclasses import dataclass + +from ..abc import Expr, Const +from ..errors import AlgebraicError + +from ..expressions import BinaryOp, Reducible +from ..expressions.matrix import MatrixExpr +from ..expression.poly import PolyRingExpr + + +# TODO: implement matrix exponential, use caley-hamilton thm magic +@dataclass(eq=False) +class MatExp(BinaryOp, MatrixExpr, Reducible): + def __init__(self): + raise NotImplementedError + + +@dataclass(eq=False) +class PolyExp(BinaryOp, PolyRingExpr, Reducible): + """ Exponentiation operator between scalar polynomials. """ + + @property # type: ignore[override] + def right(self) -> Const: # type: ignore[override] + if not isinstance(super().right, Const): + raise AlgebraicError(f"Cannot raise {self.left} to {self.right} because" + f"{self.right} is not a constant.") + + return cast(Const, super().right) + + @right.setter + def right(self, right: Expr) -> None: + # Workaround necessary because of a bug + # https://bugs.python.org/issue14965 + BinaryOp.right.fset(self, right) # type: ignore + + def reduce(self) -> Expr: + """ See :py:meth:`mdpoly.expressions.Reducible.reduce`. """ + var = self.left + if not isinstance(self.right.value, int): + raise NotImplementedError("Cannot raise to non-integer powers (yet).") + + ntimes = self.right.value - 1 + return reduce(opmul, (var for _ in range(ntimes)), var) + + def __str__(self) -> str: + return f"({self.left} ** {self.right})" diff --git a/mdpoly/operations/mul.py b/mdpoly/operations/mul.py new file mode 100644 index 0000000..a2a3cc4 --- /dev/null +++ b/mdpoly/operations/mul.py @@ -0,0 +1,186 @@ +from __future__ import annotations +from typing import TYPE_CHECKING + +from typing import Type +from itertools import product +from dataclasses import dataclass +from dataclassabc import dataclassabc + +from ..index import Shape +from ..errors import AlgebraicError, InvalidShape +from ..index import MatrixIndex, PolyIndex + +from ..expressions import BinaryOp, Reducible +from ..expressions.matrix import MatrixExpr +from ..expression.poly import PolyRingExpr + +if TYPE_CHECKING: + from ..abc import ReprT + from ..state import State + + +# ┏┳┓┏━┓╺┳╸┏━┓╻╻ ╻ ┏━┓┏━┓┏━┓╺┳┓╻ ╻┏━╸╺┳╸┏━┓ +# ┃┃┃┣━┫ ┃ ┣┳┛┃┏╋┛ ┣━┛┣┳┛┃ ┃ ┃┃┃ ┃┃ ┃ ┗━┓ +# ╹ ╹╹ ╹ ╹ ╹┗╸╹╹ ╹ ╹ ╹┗╸┗━┛╺┻┛┗━┛┗━╸ ╹ ┗━┛ + + +@dataclassabc +class MatElemMul(BinaryOp, MatrixExpr): + """ Elementwise Matrix Multiplication. """ + + @property + def shape(self) -> Shape: + """ See :py:meth:`mdpoly.abc.Expr.shape`. """ + if not self.left.shape == self.right.shape: + raise AlgebraicError("Cannot perform element-wise multiplication of matrices " + f"{self.left} and {self.right} with different shapes, " + f"{self.left.shape} and {self.right.shape}") + return self.left.shape + + def to_repr(self, repr_type: Type[ReprT], state: State) -> tuple[ReprT, State]: + """ See :py:meth:`mdpoly.abc.Expr.to_repr`. """ + r = repr_type(self.shape) + + lrepr, state = self.left.to_repr(repr_type, state) + rrepr, state = self.right.to_repr(repr_type, state) + + # Non zero entries are the intersection since if either is zero the + # result is zero + nonzero_entries = set(lrepr.entries()) & set(rrepr.entries()) + for entry in nonzero_entries: + # Compute polynomial product between non-zero entries + for lterm, rterm in product(lrepr.terms(entry), rrepr.terms(entry)): + # Compute where the results should go + term = PolyIndex.product(lterm, rterm) + + # Compute product + p = r.at(entry, term) + lrepr.at(entry, lterm) * rrepr.at(entry, rterm) + r.set(entry, term, p) + + return r, state + + def __str__(self) -> str: + return f"({self.left} .* {self.right})" + + +@dataclassabc +class MatScalarMul(BinaryOp, MatrixExpr): + """ Matrix-Scalar Multiplication. Assumes scalar is on the left and matrix + on the right. """ + + @property + def shape(self) -> Shape: + """ See :py:meth:`mdpoly.abc.Expr.shape`. """ + if not self.left.shape == Shape.scalar(): + raise InvalidShape(f"Matrix-scalar product assumes that left argumet {self.left} " + f"but it has shape {self.left.shape}") + + return self.right.shape + + + def to_repr(self, repr_type: Type[ReprT], state: State) -> tuple[ReprT, State]: + """ See :py:meth:`mdpoly.abc.Expr.to_repr`. """ + r = repr_type(self.shape) + + scalar_repr, state = self.left.to_repr(repr_type, state) + mat_repr, state = self.right.to_repr(repr_type, state) + + for entry in mat_repr.entries(): + scalar_terms = scalar_repr.terms(MatrixIndex.scalar()) + mat_terms = mat_repr.terms(entry) + for scalar_term, mat_term in product(scalar_terms, mat_terms): + term = PolyIndex.product(scalar_term, mat_term) + + p = r.at(entry, term) + scalar_repr.at(entry, scalar_term) + mat_repr.at(entry, mat_term) + r.set(entry, term, p) + + return r, state + + def __str__(self) -> str: + return f"({self.left} * {self.right})" + + +@dataclass(eq=False) +class MatMul(BinaryOp, MatrixExpr): + """ Matrix Multiplication. """ + + @property + def shape(self) -> Shape: + """ See :py:meth:`mdpoly.abc.Expr.shape`. """ + if not self.left.shape.rows == self.right.shape.cols: + raise AlgebraicError("Cannot perform matrix multiplication between " + f"{self.left} and {self.right} (shapes {self.left.shape} and {self.right.shape})") + + return Shape(self.left.shape.rows, self.right.shape.cols) + + def to_repr(self, repr_type: Type[ReprT], state: State) -> tuple[ReprT, State]: + """ See :py:meth:`mdpoly.abc.Expr.to_repr`. """ + r = repr_type(self.shape) + + lrepr, state = self.left.to_repr(repr_type, state) + rrepr, state = self.right.to_repr(repr_type, state) + + # Compute matrix product + for row in range(self.left.shape.rows): + for col in range(self.right.shape.cols): + for i in range(self.left.shape.cols): + raise NotImplementedError + + return r, state + + + def __str__(self) -> str: + return f"({self.left} @ {self.right})" + + +@dataclass(eq=False) +class MatDotProd(BinaryOp, MatrixExpr, Reducible): + """ Dot product. """ + + @property + def shape(self) -> Shape: + if not self.left.shape.is_row(): + raise AlgebraicError(f"Left operand {self.left} must be a row!") + + if not self.right.shape.is_col(): + raise AlgebraicError(f"Right operand {self.right} must be a column!") + + if self.left.shape.cols != self.right.shape.rows: + raise AlgebraicError(f"Rows of {self.right} and columns {self.left} do not match!") + + return Shape.scalar() + + +# ┏━┓┏━┓╻ ╻ ╻┏┓╻┏━┓┏┳┓╻┏━┓╻ ┏━┓┏━┓┏━┓╺┳┓╻ ╻┏━╸╺┳╸ +# ┣━┛┃ ┃┃ ┗┳┛┃┗┫┃ ┃┃┃┃┃┣━┫┃ ┣━┛┣┳┛┃ ┃ ┃┃┃ ┃┃ ┃ +# ╹ ┗━┛┗━╸ ╹ ╹ ╹┗━┛╹ ╹╹╹ ╹┗━╸ ╹ ╹┗╸┗━┛╺┻┛┗━┛┗━╸ ╹ + + +@dataclass(eq=False) +class PolyMul(BinaryOp, PolyRingExpr): + """ Multiplication operator between scalar polynomials. """ + + def to_repr(self, repr_type: Type[ReprT], state: State) -> tuple[ReprT, State]: + """ See :py:meth:`mdpoly.abc.Expr.to_repr`. """ + r = repr_type(self.shape) + + lrepr, state = self.left.to_repr(repr_type, state) + rrepr, state = self.right.to_repr(repr_type, state) + + # Non zero entries are the intersection since if either is zero the + # result is zero + nonzero_entries = set(lrepr.entries()) & set(rrepr.entries()) + for entry in nonzero_entries: + # Compute polynomial product between non-zero entries + for lterm, rterm in product(lrepr.terms(entry), rrepr.terms(entry)): + # Compute where the results should go + term = PolyIndex.product(lterm, rterm) + + # Compute product + p = r.at(entry, term) + lrepr.at(entry, lterm) * rrepr.at(entry, rterm) + r.set(entry, term, p) + + return r, state + + def __str__(self) -> str: + return f"({self.left} * {self.right})" diff --git a/mdpoly/operations/transpose.py b/mdpoly/operations/transpose.py new file mode 100644 index 0000000..750d0c0 --- /dev/null +++ b/mdpoly/operations/transpose.py @@ -0,0 +1,20 @@ +from __future__ import annotations +from typing import TYPE_CHECKING + +from dataclasses import dataclass + +from ..expressions import UnaryOp +from ..expressions.matrix import MatrixExpr + +if TYPE_CHECKING: + from ..index import Shape + + +@dataclass(eq=False) +class MatTranspose(UnaryOp, MatrixExpr): + """ Matrix transposition """ + + @property + def shape(self) -> Shape: + """ See :py:meth:`mdpoly.abc.Expr.shape`. """ + return Shape(self.inner.shape.cols, self.inner.shape.rows) diff --git a/mdpoly/representations.py b/mdpoly/representations.py index 8268c42..32c75d7 100644 --- a/mdpoly/representations.py +++ b/mdpoly/representations.py @@ -1,10 +1,14 @@ -from .abc import Repr -from .index import Number, Shape, MatrixIndex, PolyIndex +from __future__ import annotations +from typing import TYPE_CHECKING from typing import Iterable - import numpy as np -import numpy.typing as npt + +if TYPE_CHECKING: + import numpy.typing as npt + + from .abc import Repr + from .index import Number, Shape, MatrixIndex, PolyIndex class SparseRepr(Repr): diff --git a/mdpoly/state.py b/mdpoly/state.py index 0095916..092de84 100644 --- a/mdpoly/state.py +++ b/mdpoly/state.py @@ -1,5 +1,6 @@ from __future__ import annotations from typing import TYPE_CHECKING + from dataclasses import dataclass, field from .index import PolyVarIndex diff --git a/mdpoly/util.py b/mdpoly/util.py index ae726ef..52cbb87 100644 --- a/mdpoly/util.py +++ b/mdpoly/util.py @@ -1,7 +1,12 @@ -from .constants import NUMERICS_EPS - from itertools import tee, filterfalse +from typing import TYPE_CHECKING + from typing import Iterable +from .constants import NUMERICS_EPS + +if TYPE_CHECKING: + pass + def partition(pred, iterable): """Partition entries into false entries and true entries. |