/* nag::ad::e01bf Tangent Over Adjoint Example Program.
 */

#include <dco.hpp>
#include <iostream>
#include <nagad.h>

// Function which calls NAG AD routines.
template <typename T> void func(const T &xint, T &fint);

// Driver with the adjoint calls.
//  Evaluates a piecewise cubic Hermite interpolant and computes function value
//  at point xint.
// Also computes the 2nd derivative d^2(fint)/d(xint)^2.
void driver(const double &xint, double &fint, double &d);

int main()
{
  std::cout
      << "nag::ad::e01bf Tangent Over Adjoint Example Program Results\n\n";

  // Point to interpolate at
  double xint = 8.4;
  // Interpolated function value
  double fint;

  // Derivative
  double d;

  // Call driver
  driver(xint, fint, d);

  // Print outputs
  std::cout << " Derivatives calculated: Second order adjoints\n";
  std::cout << " Computational mode    : algorithmic\n";

  std::cout.setf(std::ios::scientific, std::ios::floatfield);
  std::cout.precision(6);
  std::cout << "\n fint = " << fint << std::endl;

  // Print derivative
  std::cout << "\n 2nd Derivative of fint w.r.t. interpolation point xint :\n";
  std::cout << " d^2(fint) / d(xint)^2 = " << d << std::endl;

  return 0;
}

// Driver with the adjoint calls.
//  Evaluates a piecewise cubic Hermite interpolant and computes function value
//  at point xint.
// Also computes the 2nd derivative d^2(fint)/d(xint)^2.
void driver(const double &xintv, double &fintv, double &d)
{
  using mode = dco::ga1s<dco::gt1s<double>::type>;
  using T    = mode::type;

  // Create the AD tape
  mode::global_tape = mode::tape_t::create();

  // Variable to differentiate w.r.t.
  T xint                            = xintv;
  dco::derivative(dco::value(xint)) = 1.0;

  // Register variable to differentiate w.r.t.
  mode::global_tape->register_variable(xint);

  // Interpolated function value
  T fint;

  // Call the NAG AD Lib functions
  func(xint, fint);

  // Extract the computed solutions
  fintv = dco::passive_value(fint);
  // Initialize adjoints of y
  mode::global_tape->register_output_variable(fint);
  dco::derivative(fint) = 1.0;

  // Interpret the tape
  mode::global_tape->interpret_adjoint();

  // Adjoint of xint
  d = dco::derivative(dco::derivative(xint));

  // Remove tape
  mode::tape_t::remove(mode::global_tape);
}

// Function which calls NAG AD Library routines
template <typename T> void func(const T &xint, T &fint)
{

  // Create AD configuration data object
  Integer           ifail = 0;
  nag::ad::handle_t ad_handle;

  // Initialize abscissa and ordinate data
  std::vector<T> x{7.99, 8.09, 8.19, 8.70, 9.20, 10.00, 12.00, 15.00, 20.00};
  std::vector<T> f{0.00000E+0, 0.27643E-4, 0.43750E-1, 0.16918E+0, 0.46943E+0,
                   0.94374E+0, 0.99864E+0, 0.99992E+0, 0.99999E+0};

  // Local output of e01be; approximation of derivatives at points x
  std::vector<T> d(x.size());
  ifail = 0;
  nag::ad::e01be(ad_handle, x.size(), x.data(), f.data(), d.data(), ifail);

  // Interpolate at point x
  ifail = 0;
  nag::ad::e01bf(ad_handle, x.size(), x.data(), f.data(), d.data(), 1, &xint,
                 &fint, ifail);

  ifail = 0;
}