add half newton polytope operator

4 files changed, 141 insertions, 0 deletions

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