"""
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
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.PSD):
                # TODO: move this to PSD.dual_violation
                dv = con.dual_value
                eigs = np.linalg.eigvalsh(dv)
                min_eig = np.min(eigs)
                self.tester.assertGreaterEqual(min_eig, -(10**(-places)))
            elif isinstance(con, cp.constraints.ExpCone):
                # TODO: implement this (preferably with ExpCone.dual_violation)
                raise NotImplementedError()
            elif isinstance(con, cp.constraints.SOC):
                # TODO: implement this (preferably with SOC.dual_violation)
                raise NotImplementedError()
            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)
            elif isinstance(con, cp.constraints.PowCone3D):
                raise NotImplementedError()
            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.PSD):
                dv = con.dual_value
                pv = con.args[0].value
                comp = cp.scalar_product(pv, dv).value
            elif isinstance(con, (cp.constraints.ExpCone,
                                  cp.constraints.SOC,
                                  cp.constraints.NonPos,
                                  cp.constraints.Zero)):
                comp = cp.scalar_product(con.args, con.dual_value).value
            elif isinstance(con, cp.constraints.PowCone3D):
                comp = cp.scalar_product(con.args[:3], con.dual_value).value
            elif isinstance(con, (cp.constraints.Inequality,
                                  cp.constraints.Equality)):
                comp = cp.scalar_product(con.expr, con.dual_value).value
            elif isinstance(con, cp.constraints.PowConeND):
                msg = '\nPowConeND dual variables not implemented;' \
                       + '\nSkipping complementarity check.'
                warnings.warn(msg)
            else:
                raise ValueError('Unknown constraint type %s.' % type(con))
            self.tester.assertAlmostEqual(comp, 0, places)

    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 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))
    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 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 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_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_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


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


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 = 3, 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


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.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.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.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.verify_primal_values(places)
        return sth