qp_solution_impl.h

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

    return error ("optimal objective value c^T x + c_0 incorrect");
  return true;
}

// checks whether objective function value is x^TDx + c^T x + c_0
// and fills the vector q with 2Dx
// -------------------------------------------------------------- 
template <typename ET>
template <typename Program>
bool Quadratic_program_solution<ET>::is_value_correct 
(const Program& p, std::vector<ET>& q, Tag_false /*is_linear*/)
{
  ET d = variables_common_denominator();
  if (d <= et0)
    return error ("common variable denominator is negative");
  Variable_numerator_iterator v = variable_numerators_begin();  
  int n = p.get_n();
  add_two_Dz (p, v, q);                                     // 2d * Dx
  ET et2(2);
  ET obj = et2 * d * d * ET(p.get_c0());                    // 2d^2 * c_0
  typedef typename Program::C_iterator C_column_iterator;
  C_column_iterator c = p.get_c();  
  for (int j=0; j<n; ++j, ++v, ++c) {
    obj += et2 * d * ET(*c) * *v;                           // 2d^2 * c^T x
    obj += q[j] * *v;                                       // 2d^2 * x^TDx
  }  
  if (Quotient<ET>(obj, et2*d*d) != objective_value())
    return error ("optimal objective value x^TDx + c^T x + c_0 incorrect");
  return true;
}

// checks whether lambda^TA >= 0
// -----------------------------
template <typename ET>
template <typename Program>
bool Quadratic_program_solution<ET>::is_infeasible_2 
(const Program& p, typename std::vector<ET>& lambda_a, 
 Tag_true /*is_nonnegative*/)
{
  // fill lambda_a
  add_zA (p, infeasibility_certificate_begin(), lambda_a); 
  int n = p.get_n();
  for (int j=0; j<n; ++j)
    if (lambda_a[j] < et0)
      return error ("(lambda^TA)_j < 0 for some j");
  return true;
}

// checks whether lambda^TA >= 0 (<= 0) if u_j = infty (l_j = -infty)
// ------------------------------------------------------------------
template <typename ET>
template <typename Program>
bool Quadratic_program_solution<ET>::is_infeasible_2 
(const Program& p, typename std::vector<ET>& lambda_a, 
 Tag_false /*is_nonnegative*/)
{
  // fill lambda_a
  add_zA (p, infeasibility_certificate_begin(), lambda_a); 
  int n = p.get_n();
  typedef typename Program::FL_iterator FL_column_iterator;
  typedef typename Program::FU_iterator FU_column_iterator;
  FL_column_iterator fl = p.get_fl();
  FU_column_iterator fu = p.get_fu();
  for (int j=0; j<n; ++j, ++fl, ++fu) {
    if (!*fu && lambda_a[j] < et0)
      return error ("u_j = infty but (lambda^TA)_j < 0");
    if (!*fl && lambda_a[j] > et0)
      return error ("l_j = -infty but (lambda^TA)_j > 0");
  }
  return true;
}

// checks whether lambda^T b < 0
// ----------------------------- 
template <typename ET>
template <typename Program>
bool Quadratic_program_solution<ET>::is_infeasible_3 
(const Program& p, const typename std::vector<ET>& /*lambda_a*/,
 Tag_true /*is_nonnegative*/)
{
  ET lambda_b(0);
  int m = p.get_m();
  typedef typename Program::B_iterator B_column_iterator;
  B_column_iterator b = p.get_b();
  Infeasibility_certificate_iterator inf = infeasibility_certificate_begin();
  for (int i=0; i<m; ++i, ++b, ++inf)
    lambda_b += ET(*b) * *inf;
  if (lambda_b >= et0)
    return error ("lambda_b >= 0");
  return true;
}

// checks whether lambda^T b < \sum_{j: lambda^TA_j < 0}lambda^TA_j u_j + 
//                             \sum_{j: lambda^TA_j > 0}lambda^TA_j l_j
// ---------------------------------------------------------------------- 
template <typename ET>
template <typename Program>
bool Quadratic_program_solution<ET>::is_infeasible_3 
(const Program& p, const typename std::vector<ET>& lambda_a,
 Tag_false /*is_nonnegative*/)
{
  ET lambda_b(0);
  int m = p.get_m(); 
  typedef typename Program::B_iterator B_column_iterator;
  B_column_iterator b = p.get_b();
  Infeasibility_certificate_iterator inf = infeasibility_certificate_begin();
  for (int i=0; i<m; ++i, ++b, ++inf)
    lambda_b += ET(*b) * *inf;

  ET sum(0); 
  int n = p.get_n();  
  typedef typename Program::L_iterator L_column_iterator;  
  typedef typename Program::U_iterator U_column_iterator;  
  L_column_iterator l = p.get_l();
  U_column_iterator u = p.get_u();  
  for (int j=0; j<n; ++j, ++l, ++u) {
    if (lambda_a[j] < et0) sum += lambda_a[j] * ET(*u);
    if (lambda_a[j] > et0) sum += lambda_a[j] * ET(*l);
  }
  // now compare
  if (lambda_b >= sum)
    return error("lambda^T b >= sum");
  return true;
}
 
// checks whether (Aw)_i <= (>=, ==) 0 if i-th constraint is <= (>=, ==)
// ---------------------------------------------------------------------
template <typename ET>
template <typename Program>
bool Quadratic_program_solution<ET>::is_unbounded_1 
(const Program& p)
{
  int m = p.get_m();
  std::vector<ET> aw (m, et0);
  add_Az (p, unboundedness_certificate_begin(), aw);
  typedef typename Program::R_iterator R_column_iterator;
  R_column_iterator r = p.get_r();
  for (int i=0; i<m; ++i, ++r)
    switch (*r) {
    case SMALLER:
      if (aw[i] > et0)
	return error ("i-th constraint <= but (Aw)_i > 0");
      break;
    case EQUAL:
      if (aw[i] != et0)
	return error ("i-th constraint == but (Aw)_i != 0");
      break;
    case LARGER:
      if (aw[i] < et0)
	return error ("i-th constraint >= but (Aw)_i < 0");
      break;
    }
  return true;
}

// checks whether w >= 0
// ---------------------
template <typename ET>
template <typename Program>
bool Quadratic_program_solution<ET>::is_unbounded_2 
(const Program& p, Tag_true /*is_nonnegative*/)
{
  int n = p.get_n();
  Unboundedness_certificate_iterator unb = unboundedness_certificate_begin();
  for (int j=0; j<n; ++j, ++unb)
    if (*unb < et0)
      return error ("some w_j < 0");
  return true;
}

// checks whether w_j >= (<=) 0 if l_j (u_j) is finite
// --------------------------------------------------- 
template <typename ET>
template <typename Program>
bool Quadratic_program_solution<ET>::is_unbounded_2 
(const Program& p, Tag_false /*is_nonnegative*/)
{
  int n = p.get_n();
  typedef typename Program::FL_iterator FL_column_iterator;
  typedef typename Program::FU_iterator FU_column_iterator;
  FL_column_iterator fl = p.get_fl();
  FU_column_iterator fu = p.get_fu();
  Unboundedness_certificate_iterator unb = unboundedness_certificate_begin();
  for (int j=0; j<n; ++j, ++unb, ++fl, ++fu) {
    if (*fl && *unb < et0)
      return error ("some l_j is finite but w_j < 0");
    if (*fu && *unb > et0)
      return error ("some u_j is finite but w_j > 0");
  }  
  return true;
}

// checks whether c^Tw < 0
// -----------------------
template <typename ET>
template <typename Program>
bool Quadratic_program_solution<ET>::is_unbounded_3 
(const Program& p, Tag_true /*is_linear*/)
{
  ET cw(0);
  int n = p.get_n();
  typedef typename Program::C_iterator C_column_iterator;
  C_column_iterator c = p.get_c();
  Unboundedness_certificate_iterator unb = unboundedness_certificate_begin();
  for (int j=0; j<n; ++j, ++c, ++unb) 
    cw += ET(*c) * *unb;
  if (cw >= et0)
    return error("c^Tw >= 0");
  return true;
} 

// checks whether w^TDw = 0 and (c^T + 2x^TD)w < 0
// -----------------------------------------------
template <typename ET>
template <typename Program>
bool Quadratic_program_solution<ET>::is_unbounded_3 
(const Program& p, Tag_false /*is_linear*/)
{
  int n = p.get_n();
  std::vector<ET> two_dw (n, et0);
  add_two_Dz(p, unboundedness_certificate_begin(), two_dw);        // 2Dw
  // check w^TDw = 0
  Unboundedness_certificate_iterator unb = unboundedness_certificate_begin();
  ET wdw(0);
  for (int j=0; j<n; ++j, ++unb) 
    wdw += two_dw[j] * *unb;
  if (wdw != et0)
    return error("w^TDw != 0");
  // check (c^T + 2x^TD)w = c^Tw + 2x^T Dw < 0
  typedef typename Program::C_iterator C_column_iterator;
  C_column_iterator c = p.get_c();
  unb = unboundedness_certificate_begin();
  Variable_numerator_iterator v = variable_numerators_begin();
  ET d = variables_common_denominator();
  if (d <= et0) 
    return error ("common variable denominator is negative");
  ET res(0);
  for (int j=0; j<n; ++j, ++c, ++unb, ++v) 
    res += ET(*c) * *unb * d + two_dw[j] * *v;
  if (res >= 0)
    return error("c^T + 2x^TD)w >= 0");
  return true;
}
  
// computes the product of A with the n-vector given by z
// and adds the result to v; precondition: v has size m
// ------------------------------------------------------
template <typename ET>
template <typename Program, typename Z_iterator >
void Quadratic_program_solution<ET>::add_Az 
(const Program& p, Z_iterator z, typename std::vector<ET>& v)
{
  // iterator types
  typedef typename Program::A_iterator
    A_matrix_iterator;
  typedef typename std::iterator_traits<A_matrix_iterator>::value_type
    A_column_iterator;
  // now compute the product
  A_matrix_iterator a = p.get_a();
  int n = p.get_n();
  for (int j=0; j<n; ++j, ++a, ++z) {
    if (!CGAL::is_zero(*z)) {
      // add A_j * z_j to az
      A_column_iterator a_j = *a;
      for (int i=0; i<p.get_m(); ++i, ++a_j)
  	v[i] += *z * ET(*a_j);
    }
  }
}

// computes the product of 2D with the n-vector given by z
// and adds the result to v; precondition: v has size n
// -------------------------------------------------------
template <typename ET>
template <typename Program, typename Z_iterator >
void Quadratic_program_solution<ET>::add_two_Dz 
(const Program& p, Z_iterator z, typename std::vector<ET>& v)
{
  // we compute the contribution from the enries of D on or below
  // the diagonal rowwise, and the remaining entries columnwise
  typedef typename Program::D_iterator
    D_matrix_iterator;
  typedef typename std::iterator_traits<D_matrix_iterator>::value_type
    D_row_iterator;
  int n = p.get_n();
  // make sure that we only handle the nonzero terms of z
  std::vector<int> z_indices;
  Z_iterator l = z;
  for (int i=0; i<n; ++i, ++l) 
    if (*l != 0) z_indices.push_back(i);

  // the rowwise contribution: on/below the diagonal
  D_matrix_iterator d = p.get_d();
  for (int i=0; i<n; ++i, ++d) {
    // handle row i
    D_row_iterator d_i = *d;  // row i on/below diagonal
    for (std::vector<int>::const_iterator 
	   k = z_indices.begin(); k < z_indices.end(); ++k) {
      int j = *k;
      if (j <= i) v[i] += *(z+j) * ET (*(d_i+j));
    }
  }
  
  // the columnwise contribution: above the diagonal
  d = p.get_d();
  for (std::vector<int>::const_iterator 
	 k = z_indices.begin(); k < z_indices.end(); ++k) {
    int j = *k;
    // handle column j
    D_row_iterator d_j = *(d+j); // column j above diagonal
    for (int i=0; i<j; ++i, ++d_j)
      v[i] += *(z+j) * ET (*d_j);
  }
}


// computes the product of the m-vector given by z with the matrix
// A and adds the result to v; precondition: v has size n 
// ----------------------------------------------------------------
template <typename ET>
template <typename Program, typename Z_iterator >
void  Quadratic_program_solution<ET>::add_zA 
(const Program& p, Z_iterator z, typename std::vector<ET>& v)
{
  typedef typename Program::A_iterator
    A_matrix_iterator;
  typedef typename std::iterator_traits<A_matrix_iterator>::value_type
    A_column_iterator;
  A_matrix_iterator a = p.get_a();
  int n = p.get_n();
  int m = p.get_m();
  // make sure that we only handle the nonzero terms of z
  std::vector<int> z_indices; // indices of nonzero entries
  Z_iterator l = z;
  for (int i=0; i<m; ++i, ++l) 
    if (*l != 0) z_indices.push_back(i);
  // now compute the product columnwise 
  for (int j=0; j<n; ++j, ++a) {
    // add z^T * A_j to v[j] 
    A_column_iterator a_j = *a;
    for (std::vector<int>::const_iterator 
	   k = z_indices.begin(); k < z_indices.end(); ++k)
      v[j] += *(z+*k) * ET(*(a_j+*k));
  }  
}

// adds d c^T to v; precondition: v has length n
// ---------------------------------------------
template <typename ET>
template <typename Program>
void  Quadratic_program_solution<ET>::add_c
(const Program& p, typename std::vector<ET>& v)
{
  int n = p.get_n();
  typedef typename Program::C_iterator C_column_iterator;
  C_column_iterator c = p.get_c();
  ET d = variables_common_denominator();
  for (int j=0; j<n; ++j, ++c)
    v[j] += ET(*c) * d;
}


CGAL_END_NAMESPACE

#endif // CGAL_QP_SOLUTION_IMPL_H