/* 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\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);
}
// 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;
auto lsqfun = [&](nag::ad::handle_t&,
Integer & iflag,
const Integer & m,
const Integer & n,
const T xc[],
T fvec[],
T fjac[],
const Integer & ldfjac)
{
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;
// 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] = -(t[i] * t[m + i]) / (di * di);
// dF/dx3
fjac[i + 2 * m] = -(t[i] * t[2 * m + i]) / (di * di);
}
}
};
// Create AD configuration data object
Integer ifail = 0;
nag::ad::handle_t ad_handle;
// Set computational mode
ad_handle.set_strategy(nag::ad::symbolic);
// 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);
}