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

NAG CL Interface Introduction
Example description
/* nag_fit_dim2_spline_ts_sctr (e02jdc) Example Program.
 *
 * Copyright 2024 Numerical Algorithms Group.
 *
 * Mark 30.1, 2024.
 */

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

static Integer generate_data(Integer *n, double **x, double **y, double **f,
                             Integer *lsminp, Integer *lsmaxp, Integer *nxcels,
                             Integer *nycels, Integer *lcoefs, double **coefs,
                             Integer *gsmoothness);
static Integer handle_options(Integer gsmoothness, Integer iopts[],
                              const Integer liopts, double opts[],
                              const Integer lopts);
static Integer evaluate_at_vector(const double coefs[], const Integer iopts[],
                                  const double opts[], const double pmin[],
                                  const double pmax[]);
static Integer evaluate_on_mesh(const double coefs[], const Integer iopts[],
                                const double opts[], const double pmin[],
                                const double pmax[]);

int main(void) {
  /* Scalars */
  Integer exit_status = 0;
  Integer liopts = 100, lopts = 100;
  Integer gsmoothness, i, lcoefs, lsmaxp, lsminp, n, nxcels, nycels;
  /* Arrays */
  double *coefs = 0, *f = 0, *x = 0, *y = 0;
  double opts[100], pmax[2], pmin[2];
  Integer iopts[100];
  /* Nag Types */
  NagError fail;

  INIT_FAIL(fail);

  printf("nag_fit_dim2_spline_ts_sctr (e02jdc) Example Program Results\n");

  /* Generate the data to fit and set the compulsory algorithmic control
   * parameters.
   */
  exit_status = generate_data(&n, &x, &y, &f, &lsminp, &lsmaxp, &nxcels,
                              &nycels, &lcoefs, &coefs, &gsmoothness);
  if (exit_status != 0)
    goto END;

  /* Initialize the options arrays and set/get some options. */
  exit_status = handle_options(gsmoothness, iopts, liopts, opts, lopts);
  if (exit_status != 0)
    goto END;

  /* Compute the spline coefficients using
   * nag_fit_dim2_spline_ts_sctr (e02jdc).
   */
  nag_fit_dim2_spline_ts_sctr(n, x, y, f, lsminp, lsmaxp, nxcels, nycels,
                              lcoefs, coefs, iopts, opts, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_fit_dim2_spline_ts_sctr (e02jdc).\n%s\n",
           fail.message);
    exit_status = 1;
    NAG_FREE(x);
    NAG_FREE(y);
    NAG_FREE(f);
    goto END;
  }

  /* pmin and pmax form the bounding box of the spline. We must not attempt to
   * evaluate the spline outside this box. Use nag_blast_dmin_val (f16jpc) and
   * nag_blast_dmax_val (f16jnc) to obtain the min and max values.
   */
  nag_blast_dmin_val(n, x, 1, &i, &pmin[0], &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_blast_dmin_val (f16jpc).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }
  nag_blast_dmin_val(n, y, 1, &i, &pmin[1], &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_blast_dmin_val (f16jpc).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }
  nag_blast_dmax_val(n, x, 1, &i, &pmax[0], &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_blast_dmax_val (f16jnc).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }
  nag_blast_dmax_val(n, y, 1, &i, &pmax[1], &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_blast_dmax_val (f16jnc).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }

  /* Evaluate the approximation at a vector of values. */
  exit_status = evaluate_at_vector(coefs, iopts, opts, pmin, pmax);
  if (exit_status != 0)
    goto END;

  /* Evaluate the approximation on a mesh. */
  exit_status = evaluate_on_mesh(coefs, iopts, opts, pmin, pmax);
  if (exit_status != 0)
    goto END;

END:

  NAG_FREE(x);
  NAG_FREE(y);
  NAG_FREE(f);
  NAG_FREE(coefs);

  return exit_status;
}

static Integer generate_data(Integer *n, double **x, double **y, double **f,
                             Integer *lsminp, Integer *lsmaxp, Integer *nxcels,
                             Integer *nycels, Integer *lcoefs, double **coefs,
                             Integer *gsmoothness) {
  /* Reads n from a data file and then generates an x and a y vector of n
   * pseudorandom uniformly distributed values on (0,1]. These are passed
   * to the bivariate function of R. Franke to create the data set to fit.
   * The remaining input data for the fitter are set to suitable values for
   * this problem, as discussed by Davydov and Zeilfelder.
   * Reads the global smoothing level from a data file. This value determines
   * the minimum required length of the array of spline coefficients, coefs.
   */

  /* Scalars */
  Integer exit_status = 0;
  Integer lseed = 4, mstate = 21;
  Integer i, lstate, subid;
  /* Arrays */
  Integer seed[4], state[21];
  /* Nag Types */
  NagError fail;

  INIT_FAIL(fail);

  /* Read the size of the data set to be generated and fitted.
   * (Skip the heading in the data file.)
   */
  scanf("%*[^\n]");
  scanf("%" NAG_IFMT "%*[^\n] ", n);

  if (!(*x = NAG_ALLOC(*n, double)) || !(*y = NAG_ALLOC(*n, double)) ||
      !(*f = NAG_ALLOC(*n, double))) {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }

  /* Initialize the random number generator using
   * nag_rand_init_repeat (g05kfc) and then generate the data using
   * nag_rand_dist_uniform01 (g05sac).
   */
  subid = 53;
  seed[0] = 32958;
  seed[1] = 39838;
  seed[2] = 881818;
  seed[3] = 45812;
  lstate = mstate;
  nag_rand_init_repeat(Nag_WichmannHill_I, subid, seed, lseed, state, &lstate,
                       &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_rand_init_repeat (g05kfc).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }
  nag_rand_dist_uniform01(*n, state, *x, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_rand_dist_uniform01 (g05sac).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }
  nag_rand_dist_uniform01(*n, state, *y, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_rand_dist_uniform01 (g05sac).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }

  for (i = 0; i < *n; i++)
    (*f)[i] =
        0.75 *
            exp(-(pow((9. * (*x)[i] - 2.), 2) + pow((9. * (*y)[i] - 2.), 2)) /
                4.) +
        0.75 * exp(-pow((9. * (*x)[i] + 1.), 2) / 49. -
                   (9. * (*y)[i] + 1.) / 10.) +
        0.5 * exp(-(pow((9. * (*x)[i] - 7.), 2) + pow((9. * (*y)[i] - 3.), 2)) /
                  4.) -
        0.2 * exp(-pow((9. * (*x)[i] - 4.), 2) - pow((9. * (*y)[i] - 7.), 2));

  /* Set the grid size for the approximation. */
  *nxcels = 6;
  *nycels = *nxcels;

  /* Read the required level of global smoothing. */
  scanf("%" NAG_IFMT "%*[^\n] ", gsmoothness);

  /* Identify the computation. */
  printf("\n"
         "Computing the coefficients of a C^%" NAG_IFMT
         " spline approximation to"
         " Franke's function\n"
         "Using a %" NAG_IFMT " by %" NAG_IFMT " grid\n",
         *gsmoothness, *nxcels, *nycels);

  /* Set the local-approximation control parameters. */
  *lsminp = 3;
  *lsmaxp = 100;

  /* Set up the array to hold the computed spline coefficients. */
  switch (*gsmoothness) {
  case 1:
    *lcoefs = (((*nxcels + 2) * (*nycels + 2) + 1) / 2) * 10 + 1;
    break;
  case 2:
    *lcoefs = 28 * (*nxcels + 2) * (*nycels + 2) * 4 + 1;
    break;
  default:
    *lcoefs = 0;
  }

  if (!(*coefs = NAG_ALLOC(*lcoefs, double))) {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }

END:

  return exit_status;
}

static Integer handle_options(Integer gsmoothness, Integer iopts[],
                              const Integer liopts, double opts[],
                              const Integer lopts) {
  /* Auxiliary routine for initializing the options arrays and
   * for demonstrating how to set and get optional parameters using
   * nag_fit_opt_set (e02zkc) and nag_fit_opt_get (e02zlc) respectively.
   */

  /* Scalars */
  Integer exit_status = 0;
  double rvalue;
  Integer ivalue, lcvalue;
  /* Arrays */
  char cvalue[16 + 1], optstr[80 + 1], supersmooth[9 + 1];
  /* Nag Types */
  Nag_Boolean supersmooth_enum;
  Nag_VariableType optype;
  NagError fail;

  INIT_FAIL(fail);

  nag_fit_opt_set("Initialize = e02jdc", iopts, liopts, opts, lopts, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_fit_opt_set (e02zkc).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }

  /* Configure the global approximation method. */
  sprintf(optstr, "Global Smoothing Level = %1" NAG_IFMT, gsmoothness);
  nag_fit_opt_set(optstr, iopts, liopts, opts, lopts, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_fit_opt_set (e02zkc).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }

  /* If C^2 smoothing is requested, compute the spline using additional
   * super-smoothness constraints?
   * (The default is 'No'.)
   */
  scanf("%9s%*[^\n] ", supersmooth);
  supersmooth_enum = (Nag_Boolean)nag_enum_name_to_value(supersmooth);
  if (gsmoothness == 2 && supersmooth_enum == Nag_TRUE) {
    nag_fit_opt_set("Supersmooth C2 = Yes", iopts, liopts, opts, lopts, &fail);
    if (fail.code != NE_NOERROR) {
      printf("Error from nag_fit_opt_set (e02zkc).\n%s\n", fail.message);
      exit_status = 1;
      goto END;
    }
  }

  nag_fit_opt_set("Averaged Spline = Yes", iopts, liopts, opts, lopts, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_fit_opt_set (e02zkc).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }

  /* Configure the local approximation method.
   * (The default is 'Polynomial'.)
   */
  nag_fit_opt_set("Local Method = Polynomial", iopts, liopts, opts, lopts,
                  &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_fit_opt_set (e02zkc).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }
  sprintf(optstr, "Minimum Singular Value LPA = %16.9e", 1. / 32.);
  nag_fit_opt_set(optstr, iopts, liopts, opts, lopts, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_fit_opt_set (e02zkc).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }

  switch (gsmoothness) {
  case 1:
    sprintf(optstr, "Polynomial Starting Degree = 3");
    break;
  case 2:
    if (supersmooth_enum == Nag_TRUE) {
      /* We can benefit from starting with local polynomials of greater
       * degree than with regular C^2 smoothing.
       */
      printf("Using super-smoothing\n");
      sprintf(optstr, "Polynomial Starting Degree = 6");
    } else
      sprintf(optstr, "Polynomial Starting Degree = 5");
    break;
  }
  nag_fit_opt_set(optstr, iopts, liopts, opts, lopts, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_fit_opt_set (e02zkc).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }

  /* As an example of how to get the value of an optional parameter,
   * display whether averaging of local approximations is in operation.
   */
  lcvalue = 16 + 1;
  nag_fit_opt_get("Averaged Spline", &ivalue, &rvalue, cvalue, lcvalue, &optype,
                  iopts, opts, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_fit_opt_get (e02zlc).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }

  if (!strcmp(cvalue, "YES"))
    printf("Using an averaged local approximation\n");

END:
  return exit_status;
}

static Integer evaluate_at_vector(const double coefs[], const Integer iopts[],
                                  const double opts[], const double pmin[],
                                  const double pmax[]) {
  /* Evaluates the approximation at a vector of values using
   * nag_fit_dim2_spline_ts_evalv (e02jec).
   */

  /* Scalars */
  Integer exit_status = 0;
  Integer i, nevalv;
  /* Arrays */
  double *fevalv = 0, *xevalv = 0, *yevalv = 0;
  /* Nag Types */
  NagError fail;

  INIT_FAIL(fail);
  scanf("%" NAG_IFMT "%*[^\n] ", &nevalv);

  if (!(xevalv = NAG_ALLOC(nevalv, double)) ||
      !(yevalv = NAG_ALLOC(nevalv, double)) ||
      !(fevalv = NAG_ALLOC(nevalv, double))) {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }
  for (i = 0; i < nevalv; i++)
    scanf("%lf%lf%*[^\n]", &xevalv[i], &yevalv[i]);

  /* Force the points to be within the bounding box of the spline. */
  for (i = 0; i < nevalv; i++) {
    xevalv[i] = MAX(xevalv[i], pmin[0]);
    xevalv[i] = MIN(xevalv[i], pmax[0]);
    yevalv[i] = MAX(yevalv[i], pmin[1]);
    yevalv[i] = MIN(yevalv[i], pmax[1]);
  }

  /* Evaluate using nag_fit_dim2_spline_ts_evalv (e02jec). */
  nag_fit_dim2_spline_ts_evalv(nevalv, xevalv, yevalv, coefs, fevalv, iopts,
                               opts, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_fit_dim2_spline_ts_evalv (e02jec).\n%s\n",
           fail.message);
    exit_status = 1;
    goto END;
  }

  printf("\nValues of computed spline at (x_i,y_i):\n"
         "\n%12s%12s%12s\n",
         "x_i", "y_i", "f(x_i,y_i)");

  for (i = 0; i < nevalv; i++)
    printf("%12.2f%12.2f%12.2f\n", xevalv[i], yevalv[i], fevalv[i]);

END:

  NAG_FREE(xevalv);
  NAG_FREE(yevalv);
  NAG_FREE(fevalv);

  return exit_status;
}

static Integer evaluate_on_mesh(const double coefs[], const Integer iopts[],
                                const double opts[], const double pmin[],
                                const double pmax[]) {
  /* Evaluates the approximation on a mesh of n_x * n_y values. */

  /* Scalars */
  Integer exit_status = 0;
  Integer i, j, nxeval, nyeval;
  /* Arrays */
  char print_mesh[9 + 1];
  double *fevalm = 0, *xevalm = 0, *yevalm = 0;
  double h[2], ll_corner[2], ur_corner[2];
  /* Nag Types */
  Nag_Boolean print_mesh_enum;
  NagError fail;

  INIT_FAIL(fail);
  scanf("%" NAG_IFMT "%" NAG_IFMT "%*[^\n] ", &nxeval, &nyeval);

  if (!(xevalm = NAG_ALLOC(nxeval, double)) ||
      !(yevalm = NAG_ALLOC(nyeval, double)) ||
      !(fevalm = NAG_ALLOC(nxeval * nyeval, double))) {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }

  /* Define the mesh by its lower-left and upper-right corners, which must
   * lie within the bounding box of the spline.
   */
  scanf("%lf%lf%*[^\n] ", &ll_corner[0], &ll_corner[1]);
  scanf("%lf%lf%*[^\n] ", &ur_corner[0], &ur_corner[1]);

  for (i = 0; i < 2; i++) {
    ll_corner[i] = MAX(ll_corner[i], pmin[i]);
    ur_corner[i] = MIN(ur_corner[i], pmax[i]);
  }

  /* Set the mesh spacing and the evaluation points. */
  h[0] = (ur_corner[0] - ll_corner[0]) / (double)(nxeval - 1);
  h[1] = (ur_corner[1] - ll_corner[1]) / (double)(nyeval - 1);

  for (i = 0; i < nxeval - 1; i++)
    xevalm[i] = MIN(ll_corner[0] + (double)i * h[0], ur_corner[0]);
  for (j = 0; j < nyeval - 1; j++)
    yevalm[j] = MIN(ll_corner[1] + (double)j * h[1], ur_corner[1]);
  xevalm[nxeval - 1] = ur_corner[0];
  yevalm[nxeval - 1] = ur_corner[1];

  /* Evaluate using nag_fit_dim2_spline_ts_evalm (e02jfc). */
  nag_fit_dim2_spline_ts_evalm(nxeval, nyeval, xevalm, yevalm, coefs, fevalm,
                               iopts, opts, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_fit_dim2_spline_ts_evalm (e02jfc).\n%s\n",
           fail.message);
    exit_status = 1;
    goto END;
  }

  /* Output the computed function values? */
  scanf("%9s%*[^\n] ", print_mesh);
  print_mesh_enum = (Nag_Boolean)nag_enum_name_to_value(print_mesh);

  printf("\n");
  if (!print_mesh_enum) {
    printf("Outputting of the function values on the mesh is disabled\n");
  } else {
    printf("Values of computed spline at (x_i,y_j):\n"
           "\n%12s%12s%12s\n",
           "x_i", "y_j", "f(x_i,y_j)");
    for (j = 0; j < nyeval; j++)
      for (i = 0; i < nxeval; i++)
        printf("%12.2f%12.2f%12.2f\n", xevalm[i], yevalm[j],
               fevalm[j * nxeval + i]);
  }

END:

  NAG_FREE(xevalm);
  NAG_FREE(yevalm);
  NAG_FREE(fevalm);

  return exit_status;
}