NAG Library Manual, Mark 27.2
Interfaces:  FL   CL   CPP   AD
```/* 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    : symbolic (expert mode)\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);
// Set computational mode
ifail        = 0;
Integer mode = nagad_symbolic_expert;
nag::ad::x10ac(ad_handle, mode, 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[])
{
using value_t = typename dco::mode<T>::value_t;
// Get the callback computational mode
Integer ifail = 0, cb_mode;
nag::ad::x10bd(ad_handle, cb_mode, ifail);

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

// Extract the value of our state (xc) and parameters (ruser).
// This is only needed because we are using e04gb in expert symbolic
// adjoint mode.
value_t xv1 = dco::value(x1), xv2 = dco::value(x2), xv3 = dco::value(x3);
std::vector<value_t> yv(m), tv1(m), tv2(m), tv3(m);
for (int i = 0; i < m; i++)
{
yv[i]  = dco::value(y[i]);
tv1[i] = dco::value(t1[i]);
tv2[i] = dco::value(t2[i]);
tv3[i] = dco::value(t3[i]);
}

if (cb_mode == nagad_primal)
{
// We are in the forward pass of routine, e04gb only wants
// the primal evaluation, so compute with values (but
// assign to active output, since it's the only way to pass
// value back to caller)
for (int i = 0; i < m; i++)
{
value_t di = xv2 * tv2[i] + xv3 * tv3[i];
value_t yi = xv1 + tv1[i] / di;
// The values of the residuals
fvec[i] = yi - yv[i];
// Evaluate the Jacobian
if (iflag > 0)
{
// dF/dx1
fjac[i] = 1.0;
// dF/dx2
fjac[i + m] = -(tv1[i] * tv2[i]) / (di * di);
// dF/dx3
fjac[i + 2 * m] = -(tv1[i] * tv3[i]) / (di * di);
}
}
}
else if (cb_mode == nagad_dstate)
{
// e04gb wants derivatives w.r.t. state, so use values of the
// parameter data y,t1,t2,t3; this is required to not propagate
// adjoints to parameters accidentally
for (int i = 0; i < m; i++)
{
T di = x2 * tv2[i] + x3 * tv3[i];
T yi = x1 + tv1[i] / di;
// The values of the residuals
fvec[i] = yi - yv[i];
// Evaluate the Jacobian
if (iflag > 0)
{
// dF/dx1
fjac[i] = 1.0;
// dF/dx2
fjac[i + m] = -(tv1[i] * tv2[i]) / (di * di);
// dF/dx3
fjac[i + 2 * m] = -(tv1[i] * tv3[i]) / (di * di);
}
}
}
else if (cb_mode == nagad_dparam)
{
// e04gb wants derivatives w.r.t. the parameter data, so use values
// the of state data x
for (int i = 0; i < m; i++)
{
T di = xv2 * t2[i] + xv3 * t3[i];
T yi = xv1 + 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);
}
}
}
else
{
// In this case cb_mode == nagad_dall
// e04gb wants all derivatives, so compute with active data
for (int i = 0; i < m; i++)
{
T di = x2 * t2[i] + x3 * t3[i];
T yi = xv1 + 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);
}
}
}
}
```