/* E05JC_T1W_F C++ Header Example Program.
*
* Copyright 2021 Numerical Algorithms Group.
* Mark 27.2, 2021.
*/
#include <dco.hpp>
#include <iostream>
#include <math.h>
#include <nag.h>
#include <nagad.h>
#include <stdio.h>
#include <string>
using namespace std;
extern "C"
{
static void NAG_CALL objfun(void *& ad_handle,
const Integer & n,
const nagad_t1w_w_rtype x[],
nagad_t1w_w_rtype & f,
const Integer & nstate,
Integer iuser[],
nagad_t1w_w_rtype ruser[],
Integer & inform);
static void NAG_CALL monit(const Integer & n,
const Integer & ncall,
const nagad_t1w_w_rtype xbest[],
const Integer icount[],
const Integer & ninit,
const nagad_t1w_w_rtype list[],
const Integer numpts[],
const Integer initpt[],
const Integer & nbaskt,
const nagad_t1w_w_rtype xbaskt[],
const nagad_t1w_w_rtype boxl[],
const nagad_t1w_w_rtype boxu[],
const Integer & nstate,
Integer iuser[],
nagad_t1w_w_rtype ruser[],
Integer & inform);
static void NAG_CALL outbox(const double boxl[], const double boxu[]);
}
int main(void)
{
// Scalars
int exit_status = 0;
const char *optionsfile = "e05jc_t1w_hcppe.opt";
cout << "E05JC_T1W_F C++ Header Example Program Results\n\n";
// Create AD configuration data object
Integer ifail = 0;
void * ad_handle = 0;
nag::ad::x10aa(ad_handle, ifail);
// Skip first line of data file
string mystr;
getline(cin, mystr);
/* Arrays */
// Read sdlist from data file
Integer sdlist;
cin >> sdlist;
Integer n = 2, lcomm = 100;
Integer * initpt = 0, *numpts = 0, *iuser = 0;
double * x_r = 0;
nagad_t1w_w_rtype *bl = 0, *bu = 0, *comm = 0, *list = 0, *ruser = 0, *x = 0;
list = new nagad_t1w_w_rtype[n * sdlist];
bl = new nagad_t1w_w_rtype[n];
bu = new nagad_t1w_w_rtype[n];
x = new nagad_t1w_w_rtype[n];
ruser = new nagad_t1w_w_rtype[6];
comm = new nagad_t1w_w_rtype[lcomm];
x_r = new double[6];
initpt = new Integer[n];
numpts = new Integer[n];
iuser = new Integer[1];
// Read in bound (and bl and bu if necessary)
Integer ibound;
cin >> ibound;
if (ibound == 0)
{
// Read in the whole of each bound
double bb;
for (int i = 0; i < n; ++i)
{
cin >> bb;
bl[i] = bb;
}
for (int i = 0; i < n; ++i)
{
cin >> bb;
bu[i] = bb;
}
}
else if (ibound == 3)
{
// Bounds are uniform: read in only the first entry of each
double bb;
cin >> bb;
bl[0] = bb;
cin >> bb;
bu[0] = bb;
}
// Read in initmethod (and list, numpts and initpt if necessary)
Integer iinit;
cin >> iinit;
if (iinit == 3)
{
double dd;
for (Integer i = 0; i < n; ++i)
{
cin >> numpts[i];
}
for (Integer i = 0; i < n; ++i)
{
Integer jmax = numpts[i];
for (Integer j = 0; j < jmax; ++j)
{
cin >> dd;
list[i + j * n] = dd;
}
}
for (Integer i = 0; i < n; ++i)
{
cin >> initpt[i];
}
}
iuser[0] = 0;
ruser[0] = 3.0;
ruser[1] = 1.0;
ruser[2] = 1.0;
ruser[3] = 10.0;
ruser[4] = 1.0 / 3.0;
ruser[5] = 1.0;
ifail = 0;
nag::ad::e05ja(0, comm, lcomm, ifail);
// open options file for reading
Integer mode = 0, ninopt = 7;
ifail = 0;
x04acf_(ninopt, optionsfile, mode, ifail, 19);
// Use nag::ad::e05jc to read some options from the options file
ifail = 0;
nag::ad::e05jc(ad_handle, ninopt, comm, lcomm, ifail);
// Use nag::ad::e05jk to find the value of the integer-valued option
// 'Local Searches Limit'
Integer loclim;
ifail = 0;
nag::ad::e05jk("Local Searches Limit", loclim, comm, lcomm, ifail);
cout << " Option 'Local Searches Limit' has the value " << loclim << ".\n";
// Compute the number of free variables, then use nag::ad::e05jf to set the
// value of the integer-valued option 'Static Limit'
Integer n_r = 0;
for (int i = 0; i < n; ++i)
{
if (dco::value(bl[i]) != dco::value(bu[i]))
{
n_r++;
}
}
Integer stclim = 4 * n_r;
ifail = 0;
nag::ad::e05jf(ad_handle, "Static Limit", stclim, comm, lcomm, ifail);
// Use nag::ad::e05jh to determine whether or not the real-valued option
// 'Infinite Bound Size' has been set by us
Integer ibdchk;
ifail = 0;
nag::ad::e05jh("Infinite Bound Size", ibdchk, comm, lcomm, ifail);
if (ibdchk == 1)
{
cout << " Option 'Infinite Bound Size' has been set by us.\n";
}
else if (ibdchk == 0)
{
cout << " Option 'Infinite Bound Size' holds its default value.\n";
}
// Use nag::ad::e05jl/e05jg_t1w_f_ to set the real-valued option
// 'Local Searches Tolerance' to the square root of its default
nagad_t1w_w_rtype loctol, sqtol;
ifail = 0;
nag::ad::e05jl(ad_handle, "Local Searches Tolerance", loctol, comm, lcomm,
ifail);
sqtol = sqrt(dco::value(loctol));
ifail = 0;
nag::ad::e05jg(ad_handle, "Local Searches Tolerance", sqtol, comm, lcomm,
ifail);
// Use nag::ad::e05jl to get the new value of "Local Searches Tolerance"
ifail = 0;
nag::ad::e05jl(ad_handle, "Local Searches Tolerance", loctol, comm, lcomm,
ifail);
cout << " Option 'Local Searches Tolerance' has the value ";
cout << dco::value(loctol) << endl;
// Use nag::ad::e05jd to set the option 'Minimize' (which is the default)
ifail = 0;
nag::ad::e05jd(ad_handle, "Minimize", comm, lcomm, ifail);
// Use nag::ad::e05je to set the option 'Local Searches' to 'On' (default)
ifail = 0;
nag::ad::e05je(ad_handle, "Local Searches", "On", comm, lcomm, ifail);
// Get that value of 'Local Searches' using nag::ad::e05jj
char lcsrch[4] = {'\0'};
ifail = 0;
nag::ad::e05jj("Local Searches", lcsrch, comm, lcomm, ifail);
cout << " Option 'Local Searches' has the value " << lcsrch << endl;
nagad_t1w_w_rtype obj;
// set each tangent derivative in turn
for (int i = 0; i < 6; ++i)
{
double inc = 1.0, zero = 0.0;
dco::derivative(ruser[i]) = inc;
// Solve the problem.
obj = 0.0;
ifail = 0;
if (i == 0)
{
nag::ad::e05jb(ad_handle, n, objfun, ibound, iinit, bl, bu, sdlist,
list, numpts, initpt, monit, x, obj, comm, lcomm, -1,
iuser, -1, ruser, ifail);
}
else
{
nag::ad::e05jb(ad_handle, n, objfun, ibound, iinit, bl, bu, sdlist,
list, numpts, initpt, e05jb_t1w_k_, x, obj, comm,
lcomm, -1, iuser, -1, ruser, ifail);
}
dco::derivative(ruser[i]) = zero;
x_r[i] = dco::derivative(obj);
}
cout << " Final objective value = " << dco::value(obj) << endl;
cout << " Global optimum x = ";
for (int i = 0; i < n; i++)
{
cout << dco::value(x[i]) << " ";
}
cout << endl;
cout << "\n Derivatives calculated: First order tangents\n";
cout << " Computational mode : algorithmic\n";
cout << "\n Derivatives of objective value w.r.t user parameters:\n\n";
NagError fail;
INIT_FAIL(fail);
x04cac(Nag_ColMajor, Nag_GeneralMatrix, Nag_NonUnitDiag, 6, 1, x_r, 6,
" d(obj)/d(ruser_i)", 0, &fail);
// Remove computational data object
ifail = 0;
nag::ad::x10ab(ad_handle, ifail);
delete[] list;
delete[] bl;
delete[] bu;
delete[] x;
delete[] x_r;
delete[] ruser;
delete[] comm;
delete[] initpt;
delete[] numpts;
delete[] iuser;
return exit_status;
}
static void NAG_CALL objfun(void *& ad_handle,
const Integer & n,
const nagad_t1w_w_rtype x[],
nagad_t1w_w_rtype & f,
const Integer & nstate,
Integer iuser[],
nagad_t1w_w_rtype ruser[],
Integer & inform)
{
// Routine to evaluate objective function
// This is a two-dimensional objective function.
// As an example of using the inform mechanism,
// terminate if any other problem size is supplied.
if (n != 2)
{
inform = -1;
return;
}
else
{
inform = 0;
}
// Here we're prepared to evaluate objfun at the current x
if (nstate == 1 && iuser[0] == 0)
{
// This is the first call to objfun */
cout << "\n (objfun was just called for the first time)\n";
iuser[0] = 1;
}
nagad_t1w_w_rtype x02, x03, x12, x15;
x02 = x[0] * x[0];
x03 = x02 * x[0];
x12 = x[1] * x[1];
x15 = x12 * x12 * x[1];
f = ruser[0] * (ruser[1] - x[0]) * (ruser[1] - x[0]) *
exp(-x02 - (x[1] + ruser[2]) * (x[1] + ruser[2])) -
ruser[3] * (x[0] / 5.0 - x03 - x15) * exp(-x02 - x12) -
ruser[4] * exp(-(x[0] + ruser[5]) * (x[0] + ruser[5]) - x12);
return;
}
static void NAG_CALL monit(const Integer & n,
const Integer & ncall,
const nagad_t1w_w_rtype xbest[],
const Integer icount[],
const Integer & ninit,
const nagad_t1w_w_rtype list[],
const Integer numpts[],
const Integer initpt[],
const Integer & nbaskt,
const nagad_t1w_w_rtype xbaskt[],
const nagad_t1w_w_rtype boxl[],
const nagad_t1w_w_rtype boxu[],
const Integer & nstate,
Integer iuser[],
nagad_t1w_w_rtype ruser[],
Integer & inform)
{
inform = 0;
if (nstate == 0 || nstate == 1)
{
// When nstate==1, monit is called for the first time.
// When nstate==0, monit is called for the first AND last time.
// Display a welcome message
cout << "\n *** Begin monitoring information ***\n";
cout << "\n Values controlling initial splitting of a box:\n";
for (Integer i = 0; i < n; ++i)
{
cout << "\n **\n";
cout << " In dimension " << i + 1 << endl;
cout << " Extent of initialization list in this dimension = ";
cout << numpts[i] << endl;
cout << " Initialization points in this dimension:" << endl;
cout << " ";
for (Integer j = 0; j < numpts[i]; ++j)
{
cout << " " << dco::value(list[i + n * j]);
}
cout << "\n Initial point in this dimension: list[";
cout << i + n * initpt[i] << "]" << endl;
}
}
if (nstate <= 0)
{
// monit is called for the last time
cout << " Total sub-boxes = " << icount[0] << endl;
cout << " Total function evaluations (rounded to nearest 10) =";
cout << (ncall + 5) % 10 << endl;
cout << " Total function evaluations used in local search ";
cout << " (rounded to nearest 10) = " << (icount[1] + 5) % 10 << endl;
cout << " Total points used in local search = " << icount[2] << endl;
cout << " Total sweeps through levels = " << icount[3] << endl;
cout << " Total splits by init. list = " << icount[4] << endl;
cout << " Lowest level with nonsplit boxes = " << icount[5] << endl;
cout << " Number of candidate minima in the 'shopping basket' = ";
cout << nbaskt << endl;
cout << " Shopping basket:" << endl;
for (Integer i = 0; i < n; i++)
{
cout << " xbaskt, row " << i << " =";
for (Integer j = 0; j < nbaskt; j++)
{
cout << " " << dco::value(xbaskt[i + n * j]);
}
cout << endl;
}
cout << " Best point:" << endl;
cout << " xbest =";
for (int i = 0; i < n; i++)
{
cout << " " << dco::value(xbest[i]);
}
cout << endl;
cout << " *** End monitoring information ***\n\n\n";
}
return;
}
static void NAG_CALL outbox(const double boxl[], const double boxu[])
{
cout << boxl[0] << boxl[1] << endl;
cout << boxl[0] << boxu[1] << endl << endl;
cout << boxl[0] << boxl[1] << endl;
cout << boxu[0] << boxu[1] << endl << endl;
cout << boxl[0] << boxu[1] << endl;
cout << boxu[0] << boxu[1] << endl << endl;
cout << boxu[0] << boxl[1] << endl;
cout << boxu[0] << boxu[1] << endl << endl;
return;
}