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

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

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

// 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
// gradient elements of fsumsq w.r.t. inputs y and t, and the sum of Jacobian
// 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 &                   dfdall,
            double &                   dxdall);

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

  // Problem dimensions
  const Integer       m = 15, n = 3;
  std::vector<double> yv = {0.14, 0.18, 0.22, 0.25, 0.29, 0.32, 0.35, 0.39,
                            0.37, 0.58, 0.73, 0.96, 1.34, 2.10, 4.39};
  std::vector<double> tv = {
      1.0,  2.0,  3.0,  4.0,  5.0,  6.0,  7.0,  8.0,  9.0,  10.0, 11.0, 12.0,
      13.0, 14.0, 15.0, 15.0, 14.0, 13.0, 12.0, 11.0, 10.0, 9.0,  8.0,  7.0,
      6.0,  5.0,  4.0,  3.0,  2.0,  1.0,  1.0,  2.0,  3.0,  4.0,  5.0,  6.0,
      7.0,  8.0,  7.0,  6.0,  5.0,  4.0,  3.0,  2.0,  1.0};

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

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

  // Sum of gradient elements of sum of squares fsumsqv with respect to the
  // parameters y, t1, t2, and t3
  double dfdall;
  // Sum of Jacobian elements of x with respect to the parameters y, t1, t2, and
  // t3
  double dxdall;

  // Call driver
  driver(yv, tv, xv, fvecv, fsumsqv, dfdall, dxdall);

  // 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: First order adjoints\n";
  std::cout << " Computational mode    : symbolic (expert mode)\n\n";

  // Print derivatives of fsumsq
  std::cout
      << "\n Sum of gradient elements of sum of squares fsumsq w.r.t. parameters y and t:\n";
  std::cout << " sum_i [dfsumsq/dall_i] = " << dfdall << std::endl;

  // Print derivatives of solution points
  std::cout
      << "\n Sum of Jacobian elements of solution points x w.r.t. parameters y and t:\n";
  std::cout << " sum_ij [dx/dall]_ij = " << dxdall << 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
// gradient elements of fsumsq w.r.t. inputs y and t, and the sum of Jacobian
// 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 &                   dfdall,
            double &                   dxdall)
{
  using T = dco::ga1s<double>::type;

  // Create AD tape
  dco::ga1s<double>::global_tape = dco::ga1s<double>::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(nt * m), 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];
  }

  // Register variables to differentiate w.r.t.
  dco::ga1s<double>::global_tape->register_variable(y);
  dco::ga1s<double>::global_tape->register_variable(t);

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

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

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

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

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

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

  // Remove tape
  dco::ga1s<double>::tape_t::remove(dco::ga1s<double>::global_tape);
}

template <typename T, typename U, typename V>
  void eval(int iflag, int m, int n, U const &xc, T *fvec,
            T *fjac, V const &y, V const &t) {
  const T x1 = xc[0], x2 = xc[1], x3 = xc[2];
  for (int i = 0; i < m; i++)
      {
        T di    = x2 * t[m + i] + x3 * t[2 * m + i];
        T yi    = x1 + t[i] / di;
        fvec[i] = yi - y[i];
        if (iflag > 0)
          {
            fjac[i]         = 1.0;
            fjac[i + m]     = -(t[i] * t[m + i]) / (di * di);
            fjac[i + 2 * m] = -(t[i] * t[2 * m + i]) / (di * di);
          }
      }
}

// Returns vector of values of data. The returned type is of lower
// differentiation order than T, i.e. a double type.
template <typename T>
std::vector<double> value_of(std::size_t size,  T const& data) {
  std::vector<double> out(size);
  for (std::size_t i = 0; i < size; ++i) {
    out[i] = dco::value(data[i]);
  }
  return out;
};

// 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;

  // 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;

  // Evaluates the residuals and their 1st derivatives.
  auto lsqfun = [&](nag::ad::handle_t &ad_handle, Integer &iflag,
                    const Integer &m, const Integer &n, const T xc[], T fvec[],
                    T fjac[], const Integer &ldfjac) {
                  using value_t = typename dco::mode<T>::value_t;
                  std::vector<value_t> xcp, yp, tp;
                  switch (ad_handle.active_inputs())
                    {
                    case nag::ad::none:
                      eval(iflag, m, n, value_of(n, xc), fvec, fjac, value_of(m, y), value_of(m*nt, t));
                      break;
                    case nag::ad::state:
                      eval(iflag, m, n, xc, fvec, fjac, value_of(m, y), value_of(m*nt, t));
                      break;
                    case nag::ad::params:
                      eval(iflag, m, n, value_of(n, xc), fvec, fjac, y, t);
                      break;
                    case nag::ad::all:
                      eval(iflag, m, n, xc, fvec, fjac, y, t);
                      break;
                    }
                };
  
  // Create AD configuration data object
  Integer           ifail = 0;
  nag::ad::handle_t ad_handle;

  // Set computational mode
  ad_handle.set_strategy(nag::ad::symbolic_expert);
  // 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, ifail);
}