add putinar n-satz

13 files changed, 578 insertions, 33 deletions

diff --git a/sumofsquares/__init__.py b/sumofsquares/__init__.py
index b44da40..fb40a64 100644
--- a/sumofsquares/__init__.py
+++ b/sumofsquares/__init__.py
@@ -2,5 +2,8 @@ from sumofsquares.abc.sosexpr import ParamSOSExpr, SOSExpr
 
 from sumofsquares.init.initsosexpr import init_sos_expr, init_param_expr, init_param_expr_from_reference
 from sumofsquares.init.initenergyfunction import init_energy_function
+from sumofsquares.init.initenergyfunction2 import init_energy_function2
+from sumofsquares.init.initinteriorconstraint import init_interior_constraint
+from sumofsquares.init.initputinarepsilon import init_putinar_epsilon
 
-from sumofsquares.cvxopt import solve_cone, solve_sos_problem
+from sumofsquares.cvxopt import solve_cone, solve_cone2, solve_sos_problem, solve_sos_problem2
diff --git a/sumofsquares/cvxopt.py b/sumofsquares/cvxopt.py
index ed2146d..cc7d555 100644
--- a/sumofsquares/cvxopt.py
+++ b/sumofsquares/cvxopt.py
@@ -28,7 +28,17 @@ class CVXOptConeQPResult:
 
 
 
-def solve_cone(state, variables, cost, s_inequality, l_inequality = tuple(), previous_result=None, print_info=False, print_matrix=False):
+def solve_cone(
+        state, 
+        variables, 
+        cost, 
+        s_inequality,
+        # q_inequality = tuple(), 
+        l_inequality = tuple(), 
+        previous_result=None, 
+        print_info=False, 
+        print_matrix=False,
+    ):
     state, matrix_repr = polymatrix.to_matrix_repr(
         (cost,) + l_inequality + s_inequality,
         polymatrix.v_stack(variables),
@@ -135,6 +145,164 @@ def solve_cone(state, variables, cost, s_inequality, l_inequality = tuple(), pre
 
     return state, result
 
+def solve_cone2(
+        state, 
+        variables, 
+        cost, 
+        s_inequality = tuple(),
+        q_inequality = tuple(), 
+        l_inequality = tuple(), 
+        previous_result=None, 
+        print_info=False, 
+        print_matrix=False,
+    ):
+
+    try:    
+        state, matrix_repr = polymatrix.to_matrix_repr(
+            (cost,) + l_inequality + q_inequality + s_inequality,
+            polymatrix.v_stack(variables),
+        ).apply(state)
+    except:
+        state, matrix_repr = polymatrix.to_matrix_repr(
+            (cost,),
+            polymatrix.v_stack(variables),
+        ).apply(state)
+    
+        if len(l_inequality):
+            state, matrix_repr = polymatrix.to_matrix_repr(
+                l_inequality,
+                polymatrix.v_stack(variables),
+            ).apply(state)
+    
+        if len(q_inequality):
+            state, matrix_repr = polymatrix.to_matrix_repr(
+                q_inequality,
+                polymatrix.v_stack(variables),
+            ).apply(state)
+
+        if len(s_inequality):
+            state, matrix_repr = polymatrix.to_matrix_repr(
+                s_inequality,
+                polymatrix.v_stack(variables),
+            ).apply(state)
+    
+    
+    
+    # maximum degree of cost function must be 2
+    max_degree_cost = matrix_repr.data[0].get_max_degree()
+    assert max_degree_cost <= 2, f'{max_degree_cost=}'
+
+    # maximum degree of constraint must be 1
+    max_degree_constraints = max(data.get_max_degree() for data in matrix_repr.data[1:])
+    assert max_degree_constraints == 1, f'{max_degree_constraints=}'
+
+    cost = matrix_repr.data[0]
+
+    # cost[1] is a 1 x n vector
+    q = cost[1] / 2
+        
+    Gl = matrix_repr.data[1:1+len(l_inequality)]     
+    dim_l = sum(G[0].shape[0] for G in Gl)
+
+    Gq = matrix_repr.data[1+len(l_inequality):1+len(l_inequality)+len(q_inequality)]     
+    dim_q = list(G[0].shape[0] for G in Gq)
+
+    Gs = matrix_repr.data[1+len(l_inequality)+len(q_inequality):]  
+    
+    def gen_square_matrix_dim():
+        for Gs_raw in Gs:
+            dim = np.sqrt(Gs_raw[0].shape[0])
+
+            assert math.isclose(int(dim), dim), f'{dim=}'
+            
+            yield int(dim)
+    
+    dim_s = list(gen_square_matrix_dim())
+    
+    G = np.vstack(tuple(-v[1] for v in matrix_repr.data[1:]))
+    h = np.vstack(tuple(v[0] for v in matrix_repr.data[1:]))
+
+    if print_info:
+        print(f'number of variables: {G.shape[1]}')
+        print(f'{dim_l=}, {dim_q=}, {dim_s=}')
+    
+    if print_matrix:
+        print(f'{q=}')
+        print(f'{G=}')
+        print(f'{h=}')
+
+    if previous_result is None:
+        primalstart = None
+        
+    else:
+        primalstart = {
+            'x': cvxopt.matrix(previous_result.x),
+            'y': cvxopt.matrix(previous_result.y),
+            's': cvxopt.matrix(previous_result.s),
+            'z': cvxopt.matrix(previous_result.z),
+        }
+
+    if max_degree_cost == 1:
+        return_val = cvxopt.solvers.conelp(
+            c=cvxopt.matrix(q.reshape(-1, 1)),
+            G=cvxopt.matrix(G), 
+            h=cvxopt.matrix(h),
+            dims={'l': dim_l, 'q': dim_q, 's': dim_s},
+            primalstart=primalstart,
+        )
+        
+    else:
+        # cost[2] is a 1 x n*n vector
+        nP = int(np.sqrt(cost[2].shape[1]))
+        P = cost[2].reshape(nP, nP).toarray()
+    
+        if print_matrix:
+            print(f'{P=}')
+
+        if print_matrix:
+            print(f'{dim_l=}, {dim_q=}, {dim_s=}')
+            print(f'{P=}')
+            print(f'{q=}')
+            print(f'{G=}')
+            print(f'{h=}')
+        
+        return_val = cvxopt.solvers.coneqp(
+            P=cvxopt.matrix(P), 
+            q=cvxopt.matrix(q.reshape(-1, 1)),
+            G=cvxopt.matrix(G), 
+            h=cvxopt.matrix(h),
+            dims={'l': dim_l, 'q': dim_q, 's': dim_s,},
+            initvals=primalstart,
+        )
+        
+    x = np.array(return_val['x']).reshape(-1)
+
+    # if x[0] == None:
+    #     x0 = {}
+    # else:
+    #     x0 = {var.name: matrix_repr.get_value(var, x) for var in variables}
+    
+    # result = CVXOptConeQPResult(
+    #     x=x,
+    #     y=np.array(return_val['y']).reshape(-1),
+    #     s=np.array(return_val['s']).reshape(-1),
+    #     z=np.array(return_val['z']).reshape(-1),
+    #     status=return_val['status'],
+    #     gap=return_val['gap'],
+    #     relative_gap=return_val['relative gap'],
+    #     primal_objective=return_val['primal objective'],
+    #     dual_objective=return_val['dual objective'],
+    #     primal_infeasibility=return_val['primal infeasibility'],
+    #     dual_infeasibility=return_val['dual infeasibility'],
+    #     primal_slack=return_val['primal slack'],
+    #     dual_slack=return_val['dual slack'],
+    #     iterations=return_val['iterations'],
+    #     x0=x0
+    # )
+
+    # return state, result
+    return matrix_repr.get_value(variables[0], x), return_val['status']
+
     
 def solve_sos_problem(
     cost: tuple[polymatrix.Expression], 
@@ -175,42 +343,43 @@ def solve_sos_problem(
     return state, result
 
 
-# def solve_sos_problem(
-#     cost: tuple[polymatrix.Expression], 
-#     sos_constraints: tuple[SOSConstraint],
-#     state: polymatrix.ExpressionState,
-#     subs: tuple[ParamSOSExpr] = None, 
-#     x0: dict[ParamSOSExpr, np.ndarray] = None,
-# ):
-#     if x0 is None:
-#         x0 = {}
+def solve_sos_problem2(
+    cost: tuple[polymatrix.Expression], 
+    sos_constraints: tuple[SOSConstraint],
+    state: polymatrix.ExpressionState,
+    subs: tuple[ParamSOSExpr] = None, 
+    x0: dict[ParamSOSExpr, np.ndarray] = None,
+):
+    if x0 is None:
+        x0 = {}
     
-#     if subs is None:
-#         subs = tuple()
+    if subs is None:
+        subs = tuple()
         
-#     def gen_free_param_expr():
-#         for sos_constraint in sos_constraints:
-#             for param_expr in sos_constraint.dependence:
-#                 if param_expr not in subs:
-#                     yield param_expr
+    def gen_free_param_expr():
+        for sos_constraint in sos_constraints:
+            for param_expr in sos_constraint.dependence:
+                if param_expr not in subs:
+                    yield param_expr
                     
-#     free_param_expr = tuple(set(gen_free_param_expr()))
+    free_param_expr = tuple(set(gen_free_param_expr()))
+
+    # print(f'{tuple(val.param.name for val in free_param_expr)=}')
     
-#     sub_vals = tuple((sub.param, x0[sub.param]) for sub in subs)
+    sub_vals = tuple((sub.param, x0[sub.param.name]) for sub in subs if sub is not None and sub.param.name in x0)
     
-#     def gen_inequality():
-#         for sos_constraint in sos_constraints:
-#             if any(param_expr in free_param_expr for param_expr in sos_constraint.dependence):
-#                 # print(tuple(d.name for d in sos_constraint.dependence))
-#                 yield sos_constraint.constraint.eval(sub_vals)
+    def gen_inequality():
+        for sos_constraint in sos_constraints:
+            if any(param_expr in free_param_expr for param_expr in sos_constraint.dependence):
+                yield sos_constraint.constraint.eval(sub_vals)
     
-#     inequality = tuple(gen_inequality())
+    inequality = tuple(gen_inequality())
     
-#     state, result = solve_cone(
-#         state,
-#         variables=tuple(param_expr.param for param_expr in free_param_expr),
-#         cost=sum(cost),
-#         s_inequality=inequality,
-#     )
+    state, result = solve_cone(
+        state,
+        variables=tuple(param_expr.param for param_expr in free_param_expr),
+        cost=sum(cost),
+        s_inequality=inequality,
+    )
         
-#     return state, result
+    return state, result
diff --git a/sumofsquares/impl/energyfunction2impl.py b/sumofsquares/impl/energyfunction2impl.py
new file mode 100644
index 0000000..7e6b091
--- /dev/null
+++ b/sumofsquares/impl/energyfunction2impl.py
@@ -0,0 +1,19 @@
+import dataclassabc
+import polymatrix
+
+from sumofsquares.abc.sosexpr import ParamSOSExpr, SOSExpr
+from sumofsquares.mixins.energyfunction2mixin import EnergyFunction2Mixin
+from sumofsquares.sosconstraint.abc import SOSConstraint
+
+
+@dataclassabc.dataclassabc(frozen=True)
+class EnergyFunction2Impl(EnergyFunction2Mixin):
+    name: str
+    h: SOSExpr
+    h_dot: SOSExpr
+    epsilon: ParamSOSExpr
+    gamma_b: ParamSOSExpr
+    gamma_roi: list[ParamSOSExpr]
+    # gamma_w_inner: ParamSOSExpr | None
+    # gamma_w_outer: ParamSOSExpr | None
+    sos_constraints: tuple[SOSConstraint]
diff --git a/sumofsquares/impl/interiorconstraintimpl.py b/sumofsquares/impl/interiorconstraintimpl.py
new file mode 100644
index 0000000..ac76c19
--- /dev/null
+++ b/sumofsquares/impl/interiorconstraintimpl.py
@@ -0,0 +1,14 @@
+import dataclassabc
+import polymatrix
+
+from sumofsquares.abc.sosexpr import ParamSOSExpr, SOSExpr
+from sumofsquares.mixins.energyfunction2mixin import EnergyFunction2Mixin
+from sumofsquares.mixins.interiorconstraintmixin import InteriorConstraintMixin
+from sumofsquares.sosconstraint.abc import SOSConstraint
+
+
+@dataclassabc.dataclassabc(frozen=True)
+class InteriorConstraintImpl(InteriorConstraintMixin):
+    name: str
+    gamma: ParamSOSExpr | None
+    sos_constraints: tuple[SOSConstraint]
diff --git a/sumofsquares/impl/putinarepsilonimpl.py b/sumofsquares/impl/putinarepsilonimpl.py
new file mode 100644
index 0000000..dfba3f0
--- /dev/null
+++ b/sumofsquares/impl/putinarepsilonimpl.py
@@ -0,0 +1,15 @@
+import dataclassabc
+import polymatrix
+
+from sumofsquares.abc.sosexpr import ParamSOSExpr, SOSExpr
+from sumofsquares.mixins.putinarepsilonmixin import PutinarEpsilonMixin
+from sumofsquares.sosconstraint.abc import SOSConstraint
+
+
+@dataclassabc.dataclassabc(frozen=True)
+class PutinarEpsilonImpl(PutinarEpsilonMixin):
+    name: str
+    epsilon: ParamSOSExpr
+    gamma: dict[str, ParamSOSExpr]
+    sos_constraints: tuple[SOSConstraint]
+    condition: SOSExpr
diff --git a/sumofsquares/init/initenergyfunction2.py b/sumofsquares/init/initenergyfunction2.py
new file mode 100644
index 0000000..785ef87
--- /dev/null
+++ b/sumofsquares/init/initenergyfunction2.py
@@ -0,0 +1,110 @@
+import polymatrix
+
+from sumofsquares.abc.sosexpr import SOSExpr
+from sumofsquares.impl.energyfunction2impl import EnergyFunction2Impl
+from sumofsquares.init.initsosexpr import init_param_expr, init_param_expr_from_reference
+
+
+def init_energy_function2(
+    name: str,
+    # h: SOSExpr,
+    # x_dot: SOSExpr,
+    cond: SOSExpr,
+    h: SOSExpr,
+    positive_b: bool,
+    epsilon_min: SOSExpr | None = None,
+    roi: list[SOSExpr | polymatrix.Expression] = None,
+    decrease_rate: SOSExpr | polymatrix.Expression | None = None,
+    # w_outer: SOSExpr | polymatrix.Expression | None = None,
+    # w_inner: SOSExpr | polymatrix.Expression | None = None,
+):
+    sos_constraints = tuple()
+
+    cond = cond.cache()
+    
+    # nh = h.diff().T
+    # h_dot = (nh.T @ x_dot).cache()
+
+    gamma_b = init_param_expr_from_reference(
+        name=f'gamma_b_{name}',
+        reference=cond, 
+        multiplicand=h,
+    )
+    if positive_b:
+        sos_constraints += gamma_b.sos_constraints
+
+    energy_expr = -cond
+    energy_expr -= gamma_b * h
+
+    if epsilon_min is None:
+        epsilon = None
+
+    else:
+        epsilon = init_param_expr(
+            name=f'epsilon_{name}',
+            variables=h.variables,
+        )
+
+        sos_constraints += (epsilon - epsilon_min).sos_constraints
+        energy_expr += epsilon
+
+    if roi is None:
+        gamma_roi = None,
+
+    else:
+        def gen_gamma_roi():
+            for idx, roi_i in enumerate(roi):
+                yield init_param_expr_from_reference(
+                    name=f'gamma_roi{idx}_{name}',
+                    reference=cond, 
+                    multiplicand=roi_i,
+                )
+        gamma_roi = tuple(gen_gamma_roi())
+
+        for gamma_roi_i, roi_i in zip(gamma_roi, roi):
+            sos_constraints += gamma_roi_i.sos_constraints
+
+            energy_expr += gamma_roi_i * roi_i
+
+    if decrease_rate is not None:    
+        energy_expr -= decrease_rate
+
+    sos_constraints += energy_expr.sos_constraints
+
+    # if w_outer is None:
+    #     gamma_w_outer = None
+
+    # else:
+    #     gamma_w_outer = init_param_expr_from_reference(
+    #         name=f'gamma_w_outer_{name}',
+    #         reference=h, 
+    #         multiplicand=w_outer,
+    #     )
+
+    #     sos_constraints += gamma_w_outer.sos_constraints
+    #     sos_constraints += (h - gamma_w_outer * w_outer).sos_constraints
+
+    # if w_inner is None:
+    #     gamma_w_inner = None
+
+    # else:
+    #     gamma_w_inner = init_param_expr_from_reference(
+    #         name=f'gamma_w_inner_{name}',
+    #         reference=h, 
+    #         multiplicand=w_inner,
+    #     )
+
+    #     sos_constraints += gamma_w_inner.sos_constraints
+    #     sos_constraints += (gamma_w_inner * w_inner - h).sos_constraints
+
+    return EnergyFunction2Impl(
+        name=name,
+        h=h,
+        h_dot=cond,
+        epsilon=epsilon,
+        gamma_roi=gamma_roi,
+        gamma_b=gamma_b,
+        # gamma_w_outer=gamma_w_outer,
+        # gamma_w_inner=gamma_w_inner,
+        sos_constraints=sos_constraints,
+    )
diff --git a/sumofsquares/init/initinteriorconstraint.py b/sumofsquares/init/initinteriorconstraint.py
new file mode 100644
index 0000000..3daa44e
--- /dev/null
+++ b/sumofsquares/init/initinteriorconstraint.py
@@ -0,0 +1,28 @@
+import polymatrix
+
+from sumofsquares.abc.sosexpr import SOSExpr
+from sumofsquares.impl.interiorconstraintimpl import InteriorConstraintImpl
+from sumofsquares.init.initsosexpr import init_param_expr_from_reference
+
+
+def init_interior_constraint(
+    name: str,
+    h: SOSExpr,
+    w: SOSExpr | polymatrix.Expression,
+):
+    sos_constraints = tuple()
+
+    gamma = init_param_expr_from_reference(
+        name=f'gamma_{name}',
+        reference=h, 
+        multiplicand=w,
+    )
+
+    sos_constraints += gamma.sos_constraints
+    sos_constraints += (h - gamma * w).sos_constraints
+
+    return InteriorConstraintImpl(
+        name=name,
+        gamma=gamma,
+        sos_constraints=sos_constraints,
+    )
diff --git a/sumofsquares/init/initputinarepsilon.py b/sumofsquares/init/initputinarepsilon.py
new file mode 100644
index 0000000..d3e348d
--- /dev/null
+++ b/sumofsquares/init/initputinarepsilon.py
@@ -0,0 +1,72 @@
+import polymatrix
+
+from sumofsquares.abc.sosexpr import SOSExpr
+from sumofsquares.impl.energyfunction2impl import EnergyFunction2Impl
+from sumofsquares.impl.putinarepsilonimpl import PutinarEpsilonImpl
+from sumofsquares.init.initsosexpr import init_param_expr, init_param_expr_from_reference
+
+
+def init_putinar_epsilon(
+    name: str,
+    f0: SOSExpr,
+    fi: dict[str, SOSExpr | polymatrix.Expression] = None,
+    gi: dict[str, SOSExpr | polymatrix.Expression] = None,
+    epsilon_min: SOSExpr | None = None,
+    decrease_rate: SOSExpr | polymatrix.Expression | None = None,
+):
+    sos_constraints = tuple()
+
+    condition = f0
+
+    def gen_gamma(fi):
+        for key, val in fi.items():
+            yield key, init_param_expr_from_reference(
+                name=f'gamma_{key}_{name}',
+                reference=f0,
+                # reference=condition,
+                multiplicand=val,
+            )
+
+    if fi is None:
+        gamma_fi = {}
+
+    else:
+        gamma_fi = dict(gen_gamma(fi))
+
+        for key, gamma_i in gamma_fi.items():
+            sos_constraints += gamma_i.sos_constraints
+            condition = condition - gamma_i * fi[key]
+
+    if gi is None:
+        gamma_gi = {}
+
+    else:
+        gamma_gi = dict(gen_gamma(gi))
+
+        for key, gamma_i in gamma_gi.items():
+            condition = condition - gamma_i * gi[key]
+
+    if epsilon_min is None:
+        epsilon = None
+
+    else:
+        epsilon = init_param_expr(
+            name=f'epsilon_{name}',
+            variables=f0.variables,
+        )
+
+        sos_constraints += (epsilon - epsilon_min).sos_constraints
+        condition += epsilon
+
+    if decrease_rate is not None:    
+        condition -= decrease_rate
+
+    sos_constraints += condition.sos_constraints
+
+    return PutinarEpsilonImpl(
+        name=name,
+        epsilon=epsilon,
+        gamma=gamma_fi | gamma_gi,
+        sos_constraints=sos_constraints,
+        condition=condition,
+    )
diff --git a/sumofsquares/mixins/energyfunction2mixin.py b/sumofsquares/mixins/energyfunction2mixin.py
new file mode 100644
index 0000000..f96abb6
--- /dev/null
+++ b/sumofsquares/mixins/energyfunction2mixin.py
@@ -0,0 +1,53 @@
+import abc
+import polymatrix
+
+from sumofsquares.abc.sosexpr import ParamSOSExpr, SOSExpr
+from sumofsquares.mixins.getsosconstraintmixin import GetSOSConstraintMixin
+from sumofsquares.sosconstraint.abc import SOSConstraint
+
+
+class EnergyFunction2Mixin(GetSOSConstraintMixin, abc.ABC):
+    @property
+    @abc.abstractmethod
+    def name(self) -> str:
+        ...
+
+    @property
+    @abc.abstractmethod
+    def h(self) -> SOSExpr:
+        ...
+
+    @property
+    @abc.abstractmethod
+    def h_dot(self) -> SOSExpr:
+        ...
+
+    # @property
+    # @abc.abstractmethod
+    # def variables(self) -> polymatrix.Expression:
+    #     ...
+
+    @property
+    @abc.abstractmethod
+    def epsilon(self) -> ParamSOSExpr:
+        ...
+
+    @property
+    @abc.abstractmethod
+    def gamma_b(self) -> ParamSOSExpr:
+        ...
+
+    @property
+    @abc.abstractmethod
+    def gamma_roi(self) -> list[ParamSOSExpr]:
+        ...
+
+    # @property
+    # @abc.abstractmethod
+    # def gamma_w_outer(self) -> ParamSOSExpr | None:
+    #     ...
+
+    # @property
+    # @abc.abstractmethod
+    # def gamma_w_inner(self) -> ParamSOSExpr | None:
+    #     ...
diff --git a/sumofsquares/mixins/interiorconstraintmixin.py b/sumofsquares/mixins/interiorconstraintmixin.py
new file mode 100644
index 0000000..781a3d4
--- /dev/null
+++ b/sumofsquares/mixins/interiorconstraintmixin.py
@@ -0,0 +1,21 @@
+import abc
+
+from sumofsquares.abc.sosexpr import ParamSOSExpr
+from sumofsquares.mixins.getsosconstraintmixin import GetSOSConstraintMixin
+
+
+class InteriorConstraintMixin(GetSOSConstraintMixin, abc.ABC):
+    @property
+    @abc.abstractmethod
+    def name(self) -> str:
+        ...
+
+    # @property
+    # @abc.abstractmethod
+    # def h(self) -> SOSExpr:
+    #     ...
+
+    @property
+    @abc.abstractmethod
+    def gamma(self) -> ParamSOSExpr | None:
+        ...
diff --git a/sumofsquares/mixins/putinarepsilonmixin.py b/sumofsquares/mixins/putinarepsilonmixin.py
new file mode 100644
index 0000000..6c1cd53
--- /dev/null
+++ b/sumofsquares/mixins/putinarepsilonmixin.py
@@ -0,0 +1,28 @@
+import abc
+import polymatrix
+
+from sumofsquares.abc.sosexpr import ParamSOSExpr, SOSExpr
+from sumofsquares.mixins.getsosconstraintmixin import GetSOSConstraintMixin
+from sumofsquares.sosconstraint.abc import SOSConstraint
+
+
+class PutinarEpsilonMixin(GetSOSConstraintMixin, abc.ABC):
+    @property
+    @abc.abstractmethod
+    def name(self) -> str:
+        ...
+
+    @property
+    @abc.abstractmethod
+    def epsilon(self) -> ParamSOSExpr:
+        ...
+
+    @property
+    @abc.abstractmethod
+    def gamma(self) -> dict[str, ParamSOSExpr]:
+        ...
+
+    @property
+    @abc.abstractmethod
+    def condition(self) -> SOSExpr:
+        ...
diff --git a/sumofsquares/mixins/sosexprmixin.py b/sumofsquares/mixins/sosexprmixin.py
index d0c0d17..f3d4dc9 100644
--- a/sumofsquares/mixins/sosexprmixin.py
+++ b/sumofsquares/mixins/sosexprmixin.py
@@ -30,6 +30,10 @@ class SOSExprMixin(GetSOSConstraintMixin, abc.ABC):
         return self.underlying.sos_matrix
 
     @property
+    def sos_matrix_vec(self) -> polymatrix.Expression:
+        return self.underlying.sos_matrix_as_vector
+
+    @property
     def sos_constraints(self) -> tuple[SOSConstraint]:
         return (init_sos_constraint(
             dependence=self.dependence,
diff --git a/sumofsquares/mixins/sosexpropmixin.py b/sumofsquares/mixins/sosexpropmixin.py
index 2020043..f7cf0c2 100644
--- a/sumofsquares/mixins/sosexpropmixin.py
+++ b/sumofsquares/mixins/sosexpropmixin.py
@@ -47,6 +47,10 @@ class SOSExprOPMixin(SOSExprMixin):
             )
 
         else:
+            # var_set_left = set(left.variables)
+            # var_set_right = set(right.variables)
+            # assert var_set_left.issubset(var_set_right) or var_set_right.issubset(var_set_left), f'{left.variables=}, {right.variables=}'
+
             assert left.variables == right.variables, f'{left.variables=}, {right.variables=}'
 
             underlying=init_sos_expr_base(
@@ -142,3 +146,8 @@ class SOSExprOPMixin(SOSExprMixin):
             ),
         )
 
+    def v_stack(self, other):
+        def op(left, right):
+            return polymatrix.v_stack((left, right))
+
+        return self._binary(op, self, other)