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

NAG CL Interface Introduction
Example description
/* nag_opt_handle_disable (e04tcc) 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);
#ifdef __cplusplus
}
#endif

struct usr_data {
  Integer nres;
  double *t, *y;
};

int main(void) {

  /* This example shows how to use the NAG Optimization Modeling Suite
   * to edit a nonlinear least squares problem.
   * In particular, the outlier residuals are disabled without having to
   * redefine the full problem and comparing the solutions is easy.
   *
   * We try to fit the following model with 5 parameters:
   *   f(t) = at^2 + bt + c + d sin(omega t)
   * to some noisy data (30 points) with 2 outliers (point 10 and 20)
   * The data was generated using the values:
   *   a = 0.3
   *   b = 1.0
   *   c = 0.01
   *   d = 0.2
   *   omega = 5.0 */

  Integer nvar, nres, isparse, nnzrd, idx[2];
  double x[5];
  double rinfo[100], stats[100];
  double *rx;
  void *handle = 0;
  Integer exit_status = 0, i;
  struct usr_data ud;

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

  printf("nag_opt_handle_disable (e04tcc) Example Program Results\n\n");
  fflush(stdout);

  exit_status = 0;
  scanf(" %*[^\n]"); /* Skip heading in data file */

  /* Read the data dimension */
  scanf("%" NAG_IFMT " %*[^\n]", &nres);

  /* Allocate the user data memory */
  ud.nres = nres;
  if (!(ud.t = NAG_ALLOC(nres, double)) || !(ud.y = NAG_ALLOC(nres, double)) ||
      !(rx = NAG_ALLOC(nres, double))) {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }

  /* Read the data points */
  for (i = 0; i < nres; i++) {
    scanf("%lf", &ud.t[i]);
  }
  scanf("%*[^\n]");
  for (i = 0; i < nres; i++) {
    scanf("%lf", &ud.y[i]);
  }

  /* nag_opt_handle_init (e04rac).
   * Initialize the handle
   */
  nvar = 5;
  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 printing options
   */
  nag_opt_handle_opt_set(handle, "Print Level = 1", NAGERR_DEFAULT);
  nag_opt_handle_opt_set(handle, "Print Options = No", NAGERR_DEFAULT);
  nag_opt_handle_opt_set(handle, "Print Solution = X", NAGERR_DEFAULT);

  printf(" First solve the problem with the outliers\n\n");
  printf(" --------------------------------------------------------\n");
  fflush(stdout);

  /* nag_opt_handle_solve_bxnl (e04ggc)
   * Call the solver
   */
  comm.p = &ud;
  for (i = 0; i < nvar; i++) {
    x[i] = 1.0;
  }
  INIT_FAIL(fail);
  nag_opt_handle_solve_bxnl(handle, lsqfun, lsqgrd, NULLFN, NULLFN, NULLFN,
                            nvar, x, nres, rx, rinfo, stats, &comm, &fail);

  printf(" --------------------------------------------------------\n");
  printf("\n Now remove the outlier residuals from the problem handle\n\n");
  printf(" --------------------------------------------------------\n");
  fflush(stdout);
  /* nag_opt_handle_disable (e04tcc)
   * Disable the two outlier residuals */
  idx[0] = 10;
  idx[1] = 20;
  nag_opt_handle_disable(handle, "NLS", 2, idx, NAGERR_DEFAULT);

  /* Solve the problem again */
  for (i = 0; i < nvar; i++) {
    x[i] = 1.0;
  }
  INIT_FAIL(fail);
  nag_opt_handle_solve_bxnl(handle, lsqfun, lsqgrd, NULLFN, NULLFN, NULLFN,
                            nvar, x, nres, rx, rinfo, stats, &comm, &fail);

  printf(" --------------------------------------------------------\n");
  printf("\n  Assuming the outliers points are measured again\n");
  printf("  we can enable the residuals and adjust the values\n\n");
  printf(" --------------------------------------------------------\n");
  fflush(stdout);
  /* Fix the first variable to its known value of 0.3
   * enable the residuals and adjust the values in the data */
  nag_opt_handle_set_bound(handle, "variable", 1, 0.3, 0.3, NAGERR_DEFAULT);
  nag_opt_handle_enable(handle, "NLS", 2, idx, NAGERR_DEFAULT);
  ud.y[9] = -0.515629;
  ud.y[19] = 0.54920;

  /* Solve the problem again */
  for (i = 0; i < nvar; i++) {
    x[i] = 1.0;
  }
  INIT_FAIL(fail);
  nag_opt_handle_solve_bxnl(handle, lsqfun, lsqgrd, NULLFN, NULLFN, NULLFN,
                            nvar, x, nres, rx, rinfo, stats, &comm, &fail);

END:
  NAG_FREE(ud.t);
  NAG_FREE(ud.y);
  NAG_FREE(rx);
  /* nag_opt_handle_free (e04rzc).
   * Destroy the problem handle and deallocate all the memory. */
  if (handle)
    nag_opt_handle_free(&handle, NAGERR_DEFAULT);

  return exit_status;
}

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

  Integer i;
  struct usr_data *ud;

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

  *inform = 0;
  ud = (struct usr_data *)comm->p;

  /* Fill the residuals values */
  for (i = 0; i < nres; i++) {
    rx[i] = x[0] * pow(ud->t[i], 2) + x[1] * ud->t[i] + x[2] +
            x[3] * sin(x[4] * ud->t[i]) - ud->y[i];
  }
}

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

  Integer i;
  struct usr_data *ud;

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

  *inform = 0;
  ud = (struct usr_data *)comm->p;
  for (i = 0; i < nres; i++) {
    rdx[i * nvar + 0] = pow(ud->t[i], 2);
    rdx[i * nvar + 1] = ud->t[i];
    rdx[i * nvar + 2] = 1.0;
    rdx[i * nvar + 3] = sin(x[4] * ud->t[i]);
    rdx[i * nvar + 4] = x[3] * ud->t[i] * cos(x[4] * ud->t[i]);
  }
}