/* nag_opt_handle_solve_bxnl (e04ggc) Example Program.
 *
 * Copyright 2024 Numerical Algorithms Group.
 *
 * Mark 30.0, 2024.
 */

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

#ifdef __cplusplus
extern "C" {
#endif
static void NAG_CALL lsqfun(Integer nvar, const double x[], Integer nres,
                            double rx[], Integer *inform, Nag_Comm *comm);
static void NAG_CALL lsqgrd(Integer nvar, const double x[], Integer nres,
                            Integer nnzrd, double rdx[], Integer *inform,
                            Nag_Comm *comm);
static void NAG_CALL lsqhes(Integer nvar, const double x[], Integer nres,
                            const double lambda[], double hx[], Integer *inform,
                            Nag_Comm *comm);
#ifdef __cplusplus
}
#endif

/* Model data */
const double t[24] = {0.00E+0, 5.00E-2, 1.00E-1, 1.50E-1, 2.00E-1, 2.50E-1,
                      3.00E-1, 3.50E-1, 4.00E-1, 4.50E-1, 5.00E-1, 5.50E-1,
                      6.00E-1, 6.50E-1, 7.00E-1, 7.50E-1, 8.00E-1, 8.50E-1,
                      9.00E-1, 9.50E-1, 1.00E+0, 1.05E+0, 1.10E+0, 1.15E+0};
const double y[24] = {2.5134, 2.0443, 1.6684, 1.3664, 1.1232, 0.9269,
                      0.7679, 0.6389, 0.5338, 0.4479, 0.3776, 0.3197,
                      0.2720, 0.2325, 0.1997, 0.1723, 0.1493, 0.1301,
                      0.1138, 0.1000, 0.0883, 0.0783, 0.0698, 0.0624};

int main(void) {
  const double infbnd = 1.0e20;

  Integer nvar, nres, isparse, nnzrd;
  double x[6] = {1.2, 0.3, 5.6, 5.5, 6.5, 7.6};
  double rinfo[100], stats[100];
  double *rx, *blx, *bux, *z;
  void *handle = 0;
  Integer exit_status = 0;

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

  printf("nag_opt_handle_solve_bxnl (e04ggc) Example Program Results\n\n");
  fflush(stdout);

  nvar = 6;
  nres = 24;

  /* Allocate memory */
  if (!(comm.user = NAG_ALLOC(2 * nres, double)) ||
      !(blx = NAG_ALLOC(nvar, double)) || !(bux = NAG_ALLOC(nvar, double)) ||
      !(rx = NAG_ALLOC(nres, double)) || !(z = NAG_ALLOC(nvar, double))) {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }

  /* Fill the problem data structure
   * copy t[]
   */
  memcpy(&comm.user[0], t, 24 * sizeof(double));
  /* copy y[] = */
  memcpy(&comm.user[24], y, 24 * sizeof(double));

  /* nag_opt_handle_init (e04rac).
   * Initialize the handle
   */
  nag_opt_handle_init(&handle, nvar, NAGERR_DEFAULT);

  /* nag_opt_handle_set_nlnls (e04rmc)
   * Define residuals structure, isparse=0 means the residual structure is
   * dense => irowrd and icolrd arguments can be NULL
   */
  isparse = 0;
  nnzrd = 1;
  nag_opt_handle_set_nlnls(handle, nres, isparse, nnzrd, NULL, NULL,
                           NAGERR_DEFAULT);

  /* nag_opt_handle_opt_set (e04zmc)
   * Set options
   */

  /* Use second order information to build the TR model */
  nag_opt_handle_opt_set(handle, "BXNL Use Second Derivatives = Yes",
                         NAGERR_DEFAULT);
  /* Choose the Gauss-Newton method to solve the TR problem */
  nag_opt_handle_opt_set(handle, "BXNL Model = Gauss-Newton", NAGERR_DEFAULT);
  /* Add regularization */
  nag_opt_handle_opt_set(handle, "BXNL Glob Method = Reg", NAGERR_DEFAULT);
  /* Reduce printed output info */
  nag_opt_handle_opt_set(handle, "Print Level = 1", NAGERR_DEFAULT);

  /* Define bounds */
  blx[0] = 0.0;
  bux[0] = 1.0;
  blx[1] = -1.0;
  bux[1] = infbnd;
  blx[2] = -1.0;
  bux[2] = infbnd;
  blx[3] = -1.0;
  bux[3] = infbnd;
  blx[4] = -1.0;
  bux[4] = 1.0;
  blx[5] = -1.0;
  bux[5] = 10.0;
  /* nag_opt_handle_set_simplebounds (e04rhc) */
  nag_opt_handle_set_simplebounds(handle, nvar, blx, bux, NAGERR_DEFAULT);

  /* nag_opt_handle_solve_bxnl (e04ggc)
   * Call the solver
   */
  INIT_FAIL(fail);
  nag_opt_handle_solve_bxnl(handle, lsqfun, lsqgrd, lsqhes, NULLFN, NULLFN,
                            nvar, x, nres, rx, rinfo, stats, &comm, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_opt_handle_solve_bxnl (e04ggc).\n%s\n",
           fail.message);
    exit_status = 1;
  } else {
    /* Retrieve the solution point */
    nag_opt_handle_set_get_real(handle, "X", 1, &nvar, z, NAGERR_DEFAULT);
    printf("\nSolver stored solution iterate in the handle\n");
    printf("X: %8.2e %8.2e %8.2e %8.2e %8.2e %8.2e\n", z[0], z[1], z[2], z[3],
           z[4], z[5]);
    NAG_FREE(z);
  }

END:
  /* Free allocated memory */
  if (handle)
    /* nag_opt_handle_free (e04rzc).
     * Free the problem handle and deallocate all the memory used
     */
    nag_opt_handle_free(&handle, NAGERR_DEFAULT);
  if (comm.user)
    NAG_FREE(comm.user);
  NAG_FREE(rx);
  NAG_FREE(blx);
  NAG_FREE(bux);

  return exit_status;
}

static void NAG_CALL lsqfun(Integer nvar, const double x[], Integer nres,
                            double rx[], Integer *inform, Nag_Comm *comm) {

  int i;

  /* Interrupt the solver if the comm structure is not correctly initialized */
  if (!comm || !(comm->user)) {
    *inform = -1;
    return;
  }

  /* Interrupt the solver if the data does not correspond to the problem */
  if (nvar != 6 || nres != 24) {
    *inform = -1;
    return;
  }

  *inform = 0;

  /* Fill the residuals values */
  for (i = 0; i < nres; i++) {
    rx[i] = comm->user[nres + i];
    rx[i] = rx[i] - x[0] * exp(-x[1] * comm->user[i]) -
            x[2] * exp(-x[3] * comm->user[i]) -
            x[4] * exp(-x[5] * comm->user[i]);
  }
}

static void NAG_CALL lsqgrd(Integer nvar, const double x[], Integer nres,
                            Integer nnzrd, double rdx[], Integer *inform,
                            Nag_Comm *comm) {
  int i;

  /* Interrupt the solver if the comm structure is not correctly initialized */
  if (!comm || !(comm->user)) {
    *inform = -1;
    return;
  }

  /* Interrupt the solver if the data does not correspond to the problem */
  if (nvar != 6 || nres != 24 || nnzrd != nres * nvar) {
    *inform = -1;
    return;
  }

  *inform = 0;

  for (i = 0; i < nres; i++) {
    rdx[i * nvar + 0] = -exp(-x[1] * comm->user[i]);
    rdx[i * nvar + 1] = +comm->user[i] * x[0] * exp(-x[1] * comm->user[i]);
    rdx[i * nvar + 2] = -exp(-x[3] * comm->user[i]);
    rdx[i * nvar + 3] = +comm->user[i] * x[2] * exp(-x[3] * comm->user[i]);
    rdx[i * nvar + 4] = -exp(-x[5] * comm->user[i]);
    rdx[i * nvar + 5] = +comm->user[i] * x[4] * exp(-x[5] * comm->user[i]);
  }
}

static void NAG_CALL lsqhes(Integer nvar, const double x[], Integer nres,
                            const double lambda[], double hx[], Integer *inform,
                            Nag_Comm *comm) {
  double sum21 = 0.0, sum22 = 0.0, sum43 = 0.0, sum44 = 0.0, sum65 = 0.0,
         sum66 = 0.0;
  int i;

  /* Interrupt the solver if the comm structure is not correctly initialized */
  if (!comm || !(comm->user)) {
    *inform = -1;
    return;
  }

  /* Interrupt the solver if the data does not correspond to the problem */
  if (nvar != 6 || nres != 24) {
    *inform = -1;
    return;
  }

  *inform = 0;

  for (i = 0; i < nvar * nvar; i++)
    hx[i] = 0.0;

  for (i = 0; i < nres; i++) {
    sum21 = sum21 + (lambda[i] * comm->user[i] * exp(-x[1] * comm->user[i]));
    sum22 = sum22 + (-lambda[i] * pow(comm->user[i], 2) * x[0] *
                     exp(-x[1] * comm->user[i]));
    sum43 = sum43 + (lambda[i] * comm->user[i] * exp(-x[3] * comm->user[i]));
    sum44 = sum44 + (-lambda[i] * pow(comm->user[i], 2) * x[2] *
                     exp(-x[3] * comm->user[i]));
    sum65 = sum65 + (lambda[i] * comm->user[i] * exp(-x[5] * comm->user[i]));
    sum66 = sum66 + (-lambda[i] * pow(comm->user[i], 2) * x[4] *
                     exp(-x[5] * comm->user[i]));
  }

  hx[1 + 0 * nvar] = sum21;
  hx[0 + 1 * nvar] = hx[1 + 0 * nvar];
  hx[1 + 1 * nvar] = sum22;
  hx[3 + 2 * nvar] = sum43;
  hx[2 + 3 * nvar] = hx[3 + 2 * nvar];
  hx[3 + 3 * nvar] = sum44;
  hx[5 + 4 * nvar] = sum65;
  hx[4 + 5 * nvar] = hx[5 + 4 * nvar];
  hx[5 + 5 * nvar] = sum66;
}