/* nag::ad::e04gb Tangent over 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
// 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);
int main()
{
std::cout << " nag::ad::e04gb Tangent over Adjoint Example Program Results\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 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);
}
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<dco::gt1s<double>::type> value_of(std::size_t size, T const& data) {
std::vector<dco::gt1s<double>::type> 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);
}