/**
 * A naive implementation of forward propagation of the first order derivatives
 * Written during the lecture at the school of Automatic Differentiation
 * on 25 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.
 *
 * Compile (linux): g++ forwardDiffExample.cpp -o forwardDiffExample
 * Run: ./forwardDiffExample
 *
 * Comments added on 27 January 2010, Oslo Airport
*/

#include <iostream>
#include <cmath>

// maxDim is not a number of variables but rather a number of directional derivatives
// we want to propagate.

const int maxDim = 2;

// this class will store the actual value of an expression
// in the member 'value'
// and directional derivatives with respect to maxDim directions
// in the array 'der'
class F
{
public:
  F(double _value)
    : value(_value)
  {
    for(int i=0;i<maxDim;++i)
      der[i] = 0.;
  }

  // a simple function for printing actual data stored in an object of F
  void show() const
  {
    std::cout << "value=" << value << std::endl << "derivatives:" ;
    for(int i=0;i<maxDim;++i)
      std::cout << " " << der[i];
    std::cout << std::endl;
  }

  double value;
  double der[maxDim];
};

//-------------------------------------------------------
// overloaded operator for addition
F operator+(const F& left, const F& right)
{
  // the value is just a sum of the left and right operands
  F result(left.value + right.value);

  // the same for derivatives
  for(int i=0;i<maxDim;++i)
    result.der[i] = left.der[i] + right.der[i];
  return result;
}

//-------------------------------------------------------
// same comments as for addition
F operator-(const F& left, const F& right)
{
  F result(left.value - right.value);
  for(int i=0;i<maxDim;++i)
    result.der[i] = left.der[i] - right.der[i];
  return result;
}

//-------------------------------------------------------

F operator*(const F& left, const F& right)
{
  // multiplication of values of left and right operands
  F result(left.value * right.value);

  // derivative of multiplication is given by (lr)' = lr' + l'r
  for(int i=0;i<maxDim;++i)
    result.der[i] = left.value*right.der[i] + right.value*left.der[i];
  return result;
}

//-------------------------------------------------------

F operator/(const F& left, const F& right)
{
  // division of actual values of left and right operands
  F result(left.value / right.value);

  // derivative of (l/r)' = (l' - r'* (l/r))/r
  for(int i=0;i<maxDim;++i)
    result.der[i] = (left.der[i] - right.der[i]*result.value)/right.value;
  return result;
}

//-------------------------------------------------------

F sin(const F& arg)
{
  F result(sin(arg.value));

  // we precompute cos(arg.value) . It will be used many times in the loop
  double temp = cos(arg.value);
  for(int i=0;i<maxDim;++i)
    result.der[i] = temp * arg.der[i];
  return result;
}

//-------------------------------------------------------

int main()
{
  // The function from attached pdf materials is
  //    f(x,y) = (x+y+1)/(x*y-1).
  // We compute its gradient for x=2, y=-2.
  // Feel free to overload other elemental functions
  // and test more complicated examples.

  // We replace contants by (c,0). This is highly not optimal.
  // As an excercise please overload operators +,-,*,/
  // when one of the operands is just a constant.

  F x(2.), y(-2), one(1);

  // dx/dx = 1 and dy/dy=1
  x.der[0] = 1;
  y.der[1] = 1;

  std::cout << "Printing x: \n";
  x.show();
  std::cout << "Printing y: \n";
  y.show();

  F f = (x+y+one)/(x*y-one);
  std::cout << "Printing f: \n";
  f.show();
}