Python scipy.optimize.linprog() Examples

def test_enzo_example(self):
        # http://projects.scipy.org/scipy/attachment/ticket/1252/lp2.py
        #
        # Translated from Octave code at:
        # http://www.ecs.shimane-u.ac.jp/~kyoshida/lpeng.htm
        # and placed under MIT licence by Enzo Michelangeli
        # with permission explicitly granted by the original author,
        # Prof. Kazunobu Yoshida
        c = [4, 8, 3, 0, 0, 0]
        A_eq = [
            [2, 5, 3, -1, 0, 0],
            [3, 2.5, 8, 0, -1, 0],
            [8, 10, 4, 0, 0, -1]]
        b_eq = [185, 155, 600]
        res = linprog(c=c, A_eq=A_eq, b_eq=b_eq,
                      method=self.method, options=self.options)
        _assert_success(res, desired_fun=317.5,
                        desired_x=[66.25, 0, 17.5, 0, 183.75, 0],
                        atol=6e-6, rtol=1e-7) 

def solve_linprog(self, nu):
        n_hs = len(self.hs)
        n_constraints = len(self.constraints.index)
        if self.last_linprog_n_hs == n_hs:
            return self.last_linprog_result
        c = np.concatenate((self.errors, [self.B]))
        A_ub = np.concatenate((self.gammas - self.eps, -np.ones((n_constraints, 1))), axis=1)
        b_ub = np.zeros(n_constraints)
        A_eq = np.concatenate((np.ones((1, n_hs)), np.zeros((1, 1))), axis=1)
        b_eq = np.ones(1)
        result = opt.linprog(c, A_ub=A_ub, b_ub=b_ub, A_eq=A_eq, b_eq=b_eq, method='simplex')
        Q = pd.Series(result.x[:-1], self.hs.index)
        dual_c = np.concatenate((b_ub, -b_eq))
        dual_A_ub = np.concatenate((-A_ub.transpose(), A_eq.transpose()), axis=1)
        dual_b_ub = c
        dual_bounds = [(None, None) if i == n_constraints else (0, None) for i in range(n_constraints + 1)]  # noqa: E501
        result_dual = opt.linprog(dual_c,
                                  A_ub=dual_A_ub,
                                  b_ub=dual_b_ub,
                                  bounds=dual_bounds,
                                  method='simplex')
        lambda_vec = pd.Series(result_dual.x[:-1], self.constraints.index)
        self.last_linprog_n_hs = n_hs
        self.last_linprog_result = (Q, lambda_vec, self.eval_gap(Q, lambda_vec, nu))
        return self.last_linprog_result 

def test_enzo_example_b(self):
        # rescued from https://github.com/scipy/scipy/pull/218
        c = [2.8, 6.3, 10.8, -2.8, -6.3, -10.8]
        A_eq = [[-1, -1, -1, 0, 0, 0],
                [0, 0, 0, 1, 1, 1],
                [1, 0, 0, 1, 0, 0],
                [0, 1, 0, 0, 1, 0],
                [0, 0, 1, 0, 0, 1]]
        b_eq = [-0.5, 0.4, 0.3, 0.3, 0.3]
        if self.method == "simplex":
            # Including the callback here ensures the solution can be
            # calculated correctly.
            res = linprog(c=c, A_eq=A_eq, b_eq=b_eq,
                          method=self.method, options=self.options,
                          callback=lambda x, **kwargs: None)
        else:
            with suppress_warnings() as sup:
                sup.filter(OptimizeWarning, "A_eq does not appear...")
                res = linprog(c=c, A_eq=A_eq, b_eq=b_eq,
                              method=self.method, options=self.options)
        _assert_success(res, desired_fun=-1.77,
                        desired_x=[0.3, 0.2, 0.0, 0.0, 0.1, 0.3]) 

def wasserstein_distance(p, q, D):
    """Wasserstein-distance
    p.shape=[m], q.shape=[n], D.shape=[m, n]
    p.sum()=1, q.sum()=1, p∈[0,1], q∈[0,1]
    """
    A_eq = []
    for i in range(len(p)):
        A = np.zeros_like(D)
        A[i, :] = 1
        A_eq.append(A.reshape(-1))
    for i in range(len(q)):
        A = np.zeros_like(D)
        A[:, i] = 1
        A_eq.append(A.reshape(-1))
    A_eq = np.array(A_eq)
    b_eq = np.concatenate([p, q])
    D = D.reshape(-1)
    result = linprog(D, A_eq=A_eq[:-1], b_eq=b_eq[:-1])
    return result.fun 

def linear_program_ineq(c, A, b):
    c = c.reshape((c.size,))
    b = b.reshape((b.size,))

    # make unbounded bounds
    bounds = []
    for i in range(c.size):
        bounds.append((None, None))

    A_ub, b_ub = -A, -b
    res = linprog(c, A_ub=A_ub, b_ub=b_ub, bounds=bounds, options={"disp": False}, method='simplex')
    if res.success:
        return res.x.reshape((c.size, 1))
    else:
        np.savez('bad_scipy_lp_ineq_{:010d}'.format(np.random.randint(int(1e9))),
                 c=c, A=A, b=b, res=res)
        raise Exception('Scipy did not solve the LP. Blame Scipy.') 

def test_zero_column_2(self):
        # detected in presolve?
        np.random.seed(0)
        m, n = 2, 4
        c = np.random.rand(n)
        c[1] = -1
        A_eq = np.random.rand(m, n)
        A_eq[:, 1] = 0
        b_eq = np.random.rand(m)

        A_ub = np.random.rand(m, n)
        A_ub[:, 1] = 0
        b_ub = np.random.rand(m)
        res = linprog(c, A_ub, b_ub, A_eq, b_eq, bounds=(None, None),
                      method=self.method, options=self.options)
        _assert_unbounded(res)
        assert_equal(res.nit, 0) 

def test_linprog_cyclic_bland_bug_8561(self):
        # Test that pivot row is chosen correctly when using Bland's rule
        c = np.array([7, 0, -4, 1.5, 1.5])
        A_ub = np.array([
            [4, 5.5, 1.5, 1.0, -3.5],
            [1, -2.5, -2, 2.5, 0.5],
            [3, -0.5, 4, -12.5, -7],
            [-1, 4.5, 2, -3.5, -2],
            [5.5, 2, -4.5, -1, 9.5]])
        b_ub = np.array([0, 0, 0, 0, 1])
        if self.method == "simplex":
            res = linprog(c, A_ub=A_ub, b_ub=b_ub,
                          options=dict(maxiter=100, bland=True),
                          method=self.method)
        else:
            res = linprog(c, A_ub=A_ub, b_ub=b_ub, options=dict(maxiter=100),
                          method=self.method)
        _assert_success(res, desired_x=[0, 0, 19, 16/3, 29/3]) 

def _assert_success(res, desired_fun=None, desired_x=None,
                    rtol=1e-8, atol=1e-8):
    # res: linprog result object
    # desired_fun: desired objective function value or None
    # desired_x: desired solution or None
    if not res.success:
        msg = "linprog status {0}, message: {1}".format(res.status,
                                                        res.message)
        raise AssertionError(msg)

    assert_equal(res.status, 0)
    if desired_fun is not None:
        assert_allclose(res.fun, desired_fun,
                        err_msg="converged to an unexpected objective value",
                        rtol=rtol, atol=atol)
    if desired_x is not None:
        assert_allclose(res.x, desired_x,
                        err_msg="converged to an unexpected solution",
                        rtol=rtol, atol=atol) 

def test_negative_variable(self):
        # Test linprog with a problem with one unbounded variable and
        # another with a negative lower bound.
        c = np.array([-1, 4]) * -1  # maximize
        A_ub = np.array([[-3, 1],
                         [1, 2]], dtype=np.float64)
        A_ub_orig = A_ub.copy()
        b_ub = [6, 4]
        x0_bounds = (-np.inf, np.inf)
        x1_bounds = (-3, np.inf)

        res = linprog(c, A_ub=A_ub, b_ub=b_ub, bounds=(x0_bounds, x1_bounds),
                      method=self.method, options=self.options)

        assert_equal(A_ub, A_ub_orig)   # user input not overwritten
        _assert_success(res, desired_fun=-80 / 7, desired_x=[-8 / 7, 18 / 7]) 

def test_network_flow(self):
        # A network flow problem with supply and demand at nodes
        # and with costs along directed edges.
        # https://www.princeton.edu/~rvdb/542/lectures/lec10.pdf
        c = [2, 4, 9, 11, 4, 3, 8, 7, 0, 15, 16, 18]
        n, p = -1, 1
        A_eq = [
            [n, n, p, 0, p, 0, 0, 0, 0, p, 0, 0],
            [p, 0, 0, p, 0, p, 0, 0, 0, 0, 0, 0],
            [0, 0, n, n, 0, 0, 0, 0, 0, 0, 0, 0],
            [0, 0, 0, 0, 0, 0, p, p, 0, 0, p, 0],
            [0, 0, 0, 0, n, n, n, 0, p, 0, 0, 0],
            [0, 0, 0, 0, 0, 0, 0, n, n, 0, 0, p],
            [0, 0, 0, 0, 0, 0, 0, 0, 0, n, n, n]]
        b_eq = [0, 19, -16, 33, 0, 0, -36]
        res = linprog(c=c, A_eq=A_eq, b_eq=b_eq,
                      method=self.method, options=self.options)
        _assert_success(res, desired_fun=755, atol=1e-6, rtol=1e-7) 

def test_network_flow_limited_capacity(self):
        # A network flow problem with supply and demand at nodes
        # and with costs and capacities along directed edges.
        # http://blog.sommer-forst.de/2013/04/10/
        cost = [2, 2, 1, 3, 1]
        bounds = [
            [0, 4],
            [0, 2],
            [0, 2],
            [0, 3],
            [0, 5]]
        n, p = -1, 1
        A_eq = [
            [n, n, 0, 0, 0],
            [p, 0, n, n, 0],
            [0, p, p, 0, n],
            [0, 0, 0, p, p]]
        b_eq = [-4, 0, 0, 4]

        if self.method == "simplex":
            # Including the callback here ensures the solution can be
            # calculated correctly, even when phase 1 terminated
            # with some of the artificial variables as pivots
            # (i.e. basis[:m] contains elements corresponding to
            # the artificial variables)
            res = linprog(c=cost, A_eq=A_eq, b_eq=b_eq, bounds=bounds,
                          method=self.method, options=self.options,
                          callback=lambda x, **kwargs: None)
        else:
            with suppress_warnings() as sup:
                sup.filter(RuntimeWarning, "scipy.linalg.solve\nIll...")
                sup.filter(OptimizeWarning, "A_eq does not appear...")
                sup.filter(OptimizeWarning, "Solving system with option...")
                res = linprog(c=cost, A_eq=A_eq, b_eq=b_eq, bounds=bounds,
                              method=self.method, options=self.options)
        _assert_success(res, desired_fun=14) 

def test_bug_6690(self):
        # https://github.com/scipy/scipy/issues/6690
        A_eq = np.array([[0., 0., 0., 0.93, 0., 0.65, 0., 0., 0.83, 0.]])
        b_eq = np.array([0.9626])
        A_ub = np.array([[0., 0., 0., 1.18, 0., 0., 0., -0.2, 0.,
                          -0.22],
                         [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
                         [0., 0., 0., 0.43, 0., 0., 0., 0., 0., 0.],
                         [0., -1.22, -0.25, 0., 0., 0., -2.06, 0., 0.,
                          1.37],
                         [0., 0., 0., 0., 0., 0., 0., -0.25, 0., 0.]])
        b_ub = np.array([0.615, 0., 0.172, -0.869, -0.022])
        bounds = np.array(
            [[-0.84, -0.97, 0.34, 0.4, -0.33, -0.74, 0.47, 0.09, -1.45, -0.73],
             [0.37, 0.02, 2.86, 0.86, 1.18, 0.5, 1.76, 0.17, 0.32, -0.15]]).T
        c = np.array([-1.64, 0.7, 1.8, -1.06, -1.16,
                      0.26, 2.13, 1.53, 0.66, 0.28])

        with suppress_warnings() as sup:
            sup.filter(RuntimeWarning, "scipy.linalg.solve\nIll...")
            sup.filter(OptimizeWarning, "Solving system with option...")
            sol = linprog(c, A_ub=A_ub, b_ub=b_ub, A_eq=A_eq, b_eq=b_eq,
                          bounds=bounds, method=self.method,
                          options=self.options)
        _assert_success(sol, desired_fun=-1.191, rtol=1e-6) 

def test_maxiter(self):
        # Test with a rather large problem (400 variables,
        # 40 constraints) generated by https://gist.github.com/denis-bz/8647461
        # test iteration limit
        A, b, c = lpgen_2d(20, 20)
        maxiter = np.random.randint(6) + 1  # problem takes 7 iterations
        res = linprog(c, A_ub=A, b_ub=b, method=self.method,
                      options={"maxiter": maxiter})
        # maxiter is independent of sparse/dense
        assert_equal(res.status, 1)
        assert_equal(res.nit, maxiter) 

def test_sparse_solve_options(self):
        A, b, c, N = magic_square(3)
        with suppress_warnings() as sup:
            sup.filter(OptimizeWarning, "A_eq does not appear...")
            sup.filter(OptimizeWarning, "Invalid permc_spec option")
            o = {key: self.options[key] for key in self.options}
            permc_specs = ('NATURAL', 'MMD_ATA', 'MMD_AT_PLUS_A',
                           'COLAMD', 'ekki-ekki-ekki')
            for permc_spec in permc_specs:
                o["permc_spec"] = permc_spec
                res = linprog(c, A_eq=A, b_eq=b, bounds=(0, 1),
                              method=self.method, options=o)
                _assert_success(res, desired_fun=1.730550597) 

def test_singleton_row_ub_1(self):
        # detected in presolve?
        c = [1, 1, 1, 2]
        A_ub = [[1, 0, 0, 0], [0, 2, 0, 0], [-1, 0, 0, 0], [1, 1, 1, 1]]
        b_ub = [1, 2, -2, 4]
        res = linprog(c, A_ub=A_ub, b_ub=b_ub,
                      bounds=[(None, None), (0, None), (0, None), (0, None)],
                      method=self.method, options=self.options)
        _assert_infeasible(res)
        assert_equal(res.nit, 0) 

def test_disp(self):
        # Test with a rather large problem (400 variables,
        # 40 constraints) generated by https://gist.github.com/denis-bz/8647461
        # test that display option does not break anything.
        A, b, c = lpgen_2d(20, 20)
        res = linprog(c, A_ub=A, b_ub=b, method=self.method,
                      options={"disp": True})
        # disp is independent of sparse/dense
        _assert_success(res, desired_fun=-64.049494229) 

def test_callback(self):
        def f():
            pass
        assert_raises(NotImplementedError, linprog, c=1, callback=f,
                      method=self.method) 

def _compute_probabilities(self, lists, rankings):
        '''
        Solve the optimization problem in
        (Multileaved Comparisons for Fast Online Evaluation, CIKM'14)

        lists: lists of document IDs
        rankings: a list of Ranking instances

        Return a list of probabilities for input rankings
        '''
        # probability constraints
        A_p_sum = np.array([1]*len(rankings))
        # unbiasedness constraints
        ub_cons = self._unbiasedness_constraints(lists, rankings)
        # sensitivity
        sensitivity = self._sensitivity(lists, rankings)

        # constraints
        A_eq = np.vstack((A_p_sum, ub_cons))
        b_eq = np.array([1.0] + [0.0]*ub_cons.shape[0])

        # solving the optimization problem
        res = linprog(sensitivity, # objective function
            A_eq=A_eq, b_eq=b_eq, # constraints
            bounds=[(0, 1)]*len(rankings) # 0 <= p <= 1
            )
        return res.success, res.x, res.fun 

def optimize(self):
        
        (A_in, b_in) = self.get_inequality_constraints()
        (A_equ, b_equ) = self.get_equality_constraints()
        c = self.get_cost_matrix()

        res = linprog(c, A_ub=A_in, b_ub=b_in, A_eq=A_equ, b_eq=b_equ)

        return res

    # def generate_ 

def test_alternate_initial_point(self):
        # Test with a rather large problem (400 variables,
        # 40 constraints) generated by https://gist.github.com/denis-bz/8647461
        # use "improved" initial point
        A, b, c = lpgen_2d(20, 20)
        with suppress_warnings() as sup:
            sup.filter(RuntimeWarning, "scipy.linalg.solve\nIll...")
            sup.filter(OptimizeWarning, "Solving system with option...")
            res = linprog(c, A_ub=A, b_ub=b, method=self.method,
                          options={"ip": True, "disp": True})
            # ip code is independent of sparse/dense
        _assert_success(res, desired_fun=-64.049494229) 

def test_infeasible_ub(self):
        # detected in presolve?
        c = [1]
        A_ub = [[2]]
        b_ub = 4
        bounds = (5, 6)
        res = linprog(c=c, A_ub=A_ub, b_ub=b_ub, bounds=bounds,
                      method=self.method, options=self.options)
        _assert_infeasible(res)
        assert_equal(res.nit, 0) 

def test_singleton_row_eq_1(self):
        # detected in presolve?
        c = [1, 1, 1, 2]
        A_eq = [[1, 0, 0, 0], [0, 2, 0, 0], [1, 0, 0, 0], [1, 1, 1, 1]]
        b_eq = [1, 2, 2, 4]
        res = linprog(c, A_eq=A_eq, b_eq=b_eq,
                      method=self.method, options=self.options)
        _assert_infeasible(res)
        assert_equal(res.nit, 0) 

def test_empty_constraint_1(self):
        # detected in presolve?
        res = linprog([-1, 1, -1, 1],
                      bounds=[(0, np.inf), (-np.inf, 0), (-1, 1), (-1, 1)],
                      method=self.method, options=self.options)
        _assert_unbounded(res)
        assert_equal(res.nit, 0) 

def test_magic_square_bug_7044(self):
        # test linprog with a problem with a rank-deficient A_eq matrix
        A, b, c, N = magic_square(3)
        with suppress_warnings() as sup:
            sup.filter(OptimizeWarning, "A_eq does not appear...")
            res = linprog(c, A_eq=A, b_eq=b, bounds=(0, 1),
                          method=self.method, options=self.options)
        _assert_success(res, desired_fun=1.730550597) 

def test_bounds_equal_but_infeasible2(self):
        c = [-4, 1]
        A_eq = [[7, -2], [0, 1], [2, -2]]
        b_eq = [14, 0, 3]
        bounds = [(2, 2), (0, None)]
        res = linprog(c=c, A_eq=A_eq, b_eq=b_eq, bounds=bounds,
                      method=self.method)
        _assert_infeasible(res) 

def test_bounds_equal_but_infeasible(self):
        c = [-4, 1]
        A_ub = [[7, -2], [0, 1], [2, -2]]
        b_ub = [14, 0, 3]
        bounds = [(2, 2), (0, None)]
        res = linprog(c=c, A_ub=A_ub, b_ub=b_ub, bounds=bounds,
                      method=self.method)
        _assert_infeasible(res) 

def test_callback(self):
        # Check that callback is as advertised
        callback_complete = [False]
        last_xk = []

        def cb(xk, **kwargs):
            kwargs.pop('tableau')
            assert_(isinstance(kwargs.pop('phase'), int))
            assert_(isinstance(kwargs.pop('nit'), int))

            i, j = kwargs.pop('pivot')
            assert_(np.isscalar(i))
            assert_(np.isscalar(j))

            basis = kwargs.pop('basis')
            assert_(isinstance(basis, np.ndarray))
            assert_(basis.dtype == np.int_)

            complete = kwargs.pop('complete')
            assert_(isinstance(complete, bool))
            if complete:
                last_xk.append(xk)
                callback_complete[0] = True
            else:
                assert_(not callback_complete[0])

            # no more kwargs
            assert_(not kwargs)

        c = np.array([-3, -2])
        A_ub = [[2, 1], [1, 1], [1, 0]]
        b_ub = [10, 8, 4]
        res = linprog(c, A_ub=A_ub, b_ub=b_ub, callback=cb, method=self.method)

        assert_(callback_complete[0])
        assert_allclose(last_xk[0], res.x) 

def test_unbounded_below_and_above(self):
        A = np.eye(3)
        b = np.array([1, 2, 3])
        c = np.ones(3)
        res = linprog(c, A_eq=A, b_eq=b, bounds=(-np.inf, np.inf),
                      method=self.method, options=self.options)
        _assert_success(res, desired_x=b, desired_fun=np.sum(b)) 

def test_bounded_above_only(self):
        A = np.eye(3)
        b = np.array([1, 2, 3])
        c = np.ones(3)
        res = linprog(c, A_eq=A, b_eq=b, bounds=(-np.inf, 4),
                      method=self.method, options=self.options)
        _assert_success(res, desired_x=b, desired_fun=np.sum(b)) 

def test_bounded_below_only(self):
        A = np.eye(3)
        b = np.array([1, 2, 3])
        c = np.ones(3)
        res = linprog(c, A_eq=A, b_eq=b, bounds=(0.5, np.inf),
                      method=self.method, options=self.options)
        _assert_success(res, desired_x=b, desired_fun=np.sum(b))