#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include <limits>

#include "TRandom3.h"

#include "JGeometry2D/JVector2D.hh"
#include "JGeometry2D/JCircle2D.hh"
#include "JGeometry2D/JGeometry2DTestkit.hh"

#include "Jeep/JPrint.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"


/**
 * \file
 *
 * Example program to test enclosing circle.
 * \author mdejong
 */
int main(int argc, char* argv[])
{
  using namespace std;
  using namespace JPP;

  int     numberOfEvents;
  int     numberOfPoints;
  UInt_t  seed;
  double  precision;
  int     debug;

  try {

    JParser<> zap("Example program to test enclosing circle.");

    zap['n'] = make_field(numberOfEvents)      = 100;
    zap['N'] = make_field(numberOfPoints)      =  10;
    zap['S'] = make_field(seed)                =  0;
    zap['e'] = make_field(precision)           =  1.0e-10;
    zap['d'] = make_field(debug)               =  3;

    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }

  gRandom->SetSeed(seed);


  ASSERT(numberOfPoints >= 3);

  const int number_of_duplicates = 20;

  for (int count = 0; count != numberOfEvents; ++count) {

    STATUS(setw(6) << count << '\r'); DEBUG(endl);

    vector<JVector2D> buffer;

    for (int i = 0; i != numberOfPoints; ++i) {

      const JPosition2D pos = getRandom<JPosition2D>();

      for (int n = 0; n != number_of_duplicates; ++n) {
	buffer.push_back(pos);
      }
    }

    const JCircle2D c0(buffer.begin(), buffer.end());


    JCircle2D c1(JPosition2D(0.0, 0.0), numeric_limits<double>::max());
    
    for (vector<JVector2D>::const_iterator i = buffer.begin(); i != buffer.end(); ++i) {
      for (vector<JVector2D>::const_iterator j = i; ++j != buffer.end(); ) {

	JCircle2D c2(*i, *j);

	bool status = true;

	for (vector<JVector2D>::const_iterator p = buffer.begin(); p != buffer.end() && status; ++p) {
	  status = c2.is_inside(*p, numeric_limits<double>::epsilon());
	}

	if (status && c2.getRadius() < c1.getRadius()) {
	  c1 = c2;
	}
      }
    }

    for (vector<JVector2D>::const_iterator i = buffer.begin(); i != buffer.end(); ++i) {
      for (vector<JVector2D>::const_iterator j = i; ++j != buffer.end(); ) {
	for (vector<JVector2D>::const_iterator k = j; ++k != buffer.end(); ) {

	  JCircle2D c3(*i, *j, *k);

	  bool status = true;

	  for (vector<JVector2D>::const_iterator p = buffer.begin(); p != buffer.end() && status; ++p) {
	    status = c3.is_inside(*p, numeric_limits<double>::epsilon());
	  }

	  if (status && c3.getRadius() < c1.getRadius()) {
	    c1 = c3;
	  }
	}
      }
    }

    DEBUG("Smallest    circle " << c0 << endl);
    DEBUG("Alternative circle " << c1 << endl);

    int number_of_outliers = 0;

    for (vector<JVector2D>::const_iterator i = buffer.begin(); i != buffer.end(); ++i) {
      if (c0.getDistance(*i) > c0.getRadius() + precision) {
	++number_of_outliers;
      }
    }

    ASSERT(number_of_outliers == 0);

    ASSERT(c0.getRadius() <= c1.getRadius() + precision);
    //ASSERT(c0.getRadius() >= c1.getRadius() - precision);
  }
  STATUS(endl);

  return 0;
}