#include <iostream>

#include "Sundance.h"

#define SHOW_OBJ_FUNC_ERR     // Uncomment to show this error

int main( int argc, char* argv[] )
{

  using Sundance::List;

  const int
    u_basis_dim = 1,
    a_basis_dim = 1;

  double  left    = 0.0;
  double  right   = 1.0;
  int     n       = 2;
  double  uGuess  = 1.0;
  double  aGuess  = 2.0;
  double  uRight  = 2.0*log(cosh(right/sqrt(2.0)));
  double  obj_wgt = 0.9;

  try {

    Sundance::init(&argc,(void***)&argv);

    MeshGenerator mesher = new PartitionedLineMesher(left, right, n);
    Mesh mesh_ = mesher.getMesh();

    TSFVectorType vst = new PetraVectorType();

    Expr coord_x_ =  new CoordExpr(0);

    CellSet boundary_ = new BoundaryCellSet();
    CellSet left_     = boundary_.subset( fabs(coord_x_ - left) < 1.0e-10  );
    CellSet right_    = boundary_.subset( fabs(coord_x_ - right) < 1.0e-10 );

    TSFVectorSpace discreteStateSpace
      = new SundanceVectorSpace(mesh_, new Lagrange(u_basis_dim), vst);
    TSFVectorSpace discreteControlSpace 
      = new SundanceVectorSpace(mesh_, new Lagrange(a_basis_dim), (left_+right_), vst); // was boundry_

    Expr u0 = new DiscreteFunction(discreteStateSpace, uGuess, "u0");
    Expr alpha0 = new DiscreteFunction(discreteControlSpace, aGuess, "alpha0");

    // Write the base point
    MatlabWriter("u_dump.dat").writeField(u0.name(),u0);
    MatlabWriter("alpha_dump.dat",left_+right_).writeField(alpha0.name(),alpha0);
	
    const Expr u0u = new UnknownFunction(new Lagrange(u_basis_dim), "u");
    const Expr a0u = new UnknownFunction(new Lagrange(a_basis_dim), "a");
    const Expr obj_func_expr
      = Integral(right_,0.5*obj_wgt*(u0u-uRight)*(u0u-uRight))
      + Integral(left_+right_,0.5*(1.0 - obj_wgt)*a0u*a0u);

    //		DiscreteFunction::paranoidChecking() = true;
    // Evaluate the objective
    const double obj_eval = obj_func_expr.evaluateFunctional(
                                                             mesh_, List(u0u,a0u),List(u0,alpha0)
                                                             );

    std::cerr << "obj_func_expr.evaluateFunctional(...) = " << obj_eval << std::endl;

    // Evaluate the gradient
#ifdef SHOW_OBJ_FUNC_ERR
    Expr g_u = obj_func_expr.directSensitivity(u0u,u0);
    MatlabWriter("g_u_dump.dat").writeField(g_u.name(),g_u);
    Expr g_a = obj_func_expr.directSensitivity(a0u,alpha0);
    MatlabWriter("g_alpha_dump.dat",left_+right_).writeField(g_a.name(),g_a);
#endif

    Testing::passFailCheck(__FILE__, 0.0, 1.0 );
    Testing::timeStamp(__FILE__, __DATE__, __TIME__);

    Sundance::finalize();

  }
  catch(exception& e) {
    TSFOut::println(e.what());
    Testing::crash(__FILE__);
    Testing::timeStamp(__FILE__, __DATE__, __TIME__);
  }
  return 0;

}