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

#include "TRandom3.h"

#include "JMath/JMath.hh"

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


namespace {

  using namespace JPP;

  /**
   * A random object with some arithmetic capabilities.
   */
  struct JObject :
    public JMath<JObject>
  {
  public:
    /**
     * Default constructor.
     */
    JObject() :
      x(0.0),
      y(0.0),
      z(0.0)
    {}


    /**
     * Constructor.
     *
     * \param  __x      x
     * \param  __y      y
     * \param  __z      z
     */
    JObject(const double __x,
	    const double __y,
	    const double __z) :
      x(__x),
      y(__y),
      z(__z)
    {}


    /**
     * Add object.
     *
     * \param  object   object
     * \return          this object
     */
    JObject& add(const JObject& object)
    {
      x += object.x;
      y += object.y;
      z += object.z;

      return *this;
    }


    /**
     * Scale object.
     *
     * \param  factor   multiplication factor
     * \return          this object
     */
    JObject& mul(const double factor)
    {
      x *= factor;
      y *= factor;
      z *= factor;

      return *this;
    }


    /**
     * Write object to output.
     *
     * \param  out           output stream
     * \param  object        object
     * \return               output stream
     */
    friend inline std::ostream& operator<<(std::ostream& out, const JObject& object)
    {
      using namespace std;

      out << showpos << FIXED(5,2) << object.x << ' ' 
	  << showpos << FIXED(5,2) << object.y << ' ' 
	  << showpos << FIXED(5,2) << object.z;

      return out;
    }


    double x;
    double y;
    double z;
  };
}


/**
 * \file
 *
 * Example program to test user class with arithmetic capabilities (based on class JMATH::JMath).
 * \author mdejong
 */
int main(int argc, char**argv)
{
  using namespace std;
  using namespace JPP;

  double alpha;
  double precision;
  int    debug;

  try {

    JParser<> zap("Example program to test user class with arithmetic capabilities.");

    zap['a'] = make_field(alpha);
    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(0);

  const Double_t xmin = -1.0;
  const Double_t xmax = +1.0;
 

  const JObject A(gRandom->Uniform(xmin, xmax),
		  gRandom->Uniform(xmin, xmax),
		  gRandom->Uniform(xmin, xmax));

  const JObject B(gRandom->Uniform(xmin, xmax),
		  gRandom->Uniform(xmin, xmax),
		  gRandom->Uniform(xmin, xmax));

  const JObject C = A + B;
  const JObject D = A * alpha;
  const JObject E = interpolate(A, B, alpha);

  DEBUG("A           " << A << endl);
  DEBUG("B           " << B << endl);
  DEBUG("A + B     = " << C << endl);
  DEBUG("A * alpha = " << D << endl);
  DEBUG("interpolate " << E << endl);

  ASSERT(fabs(C.x - (A.x + B.x)) <= precision, "test of addition");
  ASSERT(fabs(C.y - (A.y + B.y)) <= precision, "test of addition");
  ASSERT(fabs(C.z - (A.z + B.z)) <= precision, "test of addition");

  ASSERT(fabs(D.x - (A.x * alpha)) <= precision, "test of multiplication");
  ASSERT(fabs(D.y - (A.y * alpha)) <= precision, "test of multiplication");
  ASSERT(fabs(D.z - (A.z * alpha)) <= precision, "test of multiplication");

  ASSERT(fabs(E.x - (A.x * (1.0 - alpha) + B.x * alpha)) <= precision, "test of interpolation");
  ASSERT(fabs(E.y - (A.y * (1.0 - alpha) + B.y * alpha)) <= precision, "test of interpolation");
  ASSERT(fabs(E.z - (A.z * (1.0 - alpha) + B.z * alpha)) <= precision, "test of interpolation");

  {
    vector<JObject> buffer;

    double x = 0.0;
    double y = 0.0;
    double z = 0.0;

    const int N = 10;

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

      buffer.push_back(JObject(gRandom->Uniform(xmin, xmax),
			       gRandom->Uniform(xmin, xmax),
			       gRandom->Uniform(xmin, xmax)));
		     
      x += buffer.rbegin()->x;
      y += buffer.rbegin()->y;
      z += buffer.rbegin()->z;
    }

    x /= N;
    y /= N;
    z /= N;

    JObject U = getAverage(buffer.begin(), buffer.end());

    DEBUG("U           " << U              << endl);
    DEBUG("U'          " << JObject(x,y,z) << endl);

    ASSERT(fabs(U.x - x) <= precision, "test of averaging");
    ASSERT(fabs(U.y - y) <= precision, "test of averaging");
    ASSERT(fabs(U.z - z) <= precision, "test of averaging");
  }

  {
    double buffer[] = { 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0 };

    ASSERT(fabs(getAverage(buffer) - 5.0) <= precision, "test of averaging");
  }
  
  return 0;
}