solver_test_helpers.py

"""
Copyright 2019, the CVXPY developers.

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

    http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
"""
import warnings

import numpy as np

import cvxpy as cp
import cvxpy.settings as s
from cvxpy.tests.base_test import BaseTest


class SolverTestHelper:

    def __init__(self, obj_pair, var_pairs, con_pairs) -> None:
        self.objective = obj_pair[0]
        self.constraints = [c for c, _ in con_pairs]
        self.prob = cp.Problem(self.objective, self.constraints)
        self.variables = [x for x, _ in var_pairs]

        self.expect_val = obj_pair[1]
        self.expect_dual_vars = [dv for _, dv in con_pairs]
        self.expect_prim_vars = [pv for _, pv in var_pairs]
        self.tester = BaseTest()

    def solve(self, solver, **kwargs) -> None:
        self.prob.solve(solver=solver, **kwargs)

    def check_primal_feasibility(self, places) -> None:
        all_cons = [c for c in self.constraints]  # shallow copy
        for x in self.prob.variables():
            attrs = x.attributes
            if attrs['nonneg'] or attrs['pos']:
                all_cons.append(x >= 0)
            elif attrs['nonpos'] or attrs['neg']:
                all_cons.append(x <= 0)
            elif attrs['imag']:
                all_cons.append(x + cp.conj(x) == 0)
            elif attrs['symmetric']:
                all_cons.append(x == x.T)
            elif attrs['diag']:
                all_cons.append(x - cp.diag(cp.diag(x)) == 0)
            elif attrs['PSD']:
                all_cons.append(x >> 0)
            elif attrs['NSD']:
                all_cons.append(x << 0)
            elif attrs['hermitian']:
                all_cons.append(x == cp.conj(x.T))
            elif attrs['boolean'] or attrs['integer']:
                round_val = np.round(x.value)
                all_cons.append(x == round_val)
        for con in all_cons:
            viol = con.violation()
            if isinstance(viol, np.ndarray):
                viol = np.linalg.norm(viol, ord=2)
            self.tester.assertAlmostEqual(viol, 0, places)

    def check_dual_domains(self, places) -> None:
        # A full "dual feasibility" check would involve checking a stationary Lagrangian.
        # No such test is planned here.
        #
        # TODO: once dual variables are stored for attributes
        #   (e.g. X = Variable(shape=(n,n), PSD=True)), check
        #   domains for dual variables of the attribute constraint.
        for con in self.constraints:
            if isinstance(con, cp.constraints.Cone):
                dual_violation = con.dual_residual
                if isinstance(con, cp.constraints.SOC):
                    dual_violation = np.linalg.norm(dual_violation)
                self.tester.assertLessEqual(dual_violation, 10**(-places))
            elif isinstance(con, cp.constraints.Inequality):
                # TODO: move this to Inequality.dual_violation
                dv = con.dual_value
                min_dv = np.min(dv)
                self.tester.assertGreaterEqual(min_dv, -(10**(-places)))
            elif isinstance(con, (cp.constraints.Equality, cp.constraints.Zero)):
                dv = con.dual_value
                self.tester.assertIsNotNone(dv)
                if isinstance(dv, np.ndarray):
                    contents = dv.dtype
                    self.tester.assertEqual(contents, float)
                else:
                    self.tester.assertIsInstance(dv, float)
            else:
                raise ValueError('Unknown constraint type %s.' % type(con))

    def check_complementarity(self, places) -> None:
        # TODO: once dual variables are stored for attributes
        #   (e.g. X = Variable(shape=(n,n), PSD=True)), check
        #   complementarity against the dual variable of the
        #   attribute constraint.
        for con in self.constraints:
            if isinstance(con, (cp.constraints.Inequality,
                                cp.constraints.Equality)):
                comp = cp.vdot(con.expr, con.dual_value).value
            elif isinstance(con, (cp.constraints.ExpCone,
                                  cp.constraints.SOC,
                                  cp.constraints.NonNeg,
                                  cp.constraints.Zero,
                                  cp.constraints.PSD,
                                  cp.constraints.PowCone3D,
                                  cp.constraints.PowConeND)):
                comp = cp.vdot(con.args, con.dual_value).value
            elif isinstance(con, cp.RelEntrConeQuad) or isinstance(con, cp.OpRelEntrConeQuad):
                msg = '\nDual variables not implemented for quadrature based approximations;' \
                       + '\nSkipping complementarity check.'
                warnings.warn(msg)
            else:
                raise ValueError('Unknown constraint type %s.' % type(con))
            self.tester.assertAlmostEqual(comp, 0, places)

    def check_stationary_lagrangian(self, places) -> None:
        L = self.prob.objective.expr
        objective = self.prob.objective
        if objective.NAME == 'minimize':
            L = objective.expr
        else:
            L = -objective.expr
        for con in self.constraints:
            if isinstance(con, (cp.constraints.Inequality,
                                cp.constraints.Equality)):
                dual_var_value = con.dual_value
                prim_var_expr = con.expr
                L = L + cp.vdot(dual_var_value, prim_var_expr)
            elif isinstance(con, (cp.constraints.ExpCone,
                                  cp.constraints.SOC,
                                  cp.constraints.Zero,
                                  cp.constraints.NonNeg,
                                  cp.constraints.PSD,
                                  cp.constraints.PowCone3D,
                                  cp.constraints.PowConeND)):
                L = L - cp.vdot(con.args, con.dual_value)
            else:
                raise NotImplementedError()
        try:
            g = L.grad
        except TypeError as e:
            assert 'is not subscriptable' in str(e)
            msg = """\n
            CVXPY problems with `diag` variables are not supported for
            stationarity checks as of now
            """
            self.tester.fail(msg)
        bad_norms = []

        r"""The convention that we follow for construting the Lagrangian is: 1) Move all
        explicitly passed constraints to the problem (via Problem.constraints) into the
        Lagrangian --- dLdX == 0 for any such variables 2) Constraints that have
        implicitly been imposed on variables at the time of declaration via specific
        flags (e.g.: PSD/symmetric etc.), in such a case we check, `dLdX\in K^{*}`, where
        `K` is the convex cone corresponding to the implicit constraint on `X`
        """
        for (opt_var, v) in g.items():
            if all(not attr for attr in list(map(lambda x: x[1], opt_var.attributes.items()))):
                """Case when the variable doesn't have any special attributes"""
                norm = np.linalg.norm(v.data) / np.sqrt(opt_var.size)
                if norm > 10**(-places):
                    bad_norms.append((norm, opt_var))
            else:
                if opt_var.is_psd():
                    """The PSD cone is self-dual"""
                    g_bad_mat = cp.Constant(np.reshape(g[opt_var].toarray(), opt_var.shape))
                    tmp_con = g_bad_mat >> 0
                    dual_cone_violation = tmp_con.residual
                    if dual_cone_violation > 10**(-places):
                        bad_norms.append((dual_cone_violation, opt_var))
                elif opt_var.is_nsd():
                    """The NSD cone is also self-dual"""
                    g_bad_mat = cp.Constant(np.reshape(g[opt_var].toarray(), opt_var.shape))
                    tmp_con = g_bad_mat << 0
                    dual_cone_violation = tmp_con.residual
                    if dual_cone_violation > 10**(-places):
                        bad_norms.append((dual_cone_violation, opt_var))
                elif opt_var.is_diag():
                    """The dual cone to the set of diagonal matrices is the set of
                        'Hollow' matrices i.e. matrices with diagonal entries zero"""
                    g_bad_mat = np.reshape(g[opt_var].toarray(), opt_var.shape)
                    diag_entries = np.diag(opt_var.value)
                    dual_cone_violation = np.linalg.norm(diag_entries) / np.sqrt(opt_var.size)
                    if diag_entries > 10**(-places):
                        bad_norms.append((dual_cone_violation, opt_var))
                elif opt_var.is_symmetric():
                    r"""The dual cone to the set of symmetric matrices is the
                    set of skew-symmetric matrices, so we check if dLdX \in
                    set(skew-symmetric-matrices)
                    g[opt_var] is the problematic gradient in question"""
                    g_bad_mat = np.reshape(g[opt_var].toarray(), opt_var.shape)
                    mat = g_bad_mat + g_bad_mat.T
                    dual_cone_violation = np.linalg.norm(mat) / np.sqrt(opt_var.size)
                    if dual_cone_violation > 10**(-places):
                        bad_norms.append((dual_cone_violation, opt_var))
                elif opt_var.is_nonpos():
                    """The cone of matrices with all entries nonpos is self-dual"""
                    g_bad_mat = cp.Constant(np.reshape(g[opt_var].toarray(), opt_var.shape))
                    tmp_con = g_bad_mat <= 0
                    dual_cone_violation = np.linalg.norm(tmp_con.residual) / np.sqrt(opt_var.size)
                    if dual_cone_violation > 10**(-places):
                        bad_norms.append((dual_cone_violation, opt_var))
                elif opt_var.is_nonneg():
                    """The cone of matrices with all entries nonneg is self-dual"""
                    g_bad_mat = cp.Constant(np.reshape(g[opt_var].toarray(), opt_var.shape))
                    tmp_con = g_bad_mat >= 0
                    dual_cone_violation = np.linalg.norm(tmp_con.residual) / np.sqrt(opt_var.size)
                    if dual_cone_violation > 10**(-places):
                        bad_norms.append((dual_cone_violation, opt_var))

        if len(bad_norms):
            msg = f"""\n
        The gradient of Lagrangian with respect to the primal variables
        is above the threshold of 10^{-places}. The names of the problematic
        variables and the corresponding gradient norms are as follows:
            """
            for norm, opt_var in bad_norms:
                msg += f"\n\t\t\t{opt_var.name} : {norm}"
            msg += '\n'
            self.tester.fail(msg)
        pass

    def verify_objective(self, places) -> None:
        actual = self.prob.value
        expect = self.expect_val
        if expect is not None:
            self.tester.assertAlmostEqual(actual, expect, places)

    def verify_primal_values(self, places) -> None:
        for idx in range(len(self.variables)):
            actual = self.variables[idx].value
            expect = self.expect_prim_vars[idx]
            if expect is not None:
                self.tester.assertItemsAlmostEqual(actual, expect, places)

    def verify_dual_values(self, places) -> None:
        for idx in range(len(self.constraints)):
            actual = self.constraints[idx].dual_value
            expect = self.expect_dual_vars[idx]
            if expect is not None:
                if isinstance(actual, list):
                    for i in range(len(actual)):
                        act = actual[i]
                        exp = expect[i]
                        self.tester.assertItemsAlmostEqual(act, exp, places)
                else:
                    self.tester.assertItemsAlmostEqual(actual, expect, places)


def lp_0() -> SolverTestHelper:
    x = cp.Variable(shape=(2,))
    con_pairs = [(x == 0, None)]
    obj_pair = (cp.Minimize(cp.norm(x, 1) + 1.0), 1)
    var_pairs = [(x, np.array([0, 0]))]
    sth = SolverTestHelper(obj_pair, var_pairs, con_pairs)
    return sth


def lp_1() -> SolverTestHelper:
    # Example from
    # http://cvxopt.org/userguide/coneprog.html?highlight=solvers.lp#cvxopt.solvers.lp
    x = cp.Variable(shape=(2,), name='x')
    objective = cp.Minimize(-4 * x[0] - 5 * x[1])
    constraints = [2 * x[0] + x[1] <= 3,
                   x[0] + 2 * x[1] <= 3,
                   x[0] >= 0,
                   x[1] >= 0]
    con_pairs = [(constraints[0], 1),
                 (constraints[1], 2),
                 (constraints[2], 0),
                 (constraints[3], 0)]
    var_pairs = [(x, np.array([1, 1]))]
    obj_pair = (objective, -9)
    sth = SolverTestHelper(obj_pair, var_pairs, con_pairs)
    return sth


def lp_2() -> SolverTestHelper:
    x = cp.Variable(shape=(2,), name='x')
    objective = cp.Minimize(x[0] + 0.5 * x[1])
    constraints = [x[0] >= -100, x[0] <= -10, x[1] == 1]
    con_pairs = [(constraints[0], 1),
                 (constraints[1], 0),
                 (constraints[2], -0.5)]
    var_pairs = [(x, np.array([-100, 1]))]
    obj_pair = (objective, -99.5)
    sth = SolverTestHelper(obj_pair, var_pairs, con_pairs)
    return sth


def lp_3() -> SolverTestHelper:
    # an unbounded problem
    x = cp.Variable(5)
    objective = (cp.Minimize(cp.sum(x)), -np.inf)
    var_pairs = [(x, None)]
    con_pairs = [(x <= 1, None)]
    sth = SolverTestHelper(objective, var_pairs, con_pairs)
    return sth


def lp_4() -> SolverTestHelper:
    # an infeasible problem
    x = cp.Variable(5)
    objective = (cp.Minimize(cp.sum(x)), np.inf)
    var_pairs = [(x, None)]
    con_pairs = [(x <= 0, None),
                 (x >= 1, None)]
    sth = SolverTestHelper(objective, var_pairs, con_pairs)
    return sth


def lp_5() -> SolverTestHelper:
    # a problem with redundant equality constraints.
    #
    # 10 variables, 6 equality constraints A @ x == b (two redundant)
    x0 = np.array([0, 1, 0, 2, 0, 4, 0, 5, 6, 7])
    mu0 = np.array([-2, -1, 0, 1, 2, 3.5])
    np.random.seed(0)
    A_min = np.random.randn(4, 10)
    A_red = A_min.T @ np.random.rand(4, 2)
    A_red = A_red.T
    A = np.vstack((A_min, A_red))
    b = A @ x0  # x0 is primal feasible
    c = A.T @ mu0  # mu0 is dual-feasible
    c[[0, 2, 4, 6]] += np.random.rand(4)
    # ^  c >= A.T @ mu0 exhibits complementary slackness with respect to x0
    #    Therefore (x0, mu0) are primal-dual optimal for ...
    x = cp.Variable(10)
    objective = (cp.Minimize(c @ x), c @ x0)
    var_pairs = [(x, x0)]
    con_pairs = [(x >= 0, None), (A @ x == b, None)]
    sth = SolverTestHelper(objective, var_pairs, con_pairs)
    return sth


def lp_6() -> SolverTestHelper:
    """Test LP with no constraints"""
    x = cp.Variable()
    from cvxpy.expressions.constants import Constant
    objective = cp.Maximize(Constant(0.23) * x)
    obj_pair = (objective, np.inf)
    var_pairs = [(x, None)]
    sth = SolverTestHelper(obj_pair, var_pairs, [])
    return sth


def lp_7() -> SolverTestHelper:
    """
    An ill-posed problem to test multiprecision ability of solvers.

    This test will not pass on CVXOPT (as of v1.3.1) and on SDPA without GMP support.
    """
    n = 50
    a = cp.Variable((n+1))
    delta = cp.Variable((n))
    b = cp.Variable((n+1))
    objective = cp.Minimize(cp.sum(cp.pos(delta)))
    constraints = [
        a[1:] - a[:-1] == delta,
        a >= cp.pos(b),
    ]
    con_pairs = [(constraints[0], None),
                 (constraints[1], None)]
    var_pairs = [(a, None),
                 (delta, None),
                 (b, None)]
    obj_pair = (objective, 0.)
    sth = SolverTestHelper(obj_pair, var_pairs, con_pairs)
    return sth


def lp_bound_attr() -> SolverTestHelper:
    """An LP using the variable bounds attribute."""
    lower = np.array([-100, 1])
    upper = np.array([-10, 1])
    x = cp.Variable(shape=(2,), name='x', bounds=[lower, upper])
    objective = cp.Minimize(x[0] + 0.5 * x[1])
    var_pairs = [(x, np.array([-100, 1]))]
    obj_pair = (objective, -99.5)
    sth = SolverTestHelper(obj_pair, var_pairs, [])
    return sth


def qp_bound_attr() -> SolverTestHelper:
    """A QP using the variable bounds attribute."""
    x = cp.Variable(shape=(2,), name='x', bounds=[np.array([1, -5]), np.array([10, 5])])
    objective = cp.Minimize(cp.sum_squares(x))
    var_pairs = [(x, np.array([1, 0]))]
    obj_pair = (objective, 1)
    sth = SolverTestHelper(obj_pair, var_pairs, [])
    return sth


def qp_0() -> SolverTestHelper:
    # univariate feasible problem
    x = cp.Variable(1)
    objective = cp.Minimize(cp.square(x))
    constraints = [x[0] >= 1]
    con_pairs = [(constraints[0], 2)]
    obj_pair = (objective, 1)
    var_pairs = [(x, 1)]
    sth = SolverTestHelper(obj_pair, var_pairs, con_pairs)
    return sth


def socp_0() -> SolverTestHelper:
    x = cp.Variable(shape=(2,))
    obj_pair = (cp.Minimize(cp.norm(x, 2) + 1), 1)
    con_pairs = [(x == 0, None)]
    var_pairs = [(x, np.array([0, 0]))]
    sth = SolverTestHelper(obj_pair, var_pairs, con_pairs)
    return sth


def socp_1() -> SolverTestHelper:
    """
    min 3 * x[0] + 2 * x[1] + x[2]
    s.t. norm(x,2) <= y
         x[0] + x[1] + 3*x[2] >= 1.0
         y <= 5
    """
    x = cp.Variable(shape=(3,))
    y = cp.Variable()
    soc = cp.constraints.second_order.SOC(y, x)
    constraints = [soc,
                   x[0] + x[1] + 3 * x[2] >= 1.0,
                   y <= 5]
    obj = cp.Minimize(3 * x[0] + 2 * x[1] + x[2])
    expect_x = np.array([-3.874621860638774, -2.129788233677883, 2.33480343377204])
    expect_x = np.round(expect_x, decimals=5)
    expect_y = 5
    var_pairs = [(x, expect_x),
                 (y, expect_y)]
    expect_soc = [np.array([2.86560262]), np.array([2.22062583,  1.22062583, -1.33812252])]
    expect_ineq1 = 0.7793969212001993
    expect_ineq2 = 2.865602615049077
    con_pairs = [(constraints[0], expect_soc),
                 (constraints[1], expect_ineq1),
                 (constraints[2], expect_ineq2)]
    obj_pair = (obj, -13.548638904065102)
    sth = SolverTestHelper(obj_pair, var_pairs, con_pairs)
    return sth


def socp_2() -> SolverTestHelper:
    """
    An (unnecessarily) SOCP-based reformulation of LP_1.
    """
    x = cp.Variable(shape=(2,), name='x')
    objective = cp.Minimize(-4 * x[0] - 5 * x[1])
    expr = cp.reshape(x[0] + 2 * x[1], (1, 1), order='F')
    constraints = [2 * x[0] + x[1] <= 3,
                   cp.constraints.SOC(cp.Constant([3]), expr),
                   x[0] >= 0,
                   x[1] >= 0]
    con_pairs = [(constraints[0], 1),
                 (constraints[1], [np.array([2.]), np.array([[-2.]])]),
                 (constraints[2], 0),
                 (constraints[3], 0)]
    var_pairs = [(x, np.array([1, 1]))]
    obj_pair = (objective, -9)
    sth = SolverTestHelper(obj_pair, var_pairs, con_pairs)
    return sth


def socp_3(axis) -> SolverTestHelper:
    x = cp.Variable(shape=(2,))
    c = np.array([-1, 2])
    root2 = np.sqrt(2)
    u = np.array([[1 / root2, -1 / root2], [1 / root2, 1 / root2]])
    mat1 = np.diag([root2, 1 / root2]) @ u.T
    mat2 = np.diag([1, 1])
    mat3 = np.diag([0.2, 1.8])

    X = cp.vstack([mat1 @ x, mat2 @ x, mat3 @ x])  # stack these as rows
    t = cp.Constant(np.ones(3, ))
    objective = cp.Minimize(c @ x)
    if axis == 0:
        con = cp.constraints.SOC(t, X.T, axis=0)
        con_expect = [
            np.array([0, 1.16454469e+00, 7.67560451e-01]),
            np.array([[0, -9.74311819e-01, -1.28440860e-01],
                      [0, 6.37872081e-01, 7.56737724e-01]])
        ]
    else:
        con = cp.constraints.SOC(t, X, axis=1)
        con_expect = [
            np.array([0, 1.16454469e+00, 7.67560451e-01]),
            np.array([[0, 0],
                      [-9.74311819e-01, 6.37872081e-01],
                      [-1.28440860e-01, 7.56737724e-01]])
        ]
    obj_pair = (objective, -1.932105)
    con_pairs = [(con, con_expect)]
    var_pairs = [(x, np.array([0.83666003, -0.54772256]))]
    sth = SolverTestHelper(obj_pair, var_pairs, con_pairs)
    return sth


def socp_bounds_attr() -> SolverTestHelper:
    x = cp.Variable(bounds=[-1, 1])
    obj_pair = (cp.Minimize(x), -1)
    var_pair = (x, -1)
    con_pair = (x**2 <= 4, 0)
    sth = SolverTestHelper(obj_pair, [var_pair], [con_pair])
    return sth


def sdp_1(objective_sense) -> SolverTestHelper:
    """
    Solve "Example 8.3" from Convex Optimization by Boyd & Vandenberghe.

    Verify (1) optimal objective values, (2) that the dual variable to the PSD constraint
    belongs to the correct cone (i.e. the dual variable is itself PSD), and (3) that
    complementary slackness holds with the PSD primal variable and its dual variable.
    """
    rho = cp.Variable(shape=(4, 4), symmetric=True)
    constraints = [0.6 <= rho[0, 1], rho[0, 1] <= 0.9,
                   0.8 <= rho[0, 2], rho[0, 2] <= 0.9,
                   0.5 <= rho[1, 3], rho[1, 3] <= 0.7,
                   -0.8 <= rho[2, 3], rho[2, 3] <= -0.4,
                   rho[0, 0] == 1, rho[1, 1] == 1, rho[2, 2] == 1, rho[3, 3] == 1,
                   rho >> 0]
    if objective_sense == 'min':
        obj = cp.Minimize(rho[0, 3])
        obj_pair = (obj, -0.39)
    elif objective_sense == 'max':
        obj = cp.Maximize(rho[0, 3])
        obj_pair = (obj, 0.23)
    else:
        raise RuntimeError('Unknown objective_sense.')
    con_pairs = [(c, None) for c in constraints]
    var_pairs = [(rho, None)]
    sth = SolverTestHelper(obj_pair, var_pairs, con_pairs)
    return sth


def sdp_2() -> SolverTestHelper:
    """
    Example SDO2 from MOSEK 9.2 documentation.
    """
    X1 = cp.Variable(shape=(2, 2), symmetric=True)
    X2 = cp.Variable(shape=(4, 4), symmetric=True)
    C1 = np.array([[1, 0], [0, 6]])
    A1 = np.array([[1, 1], [1, 2]])
    C2 = np.array([[1, -3, 0, 0], [-3, 2, 0, 0], [0, 0, 1, 0], [0, 0, 0, 0]])
    A2 = np.array([[0, 1, 0, 0], [1, -1, 0, 0], [0, 0, 0, 0], [0, 0, 0, -3]])
    b = 23
    k = -3
    var_pairs = [
        (X1, np.array([[21.04711571, 4.07709873],
                       [4.07709873, 0.7897868]])),
        (X2, np.array([[5.05366214, -3., 0., 0.],
                       [-3., 1.78088676, 0., 0.],
                       [0., 0., 0., 0.],
                       [0., 0., 0., -0.]]))
    ]
    con_pairs = [
        (cp.trace(A1 @ X1) + cp.trace(A2 @ X2) == b, -0.83772234),
        (X2[0, 1] <= k, 11.04455278),
        (X1 >> 0, np.array([[21.04711571, 4.07709873],
                            [4.07709873, 0.7897868]])),
        (X2 >> 0, np.array([[1., 1.68455405, 0., 0.],
                            [1.68455405, 2.83772234, 0., 0.],
                            [0., 0., 1., 0.],
                            [0., 0., 0., 2.51316702]]))
    ]
    obj_expr = cp.Minimize(cp.trace(C1 @ X1) + cp.trace(C2 @ X2))
    obj_pair = (obj_expr, 52.40127214)
    sth = SolverTestHelper(obj_pair, var_pairs, con_pairs)
    return sth


def expcone_1() -> SolverTestHelper:
    """
    min   3 * x[0] + 2 * x[1] + x[2]
    s.t.  0.1 <= x[0] + x[1] + x[2] <= 1
          x >= 0
          x[0] >= x[1] * exp(x[2] / x[1])
    """
    x = cp.Variable(shape=(3, 1))
    cone_con = cp.constraints.ExpCone(x[2], x[1], x[0])
    constraints = [cp.sum(x) <= 1.0,
                   cp.sum(x) >= 0.1,
                   x >= 0,
                   cone_con]
    obj = cp.Minimize(3 * x[0] + 2 * x[1] + x[2])
    obj_pair = (obj, 0.23534820622420757)
    expect_exp = [np.array([-1.35348213]), np.array([-0.35348211]), np.array([0.64651792])]
    con_pairs = [(constraints[0], 0),
                 (constraints[1], 2.3534821130067614),
                 (constraints[2], np.zeros(shape=(3, 1))),
                 (constraints[3], expect_exp)]
    expect_x = np.array([[0.05462721], [0.02609378], [0.01927901]])
    var_pairs = [(x, expect_x)]
    sth = SolverTestHelper(obj_pair, var_pairs, con_pairs)
    return sth


def expcone_socp_1() -> SolverTestHelper:
    """
    A random risk-parity portfolio optimization problem.
    """
    sigma = np.array([[1.83, 1.79, 3.22],
                      [1.79, 2.18, 3.18],
                      [3.22, 3.18, 8.69]])
    L = np.linalg.cholesky(sigma)
    c = 0.75
    t = cp.Variable(name='t')
    x = cp.Variable(shape=(3,), name='x')
    s = cp.Variable(shape=(3,), name='s')
    e = cp.Constant(np.ones(3, ))
    objective = cp.Minimize(t - c * e @ s)
    con1 = cp.norm(L.T @ x, p=2) <= t
    con2 = cp.constraints.ExpCone(s, e, x)
    # SolverTestHelper data
    obj_pair = (objective, 4.0751197)
    var_pairs = [
        (x, np.array([0.576079, 0.54315, 0.28037])),
        (s, np.array([-0.55150, -0.61036, -1.27161])),
    ]
    con_pairs = [
        (con1, 1.0),
        (con2, [np.array([-0.75, -0.75, -0.75]),
                np.array([-1.16363, -1.20777, -1.70371]),
                np.array([1.30190, 1.38082, 2.67496])]
         )
    ]
    sth = SolverTestHelper(obj_pair, var_pairs, con_pairs)
    return sth


def sdp_pcp_1() -> SolverTestHelper:
    """
    Example sdp and power cone.
    """
    Sigma = np.array([[ 0.4787481 , -0.96924914],
                      [-0.96924914,  2.77788598]])

    x = cp.Variable(shape=(2,1))
    y = cp.Variable(shape=(2,1))
    X = cp.Variable(shape=(2,2), symmetric=True)
    M1 = cp.vstack([X, x.T])
    M2 = cp.vstack([x, np.ones((1, 1))])
    M3 = cp.hstack([M1, M2])

    var_pairs = [
        (x, np.array([[0.72128204],
                      [0.27871796]])),
        (y, np.array([[0.01],
                      [0.01]])),
        (X, np.array([[0.52024779, 0.20103426],
                      [0.20103426, 0.0776837 ]])),
    ]
    con_pairs = [
        (cp.sum(x) == 1, -0.1503204799112807),
        (x >= 0, np.array([[-0.],
                           [-0.]])),
        (y >= 0.01, np.array([[0.70705506],
                               [0.70715844]])),
        (M3 >> 0, np.array([[ 0.4787481 , -0.96924914, -0.07516024],
                            [-0.96924914,  2.77788598, -0.07516024],
                            [-0.07516024, -0.07516024,  0.07516094]])),
        (cp.PowCone3D(x, np.ones((2,1)), y, 0.9), [np.array([[1.17878172e-09],
                                                             [3.05162243e-09]]),
                                                   np.array([[9.20157640e-10],
                                                             [9.40823207e-10]]),
                                                   np.array([[2.41053358e-10],
                                                             [7.43432462e-10]])]),
        ]
    obj_expr = cp.Minimize(cp.trace(Sigma @ X) + cp.norm(y, p=2))
    obj_pair = (obj_expr, 0.089301671322676)
    sth = SolverTestHelper(obj_pair, var_pairs, con_pairs)
    return sth


def pcp_1() -> SolverTestHelper:
    """
    Use a 3D power cone formulation for

    min 3 * x[0] + 2 * x[1] + x[2]
    s.t. norm(x,2) <= y
         x[0] + x[1] + 3*x[2] >= 1.0
         y <= 5
    """
    x = cp.Variable(shape=(3,))
    y_square = cp.Variable()
    epis = cp.Variable(shape=(3,))
    constraints = [cp.constraints.PowCone3D(np.ones(3), epis, x, cp.Constant([0.5, 0.5, 0.5])),
                   cp.sum(epis) <= y_square,
                   x[0] + x[1] + 3 * x[2] >= 1.0,
                   y_square <= 25]
    obj = cp.Minimize(3 * x[0] + 2 * x[1] + x[2])
    expect_x = np.array([-3.874621860638774, -2.129788233677883, 2.33480343377204])
    expect_epis = expect_x ** 2
    expect_x = np.round(expect_x, decimals=5)
    expect_epis = np.round(expect_epis, decimals=5)
    expect_y_square = 25
    var_pairs = [(x, expect_x),
                 (epis, expect_epis),
                 (y_square, expect_y_square)]
    expect_ineq1 = 0.7793969212001993
    expect_ineq2 = 2.865602615049077 / 10
    expect_pc = [np.array([4.30209047, 1.29985494, 1.56211543]),
                 np.array([0.28655796, 0.28655796, 0.28655796]),
                 np.array([2.22062898, 1.22062899, -1.33811302])]
    con_pairs = [(constraints[0], expect_pc),
                 (constraints[1], expect_ineq2),
                 (constraints[2], expect_ineq1),
                 (constraints[3], expect_ineq2)]
    obj_pair = (obj, -13.548638904065102)
    sth = SolverTestHelper(obj_pair, var_pairs, con_pairs)
    return sth


def pcp_2() -> SolverTestHelper:
    """
    Reformulate

        max  (x**0.2)*(y**0.8) + z**0.4 - x
        s.t. x + y + z/2 == 2
             x, y, z >= 0
    Into

        max  x3 + x4 - x0
        s.t. x0 + x1 + x2 / 2 == 2,
             (x0, x1, x3) in Pow3D(0.2)
             (x2, 1.0, x4) in Pow3D(0.4)
    """
    x = cp.Variable(shape=(3,))
    hypos = cp.Variable(shape=(2,))
    objective = cp.Minimize(-cp.sum(hypos) + x[0])
    arg1 = cp.hstack([x[0], x[2]])
    arg2 = cp.hstack(([x[1], 1.0]))
    pc_con = cp.constraints.PowCone3D(arg1, arg2, hypos, [0.2, 0.4])
    expect_pc_con = [np.array([1.48466366, 0.24233184]),
                     np.array([0.48466367, 0.83801333]),
                     np.array([-1., -1.])]
    con_pairs = [
        (x[0] + x[1] + 0.5 * x[2] == 2, 0.4846636697795672),
        (pc_con, expect_pc_con)
    ]
    obj_pair = (objective, -1.8073406786220672)
    var_pairs = [
        (x, np.array([0.06393515, 0.78320961, 2.30571048])),
        (hypos, None)
    ]
    sth = SolverTestHelper(obj_pair, var_pairs, con_pairs)
    return sth


def pcp_3() -> SolverTestHelper:
    from scipy.optimize import Bounds, minimize
    w = cp.Variable((2, 1))
    D = np.array([
       [-1.0856306,   0.99734545],
       [0.2829785,  -1.50629471],
       [-0.57860025,  1.65143654],
       [-2.42667924, -0.42891263],
       [1.26593626, -0.8667404],
       [-0.67888615, -0.09470897],
       [1.49138963, -0.638902]])  # T-by-N
    """
    Minimize ||D @ w||_p s.t. 0 <= w, sum(w) == 1.
        Refer to https://docs.mosek.com/modeling-cookbook/powo.html#p-norm-cones
    """
    p = 1/0.4
    T = D.shape[0]
    t = cp.Variable()
    d = cp.Variable((T, 1))
    ones = np.ones((T, 1))

    powcone = cp.constraints.PowCone3D(d, t * ones, D @ w, 1/p)
    constraints = [cp.sum(w) == 1, w >= 0, powcone, cp.sum(d) == t]
    con_pairs = [
        (constraints[0], -1.51430),
        (constraints[1], np.array([0.0, 0.0])),
        (constraints[2], [
              np.array([[0.40000935],
                        [0.40000935],
                        [0.40000935],
                        [0.40000935],
                        [0.40000935],
                        [0.40000935],
                        [0.40000935]]),
              np.array([[2.84369172e-03],
                        [1.22657446e-01],
                        [1.12146997e-01],
                        [3.45802205e-01],
                        [2.76327461e-05],
                        [1.27539057e-02],
                        [3.75878155e-03]]),
              np.array([[-0.04031276],
                        [0.38577107],
                        [-0.36558292],
                        [0.71847219],
                        [0.00249992],
                        [0.09919715],
                        [-0.04765863]])]),
        (constraints[3], 0.40000935)
    ]

    def univar_obj(w0):
        return np.linalg.norm(D[:, 0] * w0 + D[:, 1] * (1 - w0), ord=p)
    univar_bounds = Bounds([0], [1])
    univar_res = minimize(univar_obj, np.array([0.4]), bounds=univar_bounds, tol=1e-16)
    w_opt = np.array([[univar_res.x], [1 - univar_res.x]])

    obj_pair = (cp.Minimize(t), univar_res.fun)
    var_pairs = [(d, np.array([
                       [7.17144981e-03],
                       [3.09557056e-01],
                       [2.83038570e-01],
                       [8.72785905e-01],
                       [6.92995408e-05],
                       [3.21904516e-02],
                       [9.48918352e-03]])),
                 (w, w_opt),
                 (t, np.array([univar_res.fun]))]
    sth = SolverTestHelper(obj_pair, var_pairs, con_pairs)
    return sth


def mi_lp_0() -> SolverTestHelper:
    x = cp.Variable(shape=(2,))
    bool_var = cp.Variable(boolean=True)
    con_pairs = [(x == bool_var, None),
                 (bool_var == 0, None)]
    obj_pair = (cp.Minimize(cp.norm(x, 1) + 1.0), 1)
    var_pairs = [(x, np.array([0, 0])),
                 (bool_var, 0)]
    sth = SolverTestHelper(obj_pair, var_pairs, con_pairs)
    return sth


def mi_lp_1() -> SolverTestHelper:
    x = cp.Variable(2, name='x')
    boolvar = cp.Variable(boolean=True)
    intvar = cp.Variable(integer=True)
    objective = cp.Minimize(-4 * x[0] - 5 * x[1])
    constraints = [2 * x[0] + x[1] <= intvar,
                   x[0] + 2 * x[1] <= 3 * boolvar,
                   x >= 0,
                   intvar == 3 * boolvar,
                   intvar == 3]
    obj_pair = (objective, -9)
    var_pairs = [(x, np.array([1, 1])),
                 (boolvar, 1),
                 (intvar,  3)]
    con_pairs = [(c, None) for c in constraints]
    sth = SolverTestHelper(obj_pair, var_pairs, con_pairs)
    return sth


def mi_lp_2() -> SolverTestHelper:
    # Instance "knapPI_1_50_1000_1" from "http://www.diku.dk/~pisinger/genhard.c"
    n = 50
    c = 995
    z = 8373
    coeffs = [[1, 94, 485, 0], [2, 506, 326, 0], [3, 416, 248, 0],
              [4, 992, 421, 0], [5, 649, 322, 0], [6, 237, 795, 0],
              [7, 457, 43, 1], [8, 815, 845, 0], [9, 446, 955, 0],
              [10, 422, 252, 0], [11, 791, 9, 1], [12, 359, 901, 0],
              [13, 667, 122, 1], [14, 598, 94, 1], [15, 7, 738, 0],
              [16, 544, 574, 0], [17, 334, 715, 0], [18, 766, 882, 0],
              [19, 994, 367, 0], [20, 893, 984, 0], [21, 633, 299, 0],
              [22, 131, 433, 0], [23, 428, 682, 0], [24, 700, 72, 1],
              [25, 617, 874, 0], [26, 874, 138, 1], [27, 720, 856, 0],
              [28, 419, 145, 0], [29, 794, 995, 0], [30, 196, 529, 0],
              [31, 997, 199, 1], [32, 116, 277, 0], [33, 908, 97, 1],
              [34, 539, 719, 0], [35, 707, 242, 0], [36, 569, 107, 0],
              [37, 537, 122, 0], [38, 931, 70, 1], [39, 726, 98, 1],
              [40, 487, 600, 0], [41, 772, 645, 0], [42, 513, 267, 0],
              [43, 81, 972, 0], [44, 943, 895, 0], [45, 58, 213, 0],
              [46, 303, 748, 0], [47, 764, 487, 0], [48, 536, 923, 0],
              [49, 724, 29, 1], [50, 789, 674, 0]]  # index, p / w / x
    X = cp.Variable(n, boolean=True)
    objective = cp.Maximize(cp.sum(cp.multiply([i[1] for i in coeffs], X)))
    constraints = [cp.sum(cp.multiply([i[2] for i in coeffs], X)) <= c]
    obj_pair = (objective, z)
    con_pairs = [(constraints[0], None)]
    var_pairs = [(X, None)]
    sth = SolverTestHelper(obj_pair, var_pairs, con_pairs)
    return sth


def mi_lp_3() -> SolverTestHelper:
    # infeasible (but relaxable) test case
    x = cp.Variable(4, boolean=True)
    from cvxpy.expressions.constants import Constant
    objective = cp.Maximize(Constant(1))
    constraints = [x[0] + x[1] + x[2] + x[3] <= 2,
                   x[0] + x[1] + x[2] + x[3] >= 2,
                   x[0] + x[1] <= 1,
                   x[0] + x[2] <= 1,
                   x[0] + x[3] <= 1,
                   x[2] + x[3] <= 1,
                   x[1] + x[3] <= 1,
                   x[1] + x[2] <= 1]
    obj_pair = (objective, -np.inf)
    con_pairs = [(c, None) for c in constraints]
    var_pairs = [(x, None)]
    sth = SolverTestHelper(obj_pair, var_pairs, con_pairs)
    return sth


def mi_lp_4() -> SolverTestHelper:
    """Test MI without constraints"""
    x = cp.Variable(boolean=True)
    from cvxpy.expressions.constants import Constant
    objective = cp.Maximize(Constant(0.23) * x)
    obj_pair = (objective, 0.23)
    var_pairs = [(x, 1)]
    sth = SolverTestHelper(obj_pair, var_pairs, [])
    return sth


def mi_lp_5() -> SolverTestHelper:
    # infeasible boolean problem - https://trac.sagemath.org/ticket/31962#comment:48
    z = cp.Variable(11, boolean=True)
    constraints = [z[2] + z[1] == 1,
                   z[4] + z[3] == 1,
                   z[6] + z[5] == 1,
                   z[8] + z[7] == 1,
                   z[10] + z[9] == 1,
                   z[4] + z[1] <= 1,
                   z[2] + z[3] <= 1,
                   z[6] + z[2] <= 1,
                   z[1] + z[5] <= 1,
                   z[8] + z[6] <= 1,
                   z[5] + z[7] <= 1,
                   z[10] + z[8] <= 1,
                   z[7] + z[9] <= 1,
                   z[9] + z[4] <= 1,
                   z[3] + z[10] <= 1]
    obj = cp.Minimize(0)
    obj_pair = (obj, np.inf)
    con_pairs = [(c, None) for c in constraints]
    var_pairs = [(z, None)]
    sth = SolverTestHelper(obj_pair, var_pairs, con_pairs)
    return sth

def mi_lp_6() -> SolverTestHelper:
    "Test MILP for timelimit and no feasible solution"
    n = 70
    m = 70
    x = cp.Variable((n,), boolean=True, name="x")
    y = cp.Variable((n,), name="y")
    z = cp.Variable((m,), pos=True, name="z")
    A = np.random.rand(m, n)
    b = np.random.rand(m)
    objective = cp.Maximize(cp.sum(y))
    constraints = [
        A @ y <= b,
        y <= 1,
        cp.sum(x) >= 10,
        cp.sum(x) <= 20,
        z[0] + z[1] + z[2] >= 10,
        z[3] + z[4] + z[5] >= 5,
        z[6] + z[7] + z[8] >= 7,
        z[9] + z[10] >= 8,
        z[11] + z[12] >= 6,
        z[13] + z[14] >= 3,
        z[15] + z[16] >= 2,
        z[17] + z[18] >= 1,
        z[19] >= 2,
        z[20] >= 1,
        z[21] >= 1,
        z[22] >= 1,
        z[23] >= 1,
        z[24] >= 1,
        z[25] >= 1,
        z[26] >= 1,
        z[27] >= 1,
        z[28] >= 1,
        z[29] >= 1,
    ]
    return SolverTestHelper(
        (objective, None),
        [(x, None), (y, None), (z, None)],
        [(con, None) for con in constraints]
    )



def mi_lp_7() -> SolverTestHelper:
    """Problem that takes significant time to solve - for testing time/iteration limits"""
    np.random.seed(0)
    n = 24 * 8
    c = cp.Variable((n,), pos=True)
    d = cp.Variable((n,), pos=True)
    c_or_d = cp.Variable((n,), boolean=True)
    big = 1e3
    s = cp.cumsum(c * 0.9 - d)
    p = np.random.random(n)
    objective = cp.Maximize(p @ (d - c))
    constraints = [
        d <= 1,
        c <= 1,
        s >= 0,
        s <= 1,
        c <= c_or_d * big,
        d <= (1 - c_or_d) * big,
    ]
    return SolverTestHelper(
        (objective, None),
        [(c, None,), (d, None), (c_or_d, None)],
        [(con, None) for con in constraints]
    )


def mi_socp_1() -> SolverTestHelper:
    """
    Formulate the following mixed-integer SOCP with cvxpy
        min 3 * x[0] + 2 * x[1] + x[2] +  y[0] + 2 * y[1]
        s.t. norm(x,2) <= y[0]
             norm(x,2) <= y[1]
             x[0] + x[1] + 3*x[2] >= 0.1
             y <= 5, y integer.
    """
    x = cp.Variable(shape=(3,))
    y = cp.Variable(shape=(2,), integer=True)
    constraints = [cp.norm(x, 2) <= y[0],
                   cp.norm(x, 2) <= y[1],
                   x[0] + x[1] + 3 * x[2] >= 0.1,
                   y <= 5]
    obj = cp.Minimize(3 * x[0] + 2 * x[1] + x[2] + y[0] + 2 * y[1])
    obj_pair = (obj, 0.21363997604807272)
    var_pairs = [(x, np.array([-0.78510265, -0.43565177,  0.44025147])),
                 (y, np.array([1, 1]))]
    con_pairs = [(c, None) for c in constraints]  # no dual values for mixed-integer problems.
    sth = SolverTestHelper(obj_pair, var_pairs, con_pairs)
    return sth


def mi_socp_2() -> SolverTestHelper:
    """
    An (unnecessarily) SOCP-based reformulation of MI_LP_1.
    Doesn't use SOC objects.
    """
    x = cp.Variable(shape=(2,))
    bool_var = cp.Variable(boolean=True)
    int_var = cp.Variable(integer=True)
    objective = cp.Minimize(-4 * x[0] - 5 * x[1])
    constraints = [2 * x[0] + x[1] <= int_var,
                   (x[0] + 2 * x[1]) ** 2 <= 9 * bool_var,
                   x >= 0,
                   int_var == 3 * bool_var,
                   int_var == 3]
    obj_pair = (objective, -9)
    var_pairs = [(x, np.array([1, 1])),
                 (bool_var, 1),
                 (int_var, 3)]
    con_pairs = [(con, None) for con in constraints]
    sth = SolverTestHelper(obj_pair, var_pairs, con_pairs)
    return sth


def mi_pcp_0() -> SolverTestHelper:
    """
    max  x3 + x4 - x0
    s.t. x0 + x1 + x2 / 2 == 2,
         (x0, x1, x3) in Pow3D(0.2)
         (x2, q, x4) in Pow3D(0.4)
         0.1 <= q <= 1.9,
         q integer
    """
    x = cp.Variable(shape=(3,))
    hypos = cp.Variable(shape=(2,))
    q = cp.Variable(integer=True)
    objective = cp.Minimize(-cp.sum(hypos) + x[0])
    arg1 = cp.hstack([x[0], x[2]])
    arg2 = cp.hstack(([x[1], q]))
    pc_con = cp.constraints.PowCone3D(arg1, arg2, hypos, [0.2, 0.4])
    con_pairs = [
        (x[0] + x[1] + 0.5 * x[2] == 2, None),
        (pc_con, None),
        (0.1 <= q, None),
        (q <= 1.9, None)
    ]
    obj_pair = (objective, -1.8073406786220672)
    var_pairs = [
        (x, np.array([0.06393515, 0.78320961, 2.30571048])),
        (hypos, None),
        (q, 1.0)
    ]
    sth = SolverTestHelper(obj_pair, var_pairs, con_pairs)
    return sth


class StandardTestLPs:

    @staticmethod
    def test_lp_0(solver, places: int = 4, duals: bool = True, **kwargs) -> SolverTestHelper:
        sth = lp_0()
        sth.solve(solver, **kwargs)
        sth.verify_primal_values(places)
        sth.verify_objective(places)
        if duals:
            sth.check_complementarity(places)
        return sth

    @staticmethod
    def test_lp_1(solver, places: int = 4, duals: bool = True, **kwargs) -> SolverTestHelper:
        sth = lp_1()
        sth.solve(solver, **kwargs)
        sth.verify_objective(places)
        sth.verify_primal_values(places)
        if duals:
            sth.verify_dual_values(places)
        return sth

    @staticmethod
    def test_lp_2(solver, places: int = 4, duals: bool = True, **kwargs) -> SolverTestHelper:
        sth = lp_2()
        sth.solve(solver, **kwargs)
        sth.verify_objective(places)
        sth.verify_primal_values(places)
        if duals:
            sth.verify_dual_values(places)
        return sth

    @staticmethod
    def test_lp_3(solver, places: int = 4, **kwargs) -> SolverTestHelper:
        sth = lp_3()
        sth.solve(solver, **kwargs)
        sth.verify_objective(places)
        return sth

    @staticmethod
    def test_lp_4(solver, places: int = 4, **kwargs) -> SolverTestHelper:
        sth = lp_4()
        sth.solve(solver, **kwargs)
        sth.verify_objective(places)
        return sth

    @staticmethod
    def test_lp_5(solver, places: int = 4, duals: bool = True,  **kwargs) -> SolverTestHelper:
        sth = lp_5()
        sth.solve(solver, **kwargs)
        sth.verify_objective(places)
        sth.check_primal_feasibility(places)
        if duals:
            sth.check_complementarity(places)
            sth.check_dual_domains(places)
        return sth

    @staticmethod
    def test_lp_6(solver, places: int = 4, duals: bool = True, **kwargs) -> SolverTestHelper:
        sth = lp_6()
        sth.solve(solver, **kwargs)
        sth.verify_objective(places)
        sth.check_primal_feasibility(places)
        if duals:
            sth.check_complementarity(places)
            sth.check_dual_domains(places)
        return sth

    @staticmethod
    def test_lp_7(solver, places: int = 4, duals: bool = True, **kwargs) -> SolverTestHelper:
        sth = lp_7()
        import sdpap
        if sdpap.sdpacall.sdpacall.get_backend_info()["gmp"]:
            sth.solve(solver, **kwargs)
            sth.verify_objective(places)
        return sth

    @staticmethod
    def test_lp_bound_attr(
            solver,
            places: int = 4,
            duals: bool = True,
            **kwargs
        ) -> SolverTestHelper:
        sth = lp_bound_attr()
        # Verify the bounded variables data path is used
        data, _, _ = sth.prob.get_problem_data(solver)
        assert data[s.LOWER_BOUNDS] is not None, \
            f"{solver}: LOWER_BOUNDS should be populated"
        assert data[s.UPPER_BOUNDS] is not None, \
            f"{solver}: UPPER_BOUNDS should be populated"
        dims = data[s.DIMS]
        assert dims.nonneg == 0, \
            f"{solver}: bounds should not be converted to {dims.nonneg} inequality constraints"
        sth.solve(solver, **kwargs)
        sth.verify_objective(places)
        sth.verify_primal_values(places)
        if duals:
            sth.verify_dual_values(places)
        return sth

    @staticmethod
    def test_mi_lp_0(solver, places: int = 4, **kwargs) -> SolverTestHelper:
        sth = mi_lp_0()
        sth.solve(solver, **kwargs)
        sth.verify_objective(places)
        sth.verify_primal_values(places)
        return sth

    @staticmethod
    def test_mi_lp_1(solver, places: int = 4, **kwargs) -> SolverTestHelper:
        sth = mi_lp_1()
        sth.solve(solver, **kwargs)
        sth.verify_objective(places)
        sth.verify_primal_values(places)
        return sth

    @staticmethod
    def test_mi_lp_2(solver, places: int = 4, **kwargs) -> SolverTestHelper:
        sth = mi_lp_2()
        sth.solve(solver, **kwargs)
        sth.verify_objective(places)
        sth.verify_primal_values(places)
        return sth

    @staticmethod
    def test_mi_lp_3(solver, places: int = 4, **kwargs) -> SolverTestHelper:
        sth = mi_lp_3()
        sth.solve(solver, **kwargs)
        sth.verify_objective(places)
        sth.verify_primal_values(places)
        return sth

    @staticmethod
    def test_mi_lp_4(solver, places: int = 4, **kwargs) -> SolverTestHelper:
        sth = mi_lp_4()
        sth.solve(solver, **kwargs)
        sth.verify_objective(places)
        sth.verify_primal_values(places)
        return sth

    @staticmethod
    def test_mi_lp_5(solver, places: int = 4, **kwargs) -> SolverTestHelper:
        sth = mi_lp_5()
        sth.solve(solver, **kwargs)
        sth.verify_objective(places)
        sth.verify_primal_values(places)
        return sth

    @staticmethod
    def test_mi_lp_6(solver, places: int = 4, **kwargs) -> SolverTestHelper:
        sth = mi_lp_6()
        sth.solve(solver, **kwargs)
        sth.verify_objective(places)
        sth.verify_primal_values(places)
        return sth

class StandardTestQPs:

    @staticmethod
    def test_qp_0(solver, places: int = 4, duals: bool = True, **kwargs) -> SolverTestHelper:
        sth = qp_0()
        sth.solve(solver, **kwargs)
        sth.verify_primal_values(places)
        sth.verify_objective(places)
        if duals:
            sth.check_complementarity(places)
            sth.verify_dual_values(places)
        return sth

    @staticmethod
    def test_qp_bound_attr(
            solver,
            places: int = 4,
            **kwargs
        ) -> SolverTestHelper:
        sth = qp_bound_attr()
        # Verify the bounded variables data path is used
        data, _, _ = sth.prob.get_problem_data(solver)
        assert data[s.LOWER_BOUNDS] is not None, \
            f"{solver}: LOWER_BOUNDS should be populated"
        assert data[s.UPPER_BOUNDS] is not None, \
            f"{solver}: UPPER_BOUNDS should be populated"
        assert data['n_ineq'] == 0, \
            f"{solver}: bounds should not be converted to {data['n_ineq']} inequality constraints"
        sth.solve(solver, **kwargs)
        sth.verify_objective(places)
        sth.verify_primal_values(places)
        return sth


class StandardTestSOCPs:

    @staticmethod
    def test_socp_0(solver, places: int = 4, duals: bool = True, **kwargs) -> SolverTestHelper:
        sth = socp_0()
        sth.solve(solver, **kwargs)
        sth.verify_objective(places)
        sth.verify_primal_values(places)
        if duals:
            sth.check_complementarity(places)
        return sth

    @staticmethod
    def test_socp_1(solver, places: int = 4, duals: bool = True, **kwargs) -> SolverTestHelper:
        sth = socp_1()
        sth.solve(solver, **kwargs)
        sth.verify_objective(places)
        sth.verify_primal_values(places)
        if duals:
            sth.check_complementarity(places)
            sth.verify_dual_values(places)
        return sth

    @staticmethod
    def test_socp_2(solver, places: int = 4, duals: bool = True, **kwargs) -> SolverTestHelper:
        sth = socp_2()
        sth.solve(solver, **kwargs)
        sth.verify_objective(places)
        sth.verify_primal_values(places)
        if duals:
            sth.check_complementarity(places)
            sth.verify_dual_values(places)
        return sth

    @staticmethod
    def test_socp_3ax0(solver, places: int = 3, duals: bool = True, **kwargs) -> SolverTestHelper:
        sth = socp_3(axis=0)
        sth.solve(solver, **kwargs)
        sth.verify_objective(places)
        sth.verify_primal_values(places)
        if duals:
            sth.check_complementarity(places)
            sth.verify_dual_values(places)
        return sth

    @staticmethod
    def test_socp_3ax1(solver, places: int = 3, duals: bool = True, **kwargs) -> SolverTestHelper:
        sth = socp_3(axis=1)
        sth.solve(solver, **kwargs)
        sth.verify_objective(places)
        sth.verify_primal_values(places)
        if duals:
            sth.check_complementarity(places)
            sth.verify_dual_values(places)
        return sth

    @staticmethod
    def test_mi_socp_1(solver, places: int = 4, **kwargs) -> SolverTestHelper:
        sth = mi_socp_1()
        sth.solve(solver, **kwargs)
        sth.verify_objective(places)
        sth.verify_primal_values(places)
        return sth

    @staticmethod
    def test_mi_socp_2(solver, places: int = 4, **kwargs) -> SolverTestHelper:
        sth = mi_socp_2()
        sth.solve(solver, **kwargs)
        sth.verify_objective(places)
        sth.verify_primal_values(places)
        return sth

    @staticmethod
    def test_socp_bounds_attr(solver, places: int = 4, **kwargs) -> SolverTestHelper:
        sth = socp_bounds_attr()
        sth.solve(solver, **kwargs)
        sth.verify_objective(places)
        sth.verify_primal_values(places)
        return sth


class StandardTestSDPs:

    @staticmethod
    def test_sdp_1min(solver, places: int = 4, duals: bool = True, **kwargs) -> SolverTestHelper:
        sth = sdp_1('min')
        sth.solve(solver, **kwargs)
        sth.verify_objective(places=2)  # only 2 digits recorded.
        sth.check_primal_feasibility(places)
        if duals:
            sth.check_complementarity(places)
            sth.check_dual_domains(places)
        return sth

    @staticmethod
    def test_sdp_1max(solver, places: int = 4, duals: bool = True, **kwargs) -> SolverTestHelper:
        sth = sdp_1('max')
        sth.solve(solver, **kwargs)
        sth.verify_objective(places=2)  # only 2 digits recorded.
        sth.check_primal_feasibility(places)
        if duals:
            sth.check_complementarity(places)
            sth.check_dual_domains(places)
        return sth

    @staticmethod
    def test_sdp_2(solver, places: int = 2, duals: bool = True, **kwargs) -> SolverTestHelper:
        # places is set to 3 rather than 4, because analytic solution isn't known.
        sth = sdp_2()
        sth.solve(solver, **kwargs)
        sth.verify_objective(places)
        sth.check_primal_feasibility(places)
        sth.verify_primal_values(places)
        if duals:
            sth.check_complementarity(places)
            sth.check_dual_domains(places)
        return sth


class StandardTestECPs:

    @staticmethod
    def test_expcone_1(solver, places: int = 4, duals: bool = True, **kwargs) -> SolverTestHelper:
        sth = expcone_1()
        sth.solve(solver, **kwargs)
        sth.verify_objective(places)
        sth.verify_primal_values(places)
        if duals:
            sth.check_complementarity(places)
            sth.verify_dual_values(places)
        return sth


class StandardTestMixedCPs:

    @staticmethod
    def test_exp_soc_1(solver, places: int = 3, duals: bool = True, **kwargs) -> SolverTestHelper:
        sth = expcone_socp_1()
        sth.solve(solver, **kwargs)
        sth.verify_objective(places)
        sth.verify_primal_values(places)
        if duals:
            sth.check_complementarity(places)
            sth.verify_dual_values(places)
        return sth

    @staticmethod
    def test_sdp_pcp_1(solver, places: int = 3, duals: bool = False, **kwargs) -> SolverTestHelper:
        sth = sdp_pcp_1()
        sth.solve(solver, **kwargs)
        sth.verify_objective(places)
        sth.check_primal_feasibility(places)
        sth.verify_primal_values(places)
        if duals:
            sth.check_complementarity(places)
            sth.check_dual_domains(places)
        return sth


class StandardTestPCPs:

    @staticmethod
    def test_pcp_1(solver, places: int = 3, duals: bool = True, **kwargs) -> SolverTestHelper:
        sth = pcp_1()
        sth.solve(solver, **kwargs)
        sth.verify_objective(places)
        sth.check_primal_feasibility(places)
        sth.verify_primal_values(places)
        if duals:
            sth.check_complementarity(places)
            sth.verify_dual_values(places)
        return sth

    @staticmethod
    def test_pcp_2(solver, places: int = 3, duals: bool = True, **kwargs) -> SolverTestHelper:
        sth = pcp_2()
        sth.solve(solver, **kwargs)
        sth.verify_objective(places)
        sth.check_primal_feasibility(places)
        sth.verify_primal_values(places)
        if duals:
            sth.check_complementarity(places)
            sth.verify_dual_values(places)
        return sth

    @staticmethod
    def test_pcp_3(solver, places: int = 3, duals: bool = True, **kwargs) -> SolverTestHelper:
        sth = pcp_3()
        sth.solve(solver, **kwargs)
        sth.verify_objective(places)
        sth.check_primal_feasibility(places)
        sth.verify_primal_values(places)
        if duals:
            sth.check_complementarity(places)
            sth.verify_dual_values(places)

    @staticmethod
    def test_mi_pcp_0(solver, places: int = 3, **kwargs) -> SolverTestHelper:
        sth = mi_pcp_0()
        sth.solve(solver, **kwargs)
        sth.verify_objective(places)
        sth.check_primal_feasibility(places)
        sth.verify_primal_values(places)
        return sth