Fix construction of PutinarPSatz multipliers

2 files changed, 8 insertions, 8 deletions

diff --git a/sumofsquares/canon.py b/sumofsquares/canon.py
index 02c9266..8c4b541 100644
--- a/sumofsquares/canon.py
+++ b/sumofsquares/canon.py
@@ -114,7 +114,7 @@ class PutinarPSatz(Canonicalization):
         new_constraints: list[Constraint] = list(constraints)
         for i, constr in enumerate(to_process):
             new_constr = constr.expression
-            for domain_poly in constr.domain.polynomials:
+            for j, domain_poly in enumerate(constr.domain.polynomials):
                 # Make a multiplier polynomial for a domain polynomial of a constraint.
                 # 
                 # Since Putinar's PSatz does not explicity tell us the order of the
@@ -140,20 +140,19 @@ class PutinarPSatz(Canonicalization):
                 if d.read(state).scalar().constant() % 2 != 0:
                     d += 1
 
-                x = poly.v_stack((1,) + tuple(
-                    init_variable_expr(v, state.get_shape(v))
-                    for v in prob.polynomial_variables))
+                x = poly.v_stack((1,) + tuple(init_variable_expr(v, state.get_shape(v))
+                                              for v in prob.polynomial_variables))
 
-                # FIXME: need to check that there is not a \gamma_{i} already!
+                # FIXME: need to check that there is not a \gamma_i,j already!
                 # TODO: deterministically generate unique names
 
                 # TODO: better documentation
                 # Generate monomials of up to degree d and take only monomials
                 # that are in half of the newton polytope
-                monomials = x.combinations(poly.arange(d + 1)).half_newton_polytope(x)
-                c = opt_variable_from_name(rf"\gamma_{i}", shape=monomials.shape)
+                monomials = x.combinations(poly.arange(d + 1)) # .half_newton_polytope(x)
+                gamma = opt_variable_from_name(rf"\gamma_{i},{j}", shape=monomials.shape)
 
-                multiplier = c.T @ monomials
+                multiplier = gamma.T @ monomials
 
                 # Subtract multiplier from expression and impose that it is also SOS
                 new_constr -= multiplier * domain_poly
diff --git a/sumofsquares/solver/cvxopt.py b/sumofsquares/solver/cvxopt.py
index 0692787..cf0dc4e 100644
--- a/sumofsquares/solver/cvxopt.py
+++ b/sumofsquares/solver/cvxopt.py
@@ -135,6 +135,7 @@ def solve_cone(prob: ConicProblem, verbose: bool = False,
             info = cvxopt.solvers.coneqp(P=P, q=q, G=G, h=h, A=A, b=b,
                                          dims=dims, *args, **kwargs)
         else:
+            assert P is None, "Solving QP with LP solver! Something has gone very wrong"
             info = cvxopt.solvers.conelp(c=q, G=G, h=h, A=A, b=b,
                                          *args, **kwargs)