Implement ExponentialCone in SOSProblem.apply, clean up comments

2 files changed, 40 insertions, 22 deletions

diff --git a/sumofsquares/problems.py b/sumofsquares/problems.py
index 4482c1e..c313419 100644
--- a/sumofsquares/problems.py
+++ b/sumofsquares/problems.py
@@ -67,13 +67,15 @@ class ConicProblem(Problem):
     types of constraints
     ::
 
-       z    zero cone                   (s = 0)
-       l    linear cone                 (s >= 0)
-       b    box cone                    (tl <= s <= tu)
-       q    second order cone           (|s|_2 <= t)
-       s    positive semidefinite cone  (s is PSD)
-       ep   exponential cone            ((x,y,z) in ExpCone)
-       ep*  dual exponential cone       ((u,v,w) in DualExpCone)
+    NAME    LONG NAME                    MEANING
+    ----    ---------------------------  ----------------------
+       z    zero cone                    s = 0
+       l    linear cone                  s >= 0
+       b    box cone                     tl <= s <= tu
+       q    second order cone            |s|_2 <= t
+       s    positive semidefinite cone   s is PSD
+       e    exponential cone             (x,y,z) in ExpCone
+       e*   dual exponential cone        (u,v,w) in DualExpCone
        p    power cone
        p*   dual power cone
 
@@ -231,6 +233,7 @@ class SOSProblem(Problem):
                 # Polynomial non-negativity contraints is converted to PSD
                 # constraint of SOS quadratic form
                 if deg.scalar().constant() > 1:
+                    # TODO: it seems to work fine even without .symmetric(). Why?
                     state, pm = c.expression.quadratic_in(x).symmetric().apply(state)
                     constraints.append(PositiveSemiDefinite(pm))
 
@@ -239,16 +242,28 @@ class SOSProblem(Problem):
                     state, pm = c.expression.apply(state)
                     constraints.append(replace(c, expression=pm))
 
-            # PSD constraint can be passed as-is
             elif isinstance(c, PositiveSemiDefinite):
                 state, pm = c.expression.apply(state)
+                nrows, ncols = pm.shape
+                if nrows != ncols:
+                    raise ValueError(f"PSD constraint cannot contain non-square matrix of shape ({nrows, ncols})!")
+
+                # PSD constraint can be passed as-is
                 constraints.append(replace(c, expression=pm))
 
             elif isinstance(c, ExponentialCone):
-                raise NotImplementedError
+                state, pm = c.expression.apply(state)
+                nrows, ncols = pm.shape
+
+                if ncols != 3:
+                    raise ValueError("Conic constraint must be a row vector [x, y, z] ",
+                                     "or for multiple constraints it must be an n x 3 "
+                                    f"matrix! Given expression has wrong shape {pm.shape}.")
+
+                constraints.append(replace(c, expresssion=pm))
 
             else:
-                raise NotImplementedError
+                raise NotImplementedError(f"Cannot process constraint of type {type(c)} (yet).")
 
         return state, InternalSOSProblem(cost, tuple(constraints),
                                          tuple(variables), polynomial_variables,
@@ -291,13 +306,13 @@ class InternalSOSProblem(Problem):
         # See docstring of ConicProblem.constraints
         dims: dict[str, list[int]] = {
             "z": [], "l": [], "b": [], "q": [], "s": [],
-            "ep": [], "ep*": [], "p": [], "p*": [],
+            "e": [], "e*": [], "p": [], "p*": [],
         }
         
         # See docstring of ConicProblem.constraints
         constraints: dict[str, list[tuple[tuple[NDArray, ...], NDArray]]] = {
             "z": [], "l": [], "b": [], "q": [], "s": [],
-            "ep": [], "ep*": [], "p": [], "p*": [],
+            "e": [], "e*": [], "p": [], "p*": [],
         }
 
         # indices of variables in the optimization problem, sorted
@@ -352,9 +367,6 @@ class InternalSOSProblem(Problem):
                 constraints["l"].append((linear, constant))
 
             elif isinstance(c, PositiveSemiDefinite):
-                if nrows != ncols:
-                    raise ValueError(f"PSD constraint cannot contain non-square matrix of shape ({nrows, ncols})!")
-
                 dims["s"].append(nrows)
                 constraints["s"].append((linear, constant))
 
@@ -362,7 +374,7 @@ class InternalSOSProblem(Problem):
                 # Interpret row-wise
                 dims["ep"].extend(1 for _ in range(nrows))
                 for i in range(nrows):
-                    constraints["ep"].append((
+                    constraints["e"].append((
                         tuple(m[i, :] for m in linear),
                         constant[i, :]))
 
diff --git a/sumofsquares/solver/cvxopt.py b/sumofsquares/solver/cvxopt.py
index 455ee5b..1c37fbe 100644
--- a/sumofsquares/solver/cvxopt.py
+++ b/sumofsquares/solver/cvxopt.py
@@ -44,10 +44,10 @@ def vectorize_matrix(m: NDArray) -> NDArray:
 def solve_cone(prob: ConicProblem, verbose: bool = False,
                *args, **kwargs) -> tuple[dict[OptVariable, NDArray | float], CVXOPTInfo]:
     r"""
-    Any `*args` and `**kwargs` other than `prob` are passed to the CVXOPT
-    solver.
+    Any `*args` and `**kwargs` other than `prob` and `vebose` are passed to the
+    CVXOPT solver.
     """
-    # CVXOPT can solve problems that have the form
+    # CVXOPT can solve problems that have the (primal) form
     #
     #    minimize      .5 * x.T @ P @ x + q.T @ x
     #
@@ -61,7 +61,7 @@ def solve_cone(prob: ConicProblem, verbose: bool = False,
         raise NotImplementedError("Automatic conversion of box constraints "
                                   "into linear constraints is not supported yet.")
 
-    if prob.constraints["ep"] or prob.constraints["ep*"]:
+    if prob.constraints["e"] or prob.constraints["e*"]:
         raise NotSupportedBySolver("CVXOPT cannot solve problems with the exponential cone.")
 
     if prob.constraints["p"] or prob.constraints["p*"]:
@@ -121,6 +121,13 @@ def solve_cone(prob: ConicProblem, verbose: bool = False,
 
     # Solve the problem
     # pprint({"A": A, "b": b, "G": G, "h": h})
+    # pprint(dims)
+
+    # print(A)
+    # print(b)
+
+    # print(G)
+    # print(h)
 
     try:
         if prob.is_qp:
@@ -128,7 +135,6 @@ 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 q is not None, "LP has no objective! Something is very wrong"
             info = cvxopt.solvers.conelp(c=q, G=G, h=h, A=A, b=b,
                                          *args, **kwargs)
 
@@ -146,7 +152,7 @@ def solve_cone(prob: ConicProblem, verbose: bool = False,
     results, i = {}, 0
     for variable in prob.variables:
         num_indices = math.prod(variable.shape)
-        values = np.array(info['x'][i:i+num_indices]).reshape(variable.shape)
+        values = np.array(info["x"][i:i+num_indices]).reshape(variable.shape)
         if values.shape == (1, 1):
             values = values[0, 0]