// Synopsis: Testsuite of package Approximate_min_ellipsoid_d
//
// Revision: $Id$
// $Date$
//
// Author: Kaspar Fischer <fischerk@inf.ethz.ch> (ETH Zurich)

#include <cstdlib>
#include <cmath>
#include <vector>
#include <iostream>
#include <string>

#include <CGAL/Cartesian.h>
#include <CGAL/Cartesian_d.h>
#include <CGAL/point_generators_d.h>
#include <CGAL/MP_Float.h>

#include <CGAL/Approximate_min_ellipsoid_d.h>
#include <CGAL/Approximate_min_ellipsoid_d_traits_2.h>
#include <CGAL/Approximate_min_ellipsoid_d_traits_3.h>
#include <CGAL/Approximate_min_ellipsoid_d_traits_d.h>

template <typename T>
std::string tostr(const T& t)
{
  std::stringstream strm;
  strm << t;
  return strm.str();
}

template <typename T>
T fromstr(const std::string& s)
{
  std::stringstream strm(s);
  T t;
  strm >> t;
  return t;
}

typedef CGAL::MP_Float                                      ET;
typedef CGAL::Cartesian<double>                             K_2;
typedef CGAL::Approximate_min_ellipsoid_d_traits_2<K_2, ET> T_2;

typedef CGAL::Cartesian<double>                             K_3;
typedef CGAL::Approximate_min_ellipsoid_d_traits_3<K_3, ET> T_3;

typedef CGAL::Cartesian_d<double>                           K_d;
typedef CGAL::Approximate_min_ellipsoid_d_traits_d<K_d, ET> T_d;

void check(bool expr, const std::string& msg)
{
  if (!expr) {
    std::cerr << "\n" << msg << ": failed\n";
    std::exit(-1);
  }
}

struct TwoD {};
struct ThreeD {};
struct DD {};

template<typename Kernel,typename Point_list>
void add_random_points(int n, int d, int multiplicity, Point_list& list,
		       const DD&)
  // Adds n random d-dimensional points to the given list.
{
  typedef typename Point_list::value_type Point;
  CGAL::Random_points_in_cube_d<Point> rpg(d,100.0);
  for (int i = 0; i < n; ++i) {
    for (int j = 0; j < multiplicity; ++j)
      list.push_back(*rpg);

    ++rpg;
  }
}

template<typename Kernel,typename Point_list>
void add_random_points(int n, int d, int multiplicity, Point_list& list,
		       const TwoD&)
  // Adds n random d-dimensional points to the given list.
{
  typedef typename Point_list::value_type Point;
  typedef CGAL::Point_d<K_d> Point_d;
  CGAL::Random_points_in_cube_d<Point_d> rpg(d,100.0);
  for (int i = 0; i < n; ++i) {
    for (int j = 0; j < multiplicity; ++j) {
      Point_d p = *rpg;
      list.push_back(Point(*p.cartesian_begin(),*(p.cartesian_begin()+1)));
    }

    ++rpg;
  }
}

template<typename Kernel,typename Point_list>
void add_random_points(int n, int d, int multiplicity, Point_list& list,
		       const ThreeD&)
  // Adds n random d-dimensional points to the given list.
{
  typedef typename Point_list::value_type Point;
  typedef CGAL::Point_d<K_d> Point_d;
  CGAL::Random_points_in_cube_d<Point_d> rpg(d,100.0);
  for (int i = 0; i < n; ++i) {
    for (int j = 0; j < multiplicity; ++j) {
      Point_d p = *rpg;
      list.push_back(Point(*p.cartesian_begin(),
			   *(p.cartesian_begin()+1),
			   *(p.cartesian_begin()+2)));
    }

    ++rpg;
  }
}

template<typename Kernel>
struct is_d_dimensional
{
  typedef DD value;
};

template<>
struct is_d_dimensional<T_2>
{
  typedef TwoD value;
};

template<>
struct is_d_dimensional<T_3>
{
  typedef ThreeD value;
};

int id = 0; // used to count files, see simple_test() below

template<typename Kernel, typename Traits>
void simple_test(int n, int d, int multiplicity, double eps)
  // Computes (1+eps)-approximations of MEL(P) for n random points in R^d.
{
  typedef typename Traits::Point                    Point;
  typedef std::vector<Point>                        Point_list;
  typedef CGAL::Approximate_min_ellipsoid_d<Traits> Mel;

  std::cerr << "n=" << n << ", d=" << d << ", mult=" << multiplicity
	    << ", eps=" << eps;
  
  // generate points
  Point_list P;
  typedef typename is_d_dimensional<Traits>::value is_d;
  add_random_points<Kernel>(n, d, multiplicity, P, is_d());
  std::random_shuffle(P.begin(), P.end());

  // compute minellipsoid:
  Traits tco;
  Mel mel(eps, P.begin(), P.end(), tco);

  // check validity
  check(mel.is_valid(true), "validity");

  // find center:
  if (mel.is_full_dimensional()) {
    mel.center_cartesian_begin(); // (Note: forces center to be computed.)
    if (d == 2 || d == 3)
      mel.axes_lengths_begin();   // (Note: forces axes to be computed.)
  }
  
  // query
  bool is_fd   = mel.is_full_dimensional();
  check(!mel.is_empty() || !is_fd, "empty but full-dimensional");
  if (is_fd) {
    // call all accessors
    mel.achieved_epsilon();
    for (int i=0; i<d; ++i) {
      mel.defining_vector(i);
      for (int j=0; j<d; ++j)
	mel.defining_matrix(i,j);
    }
    mel.defining_scalar();
    std::cerr << ", achieved_eps=" << mel.achieved_epsilon();
    if (mel.achieved_epsilon() > eps)
      std::cerr << " (desired eps not achieved)";

    // write EPS
    if (d == 2) {
      mel.write_eps(tostr(id)+".eps");
      std::cerr << " (file " << id << ")";
      ++id;
    }
  }
  std::cerr << std::endl;
}

void test(int n, int multiplicity)
{
  // 2d
  simple_test<K_d,T_d>(n, 2, multiplicity, 0.5);
  simple_test<K_d,T_d>(n, 2, multiplicity, 0.1);
  simple_test<K_d,T_d>(n, 2, multiplicity, 0.01);
  simple_test<K_d,T_d>(n, 2, multiplicity, 0.001);

  simple_test<K_2,T_2>(n, 2, multiplicity, 0.5);
  simple_test<K_2,T_2>(n, 2, multiplicity, 0.1);
  simple_test<K_2,T_2>(n, 2, multiplicity, 0.01);
  simple_test<K_2,T_2>(n, 2, multiplicity, 0.001);
  
  // 3d
  simple_test<K_d,T_d>(n, 3, multiplicity, 0.5);
  simple_test<K_d,T_d>(n, 3, multiplicity, 0.1);
  simple_test<K_d,T_d>(n, 3, multiplicity, 0.01);
  simple_test<K_d,T_d>(n, 3, multiplicity, 0.001);

  simple_test<K_3,T_3>(n, 3, multiplicity, 0.5);
  simple_test<K_3,T_3>(n, 3, multiplicity, 0.1);
  simple_test<K_3,T_3>(n, 3, multiplicity, 0.01);
  simple_test<K_3,T_3>(n, 3, multiplicity, 0.001);

  // 5d
  simple_test<K_d,T_d>(n, 5, multiplicity, 0.5);
  simple_test<K_d,T_d>(n, 5, multiplicity, 0.1);
  simple_test<K_d,T_d>(n, 5, multiplicity, 0.01);
}

int main()
{
  // multiplicity 1
  test(1,   1);
  test(2,   1);
  test(5, 1);
  test(100, 1);

  // multiplicity 2
  test(1,   2);
  test(2,   2);
  test(100, 2);

  // high multiplicity
  test(2,   2);
  test(10, 10);
  test(100,10);
}