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

#include<CGAL/QP_functions.h>

CGAL_BEGIN_NAMESPACE

// creation & initialization
// -------------------------
template < typename Q, typename ET, typename Tags >
QP_solver<Q, ET, Tags>::
QP_solver(const Q& qp, const Quadratic_program_options& options)
  : et0(0), et1(1), et2(2),
    strategyP(0),
    inv_M_B(vout4),
    d(inv_M_B.denominator()),
    m_phase(-1), is_phaseI(false), is_phaseII(false),
    is_RTS_transition(false),
    is_LP(check_tag(Is_linear())), is_QP(!is_LP),
    //no_ineq(check_tag(Has_equalities_only_and_full_rank())),
    no_ineq(QP_functions_detail::is_in_equational_form(qp)), 
    // may change after phase I
    has_ineq(!no_ineq),
    is_nonnegative(check_tag(Is_nonnegative()))
{
  // init diagnostics
  diagnostics.redundant_equations = false;

  // initialization as in the standard-form case:
  set_verbosity(options.get_verbosity());
  // only if C_entry is double, we actually get filtered strategies,
  // otherwise we fall back to the respective non-filtered ones
  set_pricing_strategy(options.get_pricing_strategy()); 

  // Note: we first set the bounds and then call set() because set()
  // accesses qp_fl, qp_l, etc.
  set_explicit_bounds(qp);
  set(qp);

  // initialize and solve immediately:
  init();
  solve();
}

// set-up of QP

template < typename Q, typename ET, typename Tags >
void QP_solver<Q, ET, Tags>::
set_D(const Q& /*qp*/, Tag_true /*is_linear*/)
{
  // dummy value, never used
  qp_D = 0;
}

template < typename Q, typename ET, typename Tags >
void QP_solver<Q, ET, Tags>::
set_D(const Q& qp, Tag_false /*is_linear*/)
{
  qp_D = qp.get_d();
}

template < typename Q, typename ET, typename Tags >
void QP_solver<Q, ET, Tags>::
set(const Q& qp)
{
  // assertions:
  CGAL_qpe_assertion(qp.get_n() >= 0);
  CGAL_qpe_assertion(qp.get_m() >= 0); 

  // store QP
  qp_n = qp.get_n(); qp_m = qp.get_m();
  qp_A = qp.get_a(); qp_b = qp.get_b(); qp_c = qp.get_c(); qp_c0 = qp.get_c0(); 

  set_D(qp, Is_linear());
  qp_r = qp.get_r();
  
  // set up slack variables and auxiliary problem
  // --------------------------------------------
  
  // reserve memory for slack and artificial part of `A':
  if (has_ineq) {
    const unsigned int eq = std::count(qp_r, qp_r+qp_m, CGAL::EQUAL);
    slack_A.reserve(qp_m - eq);
    art_A.reserve  (       eq);
    art_s.insert(art_s.end(), qp_m, A_entry(0));
  } else
    art_A.reserve( qp_m);

  // decide on which bound the variables sit initially:
  if (!check_tag(Is_nonnegative()))
    init_x_O_v_i();

  set_up_auxiliary_problem();
    
  e = qp_m-slack_A.size(); // number of equalities
  l = (std::min)(qp_n+e+1, qp_m);  // maximal size of basis in phase I
  
  // diagnostic output:
  CGAL_qpe_debug {
    if (vout.verbose()) {
      if (vout2.verbose()) {
	vout2.out() << "======" << std::endl
		    << "Set-Up" << std::endl
		    << "======" << std::endl;
      }
    }
  }
  vout    << "[ " << (is_LP ? "LP" : "QP")
	  << ", " << qp_n << " variables, " << qp_m << " constraints"
	  << " ]" << std::endl;
  CGAL_qpe_debug {   
      if (vout2.verbose() && (!slack_A.empty())) {
	vout2.out() << " (" << slack_A.size() << " inequalities)";
      }
      if (vout2.verbose()) {
	if (has_ineq)
	  vout2.out() << "flag: has inequalities or rank not full"
		      << std::endl;
	if (vout4.verbose()) print_program();
      }
  }
  
  // set up pricing strategy:
  if (strategyP != static_cast< Pricing_strategy*>(0))
    strategyP->set(*this, vout2);
  
  // set up basis inverse:
  inv_M_B.set(qp_n, qp_m, e);
  
  // set phase:
  m_phase    = 0;
  is_phaseI  = false;
  is_phaseII = false;
}

template < typename Q, typename ET, typename Tags >
void QP_solver<Q, ET, Tags>::
set_explicit_bounds(const Q& qp)
{
  set_explicit_bounds (qp, Is_nonnegative());
}

template < typename Q, typename ET, typename Tags >
void QP_solver<Q, ET, Tags>::
set_explicit_bounds(const Q& /*qp*/, Tag_true) {
  // dummy values, never used
  qp_fl = 0;
  qp_l = 0;
  qp_fu = 0;
  qp_u = 0;
}

template < typename Q, typename ET, typename Tags >
void QP_solver<Q, ET, Tags>::
set_explicit_bounds(const Q& qp, Tag_false) {
  qp_fl = qp.get_fl();
  qp_l = qp.get_l();
  qp_fu = qp.get_fu();
  qp_u = qp.get_u();
}



template < typename Q, typename ET, typename Tags >
void QP_solver<Q, ET, Tags>::
init_x_O_v_i()
{
  // allocate storage:
  x_O_v_i.reserve(qp_n);
  x_O_v_i.resize (qp_n);

  // constants for comparisions:
  const L_entry l0(0);
  const U_entry u0(0);
  
  // our initial solution will have all original variables nonbasic,
  // and so we initialize them to zero (if the bound on the variable
  // allows it), or to the variable's lower or upper bound:
  for (int i = 0; i < qp_n; ++i) {
    CGAL_qpe_assertion( !*(qp_fl+i) || !*(qp_fu+i) || qp_l[i]<=qp_u[i]);

    if (*(qp_fl+i))                    // finite lower bound?
      if (*(qp_fu+i))                  // finite lower and finite upper bound?
	if (qp_l[i] == qp_u[i])        // fixed variable?
	  x_O_v_i[i] = FIXED;
	else                           // finite lower and finite upper?
	  if (qp_l[i] <= l0 && u0 <= qp_u[i])
	    x_O_v_i[i] = ZERO;
	  else
	    x_O_v_i[i] = LOWER;
      else                             // finite lower and infinite upper?
	if (qp_l[i] <= l0)
	  x_O_v_i[i] = ZERO;
	else 
	  x_O_v_i[i] = LOWER;
    else                               // infinite lower bound?
      if (*(qp_fu+i))                  // infinite lower and finite upper?
	if (u0 <= qp_u[i])
	  x_O_v_i[i] = ZERO;
	else
	  x_O_v_i[i] = UPPER;
      else                             // infinite lower and infinite upper?
	x_O_v_i[i] = ZERO;
  }
}

template < typename Q, typename ET, typename Tags >
void QP_solver<Q, ET, Tags>::
set_up_auxiliary_problem()
{
  ET            b_max(et0);
  const C_entry c1(1);
  int           i_max = -1; // i_max-th inequality is the most infeasible one
  int           i_max_absolute = -1; // absolute index of most infeasible ineq

  for (int i = 0; i < qp_m; ++i) {
    // Note: For nonstandard form problems, our initial solution is not the
    // zero vector (but the vector with values original_variable_value(i),
    // 0<=i<qp_n), and therefore, rhs=b-Ax is not simply b as in the standard
    // form case, but Ax_init-b:
    const ET rhs = check_tag(Is_nonnegative())?
      qp_b[i] : ET(qp_b[i]) - multiply__A_ixO(i);

    if (has_ineq && (qp_r[i] != CGAL::EQUAL)) { // inequality constraint, so we
					       // add a slack variable, and (if
					       // needed) a special artificial
      if (qp_r[i] == CGAL::SMALLER) {        // '<='

	// add special artificial ('< -0') in case the inequality is
	// infeasible for our starting point (which is the origin):
	if (rhs < et0) {
	  art_s[i] = -c1;
	  if (-rhs > b_max) {
	    i_max = slack_A.size();
	    i_max_absolute = i;
	    b_max = -rhs;
	  }
	}

	// slack variable:
	slack_A.push_back(std::make_pair(i, false));
      } else {                                 // '>='

	// add special artificial ('> +0') in case the inequality is
	// infeasible for our starting point (which is the origin):
	if (rhs > et0) {
	  art_s[i] = c1;
	  if (rhs > b_max) {
	    i_max = slack_A.size();
	    i_max_absolute = i;
	    b_max = rhs;
	  }
	}
	
	// store slack column
	slack_A.push_back(std::make_pair(i, true));
      }
      
    } else                                     // equality constraint, so we
					       // add an artificial variable
      // (Note: if rhs==et0 here then the artificial variable is (at the
      // moment!) not needed. However, we nonetheless add it, for the following
      // reason. If we did and were given an equality problem with the zero
      // vector as the right-hand side then NO artificials would be added at
      // all; so our initial basis would be empty, something we do not want.)
      art_A.push_back(std::make_pair(i, rhs < et0));
  }

  // Note: in order to make our initial starting point (which is the origin) a
  // feasible point of the auxiliary problem, we need to initialize the
  // special artificial value correctly, namely to
  //
  //   max { |b_i| | i is index of an infeasible inequality constraint }. (C1)
  //
  // The index of this "most infeasible" constraint is, at this point of the
  // code, i_max (or i_max is -1 in which case all inequality constraints are
  // feasible and hence no special artififial column is needed at all).

  // prepare initialization of special artificial column:
  // Note: the real work is done in init_basis() below.
  if (i_max >= 0) {
    art_s_i = i_max;                           // Note: the actual
    art_basic = i_max_absolute;                // initialization of art_s_i
					       // will be done in init_basis()
					       // below. We misuse art_s_i to
					       // remember i_max and art_basic
                                               // to remember i_max_absolute
  } else {                                     // no special art col needed
    art_s_i = -1;
    art_s.clear();
  }
}

// initialization (phase I)
template < typename Q, typename ET, typename Tags >
void QP_solver<Q, ET, Tags>::
init()
{
  CGAL_qpe_debug {
    vout2 << std::endl
	  << "==============" << std::endl
	  << "Initialization" << std::endl
	  << "==============" << std::endl;
              
  }

  // set status:
  m_phase    = 1;
  m_status   = QP_UPDATE;
  m_pivots   = 0;
  is_phaseI  = true;
  is_phaseII = false;

  // initial basis and basis inverse