/* nag_opt_handle_solve_bounds_foas (e04kfc) Example Program.
 *
 * Copyright 2024 Numerical Algorithms Group.
 *
 * Mark 30.3, 2024.
 */

/*
 * NLP example: Nonlinear objective + box constraint
 */

#include <math.h>
#include <nag.h>
#include <stdio.h>

#ifdef __cplusplus
extern "C" {
#endif

static void NAG_CALL objfun(Integer nvar, const double x[], double *fx,
                            Integer *inform, Nag_Comm *comm);
static void NAG_CALL objgrd(Integer nvar, const double x[], Integer nnzfd,
                            double fdx[], Integer *inform, Nag_Comm *comm);
static void NAG_CALL monit(Integer nvar, const double x[], Integer *inform,
                           const double rinfo[], const double stats[],
                           Nag_Comm *comm);

#ifdef __cplusplus
}
#endif
static void NAG_CALL objfun(Integer nvar, const double x[], double *fx,
                            Integer *inform, Nag_Comm *comm) {
  *fx = pow(1.0 - x[0], 2.0) + 100.0 * pow(x[1] - pow(x[0], 2.0), 2.0);
}

static void NAG_CALL objgrd(Integer nvar, const double x[], Integer nnzfd,
                            double fdx[], Integer *inform, Nag_Comm *comm) {
  fdx[0] = 2.0 * x[0] - 400.0 * x[0] * (x[1] - pow(x[0], 2.0)) - 2.0;
  fdx[1] = 200.0 * (x[1] - pow(x[0], 2.0));
}

static void NAG_CALL monit(Integer nvar, const double x[], Integer *inform,
                           const double rinfo[], const double stats[],
                           Nag_Comm *comm) {
  return;
}

int main(void) {
/* Define problem size */
#define NVAR 2

  /* Scalars */
  Integer nvar = NVAR;
  Integer nu = 2 * nvar;

  /* Arrays */
  /* Lower bounds */
  double blx[NVAR] = {-1.0, -2.0};
  /* Upper bounds */
  double bux[NVAR] = {0.8, 2.0};
  /* Indexes for the gradient entries */
  Integer idxfd[NVAR] = {1, 2};
  /* Initial guess x0 */
  double x[NVAR] = {-1.5, 1.9};
  /* space holder for the Lagrange multiplyers */
  double u[2 * NVAR];
  double rinfo[100], stats[100];
  void *handle = NULL;
  char opt[80];

  /* Nag Types */
  Nag_Comm comm;
  NagError fail;
  Nag_FileID mon_out;

  /* nag_file_open (x04acc).
   *  Open unit number for reading, writing or appending, and
   *  associate unit with named file
   */
  nag_file_open("e04kf_cmon.txt", 1, &mon_out, NAGERR_DEFAULT);

  /* nag_file_line_write (x04bac).
   * Write formatted record to external file
   */
  Vprintf("nag_opt_handle_solve_bounds_foas (e04kfc) "
          "Example Program Results\n");
  fflush(stdout);

  /* nag_opt_handle_init (e04rac).
   * Initialize an empty problem handle with NVAR variables.
   */
  nag_opt_handle_init(&handle, nvar, NAGERR_DEFAULT);

  /* nag_opt_handle_set_simplebounds (e04rhc).
   * Define bounds on the variables
   */
  nag_opt_handle_set_simplebounds(handle, nvar, blx, bux, NAGERR_DEFAULT);

  /* nag_opt_handle_set_nlnobj (e04rgc).
   * Define nonlinear objective
   */
  nag_opt_handle_set_nlnobj(handle, nvar, idxfd, NAGERR_DEFAULT);

  /* nag_opt_handle_opt_set (e04zmc).
   * Set optional arguments of the solver
   */
  nag_opt_handle_opt_set(handle, "FOAS Print Frequency = 5", NAGERR_DEFAULT);

  nag_opt_handle_opt_set(handle, "Print Solution = yes", NAGERR_DEFAULT);

  sprintf(opt, "Monitoring File = %" NAG_IFMT, mon_out);
  nag_opt_handle_opt_set(handle, opt, NAGERR_DEFAULT);

  nag_opt_handle_opt_set(handle, "Monitoring Level = 3", NAGERR_DEFAULT);

  nag_opt_handle_opt_set(handle, "Print Level = 1", NAGERR_DEFAULT);

  /* Solver the problem
   * nag_opt_handle_solve_bounds_foas (e04kfc).
   */
  SET_FAIL(fail);
  nag_opt_handle_solve_bounds_foas(handle, objfun, objgrd, monit, nvar, x,
                                   rinfo, stats, &comm, &fail);

  /* Print objective value at solution */
  if (fail.code == NE_NOERROR || fail.code == NW_NOT_CONVERGED) {
    Vprintf("\n\n  Solution found.\n"
            "  Objective function value at solution: %8.1e\n",
            rinfo[0]);
    nag_opt_handle_set_get_real(handle, "Dual Variables", 1, &nu, u,
                                NAGERR_DEFAULT);
    /* Print gradient vector FDX */
    Vprintf("  Gradient at solution:                 %8.1e %8.1e\n\n",
            u[0] - u[1], u[2] - u[3]);
    /* Print Lagrange multipliers */
    Vprintf("  Estimated Lagrange multipliers: blx   %8.1e %8.1e\n", u[0],
            u[2]);
    Vprintf("  Estimated Lagrange multipliers: bux   %8.1e %8.1e\n\n", u[1],
            u[3]);
  }

  nag_file_close(mon_out, NAGERR_DEFAULT);

  /* Clean up */
  if (handle)
    /* nag_opt_handle_free (e04rzc).
     * Destroy the problem handle and deallocate all the memory used
     */
    nag_opt_handle_free(&handle, NAGERR_DEFAULT);

  return (0);
}