NAG Library Manual, Mark 29
Interfaces:  FL   CL   CPP   AD 

NAG AD Library Introduction
Example description
/* E05JC_T1W_F C++ Header Example Program.
 *
 * Copyright 2023 Numerical Algorithms Group.
 * Mark 29.0, 2023.
 */
#include <dco.hpp>
#include <iostream>
#include <math.h>
#include <nag.h>
#include <nagad.h>
#include <stdio.h>
#include <string>
using namespace std;

int main()
{
  // 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;
  nag::ad::handle_t ad_handle;

  // 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;

  auto objfun = [&](nag::ad::handle_t &     ad_handle,
                  const Integer &         n,
                  const nagad_t1w_w_rtype *x,
                  nagad_t1w_w_rtype &     f,
                  const Integer &         nstate,
                  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);
                };

  auto 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;
              };

  auto 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 &               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";
                  }
              };
  
  // 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, ifail);
    }
    else
    {
      nag::ad::e05jb(ad_handle, n, objfun, ibound, iinit, bl, bu, sdlist, list,
                     numpts, initpt, nullptr, x, obj, comm, lcomm, 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);

  ifail = 0;

  delete[] list;
  delete[] bl;
  delete[] bu;
  delete[] x;
  delete[] x_r;
  delete[] ruser;
  delete[] comm;
  delete[] initpt;
  delete[] numpts;
  delete[] iuser;

  return exit_status;
}