khachiyan_approximation_impl.h

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

    for (int i=0; i<d; ++i)
      for (int j=0; j<d; ++j) {
	
	// compute element (i,j) of the product of m and M(x):
	FT v(0);
	for (int k=0; k<d; ++k)
	  v += t[i+k*d] * mi[k+j*d];
	
	// check it:
	const FT exact(i == j? 1 : 0);
	max_error = (std::max)(max_error,std::abs(v-exact));
      }

    // update statistics:
    #ifdef CGAL_APPEL_STATS_MODE
    max_error_m_all = (std::max)(max_error,max_error_m_all);
    max_error_m = max_error;
    #endif
    CGAL_APPEL_LOG("appel","  The represenation error in m is: " <<
	      to_double(max_error) << (max_error == FT(0)?
              " (zero)" : "") << "." << "\n");
    
    return max_error;
  }
  #endif // CGAL_APPEL_ASSERTION_MODE

  template<bool Embed,class Traits>
  void Khachiyan_approximation<Embed,Traits>::rank_1_update(int k,
							    const FT& tau)
  {
    // check preconditions:
    CGAL_APPEL_ASSERT(!check_tag(Exact_flag()) || tau == eps/((1+eps)*d-1));
    CGAL_APPEL_ASSERT(!check_tag(Exact_flag()) || 
		      excess<ET>(tco.cartesian_begin(*P[k])) == (1+eps)*d);
    
    const FT mu = eps / ((d-1)*(1+eps));
    const FT alpha = 1 + mu;
    const FT beta = mu / (1+eps);

    // compute into tmp the vector M(x)^{-1} p_k:
    for (int i=0; i<d; ++i) { // loop to determine tmp[i]
      tmp[i] = FT(0);
      C_it pk = tco.cartesian_begin(*P[k]);
      for (int j=0; j<d_P; ++j, ++pk)
	tmp[i] += (*pk) * mi[i+j*d];
      if (Embed)
	tmp[i] += mi[i+d_P*d];
    }

    // We need to scale the current matrix m by alpha and add to it
    // the matrix (tmp tmp^T) times -beta:
    for (int i=0; i<d; ++i)
      for (int j=0; j<d; ++j) {
	mi[i+j*d] *= alpha;
	mi[i+j*d] -= beta * tmp[i]*tmp[j];
      }

    // Update ex[k]: We need to scale ex[k] by alpha and subtract from
    // it beta (p_k^T tmp)^2:
    ex_max = 0;
    for (int k=0; k<n; ++k) {
      
      // compute gamma = (p_k^T tmp)^2:
      FT gamma(0);
      C_it pk = tco.cartesian_begin(*P[k]);
      for (int i=0; i<d_P; ++i, ++pk)
	gamma += tmp[i] * (*pk);
      if (Embed)
	gamma += tmp[d_P];
      gamma *= gamma;
      
      // update ex[k]:
      ex[k] *= alpha;
      ex[k] -= beta*gamma;

      // remember the largest so far:
      if (ex[k] > ex[ex_max])
	ex_max = k;
    }

    // check postcondition:
    #ifdef CGAL_APPEL_EXP_ASSERTION_MODE
    representation_error_of_mi();
    #endif // CGAL_APPEL_EXP_ASSERTION_MODE
  }

  template<bool Embed,class Traits>
  bool Khachiyan_approximation<Embed,Traits>::improve(const double desired_eps)
  {
    // Take the smallest eps s.t. (excess(p_m) =) p_m^T M(x)^{-1} p_m <=
    // (1+eps) d holds for all m in {1,...,n}:
    eps = ex[ex_max]/d - 1;
    is_exact_eps_uptodate = false;
    
    CGAL_APPEL_LOG("appel","  The current eps is: " << to_double(eps) << "\n");
    
    // check if we have already reached an acceptable eps:
    if (eps <= desired_eps) // numerics say we're ready to stop...
      if (exact_epsilon() <= desired_eps) // ... and if it's o.k, we stop
        // Note: if FT is inexact, exact_epsilon() may return a
        // negative number here, which we will interpret as the input
        // points being degenerate.
	return true;

    // We optimize along the line
    //
    //   x' = (1 - tau) x + tau e_{ex_max},
    //
    // i.e., we replace our current solution x with x' for the value
    // tau which minimizes the objective function on this line.  It
    // turns out that the minimum is attained at tau = eps/((1+eps)d-1)
    // which then equals tau = eps/(largest_excess-1).
    const FT tau = eps / (ex[ex_max] - 1);

    #ifdef CGAL_APPEL_ASSERTION_MODE
    // replace x by x':
    for (int i=0; i<n; ++i)
      x[i] = (1-tau)*x[i];
    x[ex_max] += tau;
    CGAL_APPEL_ASSERT(!check_tag(Exact_flag()) ||
		      std::accumulate(x.begin(),x.end(),FT(0)) == FT(1));
    #endif // CGAL_APPEL_ASSERTION_MODE

    // recompute the inverse m of M(x) (where x is our new x') and
    // update the excesses in the array ex:
    rank_1_update(ex_max,tau);

    return false;
  }

  template<bool Embed,class Traits>
  double Khachiyan_approximation<Embed,Traits>::exact_epsilon(bool recompute)
  {
    CGAL_APPEL_ASSERT(!is_deg);

    // return cached result, if possible:
    if (!recompute && is_exact_eps_uptodate)
      return eps_exact;

    // find max p_i^T M(x)^{-1} p_i:
    ET max = 0;
    for (int i=0; i<n; ++i)
      max = (std::max)(max, excess<ET>(tco.cartesian_begin(*P[i])));

    // compute (using exact arithmetic) epsilon via (*):
    typedef CGAL::Quotient<ET> QET;
    QET eps_e = QET(max,d)-1;
    eps_exact = CGAL::to_interval(eps_e).second;

    // debugging output:
    CGAL_APPEL_LOG("appel",
		   "Exact epsilon is " << eps_e << " (rounded: " <<
		   eps_exact << ")." << "\n");
    
    // check whether eps is negative (which under exact arithmetic is
    // not possible, and which we will take as a sign that the input
    // points are degenerate):
    if (CGAL::is_negative(eps_e)) {
      CGAL_APPEL_LOG("appel", "Negative Exact epsilon -> degenerate!" << "\n");
      is_deg = true;
    }

    is_exact_eps_uptodate = true;
    return eps_exact;
  }

  template<bool Embed,class Traits>
  bool Khachiyan_approximation<Embed,Traits>::is_valid(bool verbose)
  {
    // debugging output:
    if (verbose) {
      CGAL_APPEL_IF_STATS(
        std::cout << "The overall error in the matrix inverse is:   " 
	          << max_error_m_all << "." << "\n");
    }

    // handle degenerate case:
    if (is_deg)
      return true;

    // check Eq. (*) for the exact epsilon:
    const double epsilon = exact_epsilon(true);
    const ET     ratio   = (ET(epsilon)+1)*d;
    for (int i=0; i<n; ++i)
      if (excess<ET>(tco.cartesian_begin(*P[i])) > ratio) {
	CGAL_APPEL_LOG("appel","ERROR: Eq. (*) violated." << "\n");
	return false;
      }

    return true;
  }

  template<bool Embed,class Traits>
  template<typename Iterator>
  bool Khachiyan_approximation<Embed,Traits>::
  compute_inverse_of_submatrix(Iterator inverse)
  {
    CGAL_APPEL_ASSERT(!is_deg);

    // copy matrix to destination:
    for (int i=0; i<d-1; ++i)
      for (int j=0; j<d-1; ++j)
	inverse[i+j*(d-1)] = mi[i+j*d];

    // solve in place:
    return Appel_impl::pd_matrix_inverse<FT>(d-1, inverse,
					     inverse, tmp.begin());
  }
  
  template<bool Embed,class Traits>
  void Khachiyan_approximation<Embed,Traits>::print(std::ostream& o)
  {
    #ifdef CGAL_APPEL_ASSERTION_MODE
    o << "xsol := [ ";
    for (int k=0; k<n; ++k) {
      o << x[k];
      if (k<n-1)
	o << ",";
    }
    o << "];" << "\n\n";
    #endif // CGAL_APPEL_ASSERTION_MODE

    o << "Mi:= matrix([" << "\n";
    for (int i=0; i<d; ++i) {
      o << " [ ";
      for (int j=0; j<d; ++j) {
	o << mi[i+j*d];
	if (j<d-1)
	  o << ",";
      }
      o << "]";
      if (i<d-1)
	o << ",";
      o << "\n";
    }
    o << "]);" << "\n";

    o << "S:= matrix([" << "\n";
    for (int i=0; i<d; ++i) {
      o << " [ ";
      for (int j=0; j<d; ++j) {
	o << sum[i+j*d];
	if (j<d-1)
	  o << ",";
      }
      o << "]";
      if (i<d-1)
	o << ",";
      o << "\n";
    }
    o << "]);" << "\n";

    o << "M:= matrix([" << "\n";
    for (int i=0; i<d; ++i) {
      o << " [ ";
      for (int j=0; j<d; ++j) {
	o << t[i+j*d];
	if (j<d-1)
	  o << ",";
      }
      o << "]";
      if (i<d-1)
	o << ",";
      o << "\n";
    }
    o << "]);" << "\n";
  }
  
  template<bool Embed,class Traits>
  int Khachiyan_approximation<Embed,Traits>::write_eps() const
  {
    CGAL_APPEL_ASSERT(d == 2 && !Embed);
    static int id = 0;

    namespace Impl = Approximate_min_ellipsoid_d_impl;
    // (Using "Impl::tostr(id)" in the following doesn't work on GCC 2.95.2.)
    Impl::Eps_export_2 epsf(Approximate_min_ellipsoid_d_impl::tostr(id)+".eps",
                            100.0);
    epsf.set_label_mode(Impl::Eps_export_2::None);

    // output the input points:
    for (int k=0; k<n; ++k) {
      C_it pk = tco.cartesian_begin(*P[k]);
      const double u = to_double(*pk++);
      const double v = to_double(*pk);
      epsf.write_circle(u,v,0.0);
    }

    // output the inscribed ellipse:
    epsf.set_stroke_mode(Impl::Eps_export_2::Dashed);
    epsf.write_cs_ellipse(to_double(mi[0+0*d]),
			  to_double(mi[1+1*d]),
			  to_double(mi[1+0*d]));

    // Output the approximation of the minellipse: Since the relaxed
    // optimality conditions (*) hold, 
    //
    //      p_i^T M(x)^{-1} p_i <= (1+eps) d,
    // 
    // we can divide by a:=(1+eps)d to get
    //
    //     p_i^T M(x)^{-1}/a p_i <= 1,
    //
    // which means that all points lie in the ellipsoid defined by the
    // matrix M(x)^{-1}/alpha.
    const FT a = (1+eps)*d;
    for (int k=0; k<n; ++k) {
      // compute M(x)^{-1}/alpha p_k into vector tmp:
      for (int i=0; i<d; ++i) { // loop to determine tmp[i]
	tmp[i] = FT(0);
	C_it pk = tco.cartesian_begin(*P[k]);
	for (int j=0; j<d; ++j, ++pk)
	  tmp[i] += (*pk) * mi[i+j*d];
      }

      // calculate p^T tmp:
      FT result(0);
      C_it pk = tco.cartesian_begin(*P[k]);
      for (int i=0; i<d; ++i, ++pk)
	result += (*pk) * tmp[i];
    }

    epsf.set_stroke_mode(Impl::Eps_export_2::Dashed);
    epsf.set_label("E2");
    epsf.write_cs_ellipse(to_double(mi[0+0*d]/a),
			  to_double(mi[1+1*d]/a),
			  to_double(mi[1+0*d]/a));
    
    return id++;
  }

}