#include "dco.hpp"
/* X10AA_A1W_F C++ Header Example Program.
 *
 * Copyright 2023 Numerical Algorithms Group.
 * Mark 29.1, 2023.
 */

#include <iostream>
#include <nag.h>
#include <nagad.h>
#include <stdio.h>
using namespace std;

int main()
{
  /* Scalars */
  Integer           exit_status = 0, ierr_not_implemented = -199;
  nagad_a1w_w_rtype x, y;
  Integer           ifail, mode;
  double            dydx, dydx_exp, y1, inc;

  /* Arrays */
  nag::ad::handle_t handle;

  cout << "X10AA_A1W_F C++ Header Example Program Results\n\n";

  /*     Create AD tape */
  dco::ga1s<double>::global_tape = dco::ga1s<double>::tape_t::create();

  x = 0.1;

  cout << "Computing ln(1+x) and its adjoint using s01ba\n\n";

  /*     Create AD configuration data object */
  ifail = 0;
  nag::ad::x10aa(handle.data(), ifail);

  /*     Set computational mode to nagad_symbolic */
  mode  = nagad_symbolic;
  ifail = 0;
  nag::ad::x10ac(handle.data(), mode, ifail);

  /*     Get and print computational mode */
  /*     Register variables to differentiate w.r.t. */
  ifail = 0;
  nag::ad::x10ad(handle.data(), mode, ifail);
  if (mode == nagad_symbolic)
  {
    cout << "Symbolic computational mode is being used\n\n";
  }
  else if (mode == nagad_algorithmic)
  {
    cout << "Algorithmic computational mode is being used\n\n";
  }

  /* Try AD computation that expects to use only the algorithmic mode */
  ifail = 1;
  nag::ad::s01ba(handle, x, y, ifail);
  if (ifail == ierr_not_implemented)
  {
    cout << "AD computation failed as expected with ifail = " << ifail
         << "\n\n";
  }
  else
  {
    cout << " Unexpected IFAIL returned in AD computation\n";
    cout << "   ifail expected: " << ierr_not_implemented << ifail << "\n";
    cout << "   ifail returned: " << ifail << "\n";
    exit_status = 1;
    goto END;
  }

  /*     Reset computational mode to nagad_algorithmic */
  mode  = nagad_algorithmic;
  ifail = 0;
  nag::ad::x10ac(handle.data(), mode, ifail);
  ifail = 0;
  nag::ad::x10ad(handle.data(), mode, ifail);
  if (mode == nagad_symbolic)
  {
    cout << "Symbolic computational mode is being used\n\n";
  }
  else if (mode == nagad_algorithmic)
  {
    cout << "Algorithmic computational mode is being used\n\n";
  }

  dco::ga1s<double>::global_tape->register_variable(x);

  ifail = 0;
  nag::ad::s01ba(handle, x, y, ifail);
  y1 = dco::value(y);

  inc = 1.0;
  dco::derivative(y) += inc;
  ifail                                              = 0;
  dco::ga1s<double>::global_tape->sparse_interpret() = true;
  dco::ga1s<double>::global_tape->interpret_adjoint();

  /*     Get derivatives */
  dydx     = dco::derivative(x);
  dydx_exp = 1.0 / (1.0 + dco::value(x));

  cout.setf(ios::scientific, ios::floatfield);
  cout.precision(4);
  cout << "Input value of x             : " << dco::value(x) << endl;
  cout << "Output value of ln(1+x)      : " << y1 << endl;
  cout << "AD evaluated derivative      : " << dydx << endl;
  cout << "Directly computed derivative : " << dydx_exp << endl;

  /*     Remove computational data object and tape */
END:
  ifail = 0;
  nag::ad::x10ab(handle.data(), ifail);
  dco::ga1s<double>::tape_t::remove(dco::ga1s<double>::global_tape);
  return exit_status;
}