qp_solver.h

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

    
            // new basis found?
            if ( B.size() == (unsigned int)qp_m) { i = -1; }

	    CGAL_optimisation_assertion( check_basis( Is_lp()));
        }
    
    // update of basis and basis inverse
    void  update_basis( Tag_false)      // QP
        {
	  if ( B.size() < max_basis) {

            // append variable to basis
            append_variable();
    
            // remove variable from basis (if any)
            if ( i >= 0) remove_variable();

	  } else {
	    remove_variable();
            inv_M_B.multiply( A_j.begin(), two_D_Bj.begin(),
                              q_lambda.begin(), q_x.begin());
            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));
            append_variable();
	  }
        }
    
    
    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[ l ] = -1;             // `i' leaves basis
            in_B[ j ] = i;              // `j' enters basis
            B[ i] = j;                  // `j' replaces `i' in basis
            if ( is_phase_I) {
                if ( j >= qp_n) ++art_basic;
                if ( 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]);
    
            // notify pricing strategy
            strategyP->leaving_basis( l);
    
            // new basis found
            i = -1;
        }
    
    
    // computation of current solution
    void   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;
                }
                 
            }
             
        }
    
    
    // iterated ratio test and update
    // ------------------------------
    void  iterated_ratio_test_update( )
        {
            loop_until_basis_found( Is_lp());
        }
    
    // loop until new basis is found
    void  loop_until_basis_found( Tag_false)    // QP
        {
            while ( 
                    i >= 0
                          ) {
    
                // ratio test (iterated)
                ratio_test_iterated();
    
                // check for unboundedness
                if ( 
                     q_i == et_0
                                ) {
                    m_phase  = 3;
                    m_status = UNBOUNDED;
                    
                    CGAL_optimisation_debug {
                        vout1 << "  ";
                        vout << "problem is UNBOUNDED" << std::endl;
                    }
                     
                    return;
                }
    
                // update (iterated)
                update_iterated();
            }
        }
    
    void  loop_until_basis_found( Tag_true)     // LP
        {
            // nop
        }
    
    
    // ratio test (iterated)
    void  ratio_test_iterated( )
        {
            
            CGAL_optimisation_debug {
                vout2 << std::endl
                      << "Ratio Test (iterated)" << std::endl
                      << "---------------------" << std::endl;
            }
             
    
            // get `q_x'
            unsigned int  l = in_B[ j];
            std::copy( inv_M_B.column_begin( qp_m+l)+qp_m,
                       inv_M_B.column_end  ( qp_m+l)     , q_x.begin());
    
            
            CGAL_optimisation_debug {
                
		vout2 << "q_x[ " << l << "]: ";
                if ( vout2.verbose()) {
                    std::copy( q_x.begin(), q_x.begin()+B.size(),
                               std::ostream_iterator<ET>( vout2.out(), " "));
                    vout2.out() << std::endl;
                }
                 
            }
             
    
            // store frequently used values
            ET            x_j = x_B[ l];
            ET            q_j = q_x[ l];
            int  sign_q_i_q_j = CGAL_NTS sign( q_j);
    
            // initialize minimum
            x_i = ( sign_q_i_q_j < 0) ? et_1 : -et_1;   // trick: initialize
            q_i = et_0;                                 // minimum with +oo
    
            // check `t_i's
            
            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 (   ( k != l)
            
                  // t_i > 0 ?
                    && (/*   (( *q_it > et_0) && ( *x_it * q_j < x_j * *q_it))
			||*/ (( *q_it < et_0) && ( *x_it * q_j > x_j * *q_it)))
            
                  // t_i < t_min ?
                    && (   (   ( CGAL_NTS sign( *q_it) * sign_q_i_q_j > 0)
                            && ( *x_it * q_i > x_i * *q_it))
                        || (   ( CGAL_NTS sign( *q_it) * sign_q_i_q_j < 0)
                            && ( *x_it * q_i < x_i * *q_it)))) {
            
                    // store new minimum
                    i = k; x_i = *x_it; q_i = *q_it;
                    sign_q_i_q_j = CGAL_NTS sign( q_j) * CGAL_NTS sign( q_i);
                }
            }
             
    
            // check `t_j'
            
            if ( ( x_j > et_0) && ( CGAL_NTS sign( x_i) * sign_q_i_q_j <= 0)) {
                i = -1; q_i = et_1;
            }
             
    
            
            CGAL_optimisation_debug {
                vout2 << std::endl;
                for ( unsigned int k = 0; k < B.size(); ++k) {
                    if ( k != l) {
                        vout2 << "t_" << k << ": "
                              << x_j * q_x[k] - x_B[k] * q_j << '/' << q_x[k]
                              << ( ( q_i != et_0) && (i == (int)k) ? " *" : "")
                              << std::endl;
            
                    }
                }
                vout2 << "t_j: " << x_j << '/' << et_1
                      << ( ( q_i > et_0) && ( i < 0) ? " *" : "")
                      << std::endl << 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;
                    }
                }
            }
             
        }
    
    
    // update (iterated)
    void  update_iterated( )
        {
            if ( i >= 0) {
    
                
                CGAL_optimisation_debug {
                    vout2 << std::endl
                          << "Update (iterated)" << std::endl
                          << "-----------------";
                }
                 
    
                // update basis and basis inverse
                remove_variable();
    
		CGAL_optimisation_assertion( check_basis( Is_lp()));

                // compute current solution
                compute_current_solution();
            }
        }

    bool check_basis( Tag_true)       // LP
        {
            return true;
        }

    bool check_basis( Tag_false)
	{
	    if ( is_phase_I) return true;

	    Values  result( qp_m+B.size());
	    Values  col_l( qp_m, et_0), col_x( B.size());
	    unsigned int i, j;

	    // first part
	    for ( i = 0; i < (unsigned int)qp_m; ++i) {

		// get source column
		for ( j = 0; j < B.size(); ++j) col_x[ j] = qp_A[ B[ j]][ i];

		// compute target column
		inv_M_B.multiply( col_l.begin(), col_x.begin(),
				  result.begin(), result.begin()+qp_m);

		// check result
		/*
		std::copy( result.begin(), result.end(),
			   std::ostream_iterator<ET>( std::cerr, " "));
		std::cerr << endl;
		*/
		for ( j = 0; j < qp_m+B.size(); ++j) {
		    if ( ( ( j == i) && result[ j] != d) ||
			 ( ( j != i) && result[ j] != et_0)) return false;
		}
	    }

	    // second part
	    for ( i = 0; i < B.size(); ++i) {
		j = B[ i];

		// get source column
                std::copy( qp_A[ j], qp_A[ j]+qp_m, col_l.begin());
                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),
                                col_x.begin(),
                                std::bind1st( std::multiplies<ET>(), et_2));

		// compute target column
		inv_M_B.multiply( col_l.begin(), col_x.begin(),
				  result.begin(), result.begin()+qp_m);

		// check result
		/*
		std::copy( result.begin(), result.end(),
			   std::ostream_iterator<ET>( std::cerr, " "));
		std::cerr << endl;
		*/
		for ( j = 0; j < qp_m+B.size(); ++j) {
		    if ( ( ( j == i+qp_m) && result[ j] != d) ||
			 ( ( j != i+qp_m) && result[ j] != et_0)) return false;
		}
	    }
	    return true;
	}
    
};
  
CGAL_END_NAMESPACE

#include <CGAL/_QP_solver/Pricing_strategy_base.h>
#include <CGAL/_QP_solver/Full_exact_pricing.h>

CGAL_BEGIN_NAMESPACE

template < class Rep_ >
QP_solver<Rep_>::
    QP_solver( int verbose, std::ostream& stream)
        : nt_0( 0), nt_1( 1), nt_minus_1( -nt_1),
          et_0( 0), et_1( 1), et_2( 2),
          vout1( verbose == 1, stream),
          vout2( verbose >= 2, stream),
          vout3( verbose == 3, stream),
          vout ( verbose >  0, stream),
          qp_n( 0), qp_m( 0), inv_M_B( vout3), m_phase( 0),
          d( inv_M_B.denominator()),
          
          strategyP( new Pricing_strategy_default)
                                       
        {
            
            CGAL_optimisation_debug {
                vout2 << "======================================" << std::endl
                      << "The CGAL Solver for Quadratic Programs" << std::endl
                      << "======================================" << std::endl;
            }
             
        }

template < class Rep_ >
const typename QP_solver<Rep_>::Pricing_strategy&
QP_solver<Rep_>::
    pricing_strategy() const { return *strategyP; }
    
template < class Rep_ >
void
QP_solver<Rep_>::
    set_pricing_strategy( 
      typename QP_solver<Rep_>::Pricing_strategy& pricing_strategy)
        {
            
            CGAL_optimisation_debug {
                vout2 << std::endl
                      << "-----------------------" << std::endl
                      << "Pricing Strategy Change" << std::endl
                      << "-----------------------" << std::endl;
            }
             
    
            strategyP = &pricing_strategy;
            strategyP->set( *this, vout2);
        }

CGAL_END_NAMESPACE

#ifdef CGAL_CFG_NO_AUTOMATIC_TEMPLATE_INCLUSION
#  include <CGAL/_QP_solver/QP_solver.C>
#endif

#endif // CGAL_QP_SOLVER_H

// ===== EOF ==================================================================