/**
 * A naive implementation of the Taylor method for ODEs
 * Written on 24 January 2010, Geilo, Norway
 * Author: Daniel Wilczak
 *
 * Remarks: the code is as bad as possible, but it shows that the basic AD
 *          can be implemented within a few minutes.
 *
 *          The code contains some errors that can be easily fixed.
 *          For example, there is no memory management, i.e.
 *          allocated memory is never released.
 *          One needs to implement a virtual destructor.
 *
 * Compile (linux): g++ taylorMethod.cpp -o taylorMethod
 * Run: ./taylorMethod
 *
 * Comments added on 27 January 2010, Oslo Airport
*/
#include <iostream>
#include <cmath>

// Order of the Taylor method. Feel free to change it.
const int maxOrder=20;

// This class is a base class for all the nodes in the DAG
// that represents a vector field.
class T
{
public:

  T()
    : left(0), right(0)
  {}

  T(T* _left, T* _right)
    : left(_left), right(_right)
  {}

  // Virtual function to be overrided in the inheriting classes.
  // By default the function does nothing
  // so the type T will be used to decribe constants or variables.
  virtual void eval(int i){}

  T* left;
  T* right;

  long double coeffs[maxOrder+1];
};

// a node for addition
class TSum : public T
{
public:
  TSum(T* _left, T* _right) : T(_left,_right) {}

  void eval(int i)
  {
    // first we compute i-th Taylor coefficients od the left and right subexpressions
    left->eval(i);
    right->eval(i);

    // then the i-th Taylor coefficient of left + right is just a sum
    this->coeffs[i] = left->coeffs[i] + right->coeffs[i];
  }
};

// see comments for TSum
class TSub : public T
{
public:
  TSub(T* _left, T* _right)
    : T(_left,_right)
  {}

  void eval(int i)
  {
    left->eval(i);
    right->eval(i);
    this->coeffs[i] = left->coeffs[i] - right->coeffs[i];
  }
};

// a node for multiplication
class TMul : public T
{
public:
  TMul(T* _left, T* _right)
    : T(_left,_right)
  {}

  void eval(int i)
  {
    // first we compute i-th Taylor coefficients od the left and right subexpressions
    left->eval(i);
    right->eval(i);

    // the i-th derivative of left*right is a convolution
    this->coeffs[i] = 0.;
    for(int n=0;n<=i;++n)
      this->coeffs[i] += left->coeffs[n] * right->coeffs[i-n];
  }
};

// this node computes -x for a given node x
// no comments necessary
class TUnaryMinus : public T
{
public:
  TUnaryMinus(T* _left)
    : T(_left,0)
  {}

  void eval(int i)
  {
    left->eval(i);
    this->coeffs[i] = -left->coeffs[i];
  }
};

// overloaded operators. The DAG of the vector field will be automatically recorded
// when evaluate an expression

T& operator+(T& left, T& right)
{
  return *(new TSum(&left,&right));
}

T& operator-(T& left,T& right)
{
  return *(new TSub(&left,&right));
}

T& operator*(T& left,T& right)
{
  return *(new TMul(&left,&right));
}

T& operator-(T& left)
{
  return *(new TUnaryMinus(&left));
}

// most important part - iterative computations of Taylor coefficients
// see attached materials in pdf
void computeCoefficients(T& x, T& y, T& Fx, T& Fy)
{
  for(int i=0;i<maxOrder;++i)
  {
    // we compute i-th coefficient of F(x(0))
    Fx.eval(i);
    Fy.eval(i);

    // then (i+1)-th coefficient of the solution is F^[i](0)/(i+1)
    x.coeffs[i+1] = Fx.coeffs[i]/(i+1);
    y.coeffs[i+1] = Fy.coeffs[i]/(i+1);
  }
}

// Given computed Taylor coefficients stored in x and y
// and a chosen time step 'h' of the Taylor method
// we compute the next point on the trajectory.
// We just sum the Taylor polynomial
// by the Horner scheme.
// The result will be stored in 0 order coefficients of x and y
// so it is perfect to start the next step of the Taylor method.

void TaylorStep(T& x, T& y, long double h)
{
  for(int i=maxOrder;i>0;--i)
  {
    x.coeffs[i-1] += x.coeffs[i]*h;
    y.coeffs[i-1] += y.coeffs[i]*h;
  }
}

int main()
{
  T x,y;
  // we initialize the vector field. The DAG will be recorded.
  T& Fx = y*(x*x+y*y);
  T& Fy = -x*(x*x+y*y);

  // we set initial condition for the integration
  x.coeffs[0] = 1;
  y.coeffs[0] = -1;

  // we fix some representable number as a time step
  long double h = 1./32;
  // compute 1000 steps and compare with the explicit solution
  for(int i=0;i<1000;++i)
  {
    computeCoefficients(x,y,Fx,Fy);
    TaylorStep(x,y,h);
    // print points on the trajectory from numerical integration
    std::cout << x.coeffs[0] << ", " << y.coeffs[0] << std::endl;

    // compare x(t) and y(t) from the Taylor method and explicit solution
    std::cout << (x.coeffs[0] - sqrt(2.)*cos(2*(i+1)*h+M_PI_4l)) << std::endl;
    std::cout << (y.coeffs[0] + sqrt(2.)*sin(2*(i+1)*h+M_PI_4l)) << std::endl;
  }
}