/* nag::ad::e04us 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> &ruser, std::vector<T> &x);
// Driver with the tangent over adjoint calls.
// Solve the sum-of-squares minimization problem
// F(x) = .5 sum_{i=0}^{44} (y_i - f_i(x))^2
// for
// f_i(x) = x_1 + (0.49 - x_1) * e^{-x_2 * (a_i - 8)}
//
// and a_i is given via ruser parameter.
//
// Also, computes the sum of all Hessian entries d2x/druser2.
void driver(std::vector<double> const &ruser,
std::vector<double> & x,
double & d2xdruser2);
int main()
{
std::cout << " nag::ad::e04us Tangent Over Adjoint Example Program Results\n";
// (Initial) Solution
std::vector<double> x{0.4, 0.0};
// Parameters a
std::vector<double> ruser{
8.0, 8.0, 10.0, 10.0, 10.0, 10.0, 12.0, 12.0, 12.0, 12.0, 14.0,
14.0, 14.0, 16.0, 16.0, 16.0, 18.0, 18.0, 20.0, 20.0, 20.0, 22.0,
22.0, 22.0, 24.0, 24.0, 24.0, 26.0, 26.0, 26.0, 28.0, 28.0, 30.0,
30.0, 30.0, 32.0, 32.0, 34.0, 36.0, 36.0, 38.0, 38.0, 40.0, 42.0};
// Sum of Hessian entries dx/druser
double d2xdruser2;
driver(ruser, x, d2xdruser2);
// Print outputs
std::cout << "\n Derivatives calculated: First order tangent over adjoints\n";
std::cout << " Computational mode : algorithmic\n";
// Print derivatives
std::cout.setf(std::ios::scientific, std::ios::floatfield);
std::cout.precision(12);
std::cout << "\n Solution:\n";
std::cout << " x = [" << x[0] << ", " << x[1] << "]" << std::endl;
std::cout << "\n Sum of all Hessian entries:\n";
std::cout << " sum_ijk [d2x/druser2]_ijk = " << d2xdruser2 << std::endl;
return 0;
}
// Driver with the tangent over adjoint calls.
// Solve the sum-of-squares minimization problem
// F(x) = .5 sum_{i=0}^{44} (y_i - f_i(x))^2
// for
// f_i(x) = x_1 + (0.49 - x_1) * e^{-x_2 * (a_i - 8)}
//
// and a_i is given via ruser parameter.
//
// Also, computes the sum of all Hessian entries d2x/druser2.
void driver(std::vector<double> const &ruserv,
std::vector<double> & xv,
double & d2xdruser2)
{
using mode = dco::ga1s<dco::gt1s<double>::type>;
using T = mode::type;
// Create AD tape
mode::global_tape = mode::tape_t::create();
std::vector<T> x(begin(xv), end(xv));
std::vector<T> ruser(begin(ruserv), end(ruserv));
dco::derivative(dco::value(ruser)) = std::vector<double>(ruser.size(), 1.0);
mode::global_tape->register_variable(ruser);
// Call the NAG AD Lib functions
func(ruser, x);
// Get output
xv = dco::passive_value(x);
// Set up evaluation of derivatives
mode::global_tape->register_output_variable(x);
dco::value(dco::derivative(x)) = std::vector<double>(x.size(), 1.0);
mode::global_tape->interpret_adjoint();
// Get sum of Hessian elements.
d2xdruser2 = 0;
for (T const &r : ruser)
{
d2xdruser2 += dco::derivative(dco::derivative(r));
}
// Remove tape
mode::tape_t::remove(mode::global_tape);
}
template <typename T> void func(std::vector<T> &ruser, std::vector<T> &x)
{
// Problem sizes
Integer m = 44, n = 2, nclin = 1, ncnln = 1, lb = n + nclin + ncnln,
liwork = 3 * n + nclin + 2 * ncnln, lda = nclin, sda = n,
ldcj = ncnln, ldfj = m, ldr = n, lwork = 20 * n + m * (n + 3);
lwork += 2 * n * n + 11 * nclin;
lwork += n * nclin + (2 * n + 21) * ncnln;
std::vector<T> c(ncnln), cjac(ncnln * n), f(m), fjac(m * n), clamda(lb),
r(ldr * n), rwsav(475), work(lwork);
std::vector<Integer> iwsav(610), iwork(liwork), istate(lb);
std::vector<logical> lwsav(120);
// Problem data
std::vector<T> a{1.0, 1.0};
std::vector<T> y{0.49, 0.49, 0.48, 0.47, 0.48, 0.47, 0.46, 0.46, 0.45,
0.43, 0.45, 0.43, 0.43, 0.44, 0.43, 0.43, 0.46, 0.45,
0.42, 0.42, 0.43, 0.41, 0.41, 0.40, 0.42, 0.40, 0.40,
0.41, 0.40, 0.41, 0.41, 0.40, 0.40, 0.40, 0.38, 0.41,
0.40, 0.40, 0.41, 0.38, 0.40, 0.40, 0.39, 0.39};
std::vector<T> bl{0.4, -4.0, 1.0, 0.0};
std::vector<T> bu{1.0E+25, 1.0E+25, 1.0E+25, 1.0E+25};
// Initialize sav arrays
Integer ifail = 0;
char cwsav;
nag::ad::e04wb("E04USA", &cwsav, 1, lwsav.data(), lwsav.size(), iwsav.data(),
iwsav.size(), rwsav.data(), rwsav.size(), ifail);
auto objfun = [&](nag::ad::handle_t &ad_handle,
Integer & mode,
const Integer & m,
const Integer & n,
const Integer & ldfj,
const Integer & needfi,
const T x[],
T f[],
T fjac[],
const Integer & nstate)
{
T x1{x[0]}, x2{x[1]};
#define FJAC(i, j) fjac[(i) + (j)*m]
if (mode == 0 && needfi > 0)
{
// Update only one element f(i)
Integer i = needfi - 1;
T ai = ruser[i] - 8.0;
f[i] = x1 + (0.49 - x1) * exp(-x2 * ai);
}
else
{
for (int i = 0; i < m; ++i)
{
T ai = ruser[i] - 8.0;
if (mode == 0 || mode == 2)
{
// Update the whole of f(i)
f[i] = x1 + (0.49 - x1) * exp(-x2 * ai);
}
if (mode == 1 || mode == 2)
{
// Update df/dx1
FJAC(i, 0) = 1.0 - exp(-x2 * ai);
// Update df/dx2
FJAC(i, 1) = -(0.49 - x1) * ai * exp(-x2 * ai);
}
}
}
#undef FJAC
};
auto confun = [&](nag::ad::handle_t &ad_handle,
Integer & mode,
const Integer & ncnln,
const Integer & n,
const Integer & ldcj,
const Integer needc[],
const T x[],
T c[],
T cjac[],
const Integer & nstate)
{
// calling confun for the first time
if (nstate == 1)
{
for (int i = 0; i < ncnln * n; ++i)
{
cjac[i] = 0.0;
}
}
// evaluate c and cjac (depending on mode) at point x
if (needc[0] > 0)
{
// c required
if (mode == 0 || mode == 2)
{
c[0] = -0.09 - x[0] * x[1] + 0.49 * x[1];
}
// cjac required
if (mode == 1 || mode == 2)
{
cjac[0] = -x[1];
cjac[ncnln] = -x[0] + 0.49;
}
}
};
// Create AD configuration data object
nag::ad::handle_t ad_handle;
ifail = 0;
// Solve the problem
Integer iter;
T objf;
ifail = 0;
nag::ad::e04us(ad_handle, m, n, nclin, ncnln, lda, ldcj, ldfj, ldr, a.data(),
bl.data(), bu.data(), y.data(), confun, objfun, iter,
istate.data(), c.data(), cjac.data(), f.data(), fjac.data(),
clamda.data(), objf, r.data(), x.data(), iwork.data(), liwork,
work.data(), lwork, lwsav.data(), iwsav.data(), rwsav.data(), ifail);
}