NAG Library Manual, Mark 27.2
Interfaces:  FL   CL   CPP   AD 

NAG AD Library Introduction
Example description
/* nag::ad::e04gb Tangent over Adjoint Example Program.
 */

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

std::stringstream filecontent("0.14  1.0 15.0  1.0\
 0.18  2.0 14.0  2.0\
 0.22  3.0 13.0  3.0\
 0.25  4.0 12.0  4.0\
 0.29  5.0 11.0  5.0\
 0.32  6.0 10.0  6.0\
 0.35  7.0  9.0  7.0\
 0.39  8.0  8.0  8.0\
 0.37  9.0  7.0  7.0\
 0.58 10.0  6.0  6.0\
 0.73 11.0  5.0  5.0\
 0.96 12.0  4.0  4.0\
 1.34 13.0  3.0  3.0\
 2.10 14.0  2.0  2.0\
 4.39 15.0  1.0  1.0");

// Function which calls NAG AD Library routines.
template <typename T>
void func(std::vector<T> &y,
          std::vector<T> &t,
          std::vector<T> &x,
          std::vector<T> &fvec,
          T &             fsumsq);

// Driver with the adjoint calls.
// Computes the minimum of the sum of squares of m nonlinear functions, the
// solution point and the corresponding residuals. Also, computes the sum of
// Hessian elements of fsumsq w.r.t. inputs y and t, and the sum of Hessian
// elements of x w.r.t. inputs y and t.
void driver(const std::vector<double> &yv,
            std::vector<double> &      tv,
            std::vector<double> &      xv,
            std::vector<double> &      fvecv,
            double &                   fsumsqv,
            double &                   d2f_dall2,
            double &                   d2x_dall2);

// Evaluates the residuals and their 1st derivatives.
template <typename T>
void NAG_CALL lsqfun(void *&        ad_handle,
                     Integer &      iflag,
                     const Integer &m,
                     const Integer &n,
                     const T        xc[],
                     T              fvec[],
                     T              fjac[],
                     const Integer &ldfjac,
                     Integer        iuser[],
                     T              ruser[]);

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

  // Problem dimensions
  const Integer       m = 15, n = 3, nt = 3;
  std::vector<double> yv(m), tv(m * nt);
  for (int i = 0; i < m; i++)
    {
      filecontent >> yv[i];
      for (int j = 0; j < nt; j++)
        {
          Integer k = j * m + i;
          filecontent >> tv[k];
        }
    }

  std::vector<double> xv(n), fvecv(m);

  // Sum of squares of the residuals at the computed point xv
  double fsumsqv;

  // Sum of Hessian elements of sum of squares fsumsqv with respect to
  // parameters y, t1, t2 and t3
  double d2f_dall2;
  // Sum of Hessian elements of x with respect to parameters y, t1, t2 and t3
  double d2x_dall2;

  // Call driver
  driver(yv, tv, xv, fvecv, fsumsqv, d2f_dall2, d2x_dall2);

  // Primal results
  std::cout.setf(std::ios::scientific, std::ios::floatfield);
  std::cout.precision(12);
  std::cout << "\n Sum of squares = ";
  std::cout.width(20);
  std::cout << fsumsqv;
  std::cout << "\n Solution point = ";
  for (int i = 0; i < n; i++)
    {
      std::cout.width(20);
      std::cout << xv[i];
    }
  std::cout << std::endl;

  std::cout << "\n Residuals :\n";
  for (int i = 0; i < m; i++)
    {
      std::cout.width(20);
      std::cout << fvecv[i] << std::endl;
    }
  std::cout << std::endl;

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

  // Print derivatives of fsumsq
  std::cout
      << "\n Sum of Hessian elements of sum of squares fsumsq w.r.t. parameters y and t:\n";
  std::cout << " sum_i [d2fsumsq/dall2]_ij = " << d2f_dall2 << std::endl;

  // Print derivatives of solution points
  std::cout
      << "\n Sum of Hessian elements of solution points x w.r.t. parameters y and t:\n";
  std::cout << " sum_ij [d2x/dall2]_ij = " << d2x_dall2 << std::endl;

  return 0;
}

// Driver with the adjoint calls.
// Computes the minimum of the sum of squares of m nonlinear functions, the
// solution point and the corresponding residuals. Also, computes the sum of
// Hessian elements of fsumsq w.r.t. inputs y and t, and the sum of Hessian
// elements of x w.r.t. inputs y and t.
void driver(const std::vector<double> &yv,
            std::vector<double> &      tv,
            std::vector<double> &      xv,
            std::vector<double> &      fvecv,
            double &                   fsumsqv,
            double &                   d2f_dall2,
            double &                   d2x_dall2)
{
  using mode = dco::ga1s<dco::gt1s<double>::type>;
  using T    = mode::type;

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

  // Problem dimensions
  const Integer m = fvecv.size(), n = xv.size();
  const Integer nt = n;

  // AD routine arguments
  std::vector<T> y(m), t(m * nt), x(n), fvec(m);
  T              fsumsq;
  for (int i = 0; i < m; i++)
    {
      y[i] = yv[i];
    }
  for (int i = 0; i < m * nt; i++)
    {
      t[i] = tv[i];
    }

  dco::derivative(dco::value(y)) = std::vector<double>(y.size(), 1.0);
  dco::derivative(dco::value(t)) = std::vector<double>(t.size(), 1.0);
  // Register variables to differentiate w.r.t.
  mode::global_tape->register_variable(y);
  mode::global_tape->register_variable(t);

  // Call the NAG AD Lib functions
  func(y, t, x, fvec, fsumsq);

  // Sum of squares
  fsumsqv = dco::passive_value(fsumsq);
  // Solution point
  xv = dco::passive_value(x);
  // Residuals
  fvecv = dco::passive_value(fvec);

  // Set up evaluation of derivatives of fsumsq via adjoints
  mode::global_tape->register_output_variable(fsumsq);
  dco::derivative(fsumsq) = 1.0;
  mode::global_tape->interpret_adjoint();

  // Get sum of Hessian elements of fsumsq w.r.t. y and t
  d2f_dall2 = 0;
  for (int i = 0; i < y.size(); i++)
    {
      d2f_dall2 += dco::derivative(dco::derivative(y[i]));
    }
  for (int i = 0; i < t.size(); i++)
    {
      d2f_dall2 += dco::derivative(dco::derivative(t[i]));
    }

  // Get sum of Hessian elements of solution points x w.r.t. y and t
  mode::global_tape->register_output_variable(x);
  mode::global_tape->zero_adjoints();
  dco::derivative(x) = std::vector<dco::gt1s<double>::type>(n, 1.0);
  mode::global_tape->interpret_adjoint();

  d2x_dall2 = 0;
  for (int i = 0; i < y.size(); i++)
    {
      d2x_dall2 += dco::derivative(dco::derivative(y[i]));
    }
  for (int i = 0; i < t.size(); i++)
    {
      d2x_dall2 += dco::derivative(dco::derivative(t[i]));
    }

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

// Function which calls NAG AD Library routines.
template <typename T>
void func(std::vector<T> &y,
          std::vector<T> &t,
          std::vector<T> &x,
          std::vector<T> &fvec,
          T &             fsumsq)
{
  // Problem dimensions
  const Integer m = fvec.size(), n = x.size();
  const Integer ldfjac = m, nt = n, ldv = n;

  // All additional data accessed in the callback MUST be in ruser.
  // Pack the parameters [y, t1, t2, t3] into the columns of ruser
  std::vector<T> ruser(m * nt + m);
  T *            _y = &ruser[0];
  T *            _t = &ruser[m];

  for (int i = 0; i < m; i++)
    {
      _y[i] = y[i];
    }

  for (int i = 0; i < m * nt; i++)
    {
      _t[i] = t[i];
    }

  // Initial guess of the position of the minimum
  dco::passive_value(x[0]) = 0.5;
  for (int i = 1; i < n; i++)
    {
      dco::passive_value(x[i]) = dco::passive_value(x[i - 1]) + 0.5;
    }

  std::vector<T> s(n), v(ldv * n), fjac(m * n);
  Integer        niter, nf;

  Integer iprint = -1;
  Integer selct  = 2;
  Integer maxcal = 200 * n;
  T       eta    = 0.5;
  T       xtol   = 10.0 * sqrt(X02AJC);
  T       stepmx = 10.0;

  // Create AD configuration data object
  Integer ifail     = 0;
  void *  ad_handle = 0;
  nag::ad::x10aa(ad_handle, ifail);
  // Routine for computing the minimum of the sum of squares of m nonlinear
  // functions.
  ifail = 0;
  nag::ad::e04gb(ad_handle, m, n, selct, lsqfun, nullptr, iprint, maxcal, eta,
                 xtol, stepmx, x.data(), fsumsq, fvec.data(), fjac.data(),
                 ldfjac, s.data(), v.data(), ldv, niter, nf, 0, nullptr,
                 ruser.size(), ruser.data(), ifail);

  // Remove computational data object
  ifail = 0;
  nag::ad::x10ab(ad_handle, ifail);
}

// Evaluates the residuals and their 1st derivatives.
template <typename T>
void NAG_CALL lsqfun(void *&        ad_handle,
                     Integer &      iflag,
                     const Integer &m,
                     const Integer &n,
                     const T        xc[],
                     T              fvec[],
                     T              fjac[],
                     const Integer &ldfjac,
                     Integer        iuser[],
                     T              ruser[])
{
  const T *y  = ruser;
  const T *t1 = ruser + m;
  const T *t2 = ruser + 2 * m;
  const T *t3 = ruser + 3 * m;
  const T  x1 = xc[0], x2 = xc[1], x3 = xc[2];

  for (int i = 0; i < m; i++)
    {
      T di = x2 * t2[i] + x3 * t3[i];
      T yi = x1 + t1[i] / di;
      // The values of the residuals
      fvec[i] = yi - y[i];
      // Evaluate the Jacobian
      if (iflag > 0)
        {
          // dF/dx1
          fjac[i] = 1.0;
          // dF/dx2
          fjac[i + m] = -(t1[i] * t2[i]) / (di * di);
          // dF/dx3
          fjac[i + 2 * m] = -(t1[i] * t3[i]) / (di * di);
        }
    }
}