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

NAG CL Interface Introduction
Example description
/* nag_opt_handle_solve_nldf (e04gnc) Example Program.
 *
 * Copyright 2023 Numerical Algorithms Group.
 *
 * Mark 29.1, 2023.
 */

#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);
#ifdef __cplusplus
}
#endif

/* Model data
 * y = f(t;A,w) := A*t*sin(-w*t)
 * Data generated using A = 0.1 and w = 0.8 and added random noise
 * Residual 4, 12, 16 and 20 are outliers
 */
const double t[24] = {1.0E+00, 2.0E+00, 3.0E+00, 4.0E+00, 5.0E+00, 6.0E+00,
                      7.0E+00, 8.0E+00, 9.0E+00, 1.0E+01, 1.1E+01, 1.2E+01,
                      1.3E+01, 1.4E+01, 1.5E+01, 1.6E+01, 1.7E+01, 1.8E+01,
                      1.9E+01, 2.0E+01, 2.1E+01, 2.2E+01, 2.3E+01, 2.4E+01};
const double y[24] = { 0.0523, -0.1442, -0.0422,  1.8106,  0.3271,  0.4684, 
                       0.4593, -0.0169, -0.7811, -1.1356, -0.5343, -3.0043, 
                       1.1832,  1.5153,  0.7120, -2.2923, -1.4871, -1.7083,
                      -0.9936, -5.2873,  1.7555,  2.0642,  0.9499, -0.6234};

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

  Integer nvar, nres, nclin, nnza, isparse, nnzrd, x_idx, idlc;
  Integer irowa[2] = {1, 1};
  Integer icola[2] = {1, 2};
  double x[2] = {0.3, 0.7};
  double rinfo[100], stats[100];
  double blx[2] = {-1.0, 0.0};
  double bux[2] = {infbnd, 1.0};
  double bla[1] = {0.2};
  double bua[1] = {1.0};
  double a[2] = {1.0, 1.0};
  double *rx = 0; 
  void *handle = 0;
  Integer exit_status = 0;

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

  printf("nag_opt_handle_solve_nldf (e04gnc) Example Program Results\n\n");
  fflush(stdout);

  nvar = 2;
  nres = 24;
  /* Allocate memory */
  if (!(comm.user = NAG_ALLOC(2 * nres, double)) || 
      !(rx = NAG_ALLOC(nres, 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 least-squares loss function */
  nag_opt_handle_opt_set(handle, "NLDF Loss Function Type = L2",
                         NAGERR_DEFAULT);
  /* Reduce printed output info */
  nag_opt_handle_opt_set(handle, "Print Level = 1", NAGERR_DEFAULT);
  nag_opt_handle_opt_set(handle, "Print Options = No", NAGERR_DEFAULT);

  /* Define bounds */
  /* nag_opt_handle_set_simplebounds (e04rhc) */
  nag_opt_handle_set_simplebounds(handle, nvar, blx, bux, NAGERR_DEFAULT);

  printf(" First solve the problem using least squares loss function\n");
  printf(" --------------------------------------------------------\n");
  fflush(stdout);

  /* nag_opt_handle_set_linconstr (e04rjc)
   * Define linear constraints */
  nclin = 1;
  nnza = 2;
  idlc = 0;
  nag_opt_handle_set_linconstr(handle, nclin, bla, bua, nnza, irowa, icola, a,
                               &idlc, NAGERR_DEFAULT);

  /* nag_opt_handle_solve_nldf (e04gnc) 
   * Call the solver
   */
  INIT_FAIL(fail);
  nag_opt_handle_solve_nldf(handle, lsqfun, lsqgrd, NULLFN, NULLFN, NULLFN,
                            nvar, x, nres, rx, rinfo, stats, &comm, &fail);
  if (fail.code != NE_NOERROR && fail.code != NW_NOT_CONVERGED) {
    printf("Error from nag_opt_handle_solve_nldf (e04gnc).\n%s\n",
           fail.message);
    exit_status = 1;
  } else {
    /* Print the solution point */
    printf("\n Optimal parameters X (Least Squares):\n");
    printf("  x_idx    Value\n");
    for (x_idx = 0; x_idx < nvar; x_idx++) {
      printf("  %5" NAG_IFMT "   %7.2E\n", x_idx + 1, x[x_idx]);
    }
  }

  printf("\n Now solve the problem using SmoothL1 loss function\n");
  printf(" --------------------------------------------------------\n");
  fflush(stdout);

  /* nag_opt_handle_opt_set (e04zmc)
   * Set options
   */

  /* Use SmoothL1 loss function */
  nag_opt_handle_opt_set(handle, "NLDF Loss Function Type = SmoothL1",
                         NAGERR_DEFAULT);

  /* nag_opt_handle_solve_nldf (e04gnc) 
   * Call the solver
   */
  INIT_FAIL(fail);
  x[0] = 0.3;
  x[1] = 0.7;
  nag_opt_handle_solve_nldf(handle, lsqfun, lsqgrd, NULLFN, NULLFN, NULLFN,
                            nvar, x, nres, rx, rinfo, stats, &comm, &fail);
  if (fail.code != NE_NOERROR && fail.code != NW_NOT_CONVERGED) {
    printf("Error from nag_opt_handle_solve_nldf (e04gnc).\n%s\n",
           fail.message);
    exit_status = 1;
  } else {
    /* Print the solution point */
    /* Print the solution point */
    printf("\n Optimal parameters X (SmoothL1):\n");
    printf("  x_idx    Value\n");
    for (x_idx = 0; x_idx < nvar; x_idx++) {
      printf("  %5" NAG_IFMT "   %7.2E\n", x_idx + 1, x[x_idx]);
    }
  }

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

  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 != 2 || 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] * comm->user[i] * sin(-x[1] * 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 != 2 || nres != 24 || nnzrd != nres * nvar) {
    *inform = -1;
    return;
  }

  *inform = 0;

  for (i = 0; i < nres; i++) {
    rdx[i * nvar + 0] = -comm->user[i] * sin(-x[1] * comm->user[i]);
    rdx[i * nvar + 1] = comm->user[i] * comm->user[i] * x[0] 
                        * cos(x[1] * comm->user[i]);
  }
}