qp_solver_impl.h

很多二维 三维几何计算算法 C++ 类库

                d_min = diff;
                q_min = q_k;
                ratio_test_bound_index = UPPER;
            }    
        }
    }
}

// test for one basic variable with mixed bounds,
// note that this function writes the member variables i, x_i, q_i and
// ratio_test_bound_index, although the second and third variable name
// are in the context of upper bounding misnomers
template < typename Q, typename ET, typename Tags >
void  QP_solver<Q, ET, Tags>::
test_mixed_bounds_dir_pos(int k, const ET& x_k, const ET& q_k, 
                                int& i_min, ET& d_min, ET& q_min)
{
    if (q_k > et0) {                                // check for lower bound
        if (k < qp_n) {                             // original variable
            if (*(qp_fl+k)) {
                ET  diff = x_k - (d * ET(qp_l[k]));
                if (diff * q_min < d_min * q_k) {
                    i_min = k;
                    d_min = diff;
                    q_min = q_k;
                    ratio_test_bound_index = LOWER;
                }
            }
        } else {                                    // artificial variable
            if (x_k * q_min < d_min * q_k) {
                i_min = k;
                d_min = x_k;
                q_min = q_k;
            }
        }
    } else {                                        // check for upper bound
        if ((q_k < et0) && (k < qp_n) && *(qp_fu+k)) {
            ET  diff = (d * ET(qp_u[k])) - x_k;
            if (diff * q_min < -(d_min * q_k)) {
                i_min = k;
                d_min = diff;
                q_min = -q_k;
                ratio_test_bound_index = UPPER;
            }
        }
    }
}

// test for one basic variable with mixed bounds,
// note that this function writes the member variables i, x_i, q_i and
// ratio_test_bound_index, although the second and third variable name
// are in the context of upper bounding misnomers
template < typename Q, typename ET, typename Tags >
void  QP_solver<Q, ET, Tags>::
test_mixed_bounds_dir_neg(int k, const ET& x_k, const ET& q_k, 
                                int& i_min, ET& d_min, ET& q_min)
{
    if (q_k < et0) {                                // check for lower bound
        if (k < qp_n) {                             // original variable
            if (*(qp_fl+k)) {
                ET  diff = x_k - (d * ET(qp_l[k]));
                if (diff * q_min < -(d_min * q_k)) {
                    i_min = k;
                    d_min = diff;
                    q_min = -q_k;
                    ratio_test_bound_index = LOWER;
                }
            }
        } else {                                    // artificial variable
            if (x_k * q_min < -(d_min * q_k)) {
                i_min = k;
                d_min = x_k;
                q_min = -q_k;
            }
        }
    } else {                                        // check for upper bound
        if ((q_k > et0) && (k < qp_n) && *(qp_fu+k)) {
            ET  diff = (d * ET(qp_u[k])) - x_k;
            if (diff * q_min < d_min * q_k) {
                i_min = k;
                d_min = diff;
                q_min = q_k;
                ratio_test_bound_index = UPPER;
            }
        }
    }
}    


template < typename Q, typename ET, typename Tags >                                         // QP case
void
QP_solver<Q, ET, Tags>::
ratio_test_2( Tag_false)
{
    // diagnostic output
    CGAL_qpe_debug {
	if ( vout2.verbose()) {
	    vout2.out() << std::endl
			<< "Ratio Test (Step 2)" << std::endl
			<< "-------------------" << std::endl;
	}
    }

    // compute `p_lambda' and `p_x' (Note: `p_...' is stored in `q_...')
    ratio_test_2__p( no_ineq);
 
    // diagnostic output
    CGAL_qpe_debug {
	if ( vout3.verbose()) {
	    vout3.out() << "p_lambda: ";
	    std::copy( q_lambda.begin(), q_lambda.begin()+C.size(),
		       std::ostream_iterator<ET>( vout3.out()," "));
	    vout3.out() << std::endl;
	    vout3.out() << "   p_x_O: ";
	    std::copy( q_x_O.begin(), q_x_O.begin()+B_O.size(),
		       std::ostream_iterator<ET>( vout3.out()," "));
	    vout3.out() << std::endl;
	    if ( has_ineq) {
		vout3.out() << "   p_x_S: ";
		std::copy( q_x_S.begin(), q_x_S.begin()+B_S.size(),
			   std::ostream_iterator<ET>( vout3.out()," "));
		vout3.out() << std::endl;
	    }
	    vout3.out() << std::endl;
	}
    }
    
    // Idea here: At this point, the goal is to increase \mu_j until either we
    // become optimal (\mu_j=0), or one of the variables in x^*_\hat{B} drops
    // down to zero.
    //
    // Let us see first how this is done in the standard-form case (where
    // Sven's thesis applies).  Eq. (2.11) in Sven's thesis holds, and by
    // multlying it by $M_\hat{B}^{-1}$ we obtain an equation for \lambda and
    // x^*_\hat{B}.  The interesting equation (the one for x^*_\hat{B}) looks
    // more or less as follows:
    //
    //    x(mu_j)      = x(0) + mu_j      q_it                          (1)
    //
    // where q_it is the vector from (2.12).  In paritcular, for
    // mu_j=mu_j(t_1) (i.e., if we plug the value of mu_j at the beginning of
    // this ratio step 2 into (1)) we have
    //
    //    x(mu_j(t_1)) = x(0) + mu_j(t_1) q_it                          (2)
    //
    // where x(mu_j(t_1)) is the current solution of the solver at this point
    // (i.e., at the beginning of ratio step 2).
    //
    // By subtracting (2) from (1) we can thus eliminate the "unkown" x(0)
    // (which is cheaper than computing it):
    //
    //    x(mu_j) = x(mu_j(t_1)) + (mu_j-mu_j(t_1)) q_it
    //                             ----------------
    //                                  := delta
    //
    // In order to compute for each variable x_k in \hat{B} the value of mu_j
    // for which x_k(mu_j) = 0, we thus evaluate
    //
    //                x(mu_j(t_1))_k
    //    delta_k:= - --------------
    //                    q_it_k
    //
    // The first variable in \hat{B} that hits zero "in the future" is then
    // the one whose delta_k equals
    //
    //    delta_min:= min {delta_k | k in \hat{B} and (q_it)_k < 0 }
    //    
    // So in order to handle the standard-form case, we would compute this
    // minimum.  Once we have delta_min, we need to check whether we get
    // optimal BEFORE a variable drops to zero.  As delta = mu_j - mu_j(t_1),
    // the latter is precisely the case if delta_min >= -mu_j(t_1).
    //
    // (Note: please forget the crap identitiy between (2.11) and (2.12); the
    // notation is misleading.)
    //
    // Now to the nonstandard-form case.
    
    // fw: By definition delta_min >= 0, such that initializing
    // delta_min with -mu_j(t_1) has the desired effect that a basic variable
    // is leaving only if 0 <= delta_min < -mu_j(t_1).
    //
    // The only initialization of delta_min as fraction x_i/q_i that works is
    // x_i=mu_j(t_1); q_i=-1; (see below).
    //
    // Since mu_j(t_1) has been computed in ratio test step 1 we can
    // reuse it.
      
    x_i = mu;                                     // initialize minimum
    q_i = -et1;                                        // with -mu_j(t_1) 

    Value_iterator  x_it = x_B_O.begin();
    Value_iterator  q_it = q_x_O.begin();
    Index_iterator  i_it;
    for ( i_it = B_O.begin(); i_it != B_O.end(); ++i_it, ++x_it, ++q_it) {
	if ( ( *q_it < et0) && ( ( *x_it * q_i) < ( x_i * *q_it))) {
	    i = *i_it; x_i = *x_it; q_i = *q_it;
	}
    }
    x_it = x_B_S.begin();
    q_it = q_x_S.begin();
    for ( i_it = B_S.begin(); i_it != B_S.end(); ++i_it, ++x_it, ++q_it) {
	if ( ( *q_it < et0) && ( ( *x_it * q_i) < ( x_i * *q_it))) {
	    i = *i_it; x_i = *x_it; q_i = *q_it;
	}
    }

    CGAL_qpe_debug {
	if ( vout2.verbose()) {
	    for ( unsigned int k = 0; k < B_O.size(); ++k) {
		vout2.out() << "mu_j_O_" << k << ": - "
			    << x_B_O[ k] << '/' << q_x_O[ k]
			    << ( ( q_i < et0) && ( i == B_O[ k]) ? " *" : "")
			    << std::endl;
	    }
	    for ( unsigned int k = 0; k < B_S.size(); ++k) {
		vout2.out() << "mu_j_S_" << k << ": - "
			    << x_B_S[ k] << '/' << q_x_S[ k]
			    << ( ( q_i < et0) && ( i == B_S[ k]) ? " *" : "")
			    << std::endl;
	    }
	    vout2.out() << std::endl;
	}
    }
    if ( i < 0) {
      vout2 << "leaving variable: none" << std::endl;
    } else {
      vout1 << ", ";
      vout  << "leaving"; vout2 << " variable"; vout << ": ";
      vout  << i;
      if ( vout2.verbose()) {
	if ( i < qp_n) {
	  vout2.out() << " (= B_O[ " << in_B[ i] << "]: original)"
		      << std::endl;
	} else {
	  vout2.out() << " (= B_S[ " << in_B[ i] << "]: slack)"
		      << std::endl;
	}
      }
    }
    
}

// update (step 1)
template < typename Q, typename ET, typename Tags >
void
QP_solver<Q, ET, Tags>::
update_1( )
{
    CGAL_qpe_debug {
	if ( vout2.verbose()) {
	    vout2.out() << std::endl;
	    if ( is_LP || is_phaseI) {
		vout2.out() << "------" << std::endl
			    << "Update" << std::endl
			    << "------" << std::endl;
	    } else {
		vout2.out() << "Update (Step 1)" << std::endl
			    << "---------------" << std::endl;
	    }
	}
    }

    // update basis & basis inverse
    update_1( Is_linear());
    CGAL_qpe_assertion(check_basis_inverse());
    
    // check the updated vectors r_C and r_S_B
    CGAL_expensive_assertion(check_r_C(Is_nonnegative()));
    CGAL_expensive_assertion(check_r_S_B(Is_nonnegative()));
    
    // check the vectors r_B_O and w in phaseII for QPs
    CGAL_qpe_debug {
        if (is_phaseII && is_QP) {
            CGAL_expensive_assertion(check_r_B_O(Is_nonnegative()));
            CGAL_expensive_assertion(check_w(Is_nonnegative()));
        }
    }

    // compute current solution
    compute_solution(Is_nonnegative());
    
    // check feasibility 
    CGAL_qpe_debug {
      if (j < 0 && !is_RTS_transition) // todo kf: is this too conservative?
                 // Note: the above condition is necessary because of the
                 // following.  In theory, it is true that the current
                 // solution is at this point in the solver always
                 // feasible. However, the solution has its x_j-entry equal to
                 // the current t from the pricing, and is not zero (in the
                 // standard-form case) or the current lower/upper bound
                 // of the variable (in the non-standard-form case), resp., as
                 // the routines is_solution_feasible() and
                 // is_solution_feasible_for_auxiliary_problem() assume.
	if (is_phaseI) {
	  CGAL_expensive_assertion(
            is_solution_feasible_for_auxiliary_problem());
	} else {
	  CGAL_expensive_assertion(is_solution_feasible());
	}
      else
	vout2 << "(feasibility not checked in intermediate step)" << std::endl;
      CGAL_expensive_assertion(check_tag(Is_nonnegative()) ||
			       r_C.size() == C.size());
    }
	 
}

// update (step 2)
template < typename Q, typename ET, typename Tags >                                         // QP case
void
QP_solver<Q, ET, Tags>::
update_2( Tag_false)
{
    CGAL_qpe_debug {
	vout2 << std::endl
	      << "Update (Step 2)" << std::endl
	      << "---------------" << std::endl;
    }

    // leave variable from basis
    leave_variable();
    CGAL_qpe_debug {
        check_basis_inverse();
    }
    
    // check the updated vectors r_C, r_S_B, r_B_O and w
    CGAL_expensive_assertion(check_r_C(Is_nonnegative()));
    CGAL_expensive_assertion(check_r_S_B(Is_nonnegative()));
    CGAL_expensive_assertion(check_r_B_O(Is_nonnegative()));
    CGAL_expensive_assertion(check_w(Is_nonnegative()));

    // compute current solution
    compute_solution(Is_nonnegative());
}
 
template < typename Q, typename ET, typename Tags >
void
QP_solver<Q, ET, Tags>::
expel_artificial_variables_from_basis( )
{
    int row_ind;
    ET r_A_Cj;
    
    CGAL_qpe_debug {
        vout2 << std::endl
	      << "---------------------------------------------" << std::endl
	      << "Expelling artificial variables from the basis" << std::endl
	      << "---------------------------------------------" << std::endl;
    }
    
    // try to pivot the artificials out of the basis
    // Note that we do not notify the pricing strategy about variables
    // leaving the basis, furthermore the pricing strategy does not
    // know about variables entering the basis.
    // The partial pricing strategies that keep the set of nonbasic vars
    // explicitly are synchronized during transition from phaseI to phaseII 
    for (unsigned int i_ = qp_n + slack_A.size(); i_ < in_B.size(); ++i_) {
        if (is_basic(i_)) { 					// is basic
	    if (has_ineq) {
	        row_ind = in_C[ art_A[i_ - qp_n - slack_A.size()].first];
	    } else {
	        row_ind = art_A[i_ - qp_n].first;
	    }
	    
	    // determine first possible entering variable,
	    // if there is any
	    for (unsigned int j_ = 0; j_ < qp_n + slack_A.size(); ++j_) {
	        if (!is_basic(j_)) {  				// is nonbasic 
		    ratio_test_init__A_Cj( A_Cj.begin(), j_, no_ineq);
		    r_A_Cj = inv_M_B.inv_M_B_row_dot_col(row_ind, A_Cj.begin());
		    if (r_A_Cj != et0) {
		        ratio_test_1__q_x_O(Is_linear());
			i = i_;
			j = j_;
			update_1(Is_linear());
			break;
		    } 
		}
	    }
	}
    }
    if ((art_basic != 0) && no_ineq) {
      // the vector in_C was not used in phase I, but now we remove redundant
      // constraints and switch to has_ineq treatment, hence we need it to
      // be correct at this stage
      for (int i=0; i<qp_m; ++i)
	in_C.push_back(i);
    }
    diagnostics.redundant_equations = (art_basic != 0);

    // now reset the no_ineq and has_ineq flags to match the situation
    no_ineq = no_ineq && !diagnostics.redundant_equations;
    has_ineq = !no_ineq;
    
    // remove the remaining ones with their corresponding equality constraints
    // Note: the special artificial variable can always be driven out of the
    // basis
    for (unsigned int i_ = qp_n + slack_A.size(); i_ < in_B.size(); ++i_) {
        if (in_B[i_] >= 0) {
	    i = i_;
	    CGAL_qpe_debug {
	        vout2 << std::endl
		      << "~~> removing artificial variable " << i
		      << " and its equality constraint" << std::endl
		      << std::endl;
	    }
	    remove_artificial_variable_and_constraint();
	}
    }
}


// replace variable in basis
template < typename Q, typename ET, typename Tags >
void
QP_solver<Q, ET, Tags>::
replace_variable( )
{
    CGAL_qpe_debug {
	vout2 <<   "<--> nonbasic (" << variable_type( j) << ") variable " << j
	      << " replaces basic (" << variable_type( i) << ") variable " << i
	      << std::endl << std::endl;