From bad4b9428956eca37bbb0eff6fb6c2abb0263e6a Mon Sep 17 00:00:00 2001 From: Michael Schneeberger Date: Wed, 29 Mar 2023 13:47:49 +0200 Subject: add half newton polytope operator --- polymatrix/expression/halfnewtonpolytopeexpr.py | 4 + .../expression/impl/halfnewtonpolytopeexprimpl.py | 10 ++ .../expression/init/inithalfnewtonpolytopeexpr.py | 14 +++ .../mixins/halfnewtonpolytopeexprmixin.py | 113 +++++++++++++++++++++ 4 files changed, 141 insertions(+) create mode 100644 polymatrix/expression/halfnewtonpolytopeexpr.py create mode 100644 polymatrix/expression/impl/halfnewtonpolytopeexprimpl.py create mode 100644 polymatrix/expression/init/inithalfnewtonpolytopeexpr.py create mode 100644 polymatrix/expression/mixins/halfnewtonpolytopeexprmixin.py diff --git a/polymatrix/expression/halfnewtonpolytopeexpr.py b/polymatrix/expression/halfnewtonpolytopeexpr.py new file mode 100644 index 0000000..8240d45 --- /dev/null +++ b/polymatrix/expression/halfnewtonpolytopeexpr.py @@ -0,0 +1,4 @@ +from polymatrix.expression.mixins.halfnewtonpolytopeexprmixin import HalfNewtonPolytopeExprMixin + +class HalfNewtonPolytopeExpr(HalfNewtonPolytopeExprMixin): + pass diff --git a/polymatrix/expression/impl/halfnewtonpolytopeexprimpl.py b/polymatrix/expression/impl/halfnewtonpolytopeexprimpl.py new file mode 100644 index 0000000..6fef1da --- /dev/null +++ b/polymatrix/expression/impl/halfnewtonpolytopeexprimpl.py @@ -0,0 +1,10 @@ +import dataclass_abc +from polymatrix.expression.halfnewtonpolytopeexpr import HalfNewtonPolytopeExpr + +from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin + +@dataclass_abc.dataclass_abc(frozen=True) +class HalfNewtonPolytopeExprImpl(HalfNewtonPolytopeExpr): + monomials: ExpressionBaseMixin + variables: ExpressionBaseMixin + filter: ExpressionBaseMixin | None diff --git a/polymatrix/expression/init/inithalfnewtonpolytopeexpr.py b/polymatrix/expression/init/inithalfnewtonpolytopeexpr.py new file mode 100644 index 0000000..e4ab6a2 --- /dev/null +++ b/polymatrix/expression/init/inithalfnewtonpolytopeexpr.py @@ -0,0 +1,14 @@ +from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin +from polymatrix.expression.impl.halfnewtonpolytopeexprimpl import HalfNewtonPolytopeExprImpl + + +def init_half_newton_polytope_expr( + monomials: ExpressionBaseMixin, + variables: ExpressionBaseMixin, + filter: ExpressionBaseMixin | None = None, +): + return HalfNewtonPolytopeExprImpl( + monomials=monomials, + variables=variables, + filter=filter +) diff --git a/polymatrix/expression/mixins/halfnewtonpolytopeexprmixin.py b/polymatrix/expression/mixins/halfnewtonpolytopeexprmixin.py new file mode 100644 index 0000000..ed1f6d9 --- /dev/null +++ b/polymatrix/expression/mixins/halfnewtonpolytopeexprmixin.py @@ -0,0 +1,113 @@ + +import abc +import dataclasses +import itertools +import math +from polymatrix.expression.utils.getmonomialindices import get_monomial_indices +from polymatrix.expression.utils.getvariableindices import get_variable_indices_from_variable + +from polymatrix.polymatrix.init.initpolymatrix import init_poly_matrix +from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin +from polymatrix.expressionstate.mixins.expressionstatemixin import ExpressionStateMixin +from polymatrix.polymatrix.mixins.polymatrixmixin import PolyMatrixMixin + + +class HalfNewtonPolytopeExprMixin(ExpressionBaseMixin): + @property + @abc.abstractclassmethod + def monomials(self) -> ExpressionBaseMixin: + ... + + @property + @abc.abstractmethod + def variables(self) -> ExpressionBaseMixin: + ... + + @property + @abc.abstractmethod + def filter(self) -> ExpressionBaseMixin | None: + ... + + # overwrites abstract method of `ExpressionBaseMixin` + def apply( + self, + state: ExpressionStateMixin, + ) -> tuple[ExpressionStateMixin, PolyMatrixMixin]: + + state, sos_monomials = get_monomial_indices(state, self.monomials) + state, variable_indices = get_variable_indices_from_variable(state, self.variables) + + if self.filter is not None: + state, filter_monomials = get_monomial_indices(state, self.filter) + else: + filter_monomials = None + + def gen_degrees(): + for monom in sos_monomials: + yield sum(d for _, d in monom) + + min_deg = min(math.floor(v / 2) for v in gen_degrees()) + max_deg = max(math.ceil(v / 2) for v in gen_degrees()) + + def gen_min_max_degree_per_variable(): + for var_index in variable_indices: + def gen_degrees_per_variable(): + for monom in sos_monomials: + for monom_var_index, count in monom: + if monom_var_index == var_index: + yield count + + min_deg = min(math.floor(v / 2) for v in gen_degrees_per_variable()) + max_deg = max(math.ceil(v / 2) for v in gen_degrees_per_variable()) + + yield var_index, (min_deg, max_deg) + + var_index_to_min_max = dict(gen_min_max_degree_per_variable()) + + def acc_combinations(acc, degree): + state, filtered_monomials = acc + + state, degree_monomials = get_monomial_indices( + state, + self.variables.combinations(degree), + ) + + if filter_monomials is not None: + def gen_filtered_monomials1(): + for monom in degree_monomials: + if monom not in filter_monomials: + yield monom + + degree_monomials_filt = tuple(gen_filtered_monomials1()) + + else: + degree_monomials_filt = degree_monomials + + def gen_filtered_monomials(): + for monom in degree_monomials_filt: + def is_candidate(): + for var_index, count in monom: + min_deg, max_deg = var_index_to_min_max[var_index] + + if not (min_deg <= count <= max_deg): + return False + + return True + + if is_candidate(): + yield monom + + return state, filtered_monomials + tuple(gen_filtered_monomials()) + + *_, (state, monomials) = itertools.accumulate( + range(min_deg, max_deg + 1), + acc_combinations, + initial=(state, tuple()), + ) + + poly_matrix = init_poly_matrix( + terms={(row, 0): {monom: 1} for row, monom in enumerate(monomials)}, + shape=(len(monomials), 1), + ) + + return state, poly_matrix -- cgit v1.2.1