qp_solution_impl.h

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

// Copyright (c) 1997-2007  ETH Zurich (Switzerland).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org); you may redistribute it under
// the terms of the Q Public License version 1.0.
// See the file LICENSE.QPL distributed with CGAL.
//
// Licensees holding a valid commercial license may use this file in
// accordance with the commercial license agreement provided with the software.
//
// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
//
// $URL: svn+ssh://[e-mail]/svn/cgal/trunk/QP_solver/include/CGAL/QP_solver/QP_solution_impl.h $
// $Id: QP_solution_impl.h 38416 2007-04-23 09:12:34Z gaertner $
// 
//
// Author(s)     : Bernd Gaertner <[e-mail]>
#ifndef CGAL_QP_SOLUTION_IMPL_H
#define CGAL_QP_SOLUTION_IMPL_H

#include<iterator>

CGAL_BEGIN_NAMESPACE
// checks whether the solution actually solves program p  
// this performs exactly the checks described in the doc
// of the class Quadratic_program_solution
// -----------------------------------------------------
template <typename ET>
template <class Program, typename Is_linear, typename Is_nonnegative>
bool Quadratic_program_solution<ET>::solves_program 
(const Program& p, 
 Is_linear is_linear, Is_nonnegative is_nonnegative)
{    
  CGAL_qpe_assertion_msg(!is_void(), "Solution not initialized");

  // first check whether the dimensions agree
  int n = variable_numerators_end() - variable_numerators_begin();
  if (n != p.get_n()) 
    return error ("wrong number of variables");
  int m = solver()->lambda_end() - solver()->lambda_begin();
  if (m != p.get_m()) 
    return error("wrong number of constraints");

  // now distinguish between the three possible cases
  switch (status()) {
  case QP_OPTIMAL: {
    std::vector<ET> ax_minus_b (m, et0); // d(Ax-b)
    std::vector<ET> two_Dx (n, et0);     // d(2Dx)
    return 
      // the following fills ax_minus_b...
      is_feasible (p, ax_minus_b, is_nonnegative) &&          // feasible?
      // the following fills two_Dx...
      is_value_correct (p, two_Dx, is_linear) &&              // obj. value?
      is_optimal_1 (p) &&                                     // condition 1?
      // ...and the following uses ax_minus_b
      is_optimal_2 (p, ax_minus_b)         &&                 // condition 2?
      // ...and the following uses two_Dx
      is_optimal_3 (p, two_Dx, is_linear, is_nonnegative);    // condition 3?
  }
  case QP_INFEASIBLE: {
    std::vector<ET> lambda_a (n, et0); // lambda^TA
    return 
      is_infeasible_1 (p) &&                                  // condition 1?
      // the following fills lambda_a...
      is_infeasible_2 (p, lambda_a, is_nonnegative) &&        // condition 2?
      // ...and the following uses lambda_a
      is_infeasible_3 (p, lambda_a, is_nonnegative);          // condition 3?
  }
  case QP_UNBOUNDED: {
    std::vector<ET> ax_minus_b (m, et0);
    return 
      is_feasible (p, ax_minus_b, is_nonnegative) &&          // feasible?
      is_unbounded_1 (p) &&                                   // condition 1?
      is_unbounded_2 (p, is_nonnegative) &&                   // condition 2?
      is_unbounded_3 (p, is_linear);                          // condition 3?
  }
  default:
    return error ("solution in undefined state");
  }
}

// tests whether Ax ~ b is satisfied and computes d(Ax-b);
// precondition: ax_minus_b has length m and is zero
// ------------------------------------------------------
template <typename ET>
template <typename Program>
bool Quadratic_program_solution<ET>::are_constraints_feasible 
(const Program& p, 
 typename std::vector<ET>& ax_minus_b)
{
  typedef typename Program::B_iterator B_column_iterator;
  typedef typename Program::R_iterator R_column_iterator;

  // fill ax_minus_b (length m)
  int m = p.get_m();
  B_column_iterator b = p.get_b();
  add_Az (p, variable_numerators_begin(), ax_minus_b);     // d(Ax)
  ET d = variables_common_denominator();
  for (int i=0; i<m; ++i, ++b)
    ax_minus_b[i] -= ET (*b) * d;                          // d(Ax-b)
  
  // now validate constraints 
  if (d <= et0) 
    return error ("common variable denominator is negative");
  R_column_iterator r = p.get_r();
  for (int i=0; i<m; ++i, ++r)
    switch (*r) {
    case SMALLER: // <=
      if (ax_minus_b[i] > et0) return error("inequality (<=) violated");
      break;
    case EQUAL:   // ==
      if (ax_minus_b[i] != et0) return error("inequality (==) violated");
      break;
    case LARGER:  // >=
      if (ax_minus_b[i] < et0) return error("inequality (>=) violated"); 
      break;
    }
  return true;
}

// tests whether Ax ~ b, x >= 0 is satisfied and computes d(Ax-b)
// precondition: ax_minus_b has length m and is zero
// --------------------------------------------------------------
template <typename ET>
template <typename Program>
bool Quadratic_program_solution<ET>::is_feasible 
(const Program& p, 
 typename std::vector<ET>& ax_minus_b,
 Tag_true /*is_nonnegative*/)
{
  // test Ax ~ b
  if (!are_constraints_feasible (p, ax_minus_b)) return false;
  // test x >= 0  
  if (variables_common_denominator() <= et0) 
    return error ("common variable denominator is negative");
  for (Variable_numerator_iterator v = variable_numerators_begin();
       v < variable_numerators_end(); ++v)
    if (*v < et0) return error("bound (>= 0) violated");
  return true;
}

// tests whether Ax ~ b, l <= x <= u is satisfied and computes d(Ax-b)
// precondition: ax_minus_b has length m and is zero
// -------------------------------------------------------------------
template <typename ET>
template <typename Program>
bool Quadratic_program_solution<ET>::is_feasible 
(const Program& p, 
 typename std::vector<ET>& ax_minus_b,
 Tag_false /*is_nonnegative*/)
{
  // test Ax ~ b
  if (!are_constraints_feasible (p, ax_minus_b)) return false; // (*)
  // test l <= x <= u   
  ET d = variables_common_denominator();
  if (d <= et0) 
    return error ("common variable denominator is negative");
  typedef typename Program::FL_iterator FL_column_iterator;
  typedef typename Program::L_iterator L_column_iterator;  
  typedef typename Program::FU_iterator FU_column_iterator;
  typedef typename Program::U_iterator U_column_iterator;
  FL_column_iterator fl = p.get_fl();
  FU_column_iterator fu = p.get_fu();
  L_column_iterator l = p.get_l();
  U_column_iterator u = p.get_u();  
  for (Variable_numerator_iterator v = variable_numerators_begin();
       v < variable_numerators_end(); ++v, ++fl, ++l, ++fu, ++u) {
    if (*fl && (*v < ET(*l) * d)) 
      return error("bound (>=l) violated");
    if (*fu && (*v > ET(*u) * d)) 
      return error("bound (<=u) violated");
  }
  return true;
}

// checks whether constraint type <= (>=) implies lambda >= (<=) 0
// --------------------------------------------------------------- 
template <typename ET>
template <typename Program>
bool Quadratic_program_solution<ET>::is_optimal_1 
(const Program& p)
{
  if (variables_common_denominator() <= et0) 
    return error ("wrong sign of optimality certificate");
  typedef typename Program::R_iterator R_column_iterator;
  int m = p.get_m();
  R_column_iterator r = p.get_r();
  Optimality_certificate_numerator_iterator opt = 
    optimality_certificate_numerators_begin();
  for (int i=0; i<m; ++i, ++r, ++opt) {
    if (*r == SMALLER &&  *opt < et0)
      return error("constraint (<=) with lambda < 0");
    if (*r == LARGER &&  *opt > et0)
      return error("constraint (>=) with lambda > 0");
  }
  return true;
}

// checks whether constraint type <= (>=) implies lambda >= (<=) 0
// --------------------------------------------------------------- 
template <typename ET>
template <typename Program>
bool Quadratic_program_solution<ET>::is_infeasible_1 
(const Program& p)
{ 
  int m = p.get_m();
  typedef typename Program::R_iterator R_column_iterator;
  R_column_iterator r = p.get_r();
  Infeasibility_certificate_iterator inf = infeasibility_certificate_begin();
  for (int i=0; i<m; ++i, ++r, ++inf) {
    if (*r == SMALLER &&  *inf < et0)
      return error("constraint (<=) with lambda < 0");
    if (*r == LARGER &&  *inf > et0)
      return error("constraint (>=) with lambda > 0");
  }
  return true;
}

// checks whether lambda and Ax-b are complementary
// ------------------------------------------------
template <typename ET>
template <typename Program>
bool Quadratic_program_solution<ET>::is_optimal_2 
(const Program& p, 
 const typename std::vector<ET>& ax_minus_b)
{
  int m = p.get_m();  
  Optimality_certificate_numerator_iterator opt = 
    optimality_certificate_numerators_begin();
  for (int i=0; i<m; ++i, ++opt)
    if (*opt != et0 && ax_minus_b[i] != et0) 
      return error("lambda and Ax-b are not complementary");
  return true;					     
}

// checks whether (c^T + lambda^TA)_j >= (==) 0 if x_j = (>) 0
// -----------------------------------------------------------
template <typename ET>
template <typename Program>
bool Quadratic_program_solution<ET>::is_optimal_3 
(const Program& p, std::vector<ET>& q,
 Tag_true /*is_linear*/, Tag_true /*is_nonnegative*/)
{
  int n = p.get_n();                                           // q = 0
  add_c (p, q);                                                // d c^T
  add_zA (p, optimality_certificate_numerators_begin(), q);    // d lambda^T A
  Variable_numerator_iterator v = variable_numerators_begin();
  for (int j=0; j<n; ++j, ++v) {
    if (q[j] < et0) 
      return error("some (c^T + lambda^TA)_j is negative");
    if (*v > et0 && q[j] > et0)
      return error("(c^T + lambda^TA) and x are not complementary");
  }
  return true;
}

// checks whether (c^T + lambda^TA + 2x^TD)_j >= (==) 0 if x_j = (>) 0
// -------------------------------------------------------------------
template <typename ET>
template <typename Program>
bool Quadratic_program_solution<ET>::is_optimal_3 
(const Program& p, std::vector<ET>& q,
 Tag_false /*is_linear*/, Tag_true /*is_nonnegative*/)
{
  int n = p.get_n();                                          // q = 2Dx
  add_c (p, q);                                               // d * c^T
  add_zA (p, optimality_certificate_numerators_begin(), q);   // d * lambda^T A
  Variable_numerator_iterator v = variable_numerators_begin();
  for (int j=0; j<n; ++j, ++v) {
    if (q[j] < et0) 
      return error("some (c^T + lambda^TA + 2Dx)_j is negative");
    if (*v > et0 && q[j] > et0)
      return error("(c^T + lambda^TA + 2Dx) and x are not complementary");
  }
  return true;
}

// checks whether (c^T + lambda^TA )_j >= (==, <=) 0 if 
// x_j = l_j < u_j (l_j < x_j < u_j,  l_j < u_j = x_j)
// ---------------------------------------------------
template <typename ET>
template <typename Program>
bool Quadratic_program_solution<ET>::is_optimal_3 
(const Program& p, std::vector<ET>& q, 
 Tag_true /*is_linear*/, Tag_false /*is_nonnegative*/)
{
  int n = p.get_n();                                          // q = 0
  add_c (p, q);                                               // d * c^T
  add_zA (p, optimality_certificate_numerators_begin(), q);   // d * lambda^T A
  typedef typename Program::FL_iterator FL_column_iterator;
  typedef typename Program::L_iterator L_column_iterator;  
  typedef typename Program::FU_iterator FU_column_iterator;
  typedef typename Program::U_iterator U_column_iterator;
  FL_column_iterator fl = p.get_fl();
  FU_column_iterator fu = p.get_fu();
  L_column_iterator l = p.get_l();
  U_column_iterator u = p.get_u();  
  Variable_numerator_iterator v = variable_numerators_begin();  
  ET d = variables_common_denominator();
  if (d <= et0)
    return error ("common variable denominator is negative");
  for (int j=0; j<n; ++j, ++v, ++fl, ++l, ++fu, ++u) {
    if (*fl && *v == ET(*l) * d && (!*fu || *l < *u) && q[j] < et0) 
      return error("x_j = l_j < u_j but (c^T + lambda^TA )_j < 0");
    if ( (!*fl || *v > ET(*l) * d) && (!*fu || *v < ET(*u) * d) && q[j] != et0)
      return error("l_j < x_j < u_j but (c^T + lambda^TA )_j != 0");
    if (*fu && *v == ET(*u) * d && (!*fl || *l < *u) && q[j] > et0) 
      return error("x_j = u_j > l_j but (c^T + lambda^TA )_j > 0");
  }
  return true;
}

// checks whether (c^T + lambda^TA +2x^*D)_j >= (==, <=) 0 if 
// x_j = l_j < u_j (l_j < x_j < u_j,  l_j < u_j = x_j)
// ----------------------------------------------------------
template <typename ET>
template <typename Program>
bool  Quadratic_program_solution<ET>::is_optimal_3 
(const Program& p, std::vector<ET>& q, 
 Tag_false /*is_linear*/, Tag_false /*is_nonnegative*/)
{
  int n = p.get_n();                                         // q = d * 2Dx
  add_c (p, q);                                              // d * c^T
  add_zA (p, optimality_certificate_numerators_begin(), q);  // d * lambda^T A
  typedef typename Program::FL_iterator FL_column_iterator;
  typedef typename Program::L_iterator L_column_iterator;  
  typedef typename Program::FU_iterator FU_column_iterator;
  typedef typename Program::U_iterator U_column_iterator;
  FL_column_iterator fl = p.get_fl();
  FU_column_iterator fu = p.get_fu();
  L_column_iterator l = p.get_l();
  U_column_iterator u = p.get_u();  
  Variable_numerator_iterator v = variable_numerators_begin();  
  ET d = variables_common_denominator();
  if (d <= et0)
    return error ("common variable denominator is negative");
  for (int j=0; j<n; ++j, ++v, ++fl, ++l, ++fu, ++u) {
    if (*fl && *v == ET(*l) * d && (!*fu || *l < *u) && q[j] < et0) 
      return error("x_j = l_j < u_j but (c^T + lambda^TA + 2Dx)_j < 0");
    if ( (!*fl || *v > ET(*l) * d) && (!*fu || *v < ET(*u) * d) && q[j] != et0)
      return error("l_j < x_j < u_j but (c^T + lambda^TA + 2Dx)_j != 0");
    if (*fu && *v == ET(*u) * d && (!*fl || *l < *u) && q[j] > et0) 
      return error("x_j = u_j > l_j but (c^T + lambda^TA + 2Dx)_j > 0");
  }
  return true;
}

// checks whether objective function value is c^T x + c_0, but
// leaves the vector q alone
// -------------------------------------------------------------- 
template <typename ET>
template <typename Program>
bool Quadratic_program_solution<ET>::is_value_correct 
(const Program& p, std::vector<ET>& /*q*/, Tag_true /*is_linear*/)
{
  // check objective value c^T x + c_0
  ET d = variables_common_denominator();
  if (d <= et0)
    return error ("common variable denominator is negative");
  Variable_numerator_iterator v = variable_numerators_begin();
  ET obj = d  * ET(p.get_c0());                             // d * c_0
  typedef typename Program::C_iterator C_column_iterator;
  C_column_iterator c = p.get_c();
  int n = p.get_n();
  for (int j=0; j<n; ++j, ++v, ++c)
    obj += ET(*c) * *v;                                     // d * c^T x
  if (Quotient<ET>(obj, d) != objective_value())