qp_solver.h

来自「CGAL is a collaborative effort of severa」· C头文件 代码 · 共 1,686 行 · 第 1/4 页

    
    
    int                      art_basic; // number of basic artificial variables
    
    
    // auxiliary problem
    A_artificial             art_A;     // artificial part of constraint matrix
    
    C_auxiliary              aux_c;     // auxiliary objective vector
    
    
    // additional variables
    Values                   b;         // exact version of `qp_b'
    Values                   minus_c_B; // exact version of `-qp_c'
                                        // restricted to basis
    
    
    bool                     is_phase_I;// flag indicating phase I
    
    
    int                      j;         // index of entering variable `x_j'
    
    
    Values                   A_j;       // exact version of j-th column of A
    Values                   two_D_Bj;  // exact version of twice the j-th
                                        // column of D restricted to basis B
    
    
    Values                   q_lambda;  // used in the ratio test
    Values                   q_x;       // --------- " ----------
    
    
    int                      i;         // index of leaving variable `x_i'
    ET                       x_i;       // numerator of leaving variable `x_i'
    ET                       q_i;       // corresponding `q_i'
    
    
    ET                       mu_j;      //   numerator of `t_j'
    ET                       nu;        // denominator of `t_j'
    
    

    
    // ------------------------------------------------------------------------
    
    // pivot step
    void  pivot_step( );
    
    
    // set-up functions
    void  set_up_auxiliary_problem( );
    
    void  set_up_basis( );
    
    void  set_up_initial_solution( );
    
    void  set_up_additional_variables( );
    
    
    // pricing
    // -------
    void  pricing( )
        {
            
            CGAL_optimisation_debug {
                vout2 << std::endl
                      << "Pricing" << std::endl
                      << "-------" << std::endl;
            }
             
    
            if ( is_phase_I && ( art_basic == 0)) {
                j = -1;                 // all artificial variables non-basic
    
                
                CGAL_optimisation_debug {
                    vout2 << "artificial variables are non-basic" << std::endl;
                }
                 
            } else {
    
                // call pricing strategy
                j = strategyP->pricing();
    
                
                CGAL_optimisation_debug {
                    if ( j < 0) {
                        vout2 << "entering variable: none" << std::endl;
                    } else {
                        vout1 << "  ";
                        vout << "entering"; vout2 << " variable"; vout << ": ";
                        vout  << j; vout2 << std::endl;
                    }
                }
                 
            }
        }
    
    /*
    // ratio test
    // ----------
    void  ratio_test( )
        {
            
            CGAL_optimisation_debug {
                vout2 << std::endl
                      << "Ratio Test" << std::endl
                      << "----------" << std::endl;
            }
             
    
            // compute `q_lambda' and `q_x'
            
            CGAL_optimisation_debug {
                
                vout2 << "     A_j: ";
                if ( vout2.verbose()) {
                    std::copy( A_j.begin(), A_j.begin()+qp_m,
                               std::ostream_iterator<ET>( vout2.out(), " "));
                    vout2.out() << std::endl;
                }
                 
            }
             
            
            compute_q( Is_lp());
            
            
            CGAL_optimisation_debug {
                
                vout2 << "     q_x: ";
                if ( vout2.verbose()) {
                    std::copy( q_x.begin(), q_x.begin()+B.size(),
                               std::ostream_iterator<ET>( vout2.out(), " "));
                    vout2.out() << std::endl;
                }
                 
            }
             
    
            // check `t_i's
            
            x_i = et_1;                 // trick: initialize
            q_i = et_0;                 // minimum with +oo
            
            
            Value_iterator  x_it = x_B.begin();
            Value_iterator  q_it = q_x.begin();
            for ( unsigned int k = 0; k < B.size(); ++k, ++x_it, ++q_it) {
                if ( ( *q_it > et_0) && ( ( *x_it * q_i) < ( x_i * *q_it))) {
                    i = k; x_i = *x_it; q_i = *q_it;
                }
            }
            
            
    
            // check `t_j'
            
            check_t_j( Is_lp());
                                
    
            
            CGAL_optimisation_debug {
                vout2 << std::endl;
                for ( unsigned int k = 0; k < B.size(); ++k) {
                    vout2 << "t_" << k << ": " << x_B[ k] << '/' << q_x[ k]
                          << ( ( q_i > et_0) && ( i == (int)k) ? " *" : "")
                          << std::endl;
                }
                if ( ! ( CGAL::check_tag( Is_lp()) || is_phase_I)) {
                    vout2 << "t_j: " << mu_j << '/' << nu
                          << ( ( q_i > et_0) && ( i < 0) ? " *" : "")
                          << std::endl;
                }
                vout2 << std::endl;
                if ( q_i > et_0) {
                    if ( i < 0) {
                        vout2 << "leaving variable: none" << std::endl;
                    } else {
                        vout1 << ", ";
                        vout  << "leaving"; vout2 << " variable"; vout << ": ";
                        vout  << B[ i];
                        vout2 << " (= B[ " << i << "])" << std::endl;
                    }
                }
            }
             
        }

    
    // initialization of ratio-test/update loop
    void init_ratio_test_update_loop( )
        {
            // store exact version of `A_j' (implicit conversion to ET)
            if ( j < qp_n) {
                // original variable
                std::copy(  qp_A[ j     ],  qp_A[ j     ]+qp_m, A_j.begin());
            } else {
                // artificial variable
                std::copy( art_A[ j-qp_n], art_A[ j-qp_n]+qp_m, A_j.begin());
            }
    
            // store exact version of `2 D_{B,j}'
            store_2_D_Bj( Is_lp());
        }

    // storing of exact version of `2 D_{B,j}'
    void  store_2_D_Bj( Tag_false)      // QP
        {
            if ( j < qp_n) {
                // original variable
                Access_D_Bj  access_D_Bj( qp_D[ j], nt_0, 0, qp_n);
                std::transform( D_Bj_iterator( B.begin(), access_D_Bj),
                                D_Bj_iterator( B.end  (), access_D_Bj),
                                two_D_Bj.begin(),
                                std::bind1st( std::multiplies<ET>(), et_2));
            } else {
                // artificial variable
                std::fill_n( two_D_Bj.begin(), B.size(), et_0);
            }
        }
    
    void  store_2_D_Bj( Tag_true)       // LP
        {
            // nop
        }
    
    
    // computation of `q_lambda' and `q_x'
    void  compute_q( Tag_false)         // QP
        {
            
            CGAL_optimisation_debug {
                
                vout2 << "  2 D_Bj: ";
                if ( vout2.verbose()) {
                    std::copy( two_D_Bj.begin(), two_D_Bj.begin()+B.size(),
                               std::ostream_iterator<ET>( vout2.out(), " "));
                    vout2.out() << std::endl;
                }
                 
            }
             
    
            inv_M_B.multiply( A_j.begin(), two_D_Bj.begin(),
                              q_lambda.begin(), q_x.begin());
    
            
            CGAL_optimisation_debug {
                vout2 << std::endl;
                
                vout2 << "q_lambda: ";
                if ( vout2.verbose()) {
                    std::copy( q_lambda.begin(), q_lambda.begin()+qp_m,
                               std::ostream_iterator<ET>( vout2.out(), " "));
                    vout2.out() << std::endl;
                }
                 
            }
             
        }
    
    void  compute_q( Tag_true)          // LP
        {
            inv_M_B.multiply_x( A_j.begin(), q_x.begin());
        }
    
    
    // computation and checking of `t_j'
    void  check_t_j( Tag_false)         // QP
        {
            // compute `nu'
            nu = std::inner_product( q_x.begin(), q_x.end(),
                                     two_D_Bj.begin(),
                 std::inner_product( q_lambda.begin(), q_lambda.end(),
                                     A_j.begin(),
                                     ( j < qp_n) ? -et_2*d*ET( qp_D[ j][ j])
                                                 : et_0));
    
            if ( ! is_phase_I) {
    
                // compute `mu_j'
                mu_j = std::inner_product( x_B.begin(), x_B.end(),
                                           two_D_Bj.begin(),
                       std::inner_product( lambda.begin(), lambda.end(),
                                           A_j.begin(),
                                           d * ( is_phase_I ? ET( aux_c[ j])
                                                            : ET( qp_c[ j]))));
    
                // check `t_j'
                if ( ( nu < et_0) && ( ( mu_j * q_i) > ( x_i * nu))) {
                    i = -1; q_i = et_1;
                }
            }
        }
    
    void  check_t_j( Tag_true)          // LP
        {
            // nop
        }

    
    // update
    // ------
    void  update( )
        {
            
            CGAL_optimisation_debug {
                vout2 << std::endl
                      << "Update" << std::endl
                      << "------";
            }
             
    
            // update basis and basis inverse
            
            update_basis( Is_lp());
                                   
    
            // compute current solution
            
            if ( is_phase_I) {
                inv_M_B.multiply_l( minus_c_B.begin(), lambda.begin());
                inv_M_B.multiply_x(         b.begin(),    x_B.begin());
            } else {
                inv_M_B.multiply( b.begin(), minus_c_B.begin(),
                                  lambda.begin(), x_B.begin());
            }
            
            
            CGAL_optimisation_debug {
                vout2 << std::endl;
                
                vout2 << "  -c_B: ";
                if ( vout2.verbose()) {
                    std::copy( minus_c_B.begin(), minus_c_B.begin()+B.size(),
                               std::ostream_iterator<ET>( vout2.out(), " "));
                    vout2.out() << std::endl;
                }
                 
                
                vout2 << "lambda: ";
                if ( vout2.verbose()) {
                    std::copy( lambda.begin(), lambda.begin()+qp_m,
                               std::ostream_iterator<ET>( vout2.out(), " "));
                    vout2.out() << std::endl;
                }
                 
                
                vout2 << "   x_B: ";
                if ( vout2.verbose()) {
                    std::copy( x_B.begin(), x_B.begin()+B.size(),
                               std::ostream_iterator<ET>( vout2.out(), " "));
                    vout2.out() << std::endl;
                }
                 
            }
             
        }

    
    // update of basis and basis inverse
    void  update_basis( Tag_false)      // QP
    {
        // append variable to basis
        if ( ( i < 0) || ( B.size() == (unsigned int)qp_m)) {
    
            
            CGAL_optimisation_debug {
                vout2 << std::endl << "--> non-basic variable "
                      << j << " enters basis" << std::endl;
            }
             
    
            // update basis
            unsigned int  k = B.size();
            B.push_back( j);            // `j' enters basis
            in_B[ j] = k;
            if ( is_phase_I && ( j >= qp_n)) ++art_basic;
    
            
            CGAL_optimisation_debug {
                vout2 << "new ";
                
                vout2 << "basis: ";
                if ( vout2.verbose()) {
                    std::copy( B.begin(), B.end(),
                               std::ostream_iterator<int>( vout2.out(), " "));
                    vout2.out() << std::endl;
                }
                 
                vout3 << "new basis-inverse:" << std::endl;
            }
             
    
            // update basis inverse
            inv_M_B.append( q_lambda.begin(), q_x.begin(), nu);
    
            // update status
            x_B.push_back( et_0);
            q_x.push_back( et_0);
            if ( k >= two_D_Bj.size()) {
                 two_D_Bj.push_back( et_0);
                minus_c_B.push_back( -( is_phase_I ? ET( aux_c[ j])
					           : ET(  qp_c[ j])));
            } else {
                minus_c_B[ k] = is_phase_I ? -ET( aux_c[ j]) : -ET( qp_c[ j]);
            }
    
            // mark variable as entered
            j = -1;
        }
    
        // remove variable from basis
        if ( i >= 0) {
    
            
            CGAL_optimisation_debug {
                vout2 << std::endl << "<-- basic variable "
                      << B[ i] << " leaves basis" << std::endl;
            }