\(\newcommand{\W}[1]{ \; #1 \; }\) \(\newcommand{\R}[1]{ {\rm #1} }\) \(\newcommand{\B}[1]{ {\bf #1} }\) \(\newcommand{\D}[2]{ \frac{\partial #1}{\partial #2} }\) \(\newcommand{\DD}[3]{ \frac{\partial^2 #1}{\partial #2 \partial #3} }\) \(\newcommand{\Dpow}[2]{ \frac{\partial^{#1}}{\partial {#2}^{#1}} }\) \(\newcommand{\dpow}[2]{ \frac{ {\rm d}^{#1}}{{\rm d}\, {#2}^{#1}} }\)

lp_box.hpp

View page source

lp_box Source Code

namespace CppAD { // BEGIN_CPPAD_NAMESPACE

// BEGIN PROTOTYPE
template <class Vector>
bool lp_box(
    size_t        level   ,
    const Vector& A       ,
    const Vector& b       ,
    const Vector& c       ,
    const Vector& d       ,
    size_t        maxitr  ,
    Vector&       xout    )
// END PROTOTYPE
{   double inf = std::numeric_limits<double>::infinity();
    //
    size_t m = b.size();
    size_t n = c.size();
    //
    CPPAD_ASSERT_KNOWN(
        level <= 3, "lp_box: level is greater than 3");
    CPPAD_ASSERT_KNOWN(
        size_t(A.size()) == m * n, "lp_box: size of A is not m * n"
    );
    CPPAD_ASSERT_KNOWN(
        size_t(d.size()) == n, "lp_box: size of d is not n"
    );
    if( level > 0 )
    {   std::cout << "start lp_box\n";
        CppAD::abs_print_mat("A", m, n, A);
        CppAD::abs_print_mat("b", m, 1, b);
        CppAD::abs_print_mat("c", n, 1, c);
        CppAD::abs_print_mat("d", n, 1, d);
    }
    //
    // count number of limits
    size_t n_limit = 0;
    for(size_t j = 0; j < n; j++)
    {   if( d[j] < inf )
            n_limit += 1;
    }
    //
    // A_simplex and b_simplex define the extended constraints
    Vector A_simplex((m + 2 * n_limit) * (2 * n) ), b_simplex(m + 2 * n_limit);
    for(size_t i = 0; i < size_t(A_simplex.size()); i++)
        A_simplex[i] = 0.0;
    //
    // put A * x + b <= 0 in A_simplex, b_simplex
    for(size_t i = 0; i < m; i++)
    {   b_simplex[i] = b[i];
        for(size_t j = 0; j < n; j++)
        {   // x_j^+ coefficient (positive component)
            A_simplex[i * (2 * n) + 2 * j]     =   A[i * n + j];
            // x_j^- coefficient (negative component)
            A_simplex[i * (2 * n) + 2 * j + 1] = - A[i * n + j];
        }
    }
    //
    // put | x_j | <= d_j in A_simplex, b_simplex
    size_t i_limit = 0;
    for(size_t j = 0; j < n; j++) if( d[j] < inf )
    {
        // x_j^+ <= d_j constraint
        b_simplex[ m + 2 * i_limit]                         = - d[j];
        A_simplex[(m + 2 * i_limit) * (2 * n) + 2 * j]      = 1.0;
        //
        // x_j^- <= d_j constraint
        b_simplex[ m + 2 * i_limit + 1]                         = - d[j];
        A_simplex[(m + 2 * i_limit + 1) * (2 * n) + 2 * j + 1]  = 1.0;
        //
        ++i_limit;
    }
    //
    // c_simples
    Vector c_simplex(2 * n);
    for(size_t j = 0; j < n; j++)
    {   // x_j+ component
        c_simplex[2 * j]     = c[j];
        // x_j^- component
        c_simplex[2 * j + 1] = - c[j];
    }
    size_t level_simplex = 0;
    if( level >= 2 )
        level_simplex = level - 1;
    //
    Vector x_simplex(2 * n);
    bool ok = CppAD::simplex_method(
        level_simplex, A_simplex, b_simplex, c_simplex, maxitr, x_simplex
    );
    for(size_t j = 0; j < n; j++)
        xout[j] = x_simplex[2 * j] - x_simplex[2 * j + 1];
    if( level > 0 )
    {   CppAD::abs_print_mat("xout", n, 1, xout);
        if( ok )
            std::cout << "end lp_box: ok = true\n";
        else
            std::cout << "end lp_box: ok = false\n";
    }
    return ok;
}

} // END_CPPAD_NAMESPACE