khachiyan_approximation.h

来自「很多二维 三维几何计算算法 C++ 类库」· C头文件 代码 · 共 632 行 · 第 1/2 页

// Copyright (c) 1997-2001  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/Approximate_min_ellipsoid_d/include/CGAL/Approximate_min_ellipsoid_d/Khachiyan_approximation.h $
// $Id: Khachiyan_approximation.h 36397 2007-02-16 15:56:11Z gaertner $
// 
//
// Author(s)     : Kaspar Fischer <[e-mail]>

// Note: whenever a comment refers to "Khachiyan's paper" then the
// paper "Rounding of polytopes in the real number model of
// computation" is ment (Mathematics of Operations Research, Vol. 21,
// No. 2, May 1996).  Nontheless, most comments refer to the
// accompanying documentation sheet (and not to the above paper), see
// the file(s) in documentation/.

#ifndef CGAL_KHACHIYAN_APPROXIMATION_H
#define CGAL_KHACHIYAN_APPROXIMATION_H

#include <CGAL/Interval_arithmetic.h>
#include <CGAL/Approximate_min_ellipsoid_d/Approximate_min_ellipsoid_d_configure.h>

// If CGAL is not being used, we need to define certain things:
#ifndef CGAL_VERSION
  namespace CGAL {
    struct Tag_true {};
    struct Tag_false {};
  }
#endif

namespace CGAL {

  // In order not to pollute the CGAL namespace with implementation
  // fuzz, we put such things into the following attic:
  namespace Appel_impl {

    template<typename FT>
    struct Selector {
      typedef Tag_true Is_exact;
    };
    
    template<>
    struct Selector<float> {
      typedef Tag_false Is_exact;
    };
  
    template<>
    struct Selector<double> {
      typedef Tag_false Is_exact;
    };

  } // namespace Appel_impl

  template<bool Embed,class Traits>
  class Khachiyan_approximation {
  public: // types:
    typedef typename Traits::FT FT;
    typedef typename Traits::ET ET;
    typedef typename Traits::Point Point;
    typedef typename Traits::Cartesian_const_iterator C_it;
    typedef typename Appel_impl::Selector<FT>::Is_exact Exact_flag;
    
  private: // member variables and invariants for them:
    Traits& tco;                    // traits class object
    std::vector<const Point *> P;   // input points
    int n;                          // number of input points, i.e., P.size()

    // This class comes in two flavours:
    //
    //  (i) When Embed is false, the input points are taken to be
    //    ordinary points in R^{d_P}, where d_P is the dimension of the
    //    input points (as obtained by Traits::dimension()).  In this
    //    case the "dimension of the ambient space" (see variable d
    //    below) simply equals d_P.
    //
    //  (ii) When Embed is true, the input points are embedded by
    //    means of p -> (p,1) into R^{d_P+1}.  In this case, the
    //    dimension of the amient space is d_P+1.
    //
    // Notice that (at almost all places in this code) whenever a
    // point "p" is mentioned in a comment then the embedded point is
    // ment in case (ii).
    const int d_P;                  // dimension of the input points
    const int d;                    // dimension of the ambient space

    // An instance of this class will compute (provided the points are
    // not degenerate, see is_degenerate() below) a solution x to
    // Khachiyan's program (D) such that the "relaxed optimality
    // conditions"
    //
    //      p_i^T M(x)^{-1} p_i <= (1+desired_eps) d             (*)
    //
    // hold for all input points p_i.  (Notice that the p_i in this
    // case are the embedded points when Embed is true.)  Khachiyan's
    // lemma then guarantees that the ellipsoid defined by the matrix
    // M(x)^{-1} is a "good" approximation in some sense (for further
    // details on "good" refer to is_valid() below).
    //
    // The current epsilon for which (*) holds is stored in the
    // following variable; the variable is only defined if is_deg is
    // false.
    FT eps;

    // In addition to the current epsilon (see variable eps) above, we
    // also maintain the exact epsilon for which (*) holds; this exact
    // epsilon is eps_exact. (Notice that (*) for the variable eps
    // only holds modulo rounding-errors if FT is an inexact type.)
    // In order not to unnecessarily recompute eps_exact, we keep the
    // flag is_exact_eps_uptodate which is true if and only if the
    // value in eps_exact is correct.
    double eps_exact;
    bool is_exact_eps_uptodate;
    
    // Khachian's algorithm is only guaranteed to work when the input
    // points linearly span the whole space (i.e., if dim(span(P)) =
    // d, or, equivalently, when the smallest enclosing ellipsoid has
    // a volume > 0).  The user of this class must call
    // is_degenerate() which returns false iff dim(span(P)) = d;
    // is_degenerate() simply looks up the following variable:
    bool is_deg;

    #ifdef CGAL_APPEL_ASSERTION_MODE
    // We store the current solution x to program (D) as a n-vector.
    // If is_deg is true, the variable x has no meaning.
    std::vector<FT> x;
    #endif // CGAL_APPEL_ASSERTION_MODE
    
    // The inverse of the (d x d)-matrix M(x) is stored in the
    // variable mi.  The element (M(x)^{-1})_{ij} is stored in
    // mi[i+j*d].  (It would actually be enough to only store the
    // upper half of the matrix since it's symmetric -- but we don't
    // care.)  If is_deg is true, the variable mi is undefined.
    //
    // Todo: optimize by only storing one half of mi.
    std::vector<FT> mi;
    
    // The following variable sum is used to build, during the
    // initialization phase, the (d x d)-matrix M(x).  The variable is
    // only defined if is_deg is true; once is_deg becomes false, the
    // variable sum isn't used any more (except in routines for
    // debugging and statistics).  If is_deg is true then sum represents
    // the martix
    //
    //    sum(p_i p_i^T,i=0..(n-1));
    //
    // the (i,j)th element of which is stored as sum[i+j*d].  (It
    // would actually be enough to only store the upper half of the
    // matrix since it's symmetric -- but we don't care.)
    //
    // Todo: optimize by only storing one half of sum.
    std::vector<FT> sum;

    // Khachiyan's algorithm heavily relies on the excess ex[i] of
    // input point P[i] w.r.t. the current solution x (see routine
    // excess() below), i in {0,...,n-1}.  For the first iteration, we
    // compute these numbers explicitly, which is quite costly, namely
    // O(d^2) per excess, for a total of O(n d^2).  However, in all
    // remaining iterations we only update the excesses which can be
    // done in O(n d).  For this we need the excesses of the previous
    // round, and that's why we store them.  Thus, if is_deg is false
    // then ex[i] is the excess of point P[i] w.r.t. the current
    // solution x:
    std::vector<FT> ex;

    // We remember the index (into ex) of a largest element in ex: If
    // is_deg is false then
    //
    //    ex_max = argmax_{0 <= i < n} ex[i].
    //
    int ex_max; 

    // We sometimes need temporary storage.  We allocate this once at
    // instance construction time:
    mutable std::vector<FT> tmp;
    mutable std::vector<FT> t;

  public: // member variables for statistics

    #ifdef CGAL_APPEL_STATS_MODE
    // The value max_error_m is the maximal absolute value in an
    // entry of (M(x) M(x)^{-1} - I), where I is the identity matrix.
    // Of course, in theory, M(x) M(x)^{-1} - I should be the zero
    // matrix, but if FT is an inexact type (i.e., double), then this
    // need not be the case anymore, due to rounding errors.
    //
    // The variable max_error_m_all is the maximum over all
    // max_error values ever encountered during the computation.
    FT max_error_m, max_error_m_all;
    #endif // CGAL_APPEL_STATS_MODE

  private: // internal helper routines:
    
    template<typename NumberType, typename InputIterator>
    NumberType excess(InputIterator p)
      // Computes p^T M(x)^{-1} p, where x is the current solution.
      // (When Embed is true, then the p in the previous sentence's
      // formula is an embedded point...)
      //
      // Complexity: O(d^2)
      //
      // Note: this routine is parametrized by a number-type argument
      // because we will use it in exact_epsilon() and is_valid()
      // below with the exact number-type ET (which may be different
      // from the possibly inexact number-type FT).
    {
      typedef NumberType NT;
      NT result(0);
      InputIterator qi(p);
      
      for (int i=0; i<d; ++i) {

	// compute i-th entry of the vector M(x)^{-1} p into tmp:
	NT tmp(0);
	InputIterator q(p);
	for (int j=0; j<d_P; ++j, ++q)
	  tmp += NT(*q) * NT(mi[i+j*d]);
	if (Embed)
	  tmp += NT(mi[i+d_P*d]);

	// add tmp*p_i to result:
	if (!Embed || i < d_P) {
	  result += NT(*qi) * NT(tmp);
	  ++qi;
	} else
	  result += NT(tmp);
      }

      return result;
    }

    void update_sum(int start)
      // Adds to the variable sum the matrices P[i] P[i]^T where i
      // ranges from start to n-1.
      //
      // Complexity: O((n-start) d^2)
      //
      // Todo: maybe use something like Kahan Summation here? See
      // <http://en.wikipedia.org/wiki/Kahan_Summation_Algorithm>.
    {
      for (int k=start; k<n; ++k) {
	C_it pi = tco.cartesian_begin(*P[k]);
	for (int i=0; i<d_P; ++i, ++pi) {
	  C_it pj = tco.cartesian_begin(*P[k]);
	  for (int j=0; j<d_P; ++j, ++pj)
	    sum[i+j*d] += (*pi) * (*pj);
	  if (Embed)
	    sum[i+d_P*d] += (*pi);
	}

	if (Embed) {
	  C_it pj = tco.cartesian_begin(*P[k]);
	  for (int j=0; j<d_P; ++j, ++pj)
	    sum[d_P+j*d] += (*pj);
	  sum[d_P+d_P*d] += 1;
	}
      }
    }

  public: // construction & destruction:

    Khachiyan_approximation(int dim,int n_est,Traits& tco) :
      // An instance of this class represents, for some given point
      // set P, a solution x to program (D) which satisfies the
      // relaxed optimality conditions (*) for some desired_epsilon
      // (to be specified by the user at a later point).  After
      // calling this constructor, the instance will represent the
      // (degenerate) solution to program (D) for P={}.  By adding
      // points to P (see routine add()), you can compute an
      // approximate solution to (D) for any point set P.
      //
      // (A solution does not always exist; the set P must linearly
      // span the whole space in order for a solution to exist.  See
      // routine is_degenerate() and the result of routine add() for
      // more information.)
      //
      // The number n_est is an estimate on the number of points you
      // are going to add; n_est need not coincide with the actual
      // number of points you will eventually have added.  It merely
      // gives a hint on how much storage is needed.
      tco(tco), n(0),
      d_P(dim), d(Embed? d_P+1 : d_P), is_deg(true),
      #ifdef CGAL_APPEL_ASSERTION_MODE 
      x(n_est),
      #endif // CGAL_APPEL_ASSERTION_MODE
      mi(d*d), sum(d*d), ex(n_est), tmp(d), t(d*d)
    {
      CGAL_APPEL_LOG("appel","Entering Khachiyan with d=" << d << " (" <<
		(Embed? "" : "not ") << "embedded)." << std::endl);
      CGAL_APPEL_TIMER_START("khachiyan");

      // In order to satisfy the invariant on m, we have to initalize
      // m with the zero matrix:
      for (int i=0; i<d; ++i)
	for (int j=0; j<d; ++j)
	  sum[i+j*d] = FT(0);
    }
    
    ~Khachiyan_approximation();

    template<typename InputIterator>
    bool add(InputIterator first,InputIterator last,double desired_eps)
      // Adds the points from the range [first,last) to the instance's
      // set P and computes a solution x to program (D) satisfying the
      // relaxed optimality conditions (*) for the given value
      // desired_eps.  (Khachiyan's lemma then guarantees that the
      // ellipsoid defined by the matrix M(x)^{-1} is a "good"