/* nag_2d_spline_fit_ts_scat (e02jdc) Example Program.
*
* Copyright 2014 Numerical Algorithms Group.
*
* Mark 25, 2014.
*/
#include <math.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <nage02.h>
#include <nagf16.h>
#include <nagg05.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_2d_spline_fit_ts_scat (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_2d_spline_fit_ts_scat (e02jdc).
*/
nag_2d_spline_fit_ts_scat(n, x, y, f, lsminp, lsmaxp, nxcels, nycels, lcoefs,
coefs, iopts, opts, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_2d_spline_fit_ts_scat (e02jdc).\n%s\n",
fail.message);
exit_status = 1;
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_dmin_val (f16jpc) and
* nag_dmax_val (f16jnc) to obtain the min and max values.
*/
nag_dmin_val(n, x, 1, &i, &pmin[0], &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_dmin_val (f16jpc).\n%s\n",
fail.message);
exit_status = 1;
goto END;
}
nag_dmin_val(n, y, 1, &i, &pmin[1], &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_dmin_val (f16jpc).\n%s\n",
fail.message);
exit_status = 1;
goto END;
}
nag_dmax_val(n, x, 1, &i, &pmax[0], &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_dmax_val (f16jnc).\n%s\n",
fail.message);
exit_status = 1;
goto END;
}
nag_dmax_val(n, y, 1, &i, &pmax[1], &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_dmax_val (f16jnc).\n%s\n",
fail.message);
exit_status = 1;
goto END;
}
NAG_FREE(x);
NAG_FREE(y);
NAG_FREE(f);
/* 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(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("%ld%*[^\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_repeatable (g05kfc) and then generate the data using
* nag_rand_basic (g05sac).
*/
subid = 53;
seed[0] = 32958;
seed[1] = 39838;
seed[2] = 881818;
seed[3] = 45812;
lstate = mstate;
nag_rand_init_repeatable(Nag_WichmannHill_I, subid, seed, lseed, state,
&lstate, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_rand_init_repeatable (g05kfc).\n%s\n",
fail.message);
exit_status = 1;
goto END;
}
nag_rand_basic(*n, state, *x, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_rand_basic (g05sac).\n%s\n",
fail.message);
exit_status = 1;
goto END;
}
nag_rand_basic(*n, state, *y, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_rand_basic (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("%ld%*[^\n] ", gsmoothness);
/* Identify the computation. */
printf("\n"
"Computing the coefficients of a C^%ld spline approximation to"
" Franke's function\n"
"Using a %ld by %ld 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_2d_spline_ts_eval (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("%ld%*[^\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_2d_spline_ts_eval (e02jec). */
nag_2d_spline_ts_eval(nevalv, xevalv, yevalv, coefs, fevalv, iopts, opts,
&fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_2d_spline_ts_eval (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("%ld%ld%*[^\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; i++) xevalm[i] = ll_corner[0] + (double) i*h[0];
for (j = 0; j < nyeval; j++) yevalm[j] = ll_corner[1] + (double) j*h[1];
/* Evaluate using nag_2d_spline_ts_eval_rect (e02jfc). */
nag_2d_spline_ts_eval_rect(nxeval, nyeval, xevalm, yevalm, coefs, fevalm,
iopts, opts, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_2d_spline_ts_eval_rect (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;
}