/* nag::ad::c05ay Vector Tangent over Vector Adjoint Example Program.
 */

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

// Function which calls NAG AD Library routines.
template <typename T> void func(T r, T &x);

// Driver with the adjoint calls.
// Finds the zero point f(xv) = 0 for the function
// f(x) = exp(-x) - r*x.
// Also, computes the Hessian d2x_dr2 = d2xv/dr2.
template <std::size_t N, std::size_t M>
void driver(const double &rv, double &xv, std::array<std::array<double,N>,M> &d2x_dr2);

int main()
{
  feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
  
  std::cout << " nag::ad::c05ay Vector Tangent over Vector Adjoint Example Program Results\n";

  // Set problem parameters
  double rv = 2.0;
  // Solution x
  double xv;
  // Number of vector elements in vector mode
  constexpr std::size_t N=3, M=4;
  // Derivative of x
  std::array<std::array<double,N>,M> d2x_dr2;

  // Call driver
  driver(rv, xv, d2x_dr2);

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

  // Print derivatives
  std::cout.setf(std::ios::scientific, std::ios::floatfield);
  std::cout.precision(12);
  std::cout << "\n Solution:\n";
  std::cout << " x = " << xv << std::endl;
  std::cout << "\n Second derivative of solution x w.r.t. parameter r:\n";
  for (std::size_t i=0; i<N; ++i) {
    for (std::size_t j=0; j<M; ++j) {
      std::cout << " d2x/dr2(x)[" << j << "][" << i << "] = " << d2x_dr2[j][i] << std::endl;
    }
  }

  return 0;
}

// Driver with the adjoint calls.
// Finds the zero point f(xv) = 0 for the function
// f(x) = exp(-x) - r*x.
// Also, computes the Hessian d2x_dr2 = d2xv/dr2.
template <std::size_t N, std::size_t M>
void driver(const double &rv, double &xv, std::array<std::array<double,N>,M> &d2x_dr2)
{
  using MODE = dco::ga1v<typename dco::gt1v<double,N>::type,M>;
  using T    = typename MODE::type;

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

  // Function parameter r
  T r = rv;

  // Register variable to differentiate w.r.t.
  dco::derivative(dco::value(r))[0] = 1.0;
  for (std::size_t i=1; i<N; ++i) {
    dco::derivative(dco::value(r))[i] = 2.0;
  }
  MODE::global_tape->register_variable(r);

  // Variable to differentiate
  T x;

  // Call the NAG AD Lib functions
  func(r, x);

  // Extract the computed solution
  xv = dco::passive_value(x);

  MODE::global_tape->register_output_variable(x);
  dco::derivative(x)[0] = 1.0;
  for (std::size_t j=1; j<M; ++j) {
    dco::derivative(x)[j] = 2.0;
  }
  MODE::global_tape->interpret_adjoint();
  // Extract the derivatives
  for (std::size_t i=0; i<N; ++i) {
    for (std::size_t j=0; j<M; ++j) {
      d2x_dr2[j][i] = dco::derivative(dco::derivative(r)[j])[i];
    }
  }

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

// Function which calls NAG AD Library routines.
template <typename T> void func(T r, T &x)
{
  // Active variables
  T a = 0.0, b = 1.0;
  T eps = 1.0e-5, eta = 0.0;
  // Create AD configuration data object
  Integer           ifail = 0;
  nag::ad::handle_t ad_handle;

  // Routine for computing an approximation to a simple zero of the function f
  ifail = 0;
  nag::ad::c05ay(ad_handle, a, b, eps, eta,
                 [&](nag::ad::handle_t const& handle, T const& x, T & z) {
                   z = exp(-x) - x * r;
                 },
                 x, ifail);
}