qp_solver.h

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

  bool  check_basis_inverse( Tag_false is_linear);

  // diagnostic output
  void  print_program ( ) const;
  void  print_basis   ( ) const;
  void  print_solution( ) const;
  void  print_ratio_1_original(int k, const ET& x_k, const ET& q_k);
  void  print_ratio_1_slack(int k, const ET& x_k, const ET& q_k);

  const char*  variable_type( int k) const;
    
  // ensure container size
  template <class Container>
  void ensure_size(Container& c, typename Container::size_type desired_size) {
    typedef typename Container::value_type Value_type;
    for (typename Container::size_type i=c.size(); i < desired_size; ++i) {
      c.push_back(Value_type());
    }
  }
    
private:

private:  // (inefficient) access to bounds of variables:
  // Given an index of an original or slack variable, returns whether
  // or not the variable has a finite lower bound.
  bool has_finite_lower_bound(int i) const;

  // Given an index of an original or slack variable, returns whether
  // or not the variable has a finite upper bound.
  bool has_finite_upper_bound(int i) const;

  // Given an index of an original or slack variable, returns its
  // lower bound.
  ET lower_bound(int i) const;

  // Given an index of an original variable, returns its upper bound.
  ET upper_bound(int i) const;

  struct Bnd { // (inefficient) utility class representing a possibly
    // infinite bound
    enum Kind { MINUS_INF=-1, FINITE=0, PLUS_INF=1 };
    const Kind kind;      // whether the bound is finite or not
    const ET value;       // bound's value in case it is finite

    Bnd(bool is_upper, bool is_finite, const ET& value) 
      : kind(is_upper? (is_finite? FINITE : PLUS_INF) :
	     (is_finite? FINITE : MINUS_INF)),
	value(value) {}
    Bnd(Kind kind, const ET& value) : kind(kind), value(value) {}
    
    bool operator==(const ET& v) const { return kind == FINITE && value == v; }
    bool operator==(const Bnd& b) const {
      return kind == b.kind && (kind != FINITE || value == b.value);
    }
    bool operator!=(const Bnd& b) const { return !(*this == b); }
    bool operator<(const ET& v) const { return kind == FINITE && value < v; }
    bool operator<(const Bnd& b) const {
      return kind < b.kind ||
	(kind == b.kind && kind == FINITE && value < b.value);
    }
    bool operator<=(const Bnd& b) const { return *this < b || *this == b; }
    bool operator>(const ET& v) const { return kind == FINITE && value > v; }
    bool operator>(const Bnd& b)  const { return !(*this <= b); }
    bool operator>=(const Bnd& b) const { return !(*this < b); }
    
    Bnd operator*(const ET& f) const { return Bnd(kind, value*f); }
  };

  // Given an index of an original, slack, or artificial variable,
  // return its lower bound.
  Bnd lower_bnd(int i) const;

  // Given an index of an original, slack, or artificial variable,
  // return its upper bound.
  Bnd upper_bnd(int i) const;

private:
  bool is_value_correct() const;
  // ----------------------------------------------------------------------------

  // ===============================
  // class implementation (template)
  // ===============================

public:

  // pricing
  // -------
  // The solver provides three methods to compute mu_j; the first
  // two below take additional information (which the pricing
  // strategy either provides in exact- or NT-form), and the third
  // simply does the exact computation. (Note: internally, we use
  // the third version, too, see ratio_test_1__t_j().)

  // computation of mu_j with standard form
  template < class RndAccIt1, class RndAccIt2, class NT >  
  NT
  mu_j( int j, RndAccIt1 lambda_it, RndAccIt2 x_it, const NT& dd) const
  {
    NT  mu_j;

    if ( j < qp_n) {                                // original variable

      // [c_j +] A_Cj^T * lambda_C
      mu_j = ( is_phaseI ? NT( 0) : dd * NT(qp_c[ j]));
      mu_j__linear_part( mu_j, j, lambda_it, no_ineq);

      // ... + 2 D_Bj^T * x_B
      mu_j__quadratic_part( mu_j, j, x_it, Is_linear());

    } else {                                        // slack or artificial

      mu_j__slack_or_artificial( mu_j, j, lambda_it, dd,
				 no_ineq);

    }

    return mu_j;
  }
    
  // computation of mu_j with upper bounding
  template < class RndAccIt1, class RndAccIt2, class NT >  
  NT
  mu_j( int j, RndAccIt1 lambda_it, RndAccIt2 x_it, const NT& w_j,
	const NT& dd) const
  {
    NT  mu_j;

    if ( j < qp_n) {                                // original variable

      // [c_j +] A_Cj^T * lambda_C
      mu_j = ( is_phaseI ? NT( 0) : dd * NT(qp_c[ j]));
      mu_j__linear_part( mu_j, j, lambda_it, no_ineq);

      // ... + 2 D_Bj^T * x_B + 2 D_Nj x_N
      mu_j__quadratic_part( mu_j, j, x_it, w_j, dd, Is_linear());

    } else {                                        // slack or artificial

      mu_j__slack_or_artificial( mu_j, j, lambda_it, dd,
				 no_ineq);

    }

    return mu_j;
  }

  // computation of mu_j (exact, both for upper bounding and standard form)
  ET
  mu_j( int j) const
  {
    CGAL_qpe_assertion(!is_basic(j));
    
    if (!check_tag(Is_nonnegative()) &&
	!check_tag(Is_linear()) &&
	!is_phaseI && is_original(j)) {
      return mu_j(j,
		  lambda.begin(),
		  basic_original_variables_numerator_begin(),
		  w_j_numerator(j),
		  variables_common_denominator());
    } else {
      return mu_j(j,
		  lambda.begin(),
		  basic_original_variables_numerator_begin(),
		  variables_common_denominator());
    }
  }

private:

  // pricing (private helper functions)
  // ----------------------------------
  template < class NT, class It > inline                      // no ineq.
  void
  mu_j__linear_part( NT& mu_j, int j, It lambda_it, Tag_true) const
  {
    mu_j += inv_M_B.inner_product_l( lambda_it, *(qp_A+ j));
  }

  template < class NT, class It > inline                      // has ineq.
  void
  mu_j__linear_part( NT& mu_j, int j, It lambda_it, Tag_false) const
  {
    mu_j += inv_M_B.inner_product_l
      ( lambda_it,
	A_by_index_iterator( C.begin(),
			     A_by_index_accessor( *(qp_A + j))));
  }

  template < class NT, class It > inline                     
  void
  mu_j__linear_part( NT& mu_j, int j, It lambda_it, 
		     bool has_no_inequalities) const {
    if (has_no_inequalities) 
      mu_j__linear_part (mu_j, j, lambda_it, Tag_true());
    else
      mu_j__linear_part (mu_j, j, lambda_it, Tag_false());     
  }

  template < class NT, class It > inline          // LP case, standard form
  void
  mu_j__quadratic_part( NT&, int, It, Tag_true) const
  {
    // nop
  }
    
  template < class NT, class It > inline          // LP case, upper bounded
  void
  mu_j__quadratic_part( NT&, int, It, const NT& /*w_j*/, const NT& /*dd*/,
			Tag_true) const
  {
    // nop
  }    

  template < class NT, class It > inline          // QP case, standard form
  void
  mu_j__quadratic_part( NT& mu_j, int j, It x_it, Tag_false) const
  {
    if ( is_phaseII) {
      // 2 D_Bj^T * x_B
      mu_j += inv_M_B.inner_product_x
	( x_it,
	  D_pairwise_iterator_input_type( B_O.begin(),
			       D_pairwise_accessor_input_type(qp_D, j)));
    }
  }

  template < class NT, class It > inline          // QP case, upper bounded
  void
  mu_j__quadratic_part( NT& mu_j, int j, It x_it, const NT& w_j,
			const NT& dd, Tag_false) const
  {
    if ( is_phaseII) {
      mu_j += dd * w_j;
      // 2 D_Bj^T * x_B
      mu_j += inv_M_B.inner_product_x
	( x_it,
	  D_pairwise_iterator_input_type( B_O.begin(),
			       D_pairwise_accessor_input_type(qp_D, j)));
    }
  }


  template < class NT, class It >  inline                     // no ineq.
  void
  mu_j__slack_or_artificial( NT& mu_j, int j, It lambda_it, const NT& dd, Tag_true) const
  {
    j -= qp_n;
    // artificial variable
    // A_j^T * lambda
    mu_j = lambda_it[ j];
    if ( art_A[ j].second) mu_j = -mu_j;

    // c_j + ...
    mu_j += dd*NT(aux_c[ j]);

  }

  template < class NT, class It >  inline                     // has ineq.
  void
  mu_j__slack_or_artificial( NT& mu_j, int j, It lambda_it, const NT& dd, Tag_false) const
  {
    j -= qp_n;

    if ( j < (int)slack_A.size()) {                 // slack variable

      // A_Cj^T * lambda_C
      mu_j = lambda_it[ in_C[ slack_A[ j].first]];
      if ( slack_A[ j].second) mu_j = -mu_j;

    } else {                                        // artificial variable
      j -= slack_A.size();

      // A_Cj^T * lambda_C
      mu_j = lambda_it[ in_C[ art_A[ j].first]];
      if ( art_A[ j].second) mu_j = -mu_j;

      // c_j + ...
      mu_j += dd*NT(aux_c[ j]);
    }
  }

  template < class NT, class It >  inline
  void
  mu_j__slack_or_artificial( NT& mu_j, int j, It lambda_it, 
			     const NT& dd, bool has_no_inequalities) const {
    if (has_no_inequalities)
      mu_j__slack_or_artificial (mu_j, j, lambda_it, dd, Tag_true());
    else
      mu_j__slack_or_artificial (mu_j, j, lambda_it, dd, Tag_false());
  }
};

// ----------------------------------------------------------------------------

// =============================
// class implementation (inline)
// =============================

// initialization
// --------------

// transition
// ----------
template < class Q, typename ET, typename Tags >  inline                                 // QP case
void  QP_solver<Q, ET, Tags>::
transition( Tag_false)
{
  typedef  Creator_2< D_iterator, int, 
    D_pairwise_accessor >  D_transition_creator_accessor;

  typedef  Creator_2< Index_iterator, D_pairwise_accessor,
    D_pairwise_iterator >  D_transition_creator_iterator;

  typedef  Join_input_iterator_1< Index_iterator, typename Bind<
    typename Compose< D_transition_creator_iterator,
    Identity< Index_iterator >, typename
    Bind< D_transition_creator_accessor, D_iterator, 1 >::Type >::Type,
    Index_iterator, 1>::Type >
    twice_D_transition_iterator;
    
  // initialization of vector w and vector r_B_O:
  if (!check_tag(Is_nonnegative())) {
    init_w();                      
    init_r_B_O();
  }

  inv_M_B.transition( twice_D_transition_iterator( B_O.begin(),
						   bind_1( compose( D_transition_creator_iterator(),
								    Identity<Index_iterator>(), bind_1(
												       D_transition_creator_accessor(), qp_D)), B_O.begin())));
}

template < typename Q, typename ET, typename Tags >  inline                                 // LP case
void  QP_solver<Q, ET, Tags>::
transition( Tag_true)
{
  inv_M_B.transition();
}

// ratio test
// ----------
template < typename Q, typename ET, typename Tags > inline                                  // LP case
void  QP_solver<Q, ET, Tags>::
ratio_test_init__2_D_Bj( Value_iterator, int, Tag_true)
{
  // nop
}

template < typename Q, typename ET, typename Tags > inline                                  // QP case
void  QP_solver<Q, ET, Tags>::
ratio_test_init__2_D_Bj( Value_iterator two_D_Bj_it, int j_, Tag_false)
{
  if ( is_phaseII) {
    ratio_test_init__2_D_Bj( two_D_Bj_it, j_,
			     Tag_false(), no_ineq);
  }
}

template < typename Q, typename ET, typename Tags > inline                                  // QP, no ineq.
void  QP_solver<Q, ET, Tags>::
ratio_test_init__2_D_Bj( Value_iterator two_D_Bj_it, int j_, Tag_false,
			 Tag_true )
{
  // store exact version of `2 D_{B_O,j}'
  D_pairwise_accessor  d_accessor( qp_D, j_);
  std::copy( D_pairwise_iterator( B_O.begin(), d_accessor),
	     D_pairwise_iterator( B_O.end  (), d_accessor),
	     two_D_Bj_it);
}

template < typename Q, typename ET, typename Tags > inline                                  // QP, has ineq
void  QP_solver<Q, ET, Tags>::
ratio_test_init__2_D_Bj( Value_iterator two_D_Bj_it, int j_, Tag_false,
			 Tag_false)
{
  // store exact version of `2 D_{B_O,j}'
  if ( j_ < qp_n) {                               // original variable
    ratio_test_init__2_D_Bj( two_D_Bj_it, j_, Tag_false(), Tag_true());
  } else {                                        // slack variable
    std::fill_n( two_D_Bj_it, B_O.size(), et0);
  }
}

template < typename Q, typename ET, typename Tags >  inline                                 // LP case
void  QP_solver<Q, ET, Tags>::