Fix error in construction of multipliers for Positivstellensatz

2 files changed, 20 insertions, 8 deletions

diff --git a/sumofsquares/error.py b/sumofsquares/error.py
index 34ea6d8..1f30b2a 100644
--- a/sumofsquares/error.py
+++ b/sumofsquares/error.py
@@ -5,3 +5,7 @@ Errors and exceptions raised by the sum of squares package.
 class NotSupportedBySolver(Exception):
     """ The chosen solver cannot solve this problem. """
 
+
+class AlgebraicError(Exception):
+    """ TODO """
+
diff --git a/sumofsquares/problems.py b/sumofsquares/problems.py
index 14a9cd9..3ec3f54 100644
--- a/sumofsquares/problems.py
+++ b/sumofsquares/problems.py
@@ -19,6 +19,7 @@ from polymatrix.variable.abc import Variable
 
 from .abc import Problem, Constraint, Solver, Result
 from .constraints import NonNegative
+from .error import AlgebraicError
 from .solver.cvxopt import solve_sos_cone as cvxopt_solve_sos_cone
 from .utils import partition
 from .variable import OptVariable, from_name as opt_variable_from_name
@@ -84,9 +85,6 @@ class PutinarSOSProblem(Problem):
         x = poly.v_stack((1,) + tuple(self.polynomial_variables))
 
         for i, constr in enumerate(to_process):
-            # For each polynomial that defines the archimedean set
-            # we create a multiplier polynomial of the same degree
-
             # FIXME: this is not optimal, look into how to determine the
             # smallest degree needed to make the positivstellensatz work.
             new_constr = constr.expression
@@ -94,9 +92,19 @@ class PutinarSOSProblem(Problem):
                 # To know the degree of the polynomial we need to evaluate the
                 # expression, hence we use a monad here
                 def make_multiplier(state):
-                    state, p = domain_poly.degree().apply(state)
-
-                    d = p.at(0, 0).constant()
+                    state, constr_pm = constr.expression.degree().apply(state)
+                    state, domain_poly_pm = domain_poly.degree().apply(state)
+                    
+                    constr_deg = constr_pm.at(0, 0).constant()
+                    domain_poly_deg = domain_poly_pm.at(0, 0).constant()
+
+                    # degree of multiplier
+                    d = constr_deg - domain_poly_deg
+                    if d < 0:
+                        raise AlgebraicError("Cannot create Positivstellensatz multiplier because "
+                                             f"constraint polynomial has degree {constr_deg} and domain "
+                                             f"polynomial {domain_poly} has degree {domain_poly_deg}")
+                    
                     n = sum(math.comb(k + q - 1, k) for k in range(0, d+1))
 
                     # TODO: check that there is not a variable with this name already
@@ -119,8 +127,8 @@ class PutinarSOSProblem(Problem):
 
         # Convert to polymatrix object to extract variables
         new_state = self.state
-        for constr in new_constraints:
-            new_state, _ = constr.expression.apply(new_state)
+        for c in new_constraints:
+            new_state, _ = c.expression.apply(new_state)
 
         def is_optvariable(v):
             return isinstance(v, OptVariable)