From 509fbe7923c67b1cd64a2c6b83372e47ac8308ef Mon Sep 17 00:00:00 2001
From: Nao Pross <np@0hm.ch>
Date: Sun, 16 Jun 2024 19:50:59 +0200
Subject: Fix vec and mat functions in SCS

---
 sumofsquares/solver/cvxopt.py |  2 +-
 sumofsquares/solver/scs.py    | 35 +++++++++++++----------------------
 2 files changed, 14 insertions(+), 23 deletions(-)

diff --git a/sumofsquares/solver/cvxopt.py b/sumofsquares/solver/cvxopt.py
index 11b7797..ef5aaaf 100644
--- a/sumofsquares/solver/cvxopt.py
+++ b/sumofsquares/solver/cvxopt.py
@@ -84,7 +84,7 @@ def solve_cone(prob: ConicProblem, verbose: bool = False,
     # pprint(f"{b_rows = }")
         
     if A_rows:
-        # b is multiplied by -1 because it is on the RHS
+        # multiplied by -1 because it is on the RHS
         b = cvxopt.matrix(np.vstack(b_rows))
         A = cvxopt.matrix(np.vstack(A_rows)) * -1
 
diff --git a/sumofsquares/solver/scs.py b/sumofsquares/solver/scs.py
index 1e474da..662f78e 100644
--- a/sumofsquares/solver/scs.py
+++ b/sumofsquares/solver/scs.py
@@ -29,38 +29,29 @@ class SCSInfo(SolverInfo, UserDict):
             raise KeyError(f"Key {key} was not found")
         return self[key]
 
-
-def vec(m: NDArray) -> NDArray:
+def vec(S: NDArray) -> NDArray:
     r"""
     Vectorize a symmetric matrix for SCS by storing only a half of the entries
     and multiplying the off-diaognal elements by :math:`\sqrt{2}`.
     """
-    n, _ = m.shape
-    # Since it is symmetric we only work with the lower triangular part
-    scaling = np.tril(np.ones(n), k=-1) * np.sqrt(2) + np.eye(n)
-    v = (m * scaling)[np.tril_indices(n)].reshape((-1, 1))
-    assert v.shape[0] == n * (n + 1) // 2, "SCS Matrix vectorization is incorrect!"
-    return v
+    n = S.shape[0]
+    S = np.copy(S)
+    S *= np.sqrt(2)
+    S[range(n), range(n)] /= np.sqrt(2)
+    return S[np.triu_indices(n)].reshape((-1,1))
 
 
-def mat(v: NDArray) -> NDArray:
+def mat(s: NDArray) -> NDArray:
     r"""
     Reconstruct a symmetric matrix from a vector that was created using
     :py:fn:`vec`. This scales the off-diagonal entries by :math:`1/\sqrt{2}`. 
     """
-    n_float = .5 * (-1 + np.sqrt(1 + 8 * v.shape[0]))
-    n = int(n_float)
-    if not math.isclose(n, n_float):
-        raise ValueError("To construct a n x n symmetric matrix, the vector "
-                         "must be an n * (n + 1) / 2 dimensional column.")
-
-    m = np.zeros((n, n))
-    m[np.tril_indices(n)] = v
-
-    scaling = (np.tril(np.ones(n), k=-1) / np.sqrt(2) + np.eye(n))
-    m = m * scaling
-    m = m + m.T - np.diag(np.diag(m))
-    return m
+    n = int((np.sqrt(8 * len(s) + 1) - 1) / 2)
+    S = np.zeros((n, n))
+    S[np.triu_indices(n)] = s / np.sqrt(2)
+    S = S + S.T
+    S[range(n), range(n)] /= np.sqrt(2)
+    return S
 
 
 def solve_cone(prob: ConicProblem, verbose: bool = False,
-- 
cgit v1.2.1