qp_solver.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.
//
// Licenseges 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://scm.gforge.inria.fr/svn/cgal/branches/CGAL-3.3-branch/QP_solver/include/CGAL/QP_solver/QP_solver.h $
// $Id: QP_solver.h 38453 2007-04-27 00:34:44Z gaertner $
// 
//
// Author(s)     : Kaspar Fischer
//               : Bernd Gaertner <[e-mail]>
//               : Sven Schoenherr 
//               : Franz Wessendorp 

#ifndef CGAL_QP_SOLVER_H
#define CGAL_QP_SOLVER_H

#include <CGAL/iterator.h>
#include <CGAL/functional.h>
#include <CGAL/QP_solver/basic.h>
#include <CGAL/QP_solver/iterator.h>
#include <CGAL/QP_solver/functors.h>
#include <CGAL/QP_options.h>
#include <CGAL/QP_solution.h>
#include <CGAL/QP_solver/QP_basis_inverse.h>
#include <CGAL/QP_solver/QP_pricing_strategy.h>
#include <CGAL/QP_solver/QP_full_exact_pricing.h>
#include <CGAL/QP_solver/QP_partial_exact_pricing.h>
#include <CGAL/QP_solver/QP_full_filtered_pricing.h>
#include <CGAL/QP_solver/QP_partial_filtered_pricing.h>
#include <CGAL/QP_solver/QP_exact_bland_pricing.h>

#include <CGAL/algorithm.h>

#include <CGAL/IO/Verbose_ostream.h>

#include <vector>
#include <numeric>
#include <algorithm>

CGAL_BEGIN_NAMESPACE

// ==================
// class declarations
// ==================
template < typename Q, typename ET, typename Tags  >
class QP_solver;

template <class ET>
class QP_solution; 

namespace QP_solver_impl {   // namespace for implemenation details
  // --------------
  // Tags generator
  // --------------
  template < typename Linear, 
	     typename Nonnegative >
  struct QP_tags {
    typedef Linear                Is_linear;
    typedef Nonnegative           Is_nonnegative;
  };

  template < class Q, class Is_linear >
  struct D_selector {};

  template <class Q>
  struct D_selector<Q, Tag_false> // quadratic
  {
    typedef typename Q::D_iterator D_iterator;
  };

  template <class Q>
  struct D_selector<Q, Tag_true> // linear
  {
    // dummy type, not used
    typedef int** D_iterator;
  };

  template < class Q, class Is_nonnegative >
  struct Bd_selector {};

  template < class Q >
  struct Bd_selector<Q, Tag_false> // nonstandard form
  {
    typedef typename Q::FL_iterator FL_iterator;
    typedef typename Q::L_iterator L_iterator;
    typedef typename Q::FU_iterator FU_iterator;
    typedef typename Q::U_iterator U_iterator;
  };

  template < class Q >
  struct Bd_selector<Q, Tag_true> // standard form
  {
    // dummy types, not used
    typedef int* FL_iterator;
    typedef int* L_iterator;
    typedef int* FU_iterator;
    typedef int* U_iterator;
  };

  // only allow filtered pricing if NT = double
  template <typename Q, typename ET, typename Tags, typename NT>
  struct Filtered_pricing_strategy_selector
  {
    typedef QP_full_exact_pricing<Q, ET, Tags> FF;
    typedef QP_partial_exact_pricing<Q, ET, Tags> PF;
  };

  template <typename Q, typename ET, typename Tags>
  struct Filtered_pricing_strategy_selector<Q, ET, Tags, double> 
  {
    typedef QP_full_filtered_pricing<Q, ET, Tags> FF;
    typedef QP_partial_filtered_pricing<Q, ET, Tags> PF;
  };

} // end of namespace for implementation details

// ================
// class interfaces
// ================


template < typename Q, typename ET, typename Tags >
class QP_solver : public QP_solver_base<ET> {

public: // public types
  typedef  QP_solver<Q, ET, Tags> Self;
  typedef  QP_solver_base<ET> Base;
  
  // types from the QP
  typedef  typename Q::A_iterator   A_iterator;
  typedef  typename Q::B_iterator   B_iterator;
  typedef  typename Q::C_iterator   C_iterator;
  typedef  CGAL::Comparison_result Row_type;
  typedef  typename Q::R_iterator Row_type_iterator;
  
  // the remaining types might not be present in the qp, so the
  // following selectors generate dummy types for them 
  typedef  typename QP_solver_impl::
  D_selector<Q, typename Tags::Is_linear>::
  D_iterator D_iterator;
  typedef typename QP_solver_impl::
  Bd_selector<Q, typename Tags::Is_nonnegative>::
  L_iterator L_iterator;
  typedef typename QP_solver_impl::
  Bd_selector<Q, typename Tags::Is_nonnegative>::
  U_iterator U_iterator;
  typedef typename QP_solver_impl::
  Bd_selector<Q, typename Tags::Is_nonnegative>::
  FL_iterator FL_iterator;
  typedef typename QP_solver_impl::
  Bd_selector<Q, typename Tags::Is_nonnegative>::
  FU_iterator FU_iterator;

  // types from the Tags
  typedef  typename Tags::Is_linear    Is_linear;
  typedef  typename Tags::Is_nonnegative Is_nonnegative;

  // friends
  template <class Q_, class ET_>
  friend bool has_linearly_independent_equations 
  (const Q_& qp, const ET_& dummy);

private: // private types

  // types of original problem:
  typedef  typename std::iterator_traits<A_iterator>::value_type  A_column;
  typedef  typename std::iterator_traits<D_iterator>::value_type  D_row;
  
  typedef  typename std::iterator_traits<A_column  >::value_type  A_entry;
  typedef  typename std::iterator_traits<B_iterator>::value_type  B_entry;
  typedef  typename std::iterator_traits<C_iterator>::value_type  C_entry;
  typedef  typename std::iterator_traits<D_row     >::value_type  D_entry;
  typedef  typename std::iterator_traits<L_iterator>::value_type  L_entry;
  typedef  typename std::iterator_traits<U_iterator>::value_type  U_entry;
  
  // slack columns:
  //
  // The following two types are used to (conceptually) add to the matrix A
  // additional columns that model the constraints "x_s>=0" for the slack
  // variables x_s.  Of course, we do not store the column (which is just
  // plus/minus a unit vector), but maintain a pair (int,bool): the first
  // entry says in which column the +-1 is and the second entry of the pair
  // says whether it is +1 (false) or -1 (true).
  typedef  std::pair<int,bool>        Slack_column;
  typedef  std::vector<Slack_column>  A_slack;

  // artificial columns
  //
  // Artificial columns that are (conceptually) added to the matrix A are
  // handled exactly like slack columns (see above).
  typedef  std::pair<int,bool>        Art_column;
  typedef  std::vector<Art_column>    A_art;

  // special artificial column:
  //
  // Also for the special artificial variable we (conceptually) add a column
  // to A. This column contains only +-1's (but it may contain several nonzero
  // entries).
  typedef  std::vector<A_entry>       S_art;
  
  // auxiliary objective vector (i.e., the objective vector for phase I):
  typedef  std::vector<C_entry>       C_aux;

public: // export some additional types:
  
  typedef  typename Base::Indices     Indices; 
  typedef  typename Base::Index_mutable_iterator   Index_iterator;
  typedef  typename Base::Index_const_iterator     Index_const_iterator;
 
  // For problems in nonstandard form we also export the following type, which
  // for an original variable will say whether it sits at is lower, upper, at
  // its lower and upper (fixed) bound, or at zero, or whether the variable is
  // basic:
  enum  Bound_index  { LOWER, ZERO, UPPER, FIXED, BASIC };

private:
  typedef  std::vector<Bound_index>    Bound_index_values;
  typedef  typename Bound_index_values::iterator
  Bound_index_value_iterator;
  typedef  typename Bound_index_values::const_iterator
  Bound_index_value_const_iterator;
  
  // values (variables' numerators):
  typedef  std::vector<ET>            Values;
  typedef  typename Values::iterator  Value_iterator;
  typedef  typename Values::const_iterator
  Value_const_iterator;
    
  // access values by basic index functor:
  typedef  CGAL::Value_by_basic_index<Value_const_iterator>
  Value_by_basic_index;

  // access to original problem by basic variable/constraint index:
  typedef  QP_vector_accessor<A_column, false, false >  A_by_index_accessor;
  typedef  Join_input_iterator_1< Index_const_iterator, A_by_index_accessor >
  A_by_index_iterator;

  // todo kf: following can be removed once we have all these (outdated)
  // accessors removed:
  typedef  QP_vector_accessor< B_iterator, false, false >
  B_by_index_accessor;
  typedef  Join_input_iterator_1< Index_const_iterator, B_by_index_accessor >
  B_by_index_iterator;

  typedef  QP_vector_accessor< C_iterator, false, false >
  C_by_index_accessor;
  typedef  Join_input_iterator_1< Index_const_iterator, C_by_index_accessor >
  C_by_index_iterator;

  typedef  QP_matrix_accessor< A_iterator, false, true, false, false>
  A_accessor;
  typedef  typename CGAL::Bind< A_accessor,
				typename A_accessor::argument2_type,2>::Type
  A_row_by_index_accessor;
  typedef  Join_input_iterator_1< Index_iterator, A_row_by_index_accessor >
  A_row_by_index_iterator;

  // Access to the matrix D sometimes converts to ET, and 
  // sometimes retruns the original input type
  typedef  QP_matrix_pairwise_accessor< D_iterator, ET >
  D_pairwise_accessor;
  typedef  Join_input_iterator_1< Index_const_iterator,
				  D_pairwise_accessor >
  D_pairwise_iterator;
  typedef  QP_matrix_pairwise_accessor< D_iterator, D_entry >
  D_pairwise_accessor_input_type;
  typedef  Join_input_iterator_1< Index_const_iterator, 
				  D_pairwise_accessor_input_type >
  D_pairwise_iterator_input_type;

  // access to special artificial column by basic constraint index:
  typedef  QP_vector_accessor< typename S_art::const_iterator, false, false>
  S_by_index_accessor;
  typedef  Join_input_iterator_1< Index_iterator, S_by_index_accessor >
  S_by_index_iterator;
  
public:
    
  typedef  typename A_slack::const_iterator
  A_slack_iterator;

  typedef  typename A_art::const_iterator
  A_artificial_iterator;
    
  typedef  typename C_aux::const_iterator
  C_auxiliary_iterator;

  typedef typename Base::Variable_numerator_iterator
  Variable_numerator_iterator;

//   typedef typename Base::Original_index_const_iterator
//   Original_index_const_iterator;

  typedef  Index_const_iterator       Basic_variable_index_iterator;
  typedef  Value_const_iterator       Basic_variable_numerator_iterator;
  typedef  Index_const_iterator       Basic_constraint_index_iterator;
        
  typedef  QP_pricing_strategy<Q, ET, Tags>  Pricing_strategy;

private:
  // compile time tag for symbolic perturbation, should be moved into traits
  // class when symbolic perturbation is to be implemented
  Tag_false                is_perturbed;
    
  // some constants
  const ET                 et0, et1, et2;

  // verbose output streams
  mutable Verbose_ostream  vout;      // used for any  diagnostic output
  mutable Verbose_ostream  vout1;     // used for some diagnostic output
  mutable Verbose_ostream  vout2;     // used for more diagnostic output
  mutable Verbose_ostream  vout3;     // used for full diagnostic output
  mutable Verbose_ostream  vout4;     // used for output of basis inverse
  mutable Verbose_ostream  vout5;     // used for output of validity tests
    
  // pricing strategy
  Pricing_strategy*        strategyP;

  // given QP
  int                      qp_n;      // number of variables
  int                      qp_m;      // number of constraints
    
  // min x^T D x + c^T x + c0
  A_iterator               qp_A;      // constraint matrix
  B_iterator               qp_b;      // right-hand-side vector
  C_iterator               qp_c;      // objective vector
  C_entry                  qp_c0;     // constant term in objective function
  // attention: qp_D represents *twice* the matrix D
  D_iterator               qp_D;      // objective matrix
  Row_type_iterator        qp_r;      // row-types of constraints
  FL_iterator              qp_fl;     // lower bound finiteness vector
  L_iterator               qp_l;      // lower bound vector
  FU_iterator              qp_fu;     // upper bound finiteness vector
  U_iterator               qp_u;      // upper bound vector

  A_slack                  slack_A;   // slack part of constraint matrix

  // auxiliary problem    
  A_art                    art_A;     // artificial part of constraint matrix
  // Note: in phase I there is an
  // additional "fake" column attached
  // to this "matrix", see init_basis()

  S_art                    art_s;     // special artificial column for slacks
  int                      art_s_i;   // art_s_i>=0  -> index of special
  //                artificial column
  // art_s_i==-1 -> no sp. art. col
  // art_s_i==-2 -> sp. art. col removed
  //                after it left basis 
  int                      art_basic; // number of basic artificial variables
  C_aux                    aux_c;     // objective function for phase I
  // initially has the same size as A_art

  Indices                  B_O;       // basis (original variables)
  // Note: the size of B_O is always
  // correct, i.e., equals the number of
  // basic original variables, plus (in
  // phase I) the number of basic
  // artificial variables.

  Indices                  B_S;       // basis (   slack variables)
    
  Indices                  C;         // basic constraints ( C = E+S_N )
  // Note: the size of C is always
  // correct, i.e., corresponds to the
  // size of the (conceptual) set
  // $E\cup S_N$.

  Indices                  S_B;       // nonbasic constraints ( S_B '=' B_S)
    
  QP_basis_inverse<ET,Is_linear>
  inv_M_B;   // inverse of basis matrix

  const ET&                d;         // reference to `inv_M_B.denominator()'
    
  Values                   x_B_O;     // basic variables (original)
  // Note: x_B_O is only enlarged,
  // so its size need not be |B|.