Major structure change to allow more general transformation of problems

7 files changed, 592 insertions, 402 deletions

diff --git a/sumofsquares/__init__.py b/sumofsquares/__init__.py
index 6ca2b39..6c02f24 100644
--- a/sumofsquares/__init__.py
+++ b/sumofsquares/__init__.py
@@ -9,6 +9,7 @@ Module Overview
 
 Class diagram: Arrows denote inheritance, diamonds indicate composition.
 ::
+                                                                                                                                       
    ┏━━sumofsquares package━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┓                                   
    ┃                                                                                               ┃                                   
    ┃ ┌─.abc──────────────┐   ┌─.problems────────────────────────────────────────────────────┐      ┃                                   
@@ -18,12 +19,12 @@ Class diagram: Arrows denote inheritance, diamonds indicate composition.
    ┃ │         │         │   │  ┌───────────────────┐  ┌────────────┐  ┌──────────────────┐ │      ┃                                   
    ┃ │         ◇         │   │  │   ConicProblem    │  │ SOSProblem │  │InternalSOSProblem│ │      ┃                                   
    ┃ │ ┌──────────────┐  │   │  └───────────────────┘  └────────────┘  └──────────────────┘ │      ┃                                   
-   ┃ │ │ Solver(Enum) │  │   │                                ▲                             │      ┃                                   
-   ┃ │ └──────────────┘  │   │                           ┌────┴────────────────┐            │      ┃                                   
-   ┃ │                   │   │                           │                     │            │      ┃                                   
-   ┃ │ ┌──────────────┐  │   │ ┌─────────────┐  ┌─────────────────┐ ┌────────────────────┐  │      ┃                                   
-   ┃ │ │    Result    │◀─┼───┼─│ ConicResult │  │PutinarSOSProblem│ │SchmuedgenSOSProblem│  │      ┃                                   
-   ┃ │ └──────────────┘  │   │ └─────────────┘  └─────────────────┘ └────────────────────┘  │      ┃                                   
+   ┃ │ │ Solver(Enum) │  │   │                                                              │      ┃                                   
+   ┃ │ └──────────────┘  │   │                                                              │      ┃                                   
+   ┃ │                   │   │                                                              │      ┃                                   
+   ┃ │ ┌──────────────┐  │   │ ┌─────────────┐                                              │      ┃                                   
+   ┃ │ │    Result    │◀─┼───┼─│ ConicResult │                                              │      ┃                                   
+   ┃ │ └──────────────┘  │   │ └─────────────┘                                              │      ┃                                   
    ┃ │         │         │   └──────────────────────────────────────────────────────────────┘      ┃                                   
    ┃ │         │         │                                                                         ┃                                   
    ┃ │         ◇         │   ┌─.solvers.cvxopt──┐  ┌─.solvers.scs─────┐  ┌─.solvers.mosek───┐      ┃                                   
@@ -50,18 +51,18 @@ Class diagram: Arrows denote inheritance, diamonds indicate composition.
    ┃ │                   │   │                                                                   │ ┃                                   
    ┃ └───────────────────┘   └───────────────────────────────────────────────────────────────────┘ ┃                                   
    ┃                                                                                               ┃                                   
-   ┃                                                                ┌─.variable─────────────────┐  ┃   ┏━━━polymatrix package━━━━━━━━━┓
-   ┃                                                                │                           │  ┃   ┃                              ┃
-   ┃                                                                │    ┌───────────┐          │  ┃   ┃  ┌──────────┐   ┌──────────┐ ┃
-   ┃                                                                │    │OptVariable│──────────┼──╋───╋─▶│ Variable │   │Expression│ ┃
-   ┃                                                                │    └───────────┘          │  ┃   ┃  └──────────┘   └──────────┘ ┃
-   ┃                                                                │          ▲                │  ┃   ┃                       │      ┃
-   ┃                                                                │          │                │  ┃   ┃             ┌─────────┘      ┃
-   ┃                                                                │          │                │  ┃   ┃             ◇                ┃
-   ┃                                                                │ ┌────────────────┐        │  ┃   ┃  ┌─────────────────────┐     ┃
-   ┃                                                                │ │OptVariableMixin│────────┼──╋───╋─▶│ ExpressionBaseMixin │     ┃
-   ┃                                                                │ └────────────────┘        │  ┃   ┃  └─────────────────────┘     ┃
-   ┃                                                                └───────────────────────────┘  ┃   ┃                              ┃
+   ┃  ┌─.canon───────────────────────────────────────────────────┐  ┌─.variable─────────────────┐  ┃   ┏━━━polymatrix package━━━━━━━━━┓
+   ┃  │                                                          │  │                           │  ┃   ┃                              ┃
+   ┃  │ ┌────────────────────┐                                   │  │    ┌───────────┐          │  ┃   ┃  ┌──────────┐   ┌──────────┐ ┃
+   ┃  │ │  Canonicalization  │                                   │  │    │OptVariable│──────────┼──╋───╋─▶│ Variable │   │Expression│ ┃
+   ┃  │ └────────────────────┘                                   │  │    └───────────┘          │  ┃   ┃  └──────────┘   └──────────┘ ┃
+   ┃  │            ▲                                             │  │          ▲                │  ┃   ┃                       │      ┃
+   ┃  │            ├────────────────────┐                        │  │          │                │  ┃   ┃             ┌─────────┘      ┃
+   ┃  │            │                    │                        │  │          │                │  ┃   ┃             ◇                ┃
+   ┃  │    ┌───────────────┐      ┌───────────┐                  │  │ ┌────────────────┐        │  ┃   ┃  ┌─────────────────────┐     ┃
+   ┃  │    │ PutinarPSatz  │      │  LogDet   │                  │  │ │OptVariableMixin│────────┼──╋───╋─▶│ ExpressionBaseMixin │     ┃
+   ┃  │    └───────────────┘      └───────────┘                  │  │ └────────────────┘        │  ┃   ┃  └─────────────────────┘     ┃
+   ┃  └──────────────────────────────────────────────────────────┘  └───────────────────────────┘  ┃   ┃                              ┃
    ┗━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┛   ┗━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┛
 
 """
@@ -72,8 +73,9 @@ from polymatrix.expression.expression import Expression
 from polymatrix.expressionstate import ExpressionState
 
 from .abc import Problem, Result, Set, Constraint, Solver
+from .canon import Canonicalization, PutinarPSatz
 from .constraints import BasicSemialgebraicSet, NonNegative
-from .problems import PutinarSOSProblem
+from .problems import SOSProblem
 from .utils import partition
 from .variable import (
         OptVariable,
@@ -97,41 +99,26 @@ def make_problem(
         cost: Expression,
         constraints: Iterable[Constraint] = (),
         solver: Solver = Solver.CVXOPT,
-        state: ExpressionState | None = None
-    ) -> Problem:
+        psatz: Canonicalization = PutinarPSatz
+    ) -> SOSProblem:
     """
     Create a sum-of-squares optimization problem.
     """
-    if not state:
-        state = poly.make_state()
-
-        # Convert to polymatrix object to extract variables
-        state, _ = cost.apply(state)
-        for constr in constraints:
-            state, _ = constr.expression.apply(state)
-
-    # Collect decision variables from state object
-    def is_optvariable(v):
-        return isinstance(v, OptVariable)
-
-    poly_variables, opt_variables = partition(is_optvariable, state.indices.keys())
-    return PutinarSOSProblem(
-            cost=cost,
-            constraints=constraints,
-            variables=tuple(opt_variables),
-            polynomial_variables=tuple(poly_variables),
-            solver=solver,
-            state=state)
+    return psatz(SOSProblem(cost, constraints, solver))
 
 
 def solve_problem(
         cost: Expression,
         constraints: Iterable[Constraint] = (),
         solver: Solver = Solver.CVXOPT,
-        verbose: bool = False
+        verbose: bool = False,
+        state: ExpressionState | None = None
     ) -> tuple[Problem, Result]:
     """
     Solve a sum-of-squares optimization problem.
     """
+    if state is None:
+        state = poly.make_state()
+
     prob = make_problem(cost, constraints, solver)
     return prob, prob.solve(verbose)
diff --git a/sumofsquares/abc.py b/sumofsquares/abc.py
index f0e95bb..20f4bd6 100644
--- a/sumofsquares/abc.py
+++ b/sumofsquares/abc.py
@@ -1,101 +1,85 @@
 """ Abstract base classes for Sum of Squares package """
 from __future__ import annotations
-from typing import Sequence, Any
 
 from abc import ABC, abstractmethod
 from enum import Enum, auto
-
-from polymatrix.expression.expression import Expression
-from polymatrix.expressionstate import ExpressionState
-from polymatrix.variable import Variable
+from typing import Any, Generic, TypeVar
 
 from sumofsquares.variable import OptVariable
 
 
+# ┏━┓┏━┓╻  ╻ ╻┏━╸┏━┓
+# ┗━┓┃ ┃┃  ┃┏┛┣╸ ┣┳┛
+# ┗━┛┗━┛┗━╸┗┛ ┗━╸╹┗╸
+
 class Solver(Enum):
     """ Enumeration for the supported solvers. """
     CVXOPT = auto()
-    # MOSEK = auto()
+    MOSEK = auto()
+    SCS = auto()
 
 
 class SolverInfo(ABC):
     """ Type that information returned by a specific solver. """
 
 
-class Result(ABC):
-    """ Result of an optimization problem. """
-    @abstractmethod
-    def value_of(self, var: OptVariable) -> float:
-        """ Retrieve value of variable. """
+# ┏━╸┏━┓┏┓╻┏━┓╺┳╸┏━┓┏━┓╻┏┓╻╺┳╸┏━┓
+# ┃  ┃ ┃┃┗┫┗━┓ ┃ ┣┳┛┣━┫┃┃┗┫ ┃ ┗━┓
+# ┗━╸┗━┛╹ ╹┗━┛ ╹ ╹┗╸╹ ╹╹╹ ╹ ╹ ┗━┛
 
 
 class Set(ABC):
     """ A set. In this context always a subset of :math:`\mathbf{R}^q` """
 
 
-class Constraint(ABC):
-    r"""
-    Optimization constraint. 
-    """
+E = TypeVar("E")
+class Constraint(ABC, Generic[E]):
+    """ Optimization constraint. """
     @property
     @abstractmethod
-    def expression(self) -> Expression:
+    def expression(self) -> E:
         """ Expression under the constraint. """
 
     def __str__(self):
         return f"{self.__class__.__qualname__}({str(self.expression)})"
 
 
-class Problem(ABC):
-    """ Optimization Problem. """
-    
-    # --- dunder methods ---
+# ┏━┓┏━┓┏━┓┏┓ ╻  ┏━╸┏┳┓┏━┓
+# ┣━┛┣┳┛┃ ┃┣┻┓┃  ┣╸ ┃┃┃┗━┓
+# ╹  ╹┗╸┗━┛┗━┛┗━╸┗━╸╹ ╹┗━┛
 
-    def __str__(self):
-        vs = ", ".join(map(str, self.variables))
-        n = len(self.constraints)
-        s = f"SumOfSquares Problem ({self.__class__.__qualname__}):\n\n  " \
-                f"find      {vs}\n  minimize  {self.cost}\n  " \
-                f"subject to ({n} constraint{'s' if n > 1 else ''})\n"
 
-        for c in self.constraints:
-            s += f"\t* {str(c)}\n"
+class Result(ABC):
+    """ Result of an optimization problem. """
+    @abstractmethod
+    def value_of(self, var: OptVariable) -> float:
+        """ Retrieve value of variable. """
+
 
-        return s
+class Problem(ABC):
+    """ Optimization Problem. """
 
     # --- Properties ---
 
     @property
     @abstractmethod
-    def cost(self) -> Expression:
-        """ Expression of cost function. """
+    def cost(self) -> Any:
+        """ Cost function of the problem. """
+        # The type here is unspecified on purpose, see
+        # programs module for why is it Any
 
     @property
     @abstractmethod
-    def constraints(self) -> Sequence[Constraint]:
+    def constraints(self) -> Any:
         """ Optimization constraints. """
-
-    @property
-    @abstractmethod
-    def variables(self) -> Sequence[OptVariable]:
-        """ Optimization variables of the problem. """
-
-    @property
-    @abstractmethod
-    def polynomial_variables(self) -> Sequence[Variable]:
-        """ Polynomial variables, or varaibles that are not optimization varaibles. """
+        # The type here is unspecified on purpose, see
+        # programs module for why is it Any
 
     @property
     @abstractmethod
     def solver(self) -> Solver:
         """ Which solver was / will be used to solve the problem. """
 
-
-    @property
-    @abstractmethod
-    def state(self) -> ExpressionState:
-        """ Expression state used to evaluate constraints and costs. """
-
     # --- Methods ---
 
     @abstractmethod
diff --git a/sumofsquares/canon.py b/sumofsquares/canon.py
index 4918758..e5960e7 100644
--- a/sumofsquares/canon.py
+++ b/sumofsquares/canon.py
@@ -1,11 +1,131 @@
 """ Canonicalization of Problems """
-from .abc import Problem
+import math
+import polymatrix as poly
 
+from abc import abstractmethod
+from dataclassabc import dataclassabc
+from dataclasses import replace
+from typing import Sequence
+from typing_extensions import override
 
-def log_canon(prob: Problem) -> Problem:
-    raise NotImplementedError
+from polymatrix.expression.from_ import Expression, from_statemonad
+from polymatrix.expressionstate import ExpressionState
+from polymatrix.statemonad import init_state_monad
+from polymatrix.polymatrix.mixins import PolyMatrixMixin
+from polymatrix.variable import Variable
 
+from polymatrix.expression.init import (
+        init_concatenate_expr,
+        init_diag_expr,
+        init_lower_triangular_expr,
+        init_transpose_expr)
+
+from .abc import Constraint, Result
+from .constraints import NonNegative, PositiveSemiDefinite, ExponentialCone
+from .error import AlgebraicError
+from .problems import SOSProblem, InternalSOSProblem
+from .variable import OptVariable, init_opt_variable_expr, from_name as opt_variable_from_name
+from .utils import partition
+
+
+class Canonicalization:
+    """ Transform a problem into another problem. """
+
+    def __str__(self):
+        s = f"Apply canonicalization procedure {self.__class__.__qualname__} to\n"
+        return s + str(self.problem)
+
+    @property
+    @abstractmethod
+    def problem(self) -> SOSProblem:
+        """ The problem that will be canonicalized """
+
+    def solve(self, verbose: bool = False, state: ExpressionState | None = None) -> Result:
+        if state is None:
+            state = poly.make_state()
+
+        state, internal_prob = self.apply(state)
+        return internal_prob.solve(verbose)
+
+    @abstractmethod
+    def apply(self, state: ExpressionState) -> tuple[ExpressionState, InternalSOSProblem]:
+        """ Apply the canonicalization procedure """
+
+
+@dataclassabc(frozen=True)
+class PutinarPSatz(Canonicalization):
+    problem: SOSProblem
+
+    @staticmethod
+    def make_multiplier(variables: Sequence[OptVariable],
+                        polynomial_variables: Sequence[Variable],
+                        constr: Expression,
+                        domain_poly: Expression,
+                        multiplier_name: str) -> Expression:
+        r"""
+        Make a multiplier polynomial for a domain polynomial of a constraint.
+
+        Since Putinar's PSatz does not explicity tell us the order of the
+        multiplier this function will make a multiplier, let's call it
+        :math:`\gamma`, such that :math:`\deg(\gamma * p) = \deg(c)`, wherein
+        :math:`p` is the domain polynomial and :math:`c` is the contraint
+        polynomial.
+        """
+        # Number of variables in the original problem
+        q = sum(math.prod(v.shape) for v in variables)
+        x = poly.v_stack((1,) + tuple(polynomial_variables))
+
+        def make_multiplier_later(state: ExpressionState) -> tuple[ExpressionState, PolyMatrixMixin]:
+            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
+            # Coefficients of the multiplier polynomial 
+            c = opt_variable_from_name(multiplier_name, shape=(1,n))
+            multiplier = c @ x.combinations(tuple(range(d +1)))
+
+            return multiplier.apply(state)
+        return poly.from_statemonad(init_state_monad(make_multiplier_later))
+
+    # TODO: rewrite this function, this is adapted from how it was done
+    # previously, and it works but it's not efficient.
+    @abstractmethod
+    def apply(self, state: ExpressionState) -> tuple[ExpressionState, InternalSOSProblem]:
+
+        def need_positivstellensatz(c: Constraint) -> bool:
+            return isinstance(c, NonNegative) and (c.domain is not None)
+
+        constraints, to_process = partition(need_positivstellensatz, self.problem.constraints)
+        constraints = tuple(constraints) # because it is a generator
+
+        state, prob = self.problem.apply(state)
+
+        new_constraints: list[Constraint] = list(constraints)
+        for i, constr in enumerate(to_process):
+            new_constr = constr.expression
+            for domain_poly in constr.domain.polynomials:
+                m = PutinarPSatz.make_multiplier(prob.variables,
+                                                 prob.polynomial_variables,
+                                                 constr, domain_poly, rf"\gamma_{i}")
+
+                # Subtract multiplier from expression and impose that it is also SOS
+                new_constr -= m * domain_poly
+                new_constraints.append(NonNegative(m))
+
+            new_constraints.append(NonNegative(new_constr))
+
+        return replace(self.problem, constraints=new_constraints).apply(state)
 
-def log_det_canon(prob: Problem) -> Problem:
-    raise NotImplementedError
 
diff --git a/sumofsquares/constraints.py b/sumofsquares/constraints.py
index 71dbb06..3d681c5 100644
--- a/sumofsquares/constraints.py
+++ b/sumofsquares/constraints.py
@@ -1,30 +1,24 @@
-import math
-
 from dataclasses import dataclass
 from dataclassabc import dataclassabc
 from typing_extensions import override
-from typing import Iterable
-
-import polymatrix as poly
+from typing import Generic, TypeVar
 
 from polymatrix.expression.expression import Expression
-from polymatrix.expressionstate import ExpressionState
-from polymatrix.polymatrix.abc import PolyMatrix
-from polymatrix.statemonad import init_state_monad
-from polymatrix.variable import Variable
+from polymatrix.polymatrix.mixins import PolyMatrixMixin
 
 from .abc import Set, Constraint
-from .error import AlgebraicError
-from .variable import OptVariable, from_name as opt_variable_from_name
+
+
+# Type variable for expression or polymatrix
+E = TypeVar("E", bound=Expression | PolyMatrixMixin)
 
 
 # ┏━┓┏━╸╺┳╸┏━┓
 # ┗━┓┣╸  ┃ ┗━┓
 # ┗━┛┗━╸ ╹ ┗━┛
 
-
 @dataclass(frozen=True)
-class BasicSemialgebraicSet(Set):
+class BasicSemialgebraicSet(Set, Generic[E]):
     r"""
     A set that is described by the intersection of the positive loci of a
     finite number of polynomials (under some conditions). Roughly speaking this
@@ -41,7 +35,7 @@ class BasicSemialgebraicSet(Set):
 
     This tuple stores the :math:`g_i(x)`.
     """
-    polynomials: tuple[Expression, ...]
+    polynomials: tuple[E, ...]
 
     def __str__(self):
         polynomials = ', '.join(map(str, self.polynomials))
@@ -54,9 +48,9 @@ class BasicSemialgebraicSet(Set):
 
 
 @dataclassabc(frozen=True)
-class EqualToZero(Constraint):
+class EqualToZero(Constraint[E]):
     """ Constrain an expression to be equal to zero. """
-    expression: Expression
+    expression: E
 
     @override
     def __str__(self):
@@ -64,10 +58,10 @@ class EqualToZero(Constraint):
 
 
 @dataclassabc(frozen=True)
-class NonNegative(Constraint):
+class NonNegative(Constraint[E]):
     """
-    Constrain an expression to be non-negative in a certain domain. If the
-    domain in set to `None` then it is interpreted as being non-negative
+    Constrain a scalar expression to be non-negative in a certain domain. If
+    the domain in set to `None` then it is interpreted as being non-negative
     everywhere.
 
     **Note:** In this package, non-negativitiy constraint will always
@@ -75,8 +69,8 @@ class NonNegative(Constraint):
     be written as a sum-of-squares, hence this constraint can also be
     understood as `expression` being SOS.
     """
-    expression: Expression 
-    domain: BasicSemialgebraicSet | None = None
+    expression: E 
+    domain: BasicSemialgebraicSet[E] | None = None
 
     @override
     def __str__(self):
@@ -87,72 +81,11 @@ class NonNegative(Constraint):
 
 
 @dataclassabc(frozen=True)
-class PositiveSemiDefinite(Constraint):
+class PositiveSemiDefinite(Constraint[E]):
     """ Constrain a matrix to be positive semidefinite. """
-    expression: Expression 
+    expression: E
 
 
 @dataclassabc(frozen=True)
-class ExponentialCone(Constraint):
-    expression: Expression
-
-
-# ┏━┓┏━┓┏━┓╻╺┳╸╻╻ ╻┏━┓╺┳╸┏━╸╻  ╻  ┏━╸┏┓╻┏━┓┏━┓┏━╸╺┳╸╺━┓┏━╸
-# ┣━┛┃ ┃┗━┓┃ ┃ ┃┃┏┛┗━┓ ┃ ┣╸ ┃  ┃  ┣╸ ┃┗┫┗━┓┣━┫┣╸  ┃ ┏━┛┣╸ 
-# ╹  ┗━┛┗━┛╹ ╹ ╹┗┛ ┗━┛ ╹ ┗━╸┗━╸┗━╸┗━╸╹ ╹┗━┛╹ ╹┗━╸ ╹ ┗━╸┗━╸
-
-def putinar_psatz(
-        constr: NonNegative,
-        variables: Iterable[OptVariable],
-        polynomial_variables: Iterable[Variable],
-        multiplier_name: str
-    ) -> Iterable[Constraint]:
-    """
-    Apply Putinar's positivstellensatz to a non-negativity constraint over an
-    archimedean domain.
-    """
-    if not constr.domain:
-        # TODO: maybe raise an exception to warn the user?
-        return (constr,)
-
-    # Number of variables in the original problem
-    q = sum(math.prod(v.shape) for v in variables)
-    x = poly.v_stack((1,) + tuple(polynomial_variables))
-
-    new_constraints = []
-    new_constr = constr.expression
-
-    for domain_poly in constr.domain.polynomials:
-        # To know the degree of the domain polynomial we need to evaluate the
-        # expression, hence we use a monad here
-        def make_multiplier(state: ExpressionState) -> tuple[ExpressionState, PolyMatrix]:
-            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
-            # Coefficients of the multiplier polynomial 
-            c = opt_variable_from_name(multiplier_name, shape=(1,n))
-            multiplier = c @ x.combinations(tuple(range(d +1)))
-
-            return multiplier.apply(state)
-
-        m = poly.from_statemonad(init_state_monad(make_multiplier))
-
-        # Subtract multiplier from expression and impose that it is also SOS
-        new_constr -= m * domain_poly
-        new_constraints.append(NonNegative(m))
-
-    new_constraints.append(NonNegative(new_constr))
-    return new_constraints
+class ExponentialCone(Constraint[E]):
+    expression: E
diff --git a/sumofsquares/expression.py b/sumofsquares/expression.py
index 08b4fbd..28f1794 100644
--- a/sumofsquares/expression.py
+++ b/sumofsquares/expression.py
@@ -6,18 +6,30 @@ from typing_extensions import override
 from dataclassabc import dataclassabc
 
 from polymatrix.expression.expression import ExpressionBaseMixin
-from polymatrix.expressionstate.abc import ExpressionState
+from polymatrix.expression.from_ import from_statemonad
+from polymatrix.expressionstate import ExpressionState
+from polymatrix.statemonad import init_state_monad
+from polymatrix.polymatrix.mixins import PolyMatrixMixin
+
+from polymatrix.expression.init import (
+        init_concatenate_expr,
+        init_diag_expr,
+        init_lower_triangular_expr,
+        init_transpose_expr,
+        init_slice_expr)
 
 from .abc import Constraint
-from .constraints import NonNegative
+from .constraints import PositiveSemidefinite, ExponentialCone
+from .error import SolverError
+from .variable import init_opt_variable_expr
 
 
 class SOSExpressionBaseMixin(ExpressionBaseMixin):
     @override
     def apply(self, state: ExpressionState) -> tuple[ExpressionState, ExpressionBaseMixin]:
-        # FIXME: think of name for this exception
-        raise RuntimeError
-
+        raise SolverError(f"Expression containing {self.__class__.__qualname__} "
+                           "cannot be used directly, they need to be rewritten "
+                           "into an equivalent form using a canonicalization function.")
 
     @abstractmethod
     def recast(self) -> tuple[ExpressionBaseMixin, Iterable[Constraint]]:
@@ -42,17 +54,79 @@ class LogDetMixin(SOSExpressionBaseMixin):
     def recast(self) -> tuple[ExpressionBaseMixin, Iterable[Constraint]]:
         # The problem
         #
-        #   minimize logdet(A)
+        #   maximize logdet(A)
         #
         # is equivalent to solving
         #
-        #   minimize      -t
+        #   maximize      t
+        #
+        #   subject to    [ A    Z       ]
+        #                 [ Z.T  diag(Z) ] >= 0
+        #
+        #                 Z lower triangular
+        #                 t <= sum_i log(Z[i,i])
+        #
+        # and the last constraint of the above is equivalent to
+        #
+        #   t <= sum_i u[i]
+        #   u_i <= log(Z[i, i])  for all i
+        #
+        # And finally to get rid of the log the latter constraint one is
+        # equivalent to
+        #
+        #   (Z[i,i], 1, u[i]) in Exponential Cone for all i
+        #
+        # Hence we can replace the original problem with
+        #
+        #   minimize      - sum_i u[i]
         #
         #   subject to    [ A    Z       ]
         #                 [ Z.T  diag(Z) ] >= 0
         #
         #                 Z lower triangular
-        #                 t <= sum_i log(Z[i, i])
+        #                 (Z[i,i], 1, u[i]) in ExpCone for all i
+
+        A = self.underlying
+
+        # FIXME: get rid of these functions, create ShapeExprMixin in polymatrix?
+        def make_u(state: ExpressionState) -> tuple[ExpressionState, PolyMatrixMixin]:
+            state, pm = A.apply(state)
+            n, m = pm.shape
+            if n != m:
+                raise ValueError(f"Matrix A of logdet(A) must be square, "
+                                 f"but it has shape {pm.shape}")
+
+            # FIXME: should check to avoid name clashes
+            u = init_opt_variable_expr("u_logdet", shape=(n, 1))
+            return u.apply(state)
+
+        def make_z(state: ExpressionState) -> tuple[ExpressionState, PolyMatrixMixin]:
+            state, pm = A.apply(state)
+            n, m = pm.shape
+            if n != m:
+                raise ValueError(f"Matrix A of logdet(A) must be square, "
+                                 f"but it has shape {pm.shape}")
+
+            # FIXME: should check to avoid name clashes
+            Z = init_lower_triangular_expr(init_opt_variable_expr("Z_logdet", shape=(n * (n + 1) // 2, 1)))
+            return Z.apply(state)
+
+        Z = from_statemonad(init_state_monad(make_z))
+        Z_T = init_transpose_expr(Z)
+        Z_diag = init_diag_expr(Z)
+
+        u = from_statemonad(init_state_monad(make_u))
+
+        # we call the new big matrix Q
+        Q = init_concatenate_expr(((A, Z), (Z_T, Z_diag)))
+
+        def make_expcones(state: ExpressionState) -> tuple[ExpressionState, tuple[Constraint]]:
+            pass
+
+
+        constraints = [PositiveSemidefinite(Q)]
+
+        raise NotImplementedError
 
 
 @dataclassabc(froze=True)
diff --git a/sumofsquares/problems.py b/sumofsquares/problems.py
index 500d06a..6161641 100644
--- a/sumofsquares/problems.py
+++ b/sumofsquares/problems.py
@@ -1,31 +1,91 @@
 """ 
 Specific optimization problems that are related to Sum of Squares.
+
+TODO: explanation for InternalSOSProblem
 """
 from __future__ import annotations
 
-import math
 import polymatrix as poly
+import numpy as np
 
 from dataclassabc import dataclassabc
 from dataclasses import replace
-from itertools import chain
+from numpy.typing import NDArray
 from typing import Any, Sequence
 from typing_extensions import override
 
 from polymatrix.expression.expression import Expression, VariableExpression, init_expression
 from polymatrix.expressionstate import ExpressionState
+from polymatrix.polymatrix.mixins import PolyMatrixMixin
+from polymatrix.polymatrix.index import MonomialIndex, VariableIndex
 from polymatrix.variable import Variable
 
 from .abc import Problem, Constraint, Solver, Result
-from .constraints import NonNegative, putinar_psatz
-from .error import AlgebraicError
+from .constraints import NonNegative, EqualToZero, PositiveSemiDefinite, ExponentialCone
 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
+from .variable import OptVariable
+
+
+# ┏━╸┏━┓┏┓╻╻┏━╸   ┏━┓┏━┓┏━┓┏┓ ╻  ┏━╸┏┳┓
+# ┃  ┃ ┃┃┗┫┃┃     ┣━┛┣┳┛┃ ┃┣┻┓┃  ┣╸ ┃┃┃
+# ┗━╸┗━┛╹ ╹╹┗━╸   ╹  ╹┗╸┗━┛┗━┛┗━╸┗━╸╹ ╹
+
+
+@dataclassabc(frozen=True)
+class ConicProblem(Problem):
+    """
+    Conic program with linear or quadratic cost functions.
+    All values stored in this class are numerical (numpy arrays).
+
+    The pproblem has the following form
+    ::
+
+        minimize      .5 * x.T @ P @ x + q.T @ x
+        
+        subject to    G @ x <= h
+                      A @ x = b
+
+    """
+    P: NDArray
+    q: NDArray
+    A: NDArray
+    b: NDArray
+    G: NDArray
+    h: NDArray
+
+    dims: dict
+    is_qp: bool
+    solver: Solver
+    variables: Sequence[OptVariable]
+
+    @property
+    @override
+    def cost(self) -> tuple[NDArray, NDArray]:
+        return self.P, self.q
+
+    @property
+    @override
+    def constraints(self) -> dict[str, tuple[NDArray, NDArray]]:
+        return {
+            "eq": (self.A, self.b),
+            "ineq": (self.G, self.h)
+        }
+
+    @override
+    def solve(self, verbose: bool = False) -> ConicResult:
+        match self.solver:
+            case Solver.CVXOPT:
+                result, info = cvxopt_solve_sos_cone(self, verbose)
+            case _:
+                raise NotImplementedError
+
+        return ConicResult(values=result, solver_info=info)
+
 
 
 @dataclassabc(frozen=True)
-class SOSResult(Result):
+class ConicResult(Result):
     values: dict[OptVariable, float]
     solver_info: Any
 
@@ -46,75 +106,234 @@ class SOSResult(Result):
         return self.values[var]
 
 
+# ┏━┓╻ ╻┏┳┓   ┏━┓┏━╸   ┏━┓┏━┓╻ ╻┏━┓┏━┓┏━╸┏━┓   ┏━┓┏━┓┏━┓┏━╸┏━┓┏━┓┏┳┓
+# ┗━┓┃ ┃┃┃┃   ┃ ┃┣╸    ┗━┓┃┓┃┃ ┃┣━┫┣┳┛┣╸ ┗━┓   ┣━┛┣┳┛┃ ┃┃╺┓┣┳┛┣━┫┃┃┃
+# ┗━┛┗━┛╹ ╹   ┗━┛╹     ┗━┛┗┻┛┗━┛╹ ╹╹┗╸┗━╸┗━┛   ╹  ╹┗╸┗━┛┗━┛╹┗╸╹ ╹╹ ╹
+
+
 @dataclassabc(frozen=True)
-class SchmuedgenSOSProblem(Problem):
-    r"""
-    SOS problem with quadratic cost, linear constraints and optinally with
-    positivity constraints over a basic closed semialgebraic set
-    :math:`\mathbf{K}`. To constrain the positivity over :math:`\mathbf{K}`
-    this problem uses Schmüdgen's Positivstellensatz.
+class SOSProblem(Problem):
+    """
+    Generic sum of squares problem.
+    This problem contains expression objects.
     """
     cost: Expression
-    constraints: Sequence[Constraint]
-    variables: Sequence[OptVariable]
-    polynomial_variables: Sequence[Variable]
+    constraints: Sequence[Constraint[Expression]]
     solver: Solver
-    state: ExpressionState
 
+    # --- dunder methods ---
+
+    def __str__(self):
+        n = len(self.constraints)
+        s = f"SumOfSquares Problem ({self.__class__.__qualname__}):\n\n  " \
+                f"minimize  {self.cost}\n  " \
+                f"subject to ({n} constraint{'s' if n > 1 else ''})\n"
+
+        for c in self.constraints:
+            s += f"\t* {str(c)}\n"
+
+        return s
+
+    @override
+    def solve(self, verbose: bool = False, state: ExpressionState | None = None) -> Result:
+        if state is None:
+            state = poly.make_state()
+
+        state, internal_prob = self.apply(state)
+        return internal_prob.solve(verbose)
+
+    def apply(self, state: ExpressionState) -> tuple[ExpressionState, InternalSOSProblem]:
+        """ Convert to internal SOS problem by applying state to the expressions. """
+
+        constraints: list[Constraint[PolyMatrixMixin]] = []
+        state, cost = self.cost.apply(state)
+
+        # Compute the polymatrix of each constraint expression. Even though the
+        # result may not be used, this is necessary to account for the case of
+        # variables that are only present in the constraints, which the state
+        # object may not contain yet.
+        # TODO: Is there a more efficient way? There has to be a better way.
+        for c in self.constraints:
+            state, _  = c.expression.apply(state)
+
+        # Collect variables
+        def is_optvariable(v):
+            return isinstance(v, OptVariable)
+
+        polynomial_variables, variables = partition(is_optvariable, state.indices.keys())
+        polynomial_variables = tuple(polynomial_variables) # because it is a generator
 
+        x = poly.v_stack(polynomial_variables)
+        for c in self.constraints:
+            # Scalar polynomial equality and non-negativity contraints need to
+            # be converted to (linear) coefficient matching constraints
+            if isinstance(c, NonNegative | EqualToZero):
+                state, pm = c.expression.quadratic_in(x).apply(state)
+
+            elif isinstance(c, PositiveSemiDefinite):
+                state, pm = c.expression.apply(state)
+
+            elif isinstance(c, ExponentialCone):
+                raise NotImplementedError
+
+            else:
+                raise NotImplementedError
+
+            constraints.append(replace(c, expression=pm))
+
+        return state, InternalSOSProblem(cost, tuple(constraints),
+                                         tuple(variables), polynomial_variables,
+                                         self.solver, state)
+
+
+# FIXME: I hate this name, and the fact that this class needs to exist
 @dataclassabc(frozen=True)
-class PutinarSOSProblem(Problem):
-    r"""
-    SOS problem with quadratic cost, linear constraints and optinally with
-    positivity constraints over an archimedean set :math:`\mathbf{K}`. To
-    constrain the positivity over :math:`\mathbf{K}` this problem uses
-    Putinar's positivstellensatz.
+class InternalSOSProblem(Problem):
     """
-    cost: Expression
-    constraints: Sequence[Constraint]
+    This is a class for internal use only. It is a SOS problem that contains
+    polymatrix objects.
+
+    It exists only because expressions are lazily evaluated, so if we want to
+    perform transformations on problems (see canon module), we may need to
+    apply the state to the expressions to convert the expressions to polymatrix
+    objects.
+    """
+    cost: PolyMatrixMixin
+    constraints: Sequence[Constraint[PolyMatrixMixin]]
     variables: Sequence[OptVariable]
     polynomial_variables: Sequence[Variable]
     solver: Solver
-    state: ExpressionState
-
-    def after_positivstellensatz(self) -> Problem:
-        """
-        Apply Putinar's positivstellensatz to convert the SOS constraints
-        into coefficient matching constraints 
-        """
-        def need_positivstellensatz(c: Constraint) -> bool:
-            return isinstance(c, NonNegative) and (c.domain is not None)
-
-        # Add langrange-multiplier-like polynomials to SOS constraints, while
-        # leaving the rest untouched.
-        constraints, to_process = partition(need_positivstellensatz, self.constraints)
-        new_constraints: list[Constraint] = []
-        for i, c in enumerate(to_process):
-            new_constraints.extend(putinar_psatz(c, self.variables, self.polynomial_variables, r"\gamma_{i}"))
-
-        # Convert to polymatrix object to extract variables created by p-satz
-        # FIXME: this is a dirty trick to make it work for now, it is not
-        # efficient, is there a better way?
-        new_state = self.state
-        for c in new_constraints:
-            new_state, _ = c.expression.apply(new_state)
+    # TODO: remove state field from this class, it is redundant
+    state: ExpressionState 
 
-        def is_optvariable(v):
-            return isinstance(v, OptVariable)
+    def to_conic_problem(self) -> ConicProblem:
+        cost = poly.to_affine(self.cost)
+        if cost.degree > 2:
+            raise ValueError("This package can solve at most quadratic conic programs, "
+                             f"but the given problem has degree {cost.degree} > 2.")
+
+        is_qp = (cost.degree == 2)
+
+        A_rows, b_rows = [], []
+        G_rows, h_rows= [], []
+        l_dims, q_dims, s_dims = [], [], []
+
+        # indices of variables in the optimization problem
+        variable_indices = sum(sorted(tuple(self.state.get_indices_as_variable_index(v))
+                                      for v in self.variables), ())
+
+        # cost linear term
+        cost_coeffs = dict(cost.affine_coefficients_by_degrees(variable_indices))
+        q = cost_coeffs.get(1) 
+        if q is None and is_qp:
+            q = np.zeros((len(variable_indices), 1))
+
+        # cost quadratic term
+        if is_qp:
+            n = len(variable_indices)
+            P = np.zeros((n, n))
+            # This works because of how monomial indices are sorted
+            # see also polymatrix.index.MonomialIndex.__lt__
+            P[np.triu_indices(n)] = cost_coeffs[2]
+
+            # Make symmetric. There is no 0.5 factor because there it is already
+            # there in the canonical form of CVXOPT
+            P = (P + P.T)
+
+        else:
+            P = None
+
+        for c in self.constraints:
+            if isinstance(c, PositiveSemiDefinite):
+                constr = poly.to_affine(c.expression)
 
-        new_variables = filter(is_optvariable, new_state.indices.keys())
+                nrows, ncols = constr.shape
+                if nrows != ncols:
+                    raise ValueError(f"PSD constraint cannot contain non-square matrix of shape ({nrows, ncols})!")
 
-        return replace(self,
-                       state=new_state,
-                       variables=tuple(new_variables),
-                       constraints=tuple(constraints) + tuple(new_constraints))
+                s_dims.append(nrows)
+
+                # constant term
+                c_coeff = constr.affine_coefficient(MonomialIndex.constant())
+                c_stacked = c_coeff.T.reshape((nrows * ncols,))
+
+                # linear terms
+                l_stacked_rows = []
+                for v in variable_indices:
+                    v_coeff = constr.affine_coefficient(v)
+                    v_stacked = v_coeff.T.reshape((nrows * ncols,))
+                    l_stacked_rows.append(v_stacked)
+
+                l_stacked = np.vstack(l_stacked_rows)
+
+                h_rows.append(c_stacked)
+                G_rows.append(-1 * l_stacked)
+
+            elif isinstance(c, ExponentialCone):
+                # TODO: implement
+                raise NotImplementedError
+
+            elif isinstance(c, EqualToZero | NonNegative):
+                constr = poly.to_affine(c.expression)
+                if constr.degree > 1:
+                    raise ValueError("To convert to conic constraints must be linear or affine "
+                                    f"but {str(c.expression)} has degree {constr.degree}.")
+
+                # For each constraint we add PSD matrix
+                nrows, ncols = constr.shape
+                s_dims.append(nrows)
+
+                # Constant term
+                c_coeff = constr.affine_coefficient(MonomialIndex.constant())
+                c_stacked = c_coeff.T.reshape((nrows * ncols,))
+
+                # Linear term
+                l_stacked_rows = []
+                for v in variable_indices:
+                    # Get the affine coefficient associated to this variable
+                    v_coeff = constr.affine_coefficient(v)
+
+                    # stack the columns of coeff matrix into a fat row vector (column major order).
+                    v_stacked = v_coeff.T.reshape((nrows * ncols,))
+                    l_stacked_rows.append(v_stacked)
+
+                l_stacked = np.vstack(l_stacked_rows)
+
+                if isinstance(c, EqualToZero):
+                    A_rows.append(l_stacked)
+                    # Factor of -1 because it is on the right hand side
+                    b_rows.append(-1 * c_stacked)
+
+                elif isinstance(c, NonNegative):
+                    if c.domain:
+                        raise ValueError("Cannot convert non-negativity constraint to conic constraint. "
+                                         "Domain restricted non-negativity must be "
+                                         "preprocessed by using a Positivstellensatz!")
+
+                    # Factor of -1 because it is on the right hand side
+                    G_rows.append(-1 * l_stacked)
+                    h_rows.append(c_stacked)
+
+            else:
+                raise NotImplementedError(f"Cannot convert constraint of type {type(c)} "
+                                          "into a conic constriant.")
+
+        G = np.hstack(G_rows).T if G_rows else None
+        h = np.hstack(h_rows) if h_rows else None
+        A = np.hstack(A_rows).T if A_rows else None
+        b = np.hstack(b_rows) if b_rows else None
+
+        dims = {'l': l_dims, 'q': q_dims, 's': s_dims}
+
+        if all((G is None, h is None, A is None, b is None)):
+            raise ValueError("Optimization problem is unconstrained!")
+
+        return ConicProblem(P=P, q=q, A=A, b=b, G=G, h=h,
+                            dims=dims, is_qp=is_qp,
+                            solver=self.solver,
+                            variables=self.variables)
 
 
     @override
     def solve(self, verbose: bool = False) -> Result:
-        if not self.solver == Solver.CVXOPT:
-            raise NotImplementedError(f"Cannot solve problem {self} using solver {str(self.solver)}.")
-
-        prob = self.after_positivstellensatz()
-        result, info = cvxopt_solve_sos_cone(prob, verbose)
-        return SOSResult(result, info)
+        return self.to_conic_problem().solve(verbose)
diff --git a/sumofsquares/solver/cvxopt.py b/sumofsquares/solver/cvxopt.py
index 24ec413..6418a2d 100644
--- a/sumofsquares/solver/cvxopt.py
+++ b/sumofsquares/solver/cvxopt.py
@@ -1,11 +1,14 @@
 """
 Solve sumofsquares problems using CVXOPT
 """
+from __future__ import annotations
+
 import cvxopt
 import numpy as np
 import math
 
 from collections import UserDict
+from typing import TYPE_CHECKING
 
 import polymatrix as poly
 from polymatrix.expression.expression import Expression
@@ -16,39 +19,28 @@ from ..constraints import NonNegative, EqualToZero, PositiveSemiDefinite, Expone
 from ..error import NotSupportedBySolver, SolverError
 from ..variable import OptVariable
 
+if TYPE_CHECKING:
+    from ..problems import ConicProblem
+
 
 class CVXOPTInfo(SolverInfo, UserDict):
     """ Dictionary returned by CVXOPT. """
 
     def __getattr__(self, attr):
+        # trick to make it behave like a dataclass
         key = " ".join(attr.split("_"))
         if key not in self.keys():
             raise KeyError(f"Key {key} was not found")
         return self[key]
 
 
-def solve_sos_cone(prob: Problem, verbose: bool = False,
+def solve_sos_cone(prob: ConicProblem, verbose: bool = False,
                    *args, **kwargs) -> tuple[dict[OptVariable, float], CVXOPTInfo]:
     r"""
-    Solve a conic problem in the cone of SOS polynomials
-    :math:`\mathbf{\Sigma}_d(x)` using CVXOPT.
-
     Any `*args` and `**kwargs` other than `prob` are passed to the CVXOPT
     solver.
     """
-    # Check that the problem can be solved by CVXOpt, i.e.
-    #  - Cost function is at most quadratic
-    #  - Constraints are at most linear 
-
-    cost = poly.to_affine(prob.cost.read(prob.state))
-    if cost.degree > 2:
-        raise ValueError("CVXOpt can solve at most quadratic conic programs, "
-                         f"but the given problem has degree {cost.degree} > 2.")
-
-    is_qp = (cost.degree == 2)
-
-    # We need to convert problem to the primal canonical form (given in the
-    # CVXOPT doc) with x = vstack(prob.variables)
+    # CVXOPT can solve problems that have the form
     #
     #    minimize      .5 * x.T @ P @ x + q.T @ x
     #
@@ -56,147 +48,28 @@ def solve_sos_cone(prob: Problem, verbose: bool = False,
     #                  A @ x = b
     #                  s >= 0
     #
-    # If the problem is an LP, the CVXOPT doc call the coefficient in the cost
-    # functions `c`, here we will use `q` for both LP and QP.
-
-    A_rows, b_rows = [], []
-    G_rows, h_rows= [], []
-    l_dims, q_dims, s_dims = [], [], []
-
-    # indices of variables in the optimization problem
-    variable_indices = sum(sorted(tuple(prob.state.get_indices_as_variable_index(v))
-                                  for v in prob.variables), ())
-
-    # cost linear term
-    cost_coeffs = dict(cost.affine_coefficients_by_degrees(variable_indices))
-    q = cost_coeffs.get(1) 
-    if q is None and is_qp:
-        q = np.zeros((len(variable_indices), 1))
-
-    # cost quadratic term
-    if is_qp:
-        n = len(variable_indices)
-        P = np.zeros((n, n))
-        # This works because of how monomial indices are sorted
-        # see also polymatrix.index.MonomialIndex.__lt__
-        P[np.triu_indices(n)] = cost_coeffs[2]
-
-        # Make symmetric. There is no 0.5 factor because there it is already
-        # there in the canonical form of CVXOPT
-        P = (P + P.T)
-
-    x = poly.v_stack(prob.polynomial_variables)
-    for c in prob.constraints:
-        if isinstance(c, PositiveSemiDefinite):
-            constr = poly.to_affine(c.expression.read(prob.state))
-
-            nrows, ncols = constr.shape
-            if nrows != ncols:
-                raise ValueError(f"PSD constraint cannot contain non-square matrix of shape ({nrows, ncols})!")
-
-            s_dims.append(nrows)
-
-            # constant term
-            c_coeff = constr.affine_coefficient(MonomialIndex.constant())
-            c_stacked = c_coeff.T.reshape((nrows * ncols,))
-
-            # linear terms
-            l_stacked_rows = []
-            for v in variable_indices:
-                v_coeff = constr.affine_coefficient(v)
-                v_stacked = v_coeff.T.reshape((nrows * ncols,))
-                l_stacked_rows.append(v_stacked)
-
-            l_stacked = np.vstack(l_stacked_rows)
-
-            h_rows.append(c_stacked)
-            G_rows.append(-1 * l_stacked)
-
-        elif isinstance(c, ExponentialCone):
-            raise NotSupportedBySolver("CVXOpt cannot solve problems with the exponential cone.")
-
-        elif isinstance(c, EqualToZero | NonNegative):
-            # Convert constraints to affine expression in decision variables
-            # TODO: Use half Newton polytope optimization to reduce size of marix problem
-            constr = poly.to_affine(c.expression.quadratic_in(x).read(prob.state))
-            if constr.degree > 1:
-                raise ValueError("CVXOpt can solve problem with constraints that are "
-                                f"at most linear, but {str(c.expression)} has degree {constr.degree}.")
-
-            # For each constraint we add PSD matrix
-            nrows, ncols = constr.shape
-            s_dims.append(nrows)
-
-            # Constant term
-            c_coeff = constr.affine_coefficient(MonomialIndex.constant())
-            c_stacked = c_coeff.T.reshape((nrows * ncols,))
-
-            # Linear term
-            l_stacked_rows = []
-            for v in variable_indices:
-                # Get the affine coefficient associated to this variable
-                v_coeff = constr.affine_coefficient(v)
-
-                # stack the columns of coeff matrix into a fat row vector (column major order).
-                v_stacked = v_coeff.T.reshape((nrows * ncols,))
-                l_stacked_rows.append(v_stacked)
-
-            l_stacked = np.vstack(l_stacked_rows)
-
-            if isinstance(c, EqualToZero):
-                A_rows.append(l_stacked)
-                # Factor of -1 because it is on the right hand side
-                b_rows.append(-1 * c_stacked)
-
-            elif isinstance(c, NonNegative):
-                if c.domain:
-                    raise ValueError("Cannot convert non-negativity constraint to CVXOPT format. "
-                                     "Domain restricted non-negativity must be "
-                                     "preprocessed by using a Positivstellensatz!")
-
-                # Factor of -1 because it is on the right hand side
-                G_rows.append(-1 * l_stacked)
-                h_rows.append(c_stacked)
-
-        else:
-            raise NotImplementedError(f"Cannot convert constraint of type {type(c)} "
-                                      "into the canonical form of CVXOPT.")
-
-    # Convert to CVXOPT matrices
-    q = cvxopt.matrix(q.T) if q is not None else None
-    G = cvxopt.matrix(np.hstack(G_rows).T) if G_rows else None
-    h = cvxopt.matrix(np.hstack(h_rows)) if h_rows else None
-    A = cvxopt.matrix(np.hstack(A_rows).T) if A_rows else None
-    b = cvxopt.matrix(np.hstack(b_rows)) if b_rows else None
-    dims = {'l': l_dims, 'q': q_dims, 's': s_dims}
-
-    # print(q)
-    # print(G)
-    # print(h)
-    # print(A)
-    # print(b)
-    # print(dims)
-
-    cvxopt.solvers.options['show_progress'] = verbose
-
-    if not any((G, h, A, b)):
-        raise ValueError("Optimization problem is unconstrained!")
-
-    # Add matrices for debugging
-    # info = {}
-    info = {"q": q, "G": G, "h": h, "A": A, "b": b, "dims": dims}
-
+    # This matches the structure of ConicProblem, so the values can be used
+    # directly.
     try:
-        if is_qp:
-            P = cvxopt.matrix(P)
-            info |= cvxopt.solvers.coneqp(P=P, q=q, G=G, h=h, A=A, b=b,
-                                         dims=dims, *args, **kwargs)
+        # Some are transposed because CVXOPT uses column major arrays
+        q = cvxopt.matrix(prob.q.T)
+        A = cvxopt.matrix(prob.A) if prob.A is not None else None
+        b = cvxopt.matrix(prob.b.T) if prob.b is not None else None
+        G = cvxopt.matrix(prob.G) if prob.G is not None else None
+        h = cvxopt.matrix(prob.h.T) if prob.h is not None else None
+
+        if prob.is_qp:
+            P = cvxopt.matrix(prob.P)
+
+            info = cvxopt.solvers.coneqp(P=P, q=q, G=G, h=h, A=A, b=b,
+                                         dims=prob.dims, *args, **kwargs)
         else:
-            info |= cvxopt.solvers.conelp(c=q, G=G, h=h, A=A, b=b,
+            info = cvxopt.solvers.conelp(c=q, G=G, h=h, A=A, b=b,
                                          *args, **kwargs)
 
     except Exception as e:
-        raise SolverError("CVXOpt can't solve this problem.") from e
+        raise SolverError("CVXOpt can't solve this problem, "
+                          "see previous exception why.") from e
 
     results = {}