qp_solver.h

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

    
            // make `i' the last index in `B'
            int  k = B.size()-1;
            int  l = B[ i];
            if ( i < k) {
                inv_M_B.swap( i, k);
                std::swap( in_B[ B[ i]], in_B[ B[ k]]);
                        B[ i] =         B[ k];
                minus_c_B[ i] = minus_c_B[ k];
                 two_D_Bj[ i] =  two_D_Bj[ k];
            }
    
            // update basis
            B.pop_back();               // `i' leaves basis
            in_B[ l] = -1;
            if ( is_phase_I && ( l >= 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.remove( k);
    
            // update status
            x_B.pop_back();
            q_x.pop_back();
    
            // notify pricing strategy
            strategyP->leaving_basis( l);
        }
    }

    
    void  update_basis( Tag_true)       // LP
    {
        
        CGAL_optimisation_debug {
            vout2 << std::endl << "<--> non-basic variable " << j
                  << " replaces basic variable " << B[ i] << std::endl;
        }
         
    
        // update basis
        int  l = B[ i];
        in_B[ B[ i] ] = -1;             // `i' leaves basis
        in_B[    j  ] = i;              // `j' enters basis
        B[ i] = j;                      // `j' replaces `i' in basis
        if ( is_phase_I && ( j >= qp_n)) ++art_basic;
        if ( is_phase_I && ( l >= 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.replace( i, q_x.begin());
    
        // update status
        minus_c_B[ i] = -ET( is_phase_I ? aux_c[ j] : qp_c[ j]);
    
        // mark variable as entered
        j = -1;
    
        // notify pricing strategy
        strategyP->leaving_basis( l);
    }
    */    

    // ratio test
    // ----------
    void  ratio_test( )
        {
            
            CGAL_optimisation_debug {
                vout2 << std::endl
                      << "Ratio Test" << std::endl
                      << "----------" << std::endl;
            }
             
    
            // initialize
            
            init_ratio_test();
                              
    
            // compute `q_lambda' and `q_x'
            
            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
    void init_ratio_test( )
        {
            // 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());
            }
    
            
            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;
                }
                 
            }
             
    
            // 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);
            }
    
            
            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;
                }
                 
            }
             
        }
    
    void  store_2_D_Bj( Tag_true)       // LP
        {
            // nop
        }
    
    
    // computation of `q_lambda' and `q_x'
    void  compute_q( Tag_false)         // QP
        {
            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
            
            compute_current_solution();
                                       
        }
    
    
    // append variable to basis
    void  append_variable( )
        {
            
            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]);
            }

	    CGAL_optimisation_assertion( check_basis( Is_lp()));
        }
    
    // remove variable from basis
    void  remove_variable( )
        {
            
            CGAL_optimisation_debug {
                vout2 << std::endl << "<-- basic variable "
                      << B[ i] << " leaves basis" << std::endl;
            }
             
    
            // make `i' the last index in `B'
            int  k = B.size()-1;
            int  l = B[ i];
            if ( i < k) {
                inv_M_B.swap( i, k);
                std::swap( in_B[ B[ i]], in_B[ B[ k]]);
                        B[ i] =         B[ k];
                minus_c_B[ i] = minus_c_B[ k];
                 two_D_Bj[ i] =  two_D_Bj[ k];
            }
    
            // update basis
            B.pop_back();               // `i' leaves basis
            in_B[ l] = -1;
            if ( is_phase_I && ( l >= 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.remove( k);
    
            // update status
            x_B.pop_back();
            q_x.pop_back();
    
            // notify pricing strategy
            strategyP->leaving_basis( l);