// suppress conversion warnings before other includes
# include <cppad/wno_conversion.hpp>
//

# include <adolc/adouble.h>
# include <adolc/taping.h>
# include <adolc/drivers/drivers.h>

// definitions not in Adolc distribution and required to use CppAD::AD<adouble>
# include <cppad/example/base_adolc.hpp>

# include <cppad/cppad.hpp>
// ==========================================================================
namespace { // BEGIN empty namespace
// define types for each level
typedef adouble           a1type;
typedef CppAD::AD<a1type> a2type;

// -------------------------------------------------------------------------
// class definition for C++ function object that defines ODE
class Ode {
private:
     // copy of a that is set by constructor and used by g(y)
     CPPAD_TESTVECTOR(a1type) a1x_;
public:
     // constructor
     Ode(const CPPAD_TESTVECTOR(a1type)& a1x) : a1x_(a1x)
     { }
     // the function g(y) is evaluated with two levels of taping
     CPPAD_TESTVECTOR(a2type) operator()
     ( const CPPAD_TESTVECTOR(a2type)& a2y) const
     {     size_t n = a2y.size();
          CPPAD_TESTVECTOR(a2type) a2g(n);
          size_t i;
          a2g[0] = a1x_[0];
          for(i = 1; i < n; i++)
               a2g[i] = a1x_[i] * a2y[i-1];

          return a2g;
     }
};

// -------------------------------------------------------------------------
// Routine that uses Taylor's method to solve ordinary differential equaitons
// and allows for algorithmic differentiation of the solution.
CPPAD_TESTVECTOR(a1type) taylor_ode_adolc(
     Ode                            G       ,  // function that defines the ODE
     size_t                         order   ,  // order of Taylor's method used
     size_t                         nstep   ,  // number of steps to take
     const a1type                   &a1dt   ,  // Delta t for each step
     const CPPAD_TESTVECTOR(a1type) &a1y_ini)  // y(t) at the initial time
{
     // some temporary indices
     size_t i, k, ell;

     // number of variables in the ODE
     size_t n = a1y_ini.size();

     // copies of x and g(y) with two levels of taping
     CPPAD_TESTVECTOR(a2type)   a2y(n), Z(n);

     // y, y^{(k)} , z^{(k)}, and y^{(k+1)}
     CPPAD_TESTVECTOR(a1type)  a1y(n), a1y_k(n), a1z_k(n), a1y_kp(n);

     // initialize x
     for(i = 0; i < n; i++)
          a1y[i] = a1y_ini[i];

     // loop with respect to each step of Taylors method
     for(ell = 0; ell < nstep; ell++)
     {     // prepare to compute derivatives using a1type
          for(i = 0; i < n; i++)
               a2y[i] = a1y[i];
          CppAD::Independent(a2y);

          // evaluate ODE using a2type
          Z = G(a2y);

          // define differentiable version of g: X -> Y
          // that computes its derivatives using a1type
          CppAD::ADFun<a1type> a1g(a2y, Z);

          // Use Taylor's method to take a step
          a1y_k            = a1y;     // initialize y^{(k)}
          a1type dt_kp = a1dt;    // initialize dt^(k+1)
          for(k = 0; k <= order; k++)
          {     // evaluate k-th order Taylor coefficient of y
               a1z_k = a1g.Forward(k, a1y_k);

               for(i = 0; i < n; i++)
               {     // convert to (k+1)-Taylor coefficient for x
                    a1y_kp[i] = a1z_k[i] / a1type(k + 1);

                    // add term for to this Taylor coefficient
                    // to solution for y(t, x)
                    a1y[i]    += a1y_kp[i] * dt_kp;
               }
               // next power of t
               dt_kp *= a1dt;
               // next Taylor coefficient
               a1y_k   = a1y_kp;
          }
     }
     return a1y;
}
} // END empty namespace
// ==========================================================================
// Routine that tests algorithmic differentiation of solutions computed
// by the routine taylor_ode.
bool mul_level_adolc_ode(void)
{     bool ok = true;
     double eps = 100. * std::numeric_limits<double>::epsilon();

     // number of components in differential equation
     size_t n = 4;

     // some temporary indices
     size_t i, j;

     // set up for thread_alloc memory allocator
     using CppAD::thread_alloc; // the allocator
     size_t capacity;           // capacity of an allocation

     // the vector x with length n (or greater) in double
     double* x = thread_alloc::create_array<double>(n, capacity);

     // the vector x with length n in a1type
     CPPAD_TESTVECTOR(a1type) a1x(n);
     for(i = 0; i < n; i++)
          a1x[i] = x[i] = double(i + 1);

     // declare the parameters as the independent variable
     int tag = 0;                     // Adolc setup
     int keep = 1;
     trace_on(tag, keep);
     for(i = 0; i < n; i++)
          a1x[i] <<= double(i + 1);  // a1x is independent for adouble type

     // arguments to taylor_ode_adolc
     Ode G(a1x);                // function that defines the ODE
     size_t   order = n;      // order of Taylor's method used
     size_t   nstep = 2;      // number of steps to take
     a1type   a1dt  = 1.;     // Delta t for each step
     // value of y(t, x) at the initial time
     CPPAD_TESTVECTOR(a1type) a1y_ini(n);
     for(i = 0; i < n; i++)
          a1y_ini[i] = 0.;

     // integrate the differential equation
     CPPAD_TESTVECTOR(a1type) a1y_final(n);
     a1y_final = taylor_ode_adolc(G, order, nstep, a1dt, a1y_ini);

     // declare the differentiable fucntion f : x -> y_final
     // (corresponding to the tape of adouble operations)
     double* y_final = thread_alloc::create_array<double>(n, capacity);
     for(i = 0; i < n; i++)
          a1y_final[i] >>= y_final[i];
     trace_off();

     // check function values
     double check = 1.;
     double t     = nstep * a1dt.value();
     for(i = 0; i < n; i++)
     {     check *= x[i] * t / double(i + 1);
          ok &= CppAD::NearEqual(y_final[i], check, eps, eps);
     }

     // memory where Jacobian will be returned
     double* jac_ = thread_alloc::create_array<double>(n * n, capacity);
     double** jac = thread_alloc::create_array<double*>(n, capacity);
     for(i = 0; i < n; i++)
          jac[i] = jac_ + i * n;

     // evaluate Jacobian of h at a
     size_t m = n;              // # dependent variables
     jacobian(tag, int(m), int(n), x, jac);

     // check Jacobian
     for(i = 0; i < n; i++)
     {     for(j = 0; j < n; j++)
          {     if( i < j )
                    check = 0.;
               else     check = y_final[i] / x[j];
               ok &= CppAD::NearEqual(jac[i][j], check, eps, eps);
          }
     }

     // make memroy avaiable for other use by this thread
     thread_alloc::delete_array(x);
     thread_alloc::delete_array(y_final);
     thread_alloc::delete_array(jac_);
     thread_alloc::delete_array(jac);
     return ok;
}