/* nag_fit_dim2_spline_ts_sctr (e02jdc) Example Program.
*
* Copyright 2022 Numerical Algorithms Group.
*
* Mark 28.6, 2022.
*/
#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;
}