qp_solver.h

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

  Values                   x_B_S;     // basic variables (slack)
  Values                   lambda;    // lambda (from KKT conditions)
  Bound_index_values       x_O_v_i;   // bounds value index vector
  // the following vectors are updated
  // with each update in order to avoid
  // evaluating a matrix vector
  // multiplication
  Values                   r_C;       // r_C = A_{C,N_O}x_{N_O}
  // Note: r_C.size() == C.size().

  Values                   r_S_B;     // r_S_B = A_{S_B,N_O}x_{N_O}

  // The following to variables are initialized (if used at all) in
  // transition().  They are not used in case Is_linear or
  // Is_nonnegative is set to Tag_true.
  Values                   r_B_O;     // r_B_O = 2D_{B_O,N_O}x_{N_O}
  Values                   w;         // w = 2D_{O, N_O}x_{N_O}
    
  int                      m_phase;   // phase of the Simplex method
  Quadratic_program_status                   m_status;  // status of last pivot step
  int                      m_pivots;  // number of pivot steps
    
  bool                     is_phaseI; // flag indicating phase I
  bool                     is_phaseII;// flag indicating phase II
  bool                     is_RTS_transition; // flag indicating transition
  // from Ratio Test Step1 to Ratio
  // Test Step2                                           
  const bool               is_LP;     // flag indicating a linear program
  const bool               is_QP;     // flag indicating a quadratic program

  // the following flag indicates whether the program is in equational form
  // AND still has all its equations; this is given in phase I for any
  // program in equational form, but it may change if redundant constraints
  // get removed from the basis. If no_ineq == true, the program is treated
  // in a more efficient manner, since in that case we need no bookkeeping 
  // for basic constraints
  bool                     no_ineq;   
  bool                     has_ineq;  // !no_ineq

  const bool               is_nonnegative; // standard form, from Tag

  // additional variables
  int                      l;         // minimum of 'qp_n+e+1' and 'qp_m'
  // Note: this is an upper bound for
  // the size of the reduced basis in
  // phase I (in phase II, the reduced
  // basis size can be arbitrarily
  // large)
    
  int 		     e;         // number of equality constraints
    
  // Given a variable number i, in_B[i] is -1 iff x_i is not in the current
  // basis.  If the number in_B[i] is >=0, it is the basis heading of x_i.
  Indices                  in_B;      // variable   in basis, -1 if non-basic

  // Given a number i in {0,...,qp_m-1} of a constraint, 
  Indices                  in_C;      // constraint in basis, -1 if non-basic
  // Note: in_C is only maintained if
  // there are inequality constraints.

  Values                   b_C;       // exact version of `qp_b'
  // restricted to basic constraints C
  Values                   minus_c_B; // exact version of `-qp_c'
  // restricted to basic variables B_O
  // Note: minus_c_B is only enlarged,
  // so its size need not be |B|.

  Values                   A_Cj;      // exact version of j-th column of A
  // restricted to basic constraints C
  Values                   two_D_Bj;  // exact version of twice the j-th
  // column of D restricted to B_O
  // Note: tmp_x_2 is only enlarged,
  // so its size need not be |B|.
    
  int                      j;         // index of entering variable `x_j'
    
  int                      i;         // index of leaving variable `x_i'
  ET                       x_i;       // numerator of leaving variable `x_i'
  ET                       q_i;       // corresponding `q_i'
  int                      direction; // indicates whether the current
  // entering variable x_j is increased
  // or decreased
  Bound_index              ratio_test_bound_index;  // indicates for leaving
  // original variables which bound
  // was hit with upper bounding

  ET                       mu;        //   numerator of `t_j'
  ET                       nu;        // denominator of `t_j'
    
  Values                   q_lambda;  // length dependent on C
  Values                   q_x_O;     // used in the ratio test & update
  // Note: q_x_O is only enlarged,
  // so its size need not be |B|.

  Values                   q_x_S;     // 

  Values                   tmp_l;     // temporary vector of size l
  Values                   tmp_x;     // temporary vector of s. >= B_O.size()
  // Note: tmp_x is only enlarged,
  // so its size need not be |B|.

  Values                   tmp_l_2;   // temporary vector of size l
  Values                   tmp_x_2;   // temporary vector of s. >= B_O.size()
  // Note: tmp_x_2 is only enlarged,
  // so its size need not be |B|.
  // Diagnostics
  struct Diagnostics {
    bool redundant_equations;
  };
    
  Diagnostics              diagnostics;

public:

  /*
   * Note: Some member functions below are suffixed with '_'.
   * They are member templates and their declaration is "hidden",
   * because they are also implemented in the class interface.
   * This is a workaround for M$-VC++, which otherwise fails to
   * instantiate them correctly.
   */

  // creation & initialization
  // -------------------------
  // creation
  QP_solver(const Q& qp, 
	    const Quadratic_program_options& options = 
	    Quadratic_program_options());

  virtual ~QP_solver()
  {
    if (strategyP != static_cast<Pricing_strategy*>(0))
      delete strategyP;
  }

	      
private:
  // set-up of QP
  void set( const Q& qp); 
  void set_D (const Q& qp, Tag_true is_linear);
  void set_D (const Q& qp, Tag_false is_linear);
	         
  // set-up of explicit bounds
  void set_explicit_bounds(const Q& qp); 
  void set_explicit_bounds(const Q& qp, Tag_true /*is_nonnegative*/); 
  void set_explicit_bounds(const Q& qp, Tag_false /*is_nonnegative*/);

  // initialization (of phase I)
  void  init( );
    
  // initialization (of phase II)
  /*
    template < class InputIterator >
    void  init( InputIterator  basic_variables_first,
    InputIterator  basic_variables_beyond);
  */
    
  // operations
  // ----------
  // pivot step
  Quadratic_program_status  pivot( )
  { CGAL_qpe_assertion( phase() > 0);
  CGAL_qpe_assertion( phase() < 3);
  pivot_step();
  return status(); }

  // solve QP
  Quadratic_program_status  solve( )
  { CGAL_qpe_assertion( phase() > 0);
  while ( phase() < 3) { pivot_step(); }
  return status(); }

public:

  // access
  // ------
  // access to QP
  int  number_of_variables  ( ) const { return qp_n; }
  int  number_of_constraints( ) const { return qp_m; }
    
  A_iterator  a_begin( ) const { return qp_A;      }
  A_iterator  a_end  ( ) const { return qp_A+qp_n; }
    
  B_iterator  b_begin( ) const { return qp_b;      }
  B_iterator  b_end  ( ) const { return qp_b+qp_m; }
    
  C_iterator  c_begin( ) const { return qp_c;      }
  C_iterator  c_end  ( ) const { return qp_c+qp_n; }

  C_entry     c_0    ( ) const { return qp_c0;}
    
  D_iterator  d_begin( ) const { return qp_D;      }
  D_iterator  d_end  ( ) const { return qp_D+qp_n; }
    
  Row_type_iterator  row_type_begin( ) const { return qp_r;      }
  Row_type_iterator  row_type_end  ( ) const { return qp_r+qp_m; }

  // access to current status
  int     phase     ( ) const { return m_phase;  }
  Quadratic_program_status  status    ( ) const { return m_status; }
  int     iterations( ) const { return m_pivots; }
    
  // access to common denominator
  const ET& variables_common_denominator( ) const 
  { 
    CGAL_qpe_assertion (d > 0);
    return d; 
  }

  // access to current solution
  ET  solution_numerator( ) const;

  // access to current solution
  ET  solution_denominator( ) const { return et2*d*d; }
    
  // access to original variables
  int  number_of_original_variables( ) const { return qp_n; }
    
  // access to slack variables
  int  number_of_slack_variables( ) const { return slack_A.size(); }

  // access to artificial variables
  int  number_of_artificial_variables( ) const { return art_A.size(); }
    
  C_auxiliary_iterator
  c_auxiliary_value_iterator_begin( ) const { return aux_c.begin(); }
  C_auxiliary_iterator
  c_auxiliary_value_iterator_end( ) const {return aux_c.end(); }

  // access to basic variables
  int  number_of_basic_variables( ) const { return B_O.size()+B_S.size(); }
  int  number_of_basic_original_variables( ) const { return B_O.size(); }
  int  number_of_basic_slack_variables( ) const { return B_S.size(); }

  Basic_variable_index_iterator
  basic_original_variable_indices_begin( ) const { return B_O.begin(); }
  Basic_variable_index_iterator
  basic_original_variable_indices_end  ( ) const { return B_O.end(); }
    
  Basic_variable_numerator_iterator
  basic_original_variables_numerator_begin( ) const { return x_B_O.begin(); }
  Basic_variable_numerator_iterator
  basic_original_variables_numerator_end  ( ) const { return x_B_O.begin()
							+ B_O.size(); }

public: // only the pricing strategies (including user-defined ones
        // need access to this) -- make them friends?

  // access to working variables
  int  number_of_working_variables( ) const { return in_B.size(); }
  
  bool is_basic( int j) const
  { 
    CGAL_qpe_assertion(j >= 0);
    CGAL_qpe_assertion(j < number_of_working_variables());
    return (in_B[ j] >= 0);
  }
  
  bool is_original(int j) const
  {
    CGAL_qpe_assertion(j >= 0);
    CGAL_qpe_assertion(j < number_of_working_variables());
    return (j < qp_n);    
  }
    
  bool phaseI( ) const {return is_phaseI;}
  
  bool is_artificial(int k) const;

  int get_l() const;

  // Returns w[j] for an original variable x_j.
  ET w_j_numerator(int j) const
  { 
    CGAL_qpe_assertion((0 <= j) && (j < qp_n) && is_phaseII);
    return w[j];
  }
  
  Bound_index nonbasic_original_variable_bound_index(int i) const
    // Returns on which bound the nonbasic variable x_i is currently
    // sitting:
    //
    // - LOWER: the variable is sitting on its lower bound.
    // - UPPER: the variable is sitting on its upper bound.
    // - FIXED: the variable is sitting on its lower and upper bound.
    // - ZERO: the variable has value zero and is sitting on its lower
    //   bound, its upper bound, or betweeen the two bounds.
    //
    // Note: in the latter case you can call state_of_zero_nonbasic_variable()
    // to find out which bound is active, if any.
  {
    CGAL_assertion(!check_tag(Is_nonnegative()) &&
		   !is_basic(i) && i < qp_n);
    if (x_O_v_i[i] == BASIC) {
      CGAL_qpe_assertion(false);
    }
    return x_O_v_i[i];  
  };
  
  int state_of_zero_nonbasic_variable(int i) const
    // Returns -1 if the original variable x_i equals its lower bound,
    // 0 if it lies strictly between its lower and upper bound, and 1 if
    // it coincides with its upper bound.
    // 
    // See also the documentation of nonbasic_original_variable_bound_index()
    // above.
  {
    CGAL_assertion(!check_tag(Is_nonnegative()) &&
		   !is_basic(i) && i < qp_n && x_O_v_i[i] == ZERO);
    if (*(qp_fl+i) && CGAL::is_zero(qp_l[i]))
      return -1;
    if (*(qp_fu+i) && CGAL::is_zero(qp_u[i]))
      return 1;
    return 0;
  }
  
private:
  // miscellaneous
  // -------------
  // setting the pricing strategy:
  void  set_pricing_strategy ( Quadratic_program_pricing_strategy strategy);

  // diagnostic output
  void  set_verbosity( int verbose = 0, std::ostream& stream = std::cout);


public:
  // access to indices of basic constraints
  int  number_of_basic_constraints( ) const { return C.size(); }

  Basic_constraint_index_iterator
  basic_constraint_indices_begin( ) const { return C.begin(); }
  Basic_constraint_index_iterator
  basic_constraint_indices_end  ( ) const { return C.end(); }

  // helper functions
  template < class RndAccIt1, class RndAccIt2, class NT >  
  NT  mu_j_( int j, RndAccIt1 lambda_it, RndAccIt2 x_it, const NT& dd) const;

  ET  dual_variable( int i)
  {
    for ( int j = 0; j < qp_m; ++j) {
      tmp_x[ j] = inv_M_B.entry( j, i);
    }
    return std::inner_product( tmp_x.begin(), tmp_x.begin()+qp_m,
			       minus_c_B.begin(), et0);
  }

public:
  // public access to compressed lambda (used in filtered base)
  Value_const_iterator get_lambda_begin() const
  {
    return lambda.begin();
  }
  Value_const_iterator get_lambda_end() const
  {
    return lambda.begin() + C.size();
  }
  

private:    

  // private member functions
  // ------------------------
  // initialization
  void  init_basis( );
  void  init_basis__slack_variables( int s_i, Tag_true  has_no_inequalities);
  void  init_basis__slack_variables( int s_i, Tag_false has_no_inequalities);
  void  init_basis__slack_variables( int s_i, bool has_no_inequalities) {
    if (has_no_inequalities)
      init_basis__slack_variables (s_i, Tag_true());
    else 
      init_basis__slack_variables (s_i, Tag_false());
  }

  void  init_basis__constraints    ( int s_i, Tag_true  has_no_inequalities);
  void  init_basis__constraints    ( int s_i, Tag_false has_no_inequalities);
  void  init_basis__constraints    ( int s_i, bool has_no_inequalities) {
    if (has_no_inequalities)
      init_basis__constraints (s_i, Tag_true());
    else 
      init_basis__constraints (s_i, Tag_false());
  }

  void  init_x_O_v_i();
  void  init_r_C(Tag_true  /*is_nonnegative*/);