khachiyan_approximation_impl.h

来自「很多二维 三维几何计算算法 C++ 类库」· C头文件 代码 · 共 666 行 · 第 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_impl.h $
// $Id: Khachiyan_approximation_impl.h 37145 2007-03-16 08:27:32Z afabri $
// 
//
// 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/.

#include <numeric>

#include <CGAL/tags.h>
#include <CGAL/Approximate_min_ellipsoid_d/Khachiyan_approximation.h>

namespace CGAL {

  namespace Appel_impl {
    
    // Computes the inverse of the positive definite (dxd)-matrix A
    // into Ai, by using the Cholesky decomposition of A. All three
    // iterator template parameters must be random access iterators of
    // value type FT. The iterator tmp must have d entries.  The
    // routine returns true in case no errors occurred; false is
    // returned in case of stability problems or when A is not
    // positive definite.
    //
    // Note: A is destroyed in this process.
    //
    // Precondition: A and Ai may point to the same matrix (i.e., might alias).
    template<typename FT,
	     typename Tmp_iterator,
	     typename A_iterator,
	     typename A_inverse_iterator>
    bool pd_matrix_inverse(const int d,
			   A_iterator A,
			   A_inverse_iterator Ai,
			   Tmp_iterator tmp)
    {
      // I use the following version (in matlab notation) from Walter
      // Gander's lecture "Ausgleichsrechnung" (ETH Zurich, around 2000, see
      // http://www.inf.ethz.ch/personal/chinella/education/ausgleich/
      // Demo2.txt):
      //
      //    # compute lower-triangular L s.t. A = LL^T:
      //    for j=1:d
      //      v = A(j:d,j) - L(j:d,1:j-1) * L(j,1:j-1)';
      //      L(j:d,j) = v/sqrt(v(1));
      //    end;
      //
      // Observe that the vector v in this pseudo-code is a
      // (d-j+1)-vector; we use (the last d-j+1 entries of) the vector
      // tmp to store v.  (Also, the above program uses one-based
      // counting, the code below is of course zero-based.)  Note also
      // that instead of computing the result into L we can as well
      // overwrite the lower part of A.
      for (int j=0; j<d; ++j) {
	// compute new column (v in above pseudo-code):
	for (int i=j; i<d; ++i) {
	  FT ll(0);
	  for (int k=0; k<j; ++k)
	    ll += A[i+k*d] * A[j+k*d];
	  tmp[i] = A[i+j*d] - ll;
	}

	// check regularity:
	if (tmp[j] <= 0) // todo: epsilon?
	  return false;

	// overwrite column:
	const FT scale = FT(1)/std::sqrt(tmp[j]);
	for (int i=j; i<d; ++i)
	  A[i+j*d] = tmp[i] * scale;
      }

      // Now that we have in the lower triangular part of A the
      // Cholesky decomposition A = LL^T of the original A, we compute
      // the inverse of A see "Numerical Recipes in C", end of Chapter
      // 2.9.
      for (int i=0; i<d; ++i) {
	A[i+i*d] = FT(1)/A[i+i*d];
	for (int j=i+1; j<d; ++j) {
	  FT sum(0);
	  for (int k=i; k<j; ++k)
	    sum -= A[j+k*d] * A[k+i*d];
	  A[j+i*d] = sum/A[j+j*d];
	}
      }
    
      // Finally, we calculate A^{-1} = (L^{-1})^T L^{-1} into Ai:
      for (int i=0; i<d; ++i)
	for (int j=0; j<=i; ++j) {

	  // compute entry (i,j) of A^{-1}:
	  FT sum(0);
	  for (int k=i; k<d; ++k)
	    sum += A[k+i*d] * A[k+j*d];
	  Ai[i+j*d] = sum;

	  // Since A^{-1} is symmetric, we set:
	  Ai[j+i*d] = sum;
	}

      return true;
    }

  } // end of namespace Appel_impl

  template<bool Embed,class Traits>
  Khachiyan_approximation<Embed,Traits>::
  ~Khachiyan_approximation()
  {
    CGAL_APPEL_ASSERT_EXPENSIVE(is_valid(false));
  }

  #ifdef CGAL_APPEL_ASSERTION_MODE
  template<bool Embed,class Traits>
  void Khachiyan_approximation<Embed,Traits>::compute_M_of_x()
  // Computes into t the matrix
  // 
  //    M(x) = sum(x_i p_i p_i^T,i=0..(n-1)).
  //
  // Complexity: O(n d^2)
  //
  // Note: we only use this routine to check assertions.
  {
    // erase m:
    for (int i=0; i<d; ++i)
      for (int j=0; j<d; ++j)
	t[i+j*d] = FT(0);

    // evaluate products:
    for (int k=0; 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)
	  t[i+j*d] += x[k] * (*pi) * (*pj);
	if (Embed)
	  t[i+d_P*d] += x[k] * (*pi);
      }
      if (Embed) {
	C_it pj = tco.cartesian_begin(*P[k]);
	for (int j=0; j<d_P; ++j, ++pj)
	  t[d_P+j*d] += x[k] * (*pj);
	t[d_P+d_P*d] += x[k];
      }
    }
  }
  #endif // CGAL_APPEL_ASSERTION_MODE

  template<bool Embed,class Traits>
  bool Khachiyan_approximation<Embed,Traits>::
    compute_inverse_of_t_into_mi(const Tag_true exact)
  // Note: this routine is not used in CGAL; it turned out that using
  // exact arithmetic, the algorithm is very slow.
  {
    // We need to compute into mi the inverse of the matrix t.  We use
    // Gauss-Jordan elimination to do this.  So we write t and the
    // identity matrix I as [I|t] and transform this (by means of row
    // operations) into a system of the form [N|I].  Then N is the
    // inverse of t.  Why?  This is like solving a linear system with
    // several right-hand sides simultaneously by Gauss elimination:
    // Since the transformations we apply do not change the solution
    // space of the intermediate systems, we can say: The system t x =
    // e_j has, for any i in {1,...,d}, the same solution space as I x
    // = n_i (with n_i being the i-th colum of N); it follows that
    // x=n_i.

    // store the identity matrix in mi:
    for (int i=0; i<d; ++i)
      for (int j=0; j<d; ++j)
	mi[i+j*d] = FT(i==j? 1 : 0);

    // Since it is not always possible to achieve a final form [*|I]
    // without row exchanges, we try to get the form [*|P(I)] where
    // P(I) stands for a permutation of the rows of the matrix I ---
    // in other words: we want a "1" in every row.  Therefore, we
    // maintain a integer for every row with the number of the column
    // into which we placed a "1", or a -1 in case we haven't yet
    // place a "1" in this row.  Also, we maintain for every column
    // the number of the row into which we placed the one.
    std::vector<int> col_with_one(d,-1);
    std::vector<int> row_with_one(d);

    for (int j=d-1; j>=0; --j) {
      // In this iteration of the loop we try to make the column j of
      // m a zero vector with exactly one 1 in an unused place (i.e.,
      // in a place k for which one_in_row(k) isn't yet true).

      // search for a suitable place k:
      bool found = false;
      int k = d-1;
      for (; k>=0; --k)
	if (!is_zero(t[k+j*d]) && col_with_one[k]<0) {
	  found = true;
	  break;
	}
      if (!found)
	return false;
      col_with_one[k] = j;
      row_with_one[j] = k;

      // zero out the entries above and below entry k:
      for (int i=0; i<d; ++i)
	if (i != k) {
	  // we add l times row k to row i:
	  const FT l = -t[i+j*d]/t[k+j*d];
	  for (int jj=0; jj<d; ++jj)
	    mi[i+jj*d] += l*mi[k+jj*d];
	  for (int jj=0; jj<=j; ++jj)
	    t[i+jj*d]  += l*t[k+jj*d];
	}

      // finally, we scale row k to get a one at (k,j):
      for (int jj=0; jj<d; ++jj)
	mi[k+jj*d] /= t[k+j*d];
      for (int jj=0; jj<=j; ++jj)
	t[k+jj*d] /= t[k+j*d];
    }

    // At this point we know that for any i in {1,...,d} the system m
    // x = e_i has the some solution as P(I) x = n_i.  So x =
    // P(I)^{-1} n_i and it thus suffices to undo the permutation:
    for (int i=0; i<d; ++i)
      if (row_with_one[i] != i) {
	const int repl_row = row_with_one[i];
	const int col = col_with_one[i];
	for (int j=0; j<d; ++j)
	  std::swap(mi[i+j*d],mi[repl_row+j*d]);
	row_with_one[col] = repl_row;
	col_with_one[repl_row] = col;
	row_with_one[i] = col_with_one[i] = i;
    }
    return true;
  }

  template<bool Embed,class Traits>
  bool Khachiyan_approximation<Embed,Traits>::
  compute_inverse_of_t_into_mi(const Tag_false /* exact */)
  {
    // handle the obvious case when the points cannot span \R^d:
    if (P.size() <= static_cast<unsigned int>(d))
      return false;

    return Appel_impl::pd_matrix_inverse<FT>(d,
					     t.begin(),
					     mi.begin(),
					     tmp.begin());
  }

  template<bool Embed,class Traits>
  bool Khachiyan_approximation<Embed,Traits>::
    compute_initial_inverse_from_sum()
  {
    // assertions:
    CGAL_APPEL_ASSERT(is_deg);

    // check number of points:
    if (n == 0)
      return false;

    // When this routine is called, the variable sum contains the
    // matrix sum(p_i p_i^T,i=0...(n-1)), which coincides with n M(x)
    // for x = (1/n,...,1/n).  Our aim is to compute M(x)^{-1} into
    // the variable mi and, if the latter matrix exits, to set x to
    // (1/n,...,1/n).  For this, we first compute M(x) into variable t:
    const FT invOfn = 1/FT(n);
    for (int i=0; i<d; ++i)
      for (int j=0; j<d; ++j)
	t[i+j*d] = sum[i+j*d] * invOfn;

    if (!compute_inverse_of_t_into_mi(Exact_flag()))
      return false;

    #ifdef CGAL_APPEL_ASSERTION_MODE
    // We start with the solution (1/n,...,1/n):
    for (int k=0; k<n; ++k)
      x[k] = invOfn;
    #endif // CGAL_APPEL_ASSERTION_MODE

    // Finally, we compute the excess of P[k] (w.r.t. x) into ex[k]
    // for all k:
    ex_max = 0;
    for (int k=0; k<n; ++k) 
      if ((ex[k] = excess<FT>(tco.cartesian_begin(*P[k]))) > ex[ex_max])
	ex_max = k;
    CGAL_APPEL_LOG("appel","  Largest excess after initialization is " <<
	      to_double(ex[ex_max]) << "." << "\n");
    
    // Accoding to Khachiyam (Lemma 3, eq. (2.20) in "Rounding of
    // polytopes in the real number model of computation"), the
    // following eps makes (*) hold:
    eps = n-1;
    is_exact_eps_uptodate = false;

    return true;
  }

  #ifdef CGAL_APPEL_ASSERTION_MODE
  template<bool Embed,class Traits>
  typename Traits::FT Khachiyan_approximation<Embed,Traits>::
    representation_error_of_mi()
  {
    using std::abs;

    // If the points are degenerate then the inverse of M(x) doesn't
    // exist, and so we exit immediately:
    if (is_deg)
      return FT(0);

    // compute M(x) into the matrix represented by t:
    compute_M_of_x();

    // compute mi times the matrix M(x) (which should give the
    // identity matrix):
    FT max_error(0);