/* E05JC_A1W_F C++ Header Example Program.
*
* Copyright 2020 Numerical Algorithms Group.
* Mark 27.1, 2020.
*/
#include <nag.h>
#include <dco.hpp>
#include <nagad.h>
#include <stdio.h>
#include <math.h>
#include <iostream>
#include <string>
using namespace std;
extern "C" {
static void NAG_CALL objfun(void* &ad_handle,
const Integer &n,
const nagad_a1w_w_rtype x[],
nagad_a1w_w_rtype &f,
const Integer &nstate,
Integer iuser[],
nagad_a1w_w_rtype ruser[],
Integer &inform);
static void NAG_CALL monit(const Integer &n,
const Integer &ncall,
const nagad_a1w_w_rtype xbest[],
const Integer icount[],
const Integer &ninit,
const nagad_a1w_w_rtype list[],
const Integer numpts[],
const Integer initpt[],
const Integer &nbaskt,
const nagad_a1w_w_rtype xbaskt[],
const nagad_a1w_w_rtype boxl[],
const nagad_a1w_w_rtype boxu[],
const Integer &nstate,
Integer iuser[],
nagad_a1w_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_a1w_hcppe.opt";
cout << "E05JC_A1W_F C++ Header Example Program Results\n\n";
// Create AD tape
nagad_a1w_ir_create();
// Create AD configuration data object
Integer ifail = 0;
void *ad_handle = 0;
x10aa_a1w_f_(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;
nagad_a1w_w_rtype *bl=0, *bu=0, *comm=0, *list=0, *ruser=0, *x=0;
list = new nagad_a1w_w_rtype [n*sdlist];
bl = new nagad_a1w_w_rtype [n];
bu = new nagad_a1w_w_rtype [n];
x = new nagad_a1w_w_rtype [n];
ruser = new nagad_a1w_w_rtype [6];
comm = new nagad_a1w_w_rtype [lcomm];
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];
}
}
// Plot determines whether monit displays information on
// the current search box
Integer plot;
cin >> plot;
// Communicate plot through to monit
iuser[0] = plot;
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;
e05ja_a1w_f_(0,comm,lcomm,ifail);
// open options file for reading
Integer mode = 0, ninopt = 7;
ifail = 0;
x04acf_(ninopt,optionsfile,mode,ifail,19);
// Use e05jc_a1w_f_ to read some options from the options file
ifail = 0;
e05jc_a1w_f_(ad_handle,ninopt,comm,lcomm,ifail);
// Use e05jk_a1w_f_ to find the value of the integer-valued option
// 'Local Searches Limit'
Integer loclim;
ifail = 0;
e05jk_a1w_f_("Local Searches Limit",loclim,comm,lcomm,ifail,20);
cout << " Option 'Local Searches Limit' has the value " << loclim << ".\n";
// Compute the number of free variables, then use e05jf_a1w_f_ to set the
// value of the integer-valued option 'Static Limit'
Integer n_r = 0;
for (int i=0; i<n; ++i) {
if (nagad_a1w_get_value(bl[i]) != nagad_a1w_get_value(bu[i])) {
n_r++;
}
}
Integer stclim = 4*n_r;
ifail = 0;
e05jf_a1w_f_(ad_handle,"Static Limit",stclim,comm,lcomm,ifail,12);
// Use e05jh_a1w_f_ to determine whether or not the real-valued option
// 'Infinite Bound Size' has been set by us
Integer ibdchk;
ifail = 0;
e05jh_a1w_f_("Infinite Bound Size",ibdchk,comm,lcomm,ifail,19);
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 e05jl_a1w_f_/e05jg_a1w_f_ to set the real-valued option
// 'Local Searches Tolerance' to the square root of its default
nagad_a1w_w_rtype loctol, sqtol;
ifail = 0;
e05jl_a1w_f_(ad_handle,"Local Searches Tolerance",loctol,comm,lcomm,ifail,24);
sqtol = sqrt(nagad_a1w_get_value(loctol));
ifail = 0;
e05jg_a1w_f_(ad_handle,"Local Searches Tolerance",sqtol,comm,lcomm,ifail,24);
// Use e05jl_a1w_f_ to get the new value of "Local Searches Tolerance"
ifail = 0;
e05jl_a1w_f_(ad_handle,"Local Searches Tolerance",loctol,comm,lcomm,ifail,24);
cout << " Option 'Local Searches Tolerance' has the value ";
cout << nagad_a1w_get_value(loctol) << endl;
// Use e05jd_a1w_f_ to set the option 'Minimize' (which is the default)
ifail = 0;
e05jd_a1w_f_(ad_handle,"Minimize",comm,lcomm,ifail,8);
// Use e05je_a1w_f_ to set the option 'Local Searches' to 'On' (default)
ifail = 0;
e05je_a1w_f_(ad_handle,"Local Searches","On",comm,lcomm,ifail,14,2);
// Get that value of 'Local Searches' using e05jj_a1w_f_
char lcsrch[4] = {'\0'};
ifail = 0;
e05jj_a1w_f_("Local Searches",lcsrch,comm,lcomm,ifail,14,3);
cout << " Option 'Local Searches' has the value " << lcsrch << endl;
// Register variables to differentiate w.r.t.
for (int i=0; i<6; ++i) {
nagad_a1w_ir_register_variable(&ruser[i]);
}
// Solve the problem.
nagad_a1w_w_rtype obj;
ifail = 0;
e05jb_a1w_f_(ad_handle,n,objfun,ibound,iinit,bl,bu,sdlist,list,numpts,
initpt,monit,x,obj,comm,lcomm,iuser,ruser,ifail);
cout << " Final objective value = " << nagad_a1w_get_value(obj) << endl;
cout << " Global optimum x = ";
for (int i=0; i<n; i++) {
cout << nagad_a1w_get_value(x[i]) << " ";
}
cout << endl;
cout << "\n Derivatives calculated: First order adjoints\n";
cout << " Computational mode : algorithmic\n\n";
cout << " Derivatives:\n\n";
// Setup evaluation of derivatives of objf via adjoints.
Integer inc = 1.0;
nagad_a1w_inc_derivative(&obj,inc);
ifail = 0;
nagad_a1w_ir_interpret_adjoint(ifail);
cout.setf(ios::fixed);
cout.setf(ios::right);
cout.precision(4);
// Get derivatives of objf w.r.t. ruser
cout << " derivatives of obj w.r.t ruser[0:5]:\n";
for (int i=0; i<6; i++) {
double d = nagad_a1w_get_derivative(ruser[i]);
cout.width(4); cout << i << " ";
cout.width(12); cout << d << endl;
}
cout << endl;
// Remove computational data object and tape
ifail = 0;
x10ab_a1w_f_(ad_handle,ifail);
nagad_a1w_ir_remove();
delete [] list;
delete [] bl;
delete [] bu;
delete [] x;
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_a1w_w_rtype x[],
nagad_a1w_w_rtype &f,
const Integer &nstate,
Integer iuser[],
nagad_a1w_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) {
// This is the first call to objfun */
cout << "\n (objfun was just called for the first time)\n";
}
nagad_a1w_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_a1w_w_rtype xbest[],
const Integer icount[],
const Integer &ninit,
const nagad_a1w_w_rtype list[],
const Integer numpts[],
const Integer initpt[],
const Integer &nbaskt,
const nagad_a1w_w_rtype xbaskt[],
const nagad_a1w_w_rtype boxl[],
const nagad_a1w_w_rtype boxu[],
const Integer &nstate,
Integer iuser[],
nagad_a1w_w_rtype ruser[],
Integer &inform)
{
logical plot = (iuser[0] == 1);
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 << " " << nagad_a1w_get_value(list[i+n*j]);
}
cout << "\n Initial point in this dimension: list[";
cout << i + n*initpt[i] << "]" << endl;
}
if (plot && n==2) {
cout << "\n <Begin displaying search boxes>\n";
}
}
if (plot && n==2) {
// Display the coordinates of the edges of the current search box
double boxl_r[2], boxu_r[2];
for (int i=0; i<2; ++i) {
boxl_r[i] = nagad_a1w_get_value(boxl[i]);
boxu_r[i] = nagad_a1w_get_value(boxu[i]);
}
outbox(boxl_r,boxu_r);
}
if (nstate<=0) {
// monit is called for the last time
if (plot && n==2) {
cout << " <End displaying search boxes>\n\n";
}
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 << " " << nagad_a1w_get_value(xbaskt[i + n*j]);
}
cout << endl;
}
cout << " Best point:" << endl;
cout << " xbest =";
for (int i = 0; i < n; i++) {
cout << " " << nagad_a1w_get_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;
}