#include <string>
#include <iostream>
#include <iomanip>
#include <initializer_list>

#include "JMath/JMathlib3D.hh"

#include "JLang/JManip.hh"

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


int debug;

namespace {

  const double epsilon    =  1.0e-7;      //!< step size
  const double precision  =  1.0e-4;      //!< precision


  /**
   * Test function.
   *
   * \param  title          title
   * \param  f1             function
   * \param  X              abscissa values
   * \param  Y              abscissa values
   * \param  Z              abscissa values
   */
  template<class JF3_t, class T>
  inline void test(const char* const title,
		   const JF3_t& f1,
		   const std::initializer_list<T>& X,
		   const std::initializer_list<T>& Y,
		   const std::initializer_list<T>& Z)
  {
    using namespace std;
    using namespace JPP;

    // test gradient

    JF3_t fa = f1;
    JF3_t fb = f1;

    for (size_t i = 0; i != JF3_t::parameters.size(); ++i) {

      double JF3_t::* p = JF3_t::parameters[i];

      fa.*p -= 0.5*epsilon;
      fb.*p += 0.5*epsilon;

      for (double x : X) {
	for (double y : Y) {
	  for (double z : Z) {

	    if (debug >= debug_t) {

	      cout << setw(2)    <<  i                                                           << ' '
		   << FIXED(9,5) <<  x                                                           << ' '
		   << FIXED(9,5) <<  y                                                           << ' '
		   << FIXED(9,5) <<  z                                                           << ' '
		   << FIXED(9,5) << (fb(x,y,z) - fa(x,y,z)) / epsilon                            << ' '
		   << FIXED(9,5) <<  f1.getGradient(x,y,z).*p                                    << ' '
		   << FIXED(9,5) << (fb(x,y,z) - fa(x,y,z)) / epsilon - f1.getGradient(x,y,z).*p << endl;
	    }

	    ASSERT(fabs((fb(x,y,z) - fa(x,y,z)) / epsilon  -  f1.getGradient(x,y,z).*p) <= precision, title << "; gradient[" << i << "](" << FIXED(7,3) << x << "," << FIXED(7,3) << "," << FIXED(7,3) << z << ")");
	  }
	}
      }

      fa.*p += 0.5*epsilon;
      fb.*p -= 0.5*epsilon;
    }
  }
}

/**
 * \file
 *
 * Test application for JMathlib3D.hh.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;

  try {

    JParser<> zap;

    zap['d'] = make_field(debug)               = 1;

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

  const auto X = { -3.0, -2.0, -1.0, 0.0, +1.0, +2.0, +3.0 };
  const auto Y = { -3.0, -2.0, -1.0, 0.0, +1.0, +2.0, +3.0 };
  const auto Z = { -3.0, -2.0, -1.0, 0.0, +1.0, +2.0, +3.0 };

  test("polynome3x", JPolynome3X<1, 2>(1.0, 1.0, 1.0), X, Y, Z);
  test("polynome3y", JPolynome3Y<1, 2>(1.0, 1.0, 1.0), X, Y, Z);
  test("polynome3z", JPolynome3Z<1, 2>(1.0, 1.0, 1.0), X, Y, Z);
  test("polynome3x * polynome3y * polynome3z", JPolynome3X<1, 1>(1.0, 1.0) * JPolynome3Y<2, 1>(1.0, 1.0) * JPolynome3Z<3, 1>(1.0, 1.0), X, Y, Z);
  
  test("gauss3x", JGauss3X<1>(0.0, 1.0), X, Y, Z);
  test("gauss3y", JGauss3Y<1>(0.0, 1.0), X, Y, Z);
  test("gauss3z", JGauss3Z<1>(0.0, 1.0), X, Y, Z);
  test("gauss3d", JGauss3D<1>(0.0, 1.0), X, Y, Z);

  test("gauss3x * gaus3y * gauss3z * p0 + p0", JGauss3X<1>(0.0, 1.0) * JGauss3Y<2>(0.0, 1.0) * JGauss3Z<3>(0.0, 1.0) * JP0<1>(10.0)  +  JP0<2>(1.0), X, Y, Z);

  test("(polynome3x * polynome3y * polynome3z)^2", (JPolynome3X<1, 1>(1.0, 1.0) * JPolynome3Y<2, 1>(1.0, 1.0) * JPolynome3Z<3, 1>(1.0, 1.0))^2, X, Y, Z);
}