/* nag_2d_spline_fit_ts_scat (e02jdc) Example Program. * * Copyright 2013 Numerical Algorithms Group. * * Mark 24, 2013. */ #include #include #include #include #include #include #include static Integer generate_data(Integer *n, double **x, double **y, double **f, Integer *lsminp, Integer *lsmaxp, Integer *nxcels, Integer *nycels, Integer *lcoefs, double **coefs); static Integer handle_options(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 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. */ exit_status = generate_data(&n, &x, &y, &f, &lsminp, &lsmaxp, &nxcels, &nycels, &lcoefs, &coefs); if (exit_status != 0) goto END; /* Initialize the options arrays and set/get some options. */ exit_status = handle_options(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 mi and max values. */ nag_dmin_val(n, x, 1, &i, &pmin[0], &fail); if (fail.code == NE_NOERROR) nag_dmin_val(n, y, 1, &i, &pmin[1], &fail); if (fail.code == NE_NOERROR) nag_dmax_val(n, x, 1, &i, &pmax[0], &fail); if (fail.code == NE_NOERROR) nag_dmax_val(n, y, 1, &i, &pmax[1], &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_dmin_val (f16jpc) or 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) { /* 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. */ /* 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) nag_rand_basic(*n, state, *x, &fail); if (fail.code == NE_NOERROR) nag_rand_basic(*n, state, *y, &fail); if (fail.code != NE_NOERROR) { printf("Error in obtaining repeatable data from\nnag_rand_init_repeatable " "(g05kfc) and nag_rand_basic (g05sac).\nError message is:\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)); /* Grid size for the approximation. */ *nxcels = 6; *nycels = *nxcels; /* Identify the computation. */ printf("\n" "Computing the coefficients of a C^1 spline approximation to " "Franke's function\n" "Using a %ld by %ld grid\n", *nxcels, *nycels); /* Local-approximation control parameters. */ *lsminp = 3; *lsmaxp = 100; /* Set up the array to hold the computed spline coefficients. */ *lcoefs = (((*nxcels+2)*(*nycels+2)+1)/2)*10 + 1; if (!(*coefs = NAG_ALLOC(*lcoefs, double))) { printf("Allocation failure\n"); exit_status = -1; goto END; } END: return exit_status; } static Integer handle_options(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[3+1], optstr[80+1]; /* Nag Types */ Nag_VariableType optype; NagError fail; INIT_FAIL(fail); /* Set some non-default parameters for the local approximation method. */ sprintf(optstr, "Minimum Singular Value LPA = %16.9e", 1./32.); nag_fit_opt_set("Initialize = e02jdc", iopts, liopts, opts, lopts, &fail); if (fail.code == NE_NOERROR) nag_fit_opt_set(optstr, iopts, liopts, opts, lopts, &fail); if (fail.code == NE_NOERROR) nag_fit_opt_set("Polynomial Starting Degree = 3", iopts, liopts, opts, lopts, &fail); /* Set some non-default parameters for the global approximation method. */ if (fail.code == NE_NOERROR) nag_fit_opt_set("Averaged Spline = Yes", iopts, liopts, opts, lopts, &fail); /* As an example of how to get the value of an optional parameter, * display whether averaging of local approximations is in operation. */ if (fail.code != NE_NOERROR) { printf("Error from nag_fit_opt_set (e02zkc).\n%s\n", fail.message); exit_status = 1; goto END; } lcvalue = 3 + 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; }