diff options
| -rw-r--r-- | mdpoly/__init__.py | 5 | ||||
| -rw-r--r-- | mdpoly/abc.py | 48 | ||||
| -rw-r--r-- | mdpoly/algebra.py | 61 | ||||
| -rw-r--r-- | mdpoly/constants.py | 3 | ||||
| -rw-r--r-- | mdpoly/errors.py | 14 | ||||
| -rw-r--r-- | mdpoly/leaves.py | 2 | ||||
| -rw-r--r-- | mdpoly/representations.py | 56 | ||||
| -rw-r--r-- | mdpoly/state.py | 6 | ||||
| -rw-r--r-- | mdpoly/types.py | 68 |
9 files changed, 135 insertions, 128 deletions
diff --git a/mdpoly/__init__.py b/mdpoly/__init__.py index 3f56b27..42a800d 100644 --- a/mdpoly/__init__.py +++ b/mdpoly/__init__.py @@ -1,7 +1,5 @@ -""" -A Library to represent multivariate polynomials +""" A Library to represent multivariate polynomials """ -""" # internal classes imported with underscore because # they should not be exposed to the end users @@ -31,6 +29,7 @@ State = _State # ┃╻┃┣┳┛┣━┫┣━┛┣━┛┣╸ ┣┳┛┗━┓ # ┗┻┛╹┗╸╹ ╹╹ ╹ ┗━╸╹┗╸┗━┛ +# FIXME: move out of this file class _FromHelpers: @classmethod def from_names(cls, comma_separated_names: str, strip: bool =True) -> Sequence[Self]: diff --git a/mdpoly/abc.py b/mdpoly/abc.py index 77285aa..a38a69d 100644 --- a/mdpoly/abc.py +++ b/mdpoly/abc.py @@ -1,9 +1,5 @@ -""" Abstract Base Classes of MDPoly +""" Abstract Base Classes of MDPoly """ -Defines: - - Tree data structure to represent expressions - - Representations of polynomials -""" from .types import Number, Shape, MatrixIndex, PolyIndex from .constants import NUMERICS_EPS @@ -13,18 +9,14 @@ from abc import abstractmethod @runtime_checkable class Leaf(Protocol): - """ - Leaf of the binary tree. - """ + """ Leaf of the binary tree. """ name: str shape: Shape @runtime_checkable class Expr(Protocol): - """ - Binary tree to represent a mathematical expression. - """ + """ Binary tree to represent a mathematical expression. """ left: Self | Leaf right: Self | Leaf shape: Shape @@ -34,9 +26,12 @@ class Expr(Protocol): return f"{name}(left={self.left}, right={self.right})" def nodes(self) -> Sequence[Self]: + """ Iterate over the two nodes """ return self.left, self.right def leaves(self) -> Sequence[Leaf]: + """ Returns the leaves of the tree. This is done recursively and in + :math:`\mathcal{O}(n)`.""" if isinstance(self.left, Leaf): yield self.left else: @@ -52,44 +47,55 @@ class Expr(Protocol): class Repr(Protocol): - """ - Representation of a polynomial expression. + r""" Representation of a multivariate matrix polynomial expression. + + Concretely, a representation must be able to store a mathematical object + like :math:`Q \in \mathbb{R}^2[x, y]`, for example + + .. math:: + Q = \begin{bmatrix} + x^2 + 1 & 5y \\ + 0 & y^2 + \end{bmatrix} + + See also :py:mod:`mdpoly.representations`. """ @abstractmethod def at(self, entry: MatrixIndex, term: PolyIndex) -> Number: - """ Access polynomial entry """ + """ Access polynomial coefficient. """ @abstractmethod def set(self, entry: MatrixIndex, term: PolyIndex, value: Number) -> None: - """ Set value of polynomial entry """ + """ Set value of polynomial coefficient. """ def is_zero(self, entry: MatrixIndex, term: PolyIndex, eps: Number = NUMERICS_EPS) -> bool: - """ Check if a polynomial entry is zero """ + """ Check if a polynomial coefficient is zero. """ return self.at(entry, term) < eps @abstractmethod def set_zero(self, entry: MatrixIndex, term: PolyIndex) -> None: - """ Set an entry to zero """ + """ Set a coefficient to zero (delete it). """ @abstractmethod def entries(self) -> Sequence[MatrixIndex]: - """ Return indices to non-zero entries of the matrix """ + """ Return indices to non-zero entries of the matrix. """ @abstractmethod def terms(self, entry: MatrixIndex) -> Sequence[PolyIndex]: """ Return indices to non-zero terms in the polynomial at the given - matrix entry """ + matrix entry. """ def zero_out_small(self, eps: Number = NUMERICS_EPS) -> None: - """ Zero out (set to zero) values that are smaller than EPS """ + """ Zero out (set to zero) values that are smaller than ``eps`` + See also :py:data:`mdpoly.constants.NUMERICS_EPS`. """ for entry in self.entries(): for term in self.terms(entry): if self.at(entry, term) < eps: self.set_zero(entry, term) def __iter__(self) -> Sequence[tuple[MatrixIndex, PolyIndex, Number]]: - """ Iterate over non-zero entries of the representations """ + """ Iterate over non-zero entries of the representations. """ for entry in self.entries(): for term in self.terms(entry): yield entry, term, self.at(entry, term) diff --git a/mdpoly/algebra.py b/mdpoly/algebra.py index 4d8e47e..12c9d10 100644 --- a/mdpoly/algebra.py +++ b/mdpoly/algebra.py @@ -1,7 +1,5 @@ -""" -Algebraic Structures for Expressions +""" Algebraic Structures for Expressions """ -""" from .abc import Expr, Repr from .leaves import Nothing, Const, Param from .errors import AlgebraicError @@ -19,7 +17,7 @@ import operator def binary_operator(cls: Expr) -> Expr: - """ Class decorator to modify constructor of binary operators """ + """ Class decorator to modify constructor of binary operators. """ __init_cls__ = cls.__init__ @wraps(cls.__init__) def __init_binary__(self, left, right, *args, **kwargs): @@ -31,7 +29,7 @@ def binary_operator(cls: Expr) -> Expr: def unary_operator(cls: Expr) -> Expr: - """ Class decorator to modify constructor of unary operators """ + """ Class decorator to modify constructor of unary operators. """ __init_cls__ = cls.__init__ @wraps(cls.__init__) def __init_unary__(self, left, *args, **kwargs): @@ -46,11 +44,15 @@ T = TypeVar("T") @runtime_checkable class AlgebraicStructure(Protocol): + """ Provides methods to enforce algebraic closure of operations. + + Internally, this is done by adding a class variable ``_algebra`` of type + :py:class:`mdpoly.algebra.AlgebraicStructure.Algebra`, and checking that + they match in the methods for operator overloading. """ - Provides methods to enforce algebraic closure of operations - """ + class Algebra(Enum): - """ Types of algebras """ + """ Types of algebras. """ poly_ring = 0 matrix_ring = 1 @@ -99,7 +101,8 @@ class AlgebraicStructure(Protocol): class ReducibleExpr(HasRepr, Protocol): - """ + """ 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 @@ -107,10 +110,12 @@ class ReducibleExpr(HasRepr, Protocol): """ def to_repr(self, repr_type: type, state: State) -> Repr: + """ See :py:meth:`HasRepr.to_repr` """ return self.reduce().to_repr(repr_type, state) @abstractmethod - def reduce(self) -> Expr: ... + def reduce(self) -> Expr: + """ Reduce the expression to its basic elements """ # ┏━┓┏━┓╻ ╻ ╻┏┓╻┏━┓┏┳┓╻┏━┓╻ ┏━┓╻ ┏━╸┏━╸┏┓ ┏━┓┏━┓ @@ -118,9 +123,8 @@ class ReducibleExpr(HasRepr, Protocol): # ╹ ┗━┛┗━╸ ╹ ╹ ╹┗━┛╹ ╹╹╹ ╹┗━╸ ╹ ╹┗━╸┗━┛┗━╸┗━┛╹┗╸╹ ╹ class PolyRingAlgebra(AlgebraicStructure, Protocol): - """ - Endows with the algebraic structure of a polynomial ring. - """ + """ Endows with the algebraic structure of a polynomial ring. """ + _algebra = AlgebraicStructure.Algebra.poly_ring _parameter = Param _constant = Const @@ -169,9 +173,10 @@ class PolyRingAlgebra(AlgebraicStructure, Protocol): @binary_operator class Add(Expr, HasRepr, PolyRingAlgebra): - """ Sum operator between scalar polynomials """ + """ Addition operator between scalar polynomials. """ def to_repr(self, repr_type: type, state: State) -> Repr: + """ See :py:meth:`mdpoly.representations.HasRepr.to_repr`. """ # Make a new empty representation r = repr_type() entry = MatrixIndex.scalar @@ -193,9 +198,10 @@ class Add(Expr, HasRepr, PolyRingAlgebra): @binary_operator class Subtract(Expr, ReducibleExpr, PolyRingAlgebra): - """ Subtraction operation between scalar polynomials """ + """ Subtraction operator between scalar polynomials. """ def reduce(self) -> Expr: + """ See :py:meth:`mdpoly.algebra.ReducibleExpr.reduce`. """ return Add(self.left, Multiply(self._constant(value=-1), self.right)) def __repr__(self) -> str: @@ -204,9 +210,10 @@ class Subtract(Expr, ReducibleExpr, PolyRingAlgebra): @binary_operator class Multiply(Expr, HasRepr, PolyRingAlgebra): - """ Multiplication operator between scalar polynomials """ + """ Multiplication operator between scalar polynomials. """ def to_repr(self, repr_type: type, state: State) -> Repr: + """ See :py:meth:`mdpoly.representations.HasRepr.to_repr`. """ r = repr_type() lrepr, state = self.left.to_repr(repr_type, state) @@ -229,9 +236,10 @@ class Multiply(Expr, HasRepr, PolyRingAlgebra): @binary_operator class Divide(Expr, ReducibleExpr, PolyRingAlgebra): - """ Division operator between scalar polynomials """ + """ Division operator between scalar polynomials. """ def reduce(self) -> Expr: + """ See :py:meth:`mdpoly.algebra.ReducibleExpr.reduce`. """ inverse = self._constant(value=1/self.right.value) return Multiply(inverse, self.left) @@ -241,9 +249,10 @@ class Divide(Expr, ReducibleExpr, PolyRingAlgebra): @binary_operator class Exponentiate(Expr, ReducibleExpr, PolyRingAlgebra): - """ Exponentiation operator of scalar polynomials """ + """ Exponentiation operator of scalar polynomials. """ def reduce(self) -> Expr: + """ See :py:meth:`mdpoly.algebra.ReducibleExpr.reduce`. """ var = self.left ntimes = self.right.value - 1 return reduce(operator.mul, (var for _ in range(ntimes)), var) @@ -254,7 +263,7 @@ class Exponentiate(Expr, ReducibleExpr, PolyRingAlgebra): @unary_operator class Diff(Expr, PolyRingAlgebra): - """ Total differentiation """ + """ Total differentiation. """ # ┏┳┓┏━┓╺┳╸┏━┓╻╻ ╻ ┏━┓╻ ┏━╸┏━╸┏┓ ┏━┓┏━┓ @@ -262,9 +271,8 @@ class Diff(Expr, PolyRingAlgebra): # ╹ ╹╹ ╹ ╹ ╹┗╸╹╹ ╹ ╹ ╹┗━╸┗━┛┗━╸┗━┛╹┗╸╹ ╹ class MatrixAlgebra(AlgebraicStructure, Protocol): - """ - Endows with the algebraic structure of a matrix ring. - """ + """ Endows with the algebraic structure of a matrix ring. """ + _algebra = AlgebraicStructure.Algebra.matrix_ring # FIXME: consider MatParam or something like that? _parameter = Param @@ -299,26 +307,29 @@ class MatrixAlgebra(AlgebraicStructure, Protocol): raise NotImplementedError def transpose(self): + """ Matrix transposition. """ raise NotImplementedError @property def T(self): + """ Shorthand for :py:meth:`mdpoly.algebra.MatrixAlgebra.transpose`. """ return self.transpose() def as_scalar(self) -> PolyRingAlgebra: + """ Convert to a scalar expression. """ raise NotImplementedError @binary_operator class MatAdd(Expr, MatrixAlgebra): - """ Matrix Addition """ + """ Matrix Addition. """ @binary_operator class ScalarMul(Expr, MatrixAlgebra): - """ Matrix-Scalar Multiplication """ + """ Matrix-Scalar Multiplication. """ @binary_operator class MatMul(Expr, MatrixAlgebra): - """ Matrix Multiplication """ + """ Matrix Multiplication. """ diff --git a/mdpoly/constants.py b/mdpoly/constants.py index 2d62186..4ea75f0 100644 --- a/mdpoly/constants.py +++ b/mdpoly/constants.py @@ -1 +1,4 @@ NUMERICS_EPS = 1e-12 +r""" Value of :math:`\varepsilon` used for floating point computations. + Floating point numbers smaller than :math:`\varepsilon` are considered + as being equal to zero. """ diff --git a/mdpoly/errors.py b/mdpoly/errors.py index 578c5c8..152d6df 100644 --- a/mdpoly/errors.py +++ b/mdpoly/errors.py @@ -1,16 +1,12 @@ class AlgebraicError(Exception): - """ - This is raised whenever an invalid mathematical operation is performed. - See more in mdpoly/ops.py + """ This is raised whenever an invalid mathematical operation is performed. + See also :py:mod:`mdpoly.algebra` and :py:class:`mdpoly.algebra.AlgebraicStructure` """ class InvalidShape(Exception): - """ - This is raised whenever there is an invalid shape - """ + """ This is raised whenever an operation cannot be perfomed because the + shapes of the inputs do not match. """ class MissingParameters(Exception): - """ - This is raised whenever - """ + """ This is raised whenever a parameter was not given. """ diff --git a/mdpoly/leaves.py b/mdpoly/leaves.py index a9d2005..14562a9 100644 --- a/mdpoly/leaves.py +++ b/mdpoly/leaves.py @@ -50,7 +50,7 @@ class Var(Leaf, HasRepr): def to_repr(self, repr_type: type[T], state: State) -> T: r = repr_type() - idx = PolyVarIndex.of_var(self, state), + idx = PolyVarIndex.from_var(self, state), r.set(MatrixIndex.scalar, PolyIndex(idx), 1) return r, state diff --git a/mdpoly/representations.py b/mdpoly/representations.py index 5e58cb4..802167b 100644 --- a/mdpoly/representations.py +++ b/mdpoly/representations.py @@ -10,21 +10,26 @@ import numpy.typing as npt class HasRepr(Protocol): + """ Has a representation. """ @abstractmethod - def to_repr(self, repr_type: type, state: State) -> tuple[Repr, State]: ... + def to_repr(self, repr_type: type, state: State) -> tuple[Repr, State]: + """ Convert to a concrete representation + To convert an abstract object from the expression tree, we ... TODO: finish. + + See also :py:class:`mdpoly.state.State`. + """ class SparseRepr(Repr): - """ - Sparse representation of polynomial - """ + """ Sparse representation of polynomial. """ + data: dict[MatrixIndex, dict[PolyIndex, Number]] def __init__(self): self.data = {} def at(self, entry: MatrixIndex, term: PolyIndex) -> Number: - """ Access polynomial entry """ + """ Access polynomial entry. """ if entry in self.data.keys(): if term in self.data[entry].keys(): return self.data[entry][term] @@ -32,42 +37,37 @@ class SparseRepr(Repr): return 0 def set(self, entry: MatrixIndex, term: PolyIndex, value: Number) -> None: - """ Set value of polynomial entry """ + """ Set value of polynomial entry. """ if entry not in self.data.keys(): self.data[entry] = {} self.data[entry][term] = value def set_zero(self, entry: MatrixIndex, term: PolyIndex): - """ Set an entry to zero """ + """ Set an entry to zero. """ del self.data[entry][term] if not self.data[entry]: del self.data[entry] def entries(self) -> Sequence[MatrixIndex]: - """ Return indices to non-zero entries of the matrix """ + """ Return indices to non-zero entries of the matrix. """ yield from self.data.keys() def terms(self, entry: MatrixIndex) -> Sequence[PolyIndex]: - """ Return indices to non-zero terms in the polynomial at the given - matrix entry """ + """ Return indices to terms with a non-zero coefficient in the + polynomial at the given matrix entry. """ if entry in self.data.keys(): yield from self.data[entry].keys() return tuple() -class HasSparseRepr(Protocol): - @abstractmethod - def to_sparse(self) -> SparseRepr: ... - - class SparseMatrixRepr(Repr): - """ - Sparse matrix representation of a polynomial. + """ Sparse matrix representation of a polynomial. - Collect the monomials in a vector x (the "basis") and store a matrix Q and - such that the quadratic form x^T Q x gives the polynomial. + Collect the monomials in a vector :math:`x` (the "basis") and store a + matrix :math:`Q` and such that the quadratic form :math:`x^T Q x` gives the + polynomial. """ def at(self, entry: MatrixIndex, term: PolyIndex) -> Number: @@ -94,15 +94,8 @@ class SparseMatrixRepr(Repr): raise NotImplementedError -class HasSparseMatrixRepr(Protocol): - @abstractmethod - def to_sparse_matrix(self) -> SparseMatrixRepr: ... - - class DenseRepr(Repr): - """ - Dense representation of a polynomial - """ + """ Dense representation of a polynomial that uses Numpy arrays. """ data: npt.NDArray[Number] def __init__(self, shape: Shape, dtype:type =float): @@ -126,12 +119,5 @@ class DenseRepr(Repr): raise NotImplementedError -class HasDenseRepr(Protocol): - @abstractmethod - def to_dense(self) -> DenseRepr: ... - - class DenseMatrixRepr(Repr): - """ - Dense matrix representation of a polynomials - """ + """ Dense matrix representation of a polynomials """ diff --git a/mdpoly/state.py b/mdpoly/state.py index 6f47964..001d9b6 100644 --- a/mdpoly/state.py +++ b/mdpoly/state.py @@ -17,12 +17,11 @@ class State: _last_index: Index = -1 def _make_index(self) -> Index: - """ Make a new index """ self._last_index += 1 return self._last_index def index(self, var: Var) -> Index: - """ Get the index for a variable """ + """ Get the index for a variable. """ if var not in self.variables.keys(): new_index = self._make_index() self.variables[var] = new_index @@ -31,6 +30,9 @@ class State: return self.variables[var] def from_index(self, index: Index) -> Var | Const: + """ Get a variable object from the index. + This is a reverse lookup operation and it is **not** + :math:`\mathcal{O}(1)`!""" if index == -1: from .leaves import Const return Const(value=1, name="1") diff --git a/mdpoly/types.py b/mdpoly/types.py index a8ffd90..1132260 100644 --- a/mdpoly/types.py +++ b/mdpoly/types.py @@ -14,13 +14,12 @@ Number = int | float def filtertype(t: type, seq: Iterable): + """ Filter based on type. """ yield from filter(lambda x: isinstance(x, t), seq) class Shape(NamedTuple): - """ - Describes the shape of a mathematical object. - """ + """ Describes the shape of a mathematical object. """ rows: int cols: int @@ -48,9 +47,7 @@ class Shape(NamedTuple): class MatrixIndex(NamedTuple): - """ - Tuple to index an element of a matrix or vector. - """ + """ Tuple to index an element of a matrix or vector. """ row: int col: int @@ -67,14 +64,21 @@ class MatrixIndex(NamedTuple): class PolyVarIndex(NamedTuple): - """ - Tuple to index a power of a variable. + """ Tuple to represent a variable and its exponent. They are totally + ordered with respect to the variable index stored in the ``state`` object + (see :py:attr:`mdpoly.state.State.variables`). - Concretely this represents things like x^2, y^4 but not xy + Concretely this represents things like :math:`x^2`, :math:`y^4` but not + :math:`xy`. """ var_idx: int # Index in State.variables, not variable! power: Number + @classmethod + def from_var(cls, variable: Var, state: State, power: Number =1) -> Self: + """ Make an index from a variable object. """ + return cls(var_idx=state.index(variable), power=power) + def __eq__(self, other: Self): if self.var_idx != other.var_idx: return False @@ -90,41 +94,38 @@ class PolyVarIndex(NamedTuple): @classmethod @property def constant(cls): - """ - Special index for constants. + """ Special index for constants. - Constants do not have an associated variable, - and the field power is ignored. + Constants do not have an associated variable, and the field power is + ignored. """ return cls(var_idx=-1, power=0) @staticmethod def is_constant(index: Self) -> bool: - """ Check if is index of a constant term """ + """ Check if is index of a constant term. """ return index.var_idx == -1 - @classmethod - def of_var(cls, variable: Var, state: State, power: Number =1) -> Self: - return cls(var_idx=state.index(variable), power=power) - class PolyIndex(tuple[PolyVarIndex]): """ Tuple to index a coefficient of a monomial in a polynomial. - Example: + For example, suppose there are two variables :math:`x, y` with indices 0, 1 + respectively (in the State object, see :py:attr:`mdpoly.state.State.variables`). + Then the monomal :math:`xy^2` has an index - Suppose there are two variables x, y with indices 0, 1 respectively (in - the State object). Then the monomal xy^2 has PolyIndex + .. code:: py - (PolyVarIndex(var_idx=0, power=1), PolyVarIndex(var_idx=1, power=2) + PolyIndex(PolyVarIndex(var_idx=0, power=1), PolyVarIndex(var_idx=1, power=2) - Then given the polynomial + Then given the polynomial - p(x, y) = 3x^2 + 5xy^2 + .. math:: + p(x, y) = 3x^2 + 5xy^2 - The with the PolyIndex above we can retrieve the coefficient 5 in front - of xy^2 from the Repr object. + The with the ``PolyIndex`` above we can retrieve the coefficient 5 in front + of :math:`xy^2` from the ``Repr`` object (see also :py:meth:`mdpoly.abc.Repr.at`). """ def __repr__(self) -> str: @@ -134,17 +135,20 @@ class PolyIndex(tuple[PolyVarIndex]): @classmethod def from_dict(cls, d) -> Self: + """ Construct an index froma dictionary, where the keys are the + variable index from the state object and the values are the exponent + (power). """ return cls(PolyVarIndex(k, v) for k, v in d.items()) @classmethod @property def constant(cls) -> Self: - """ Index of the constant term """ + """ Index of the constant term. """ return cls((PolyVarIndex.constant,)) @staticmethod def is_constant(index: Self) -> bool: - """ Check if is index of a constant term """ + """ Check if is index of a constant term. """ if len(index) != 1: return False @@ -152,16 +156,16 @@ class PolyIndex(tuple[PolyVarIndex]): @classmethod def sort(cls, index: Self) -> Self: - """ Sort a tuple of indices """ + """ Sort a tuple of indices. """ return cls(sorted(index)) @classmethod def product(cls, left: Self, right: Self) -> Self: """ - Compute index of product + Compute index of product. - For example if left is the index of xy and right is the index of y^2 - this functions returns the index of xy^3 + For example if left is the index of :math:`xy` and right is the index + of :math:`y^2` this functions returns the index of :math:`xy^3`. """ if cls.is_constant(left): return right |