/* nag::ad::e04gb Tangent 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 tangent 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);
// 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 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 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 tangents\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 tangent 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::gt1s<double>::type;
// 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];
}
dco::derivative(y) = std::vector<double>(y.size(), 1.0);
dco::derivative(t) = std::vector<double>(t.size(), 1.0);
// 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);
// Get sum of gradient elements of fsumsq w.r.t. y and t
dfdall = dco::derivative(fsumsq);
// Get sum of Jacobian elements of solution points x w.r.t. y and t
dxdall = 0;
for (int i = 0; i < n; i++)
{
dxdall += dco::derivative(x[i]);
}
}
// 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);
// 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[])
{
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];
for (int i = 0; i < m; i++)
{
T di = x2 * t2[i] + x3 * t3[i];
T yi = x1 + 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);
}
}
}