User's Guide
|
Contents
Example 5: Constrained Surface Fit
RBF and LSD data structure declarations
This document describes functions that compute and manipulate Radial Basis Functions (RBFs). For our purposes a RBF is a function f mapping Rn to Rm. It takes the form:
f (x) = a0 + å {i=1...k} ai j(pi ,||S (x - ci)||)
where ai Î Rm, x, ci Î Rn, S is a diagonal matrix, and the double bars denote the Euclidean norm. Several support functions are available including those which read least squares data and those which print information concerning the radial basis functions as well as relevant information on least squares fits.
The functions in RBFpack are designed to solve the least squares approximation problem:
min å {j=1... nobs} wj ||{yj - [a0 + å {i=1...k} ai j (pi ,||S (x - ci)||)]}||2
where yj Î Rm and wj Î R+ and the double bars denote the Euclidean norm. The upper limit on the first sum is nobs which stands for number of observations. We use the CLAPACK function dgelsx to solve the least squares problems. When m = 1 this has the form
min å {j=1... nobs} wj {yj - [a0 + å {i=1... k} ai j(pi ,||S (x - ci)||)] }2
If the minimization takes place only with respect to the linear parameters ai Î Rm, then this is a standard linear least squares problem and it can be solved very efficiently.
Standard choices for the radial function j(p, r) are:
exp(-{pr}2), sqrt(p2 + r2), 1/sqrt(p2 + r2), {(1 - pr)+}4 (1 + 4 pr)
These radial functions are referred to as the Gaussian, Hardy multiquadric, the inverse Hardy multiquadric, and the radial polynomial (of degree 5). For convenience, we store the information for a given RBF in one structure. The user can add a favorite radial function by following the directions in the section Built-in Radial Functions.
The Matlab user should go directly to the Matlab section.
Using RBFpack is easy provided the user understands the Surface Fitting Paradigm. This paradigm is:
|
Data (LSD *) ® |
Surface (RBF *) ® |
Evaluate / Operate |
The easiest way to envision an appropriate data set for RBFpack is to think of a table or spreadsheet of data arranged in a rectangle. The first M columns represent independent values, while the last N columns contain the dependent values. We illustrate this with surface values drawn from the function
(w, x, y, z) ® (w + x + y + z, 2w + x + y + z)
which has a 4 dimensional domain and a 2 dimensional range. Some data from this surface can be found in the table below:
|
Observation # |
Ind Var 1 |
Ind Var 2 |
Ind Var 3 |
Ind Var 4 |
Dep Var 1 |
Dep Var 2 |
|
1 |
.3 |
.2 |
.4 |
.5 |
1.4 |
1.7 |
|
2 |
.5 |
.5 |
.4 |
.1 |
1.5 |
2.0 |
|
3 |
.1 |
.1 |
.2 |
.4 |
.7 |
0.8 |
Since C does not naturally support two-dimensional arrays, RBFpack has a data structure (LSD) to easily store and retrieve these values. Thus, the user must first initialize an LSD data structure of the correct size with lsd_initialize and then stuff the problem data into the structure, or the user can use lsd_read to read this data from a file. After the data is in the LSD data structure, the user can view this data by calling lsd_view.
Once having stored the data in an LSD data structure, the next operation is to produce an RBF approximation to this data. There are three main methods for doing this. The function, rbf_fit, is recommended if a least squares fit is desired. This function attempts to find a global, nonlinear least-squares fit to the data by varying the centers and the center parameters. A sophisticated local minimizer (LSGRG2) is used to complement the global techniques. It is possible that rbf_fit could take too much time to arrive at an answer. In this case, we recommend rbf_lsq_fit. This will produce a least squares fit using linear techniques. It requires the user to supply all of the parameters for the RBF except the coefficients (rbf->rad_coef). The last option is to use the interpolation routine rbf_interpolate. This function produces an interpolating RBF to the data.
At this point, it is assumed that an RBF has been computed and stored in an RBF data structure. Having a mathematical representation of the surface is nice, but it is not very helpful if there are no evaluation or differentiation routines. Thus, one uses rbf_value to get function values and derivatives. Other functions are available to compute the least squares error (rbf_error), to view the RBF parameters (rbf_view), or to compute statistical information on the fit (rbf_error_stat).
RBFpack is designed to produce a surface from discrete data. Once produced, we can evaluate, differentiate, or view the error between the surface and data. The general methodology is to minimize the least squares error between the data and a candidate surface (using the RBF representations mentioned above). Thus the general sequence of events is:
|
Data (LSD *) ® |
Surface (RBF *) ® |
Evaluate / Operate |
Consider the problem of approximating the surface (x, y, sin(x+y^2)) from the data
|
X values |
Y values |
Sin(x+y^2) |
|
0 |
0 |
0 |
|
0.5 |
0 |
0.479426 |
|
1 |
0 |
0.841471 |
|
0 |
0.5 |
0.247404 |
|
0.5 |
0.5 |
0.681639 |
|
1 |
0.5 |
0.948985 |
|
0 |
1 |
0.841471 |
|
0.5 |
1 |
0.997495 |
|
1 |
1 |
0.909297 |
This data has nine observations. The domain dimension is two and the range dimension is one. We present this data to the RBFpack functions through a Least Squares Data structure called LSD. It is a C data structure and it stores all the relevant data for the general surface-fitting problem as follows:
struct LSD *lsd; int error_flag; /* Allocate space for the LSD data structure for 9 observations, domain dimension = 2 and range dimension = 1 */ lsd =lsd_initialize(9 /* nobs */, 2 /* domain dimension */, 1 /* range dimension */, &error_flag); lsd->xdata[0][0] = 0., lsd->xdata[0][1] = 0., lsd->ydata[0][0] = 0.; lsd->xdata[1][0] = .5, lsd->xdata[1][1] = 0., lsd->ydata[1][0] = .479426; lsd->xdata[2][0] = 1., lsd->xdata[2][1] = 0., lsd->ydata[2][0] = .841471; lsd->xdata[3][0] = 0., lsd->xdata[3][1] = .5, lsd->ydata[3][0] = .247404; lsd->xdata[4][0] = .5, lsd->xdata[4][1] = .5, lsd->ydata[4][0] = .681639; lsd->xdata[5][0] = 1., lsd->xdata[5][1] = .5, lsd->ydata[5][0] = .948985; lsd->xdata[6][0] = 0., lsd->xdata[6][1] = 1., lsd->ydata[6][0] = .841471; lsd->xdata[7][0] = .5, lsd->xdata[7][1] = 1., lsd->ydata[7][0] = .997495; lsd->xdata[8][0] = 1., lsd->xdata[8][1] = 1., lsd->ydata[8][0] = .909297;
The function lsd_initialize defaults the weights (in the least squares problem) to one. A more compact representation for the last 9 lines is:
for(j=0; j<3; j++)
for(k=0; k<3; k++)
lsd->xdata[k+3*j][0] = k/2. , lsd->xdata[k+3*j][1] = j/2.,
lsd->ydata[k+3*j][0] = sin( k/2. + (j/2.)*(j/2.) );
The data has been moved into the LSD data structure. Now we compute the approximating surface (whose parameters are stored in the RBF data structure). There are several functions which return RBF surface fits. We will discuss two functions in the overview. The first function is rbf_fit. This function uses a rather complex nonlinear algorithm to find a good fit to the data with as few parameters as possible. The function rbf_fit has three arguments. The first argument is the pointer to the LSD data structure, the second argument is the maximum number of centers that will be used in the fit, and the third argument is an error flag. Below is a code fragment which returns the surface fit in the RBF data structure rbf.
struct RBF * rbf; rbf = rbf_fit(lsd, 3, &error_flag);
A second function, rbf_interpolate, produces an interpolating surface to the given data. One generally uses the interpolation function when there is no noise in the data as in this simple example. The four arguments are the data (in an LSD data structure), a parameter to control the linear algebra solution, a pointer to an integer which returns the numerical rank, and an error flag.
struct RBF * rbf; int rank; rbf = rbf_interpolate(lsd, 1.e-6, &rank, &error_flag);
From the previous discussion, we assume that a surface has been computed and its parameters are stored in the RBF data structure. We now want to evaluate or operate on the computed surface (perhaps for the purpose of obtaining a graph). The following functions are useful in this regard:
|
Function name |
Function result |
|
Returns zero, first, or second derivatives |
|
|
Returns least squares error |
|
|
Prints information on the RBF parameters |
|
|
Computes and/or prints error statistics |
For example, if one wants the values of the surface on a 10 by 10 grid the following code could be used.
double x[2], Value[10][10];
for(j=0; j<10; j++)
for(k=0; k<10; k++)
x[0] = j/10. , x[1] = k/10.,
rbf_value(rbf, x, 0, 0, &Value[j][k]);
From a user perspective there are only a few functions that one needs to understand. The program listed below demonstrates an elementary use of RBFpack based on the discussion above. The fit is stored in the Radial Basis Function labeled rbf. The call to rbf_fit produces the least squares fit. The user should study this code to understand the flow of data from one call to the next, since this is a standard method of producing a least squares fit.
#include "rbf.h"
void main()
{ /* Example 1 */
struct LSD *lsd;
struct RBF *rbf;
int inform, j, k, nout=9, npts = 5;
double error, Value[10][10], x[2];
lsd = lsd_initialize(npts*npts /* nobs */, 2 /* domain dimension */,
1 /* range dimension */, &inform);
for(j=0; j<npts; j++)
for(k=0; k<npts; k++)
lsd->xdata[k+npts*j][0] = k/(npts-1.), lsd->xdata[k+npts*j][1] = j/(npts -1.),
lsd->ydata[k+npts*j][0] = sin( k/(npts -1.) + (j/(npts -1.))*(j/(npts -1.)) );
rbf = rbf_fit(lsd, 3, &inform); /* fit data */
error = rbf_error(rbf, lsd, &inform); /* compute error */
printf("The error is %5.6f \n", error);
for(j=0; j<nout; j++){ /* Print surface values */
printf(" \n");
for(k=0; k<nout; k++){
x[0] = j/(nout -1.) , x[1] = k/(nout-1.),
rbf_value(rbf, x, 0, 0, &Value[j][k]);
printf(" %4.2f ", Value[j][k] );}
}
}
The output on the screen should be "The error is .000794" followed by a 10 by 10 matrix of surface values.
We continue the discussion of Example 1. In this example, the same computation is made, but we use the functions lsd_view, rbf_view, and rbf_error_stat to view the data, the RBF fit, and the error characteristics. The function set_io_file opens the output file whose name is "example2.txt".
#include "rbf.h"
void main()
{ /* Example 2 */
struct LSD *lsd;
struct RBF *rbf;
int inform, j, k, nout=9, npts = 5;
double error, stat_out[5];
lsd = lsd_initialize(npts*npts /* nobs */, 2 /* domain dimension */,
1 /* range dimension */, &inform);
/* get data */
for(j=0; j<npts; j++)
for(k=0; k<npts; k++)
lsd->xdata[k+npts*j][0] = k/(npts-1.), lsd->xdata[k+npts*j][1] = j/(npts -1.),
lsd->ydata[k+npts*j][0] = sin( k/(npts -1.) + (j/(npts -1.))*(j/(npts -1.)) );
set_io_file("example2.txt", "open output file", &inform);
lsd_view(lsd, 10, &inform); /* view first 10 data points */
rbf = rbf_fit(lsd, 3, &inform); /* fit data */
rbf_view(rbf, 2, &inform); /* view RBF coefs and centers */
error = rbf_error(rbf, lsd, &inform); /* compute error */
printf("The error is %5.3f \n", error);
rbf_error_stat(rbf, lsd, 0, stat_out, &inform); /* view error stats */
set_io_file("example2.txt", "close output file", &inform);
}
The data written into example2.txt follows:
The following is the output fromlsd_view at level 10 The number of observations is 25 The domain dimension is 2 The range dimension is 1 0 xdata 0.000000 0.000000 0 ydata 0.000000 0 weight 1.000000 1 xdata 0.250000 0.000000 1 ydata 0.247404 1 weight 1.000000 2 xdata 0.500000 0.000000 2 ydata 0.479426 2 weight 1.000000 3 xdata 0.750000 0.000000 3 ydata 0.681639 3 weight 1.000000 4 xdata 1.000000 0.000000 4 ydata 0.841471 4 weight 1.000000 5 xdata 0.000000 0.250000 5 ydata 0.062459 5 weight 1.000000 6 xdata 0.250000 0.250000 6 ydata 0.307439 6 weight 1.000000 7 xdata 0.500000 0.250000 7 ydata 0.533303 7 weight 1.000000 8 xdata 0.750000 0.250000 8 ydata 0.726009 8 weight 1.000000 9 xdata 1.000000 0.250000 9 ydata 0.873575 9 weight 1.000000 The following is the output from rbf_view at level 2 The number of centers is 3 The domain dimension is 2 The range dimension is 1 The function switch is gaussian. coefficients for center 0 -41.360094 coefficients for center 1 40.127099 coefficients for center 2 0.098018 coefficients for constant 0.926057 Center 0 -0.098248 -0.101000 Center 1 -0.089618 -0.097606 Center 2 0.346511 0.807048 This is output from rbf_error_stat at level 0 The number of observations is 25 The R-square value is 0.991015 The SSE value is 0.000794 The SST value is 0.088401 The Max absolute error is 0.072847 The Ave absolute error is 0.020636
This example features the function rbf_lsq_fit. This function solves a linear least squares problem. The RBF data structure must be initialized. Then rbf_lsq_fit computes the RBF coefficients stored in rbf->rad_coef and rbf->const_coef.
#include "rbf.h"
void main()
{ /* Example 3 */
struct LSD *lsd;
struct RBF *rbf;
int inform, j, k, rank, npts = 5;
double error, stat_out[5], rcond = 1.e-6;
lsd = lsd_initialize(npts*npts /* nobs */, 2 /* domain dimension */,
1 /* range dimension */, &inform);
/* get data */
for(j=0; j<npts; j++)
for(k=0; k<npts; k++)
lsd->xdata[k+npts*j][0] = k/(npts-1.), lsd->xdata[k+npts*j][1] = j/(npts -1.),
lsd->ydata[k+npts*j][0] = sin( k/(npts -1.) + (j/(npts -1.))*(j/(npts -1.)) );
set_io_file("example3.txt", "open output file", &inform);
rbf = rbf_initialize(9, lsd->domain_dim, lsd->range_dim,
gaussian, 1, .77, &inform); /* initialize rbf */
rbf_set_values(rbf, lsd, .777, &inform); /* use lsd to set rbf parameters */
rbf_lsq_fit(rbf, lsd, rcond, &rank, &inform); /* fit data */
rbf_view(rbf, 2, &inform); /* view RBF coefs and centers */
error = rbf_error(rbf, lsd, &inform); /* compute error */
printf("The error is %5.3f \n", error);
rbf_error_stat(rbf, lsd, 0, stat_out, &inform); /* view error stats */
set_io_file("example3.txt", "close output file", &inform);
}
The output follows: Notice that although the RBF computed in this example has nine centers while Examples 1 and 2 have three centers. The centers have not been optimized and consequently the least squares error is much larger than in the earlier examples (4.156e-3 versus 7.94e-4). Of course, the computational time for Example 3 is somewhat less.
The following is the output from rbf_view at level 2 The number of centers is 9 The domain dimension is 2 The range dimension is 1 The function switch is gaussian. coefficients for center 0 0.089003 coefficients for center 1 0.481678 coefficients for center 2 1.377656 coefficients for center 3 0.127444 coefficients for center 4 6.003829 coefficients for center 5 -6.559107 coefficients for center 6 0.191392 coefficients for center 7 -1.891668 coefficients for center 8 -0.055148 coefficients for constant 0.754600 Center 0 0.102205 0.757105 Center 1 0.669371 0.121016 Center 2 0.918494 0.125301 Center 3 0.942194 0.459643 Center 4 0.223915 0.332254 Center 5 0.194164 0.318733 Center 6 0.953705 0.917702 Center 7 0.817704 0.151333 Center 8 0.457494 0.098142 This is output fromrbf_error_stat at level 0 The number of observations is 25 The R-square value is 0.952989 The SSE value is 0.004156 The SST value is 0.088401 The Max absolute error is 0.163329 The Ave absolute error is 0.046738
The next example produces an interpolant to the data. The interpolant uses the radial_poly radial basis function.
#include "rbf.h"
void main()
{ /* Example 4 */
struct LSD *lsd;
struct RBF *rbf;
int inform, j, k, rank, npts = 5;
double error, stat_out[5], rcond = 1.e-6;
lsd = lsd_initialize(npts*npts /* nobs */, 2 /* domain dimension */,
1 /* range dimension */, &inform);
/* get data */
for(j=0; j<npts; j++)
for(k=0; k<npts; k++)
lsd->xdata[k+npts*j][0] = k/(npts-1.), lsd->xdata[k+npts*j][1] = j/(npts -1.),
lsd->ydata[k+npts*j][0] = sin( k/(npts -1.) + (j/(npts -1.))*(j/(npts -1.)) );
set_io_file("example4.txt", "open output file", &inform);
rbf_fit_opt_radfcn(radial_poly); /* change default radial function */
rbf = rbf_interpolate(lsd, rcond, &rank, &inform); /* interpolate data */
rbf_view(rbf, 0, &inform); /* view RBF coefs and centers */
error = rbf_error(rbf, lsd, &inform); /* compute error */
printf("The error is %5.3f \n", error);
rbf_error_stat(rbf, lsd, 0, stat_out, &inform); /* view error stats */
set_io_file("example4.txt", "close output file", &inform);
}
The output from this example is listed below. Note that the R-square is one and the error is zero to six places. This is to be expected for interpolation. Notice that the max error is not zero since we are allowing the possibility of errors in the third decimal to creep into the calculation. This is a matter of conditioning of the linear interpolation problem. If this is not acceptable then one can adjust the rcond parameter.
The following is the output fromrbf_view at level 0 The number of centers is 25 The domain dimension is 2 The range dimension is 1 The function switch is radial_poly. This is output from rbf_error_stat at level 0 The number of observations is 25 The R-square value is 1.000000 The SSE value is 0.000000 The SST value is 0.088401 The Max absolute error is 0.000474 The Ave absolute error is 0.000075
This last example is the most realistic use of RBFpack. It is also the most complicated. The problem is to find an approximating surface to two-dimensional histogram data. This problem is prototypical of many problems which arise in engineering and finance. That is, a surface is to be erected over a domain, but in addition to the raw data used to determine the fit there are certain qualitative features which must be reflected by the surface. For the example below where we are trying to produce a smooth surface from histogram data these qualitative features are A) the integral of the surface should be 1 and B) the surface should be nonnegative. We use lsgrg2 to produce an approximating surface satisfying these "qualitative" constraints. Lsgrg2 is an optimizer which allows equality and inequality constraints to be applied while optimizing some quantity. Please consult our web pages at web.wt.net/~wti for more information on lsgrg2.
The data is generated in the function get_density_data. Now a histogram surface should have several characteristics. It should be non-negative and the integral of the surface over its domain should be one. These characteristics are (partially) enforced by the constraints generated in the functions rbf_optimize_const (where the upper and lower bounds are set for these constraints) and in lsq_error_const where the corresponding surface values are computed. We would like to enforce non-negativity everywhere, however since that is not possible we require that the function be non-negative on the 7 by 7 grid where the histogram data is taken. This is just a convenience and a larger grid could be used for this purpose if necessary.
The constrained fit looks like the surface below, except we have graphed the max(0, s(x,y)) where s is the surface computed. You might wonder why the surface is not nonnegative already since we are enforcing s >= 0 on a 7 by 7 grid. The problem is that we are using a finer grid to obtain more detail and on this finer grid a few of the values are negative. We could have enforced nonnegativity on this finer grid but that would have increased the computing time:
The code for this example follows: Note that global variables centers_1, …, const_coef_1 are used to specify which of the RBF parameters are to be used in the optimization. For this example, we allow only the radial coefficient parameters to vary (ie rad_coef_1 =1). This example requires 3 minutes to run on a 166 MH pentium. This is a rather complicated example and it is suggested that the user spend some time studying it if a similar problem is to be attempted.
#include "rbf.h"
struct RBF * glob_rbf_1;
struct LSD * glob_lsd_1;
int NPTS = 7;
int NUM_SIMS = 10000;
int FIRST = 1;
int NVARS, NFUNS;
double x[10000][2], y[20][20];
/* These are the switches to determine which
variables are in play when the optimizer is used */
int centers_1 = 0; /* InActive */
int rad_coef_1 = 1; /* Active */
int center_parameter_1 = 0; /* InActive */
int domain_scale_1 = 0; /* InActive */
int const_coef_1 = 0; /* InActive */
TYPE_export(struct RBF *) rbf_optimize_const(struct RBF * rbf_in, struct LSD * lsd,
int itmax, int * inform);
TYPE_export(void) lsq_error_const(double g[], double x[]);
TYPE_export(struct LSD *) get_density_data();
void print_graph_data(struct RBF * rbf );
double lsq_nlin(double *x);
TYPE_export(void) const_P (double x[], double g[],
long int igrad[], long int jgrad[],
double grad[], long int nzg);
TYPE_export(void) get_mat (int nfuns, int nvars, double x[], double g[],
long int igrad[], long int jgrad[],
double grad[]) ;
double f(double x,double y);
double g(double x, double y);
double drand(double);
void main()
{/* Example 5 */
struct LSD *lsd;
struct RBF *rbf, *rbf1;
int inform, k, count, npts1, i, j;
double sum=0., xmin, xmax, ymin, ymax, x[2];
lsd = get_density_data();
for(k=0; k<lsd->nobs; k++)
sum += (lsd->ydata[k][0]);
sum = sum/lsd->nobs;
printf("sum = %f \n", sum);
set_io_file("example5.txt", "open output file", &inform);
rbf = rbf_initialize(16+36, lsd->domain_dim, lsd->range_dim,
radial_poly, 1, .677, &inform); /* initialize rbf */
rbf_set_values(rbf, lsd, .5, &inform);
rbf->domain_scale[0]=1. , rbf->domain_scale[1]=1.;
rbf->const_coef[0] = .0;
count = 0;
npts1 = 4;
xmin = .35 , xmax = .65;
ymin=.35 , ymax = .65;
for(i=0; i<npts1; i++){
for(j=0; j<npts1; j++){
x[0] = xmin+i*(xmax-xmin)/(npts1-1.);
x[1] = ymin+j*(ymax-ymin)/(npts1-1.);
rbf->center_parameter[count] = 1.;
rbf->centers[count][0] = x[0];
rbf->centers[count++][1] = x[1];
}
}
npts1 = 6;
xmin = .3 , xmax = .7;
ymin=.3 , ymax = .7;
for(i=0; i<npts1; i++){
for(j=0; j<npts1; j++){
x[0] = xmin+i*(xmax-xmin)/(npts1-1.);
x[1] = ymin+j*(ymax-ymin)/(npts1-1.);
rbf->center_parameter[count] = 5.;
rbf->centers[count][0] = x[0];
rbf->centers[count++][1] = x[1];
}
}
rbf->num_centers = 16+36;
for(i=0; i<rbf->num_centers; i++)
rbf->rad_coef[i][0] = .17;
rbf_view(rbf, 5, &inform);
rbf1 = rbf_optimize_const(rbf, lsd, 500, &inform);
printf("the error is %f\n", sqrt(rbf_error(rbf1, lsd, &inform)));
rbf1->const_coef[0] = .001;
print_graph_data(rbf1);
rbf_view(rbf1, 5, &inform);
lsd_view(lsd, 225, &inform);
set_io_file("example5.txt", "close output file", &inform);
}
TYPE_export(struct RBF *) rbf_optimize_const(struct RBF * rbf_in, struct LSD * lsd,
int itmax, int * inform)
{
/* Purpose: This function returns the result
of applying optimization to the standard regression
problem. The starting point for the nonlinear
regression is contained in rbf_in and the data for
the regression is in lsd.
Return Value: A pointer to the solution of the NL regression
problem.
Arguments:
rbf_in The initial guess for the Nonlinear
regression proble. (input)
lsd The regression data (input)
itmax An upperbound on the number of line searches (input)
in the optimization algorithm (input)
inform Error code (output)
Fatal Errors:
inform Error Condition
100 rbf_in has some int members <= 0
200 rbf_copy failed
300 allocation error
J return value from lsgrg2 whose value is J
*/
struct RBF * rbf_out;
int i, j, info, count = 0;
double *max_box, *min_box, big = 1.e30, temp;
int nvars, nfuns, nobj;
double *xx, *xlb, *xub, *glb, *gub;
char *title, *report;
double *gg;
/* Check initial input */
if(rbf_in->num_centers <= 0 || rbf_in->domain_dim <= 0 ||
rbf_in->range_dim <= 0){
*inform = 100;
return NULL; /* Either rbf_in int members are mis-specified
or rbf_in and rbf_out have different int data */
}
/* Set pointers to global data structures
and copy the input info from rbf_in into glob_rbf_1 */
glob_rbf_1 = rbf_copy(rbf_in, &info);
if(info != 0){
*inform = 200;
return NULL; /* Failure from rbf_copy */
}
glob_lsd_1 = lsd;
/* Set optimization parameters */
nfuns = 2+glob_lsd_1->nobs+1;
nobj = -1;
nvars = centers_1 *(rbf_in->num_centers * (rbf_in->domain_dim)) +
rad_coef_1 *(rbf_in->num_centers * (rbf_in->range_dim)) +
center_parameter_1 *(rbf_in->num_centers) +
domain_scale_1 *(rbf_in->domain_dim) +
const_coef_1 *(rbf_in->range_dim)+1;
NVARS = nvars, NFUNS = nfuns;
title = "Output from rbf_optimize";
report = "rbopt.txt";
/* Allocate space for the extreme corners of the
smallest box that contains the points in lsd->xdata
and store these in min_box and max_box */
max_box = (double *) CALLOC(rbf_in->domain_dim, sizeof(double));
min_box = (double *) CALLOC(rbf_in->domain_dim, sizeof(double));
xlb = (double *) CALLOC(nvars+1, sizeof(double));
xub = (double *) CALLOC(nvars+1, sizeof(double));
xx = (double *) CALLOC(nvars+1, sizeof(double));
glb = (double *) CALLOC(nfuns+1, sizeof(double));
gub = (double *) CALLOC(nfuns+1, sizeof(double));
gg = (double *) CALLOC(nfuns+1, sizeof(double));
/* Error checking for allocation failure */
if((void *) max_box == NULL || (void *) min_box == NULL ||
(void *) xlb == NULL || (void *) xub == NULL ||
(void *) glb == NULL || (void *) gub == NULL ||
(void *) gg == NULL ) {
*inform = 300;
return NULL;
}
/* Set function constraints */
glb[2] = 1.0, gub[2] = 1.0;
glb[0] = -big, gub[0] = big;
for(i=0; i<glob_lsd_1->nobs; i++)
glb[3+i] = 0, gub[3+i] = big;
get_box(lsd, min_box, max_box);
/* Set the upper and lower bounds for the variables */
for(i=0; i<=nvars; i++)
xlb[i] = -big, xub[i] = big;
count = 0;
if(centers_1 != 0)
for(i=0; i<rbf_in->num_centers; i++)
for(j=0; j< rbf_in->domain_dim; j++){
temp = .1*(.01+ max_box[j] - min_box[j]);
xlb[++count] = min_box[j] - temp;
xub[ count] = max_box[j] + temp;
}
if(rad_coef_1 != 0) /* No upper bound, lower bound = 0 */
for(i=0; i<rbf_in->num_centers; i++)
for(j=0; j<rbf_in->range_dim; j++){
xlb[++count] = -big;
xub[ count] = big;
}
if(center_parameter_1 != 0)
for(i=0; i<rbf_in->num_centers; i++){
xlb[++count] = .3/sqrt((double) rbf_in->num_centers);
xub[ count] = sqrt((double) rbf_in->num_centers);
}
if(domain_scale_1 != 0)
for(i=0; i<rbf_in->domain_dim; i++){
temp = 1.1*(.01+ max_box[j] - min_box[j]);
xlb[++count] = .1/temp;
xub[ count] = 10./temp;
}
if(const_coef_1 != 0) /* No upper bound, lower bound = 0 */
for(i=0; i<rbf_in->range_dim; i++){
xlb[++count] = -big;
xub[ count] = big;
}
/* Set optimization parameters and call the optimization routine */
rbf_to_x(glob_rbf_1, xx, centers_1, rad_coef_1, center_parameter_1,
domain_scale_1, const_coef_1);
lgoptD( "limser", itmax );
lgoptD("report", 1.); /* control printing */
/* Must use "fullj" option if centers or center parmeters are
variables since the structure of the jacobion might change */
if(centers_1 == 1 || center_parameter_1 == 1) lgoptD ("fullj", 1);
lgoptG( lsq_error_const ) ;
/* if centers and center parameters are not used then we can
take advantage of the fact that all the functions are linear
with the exception of the objective */
if(centers_1 ==0 && center_parameter_1 ==0) lgoptP(const_P);
lsgrg2 ( nvars, xlb, xub, nfuns, nobj, glb, gub, title, report,
xx, &info ) ;
*inform = info;
if(info <0) return NULL; /* Failure from rbf_optimize */
/* Put the final results into rbf_out */
lsq_error_const(gg, xx);
rbf_out = rbf_copy(glob_rbf_1, &info);
free_RBF(glob_rbf_1), free(xx), free(xlb), free(xub);
free(min_box), free(max_box), free(glb), free(gub), free(gg);
return rbf_out;
}
TYPE_export(void) lsq_error_const(double g[], double x[])
{
/*
Purpose: This function computes the square root of the
least squares errorgiven the values x and
the global variables centers_1, rad_coef_1, center_parameter_1,
domain_scale_1, const_coef_1. This version assumes that
rad_coef_1 = 0 and it computes the linear variables by calling
rbf_lsq_fit. This routine is designed to be called
by grg2.
Return Value: void
Arguments:
g The objective value for the NL regression
is stored in g[1] (output)
x The input independent variables for the
grg2 function (input)
Global Variables Used:
glob_lsd_1 Pointer to global Least Squares Data set
glob_rbf_1 Pointer to global rbf
centers_1 = 0 centers_1 do not enter the optimizaton
rad_coef_1 = 1 coeficients do explicitly enter
center_parameter_1 = 0 center_parameter_1 do not enter the optimization
domain_scale_1 = 0 domain_scale_1 is not optimized
const_coef_1 = 0 const_coef_1 does not explicitly enter the
optimization problem.
*/
int i, j, count,info;
double v[1], sum;
x_to_rbf(x, glob_rbf_1, centers_1, rad_coef_1,
center_parameter_1, domain_scale_1, const_coef_1);
g[1] = sqrt(rbf_error(glob_rbf_1, glob_lsd_1, &info));
sum = 0.;
for(i=0; i< glob_lsd_1->nobs; i++){
rbf_value(glob_rbf_1, glob_lsd_1->xdata[i], 0, 0, v);
sum = sum + v[0];
}
g[2] = sum/(glob_lsd_1->nobs);
/* for(i=0; i< glob_lsd_1->nobs; i++)
{rbf_value(glob_rbf_1, glob_lsd_1->xdata[i], 0, 0, v);
g[3+i] = v[0];}
*/
count = 0;
for(i=0; i<NPTS; i++)
for(j=0; j<NPTS; j++){
rbf_value(glob_rbf_1, glob_lsd_1->xdata[j+NPTS*i], 0, 0, v);
g[3+count++] = v[0];
}
g[0] =0.; /* define g[0] so that the optimization can "start" at 0 */
return;
}
double f(double x, double y)
{double temp;
temp = .1+(2*x+y)/4.;
return temp;
}
double g(double x, double y)
{double temp;
temp = .3+ y/2;
return temp;
}
TYPE_export(struct LSD *) get_density_data()
{struct LSD * lsd;
int i , j, k, inform;
double xsum, ysum;
/* fill simulation data */
drand(.777);
for(j=0; j<NUM_SIMS; j++){
xsum = (drand(-1.)+drand(-1.)+drand(-1.))/3.;
ysum = (drand(-1.)+drand(-1.)+drand(-1.))/3.;
x[j][0] = f(xsum, ysum);
x[j][1] = g(xsum, ysum);
}
for(j = 0; j<NPTS; j++)
for(k = 0; k<NPTS; k++)
y[j][k] = 0.;
/* first pass, find out how many points are in each bin */
for(i=0; i < NUM_SIMS; i++)
for(j = 0; j<NPTS; j++)
for(k = 0; k<NPTS; k++)
if(j/(NPTS-0.) <= x[i][0] &&
k/(NPTS-0.) <= x[i][1] &&
(j+1)/(NPTS-0.) > x[i][0] &&
(k+1)/(NPTS-0.) > x[i][1]) y[j][k] += 1.;
for(j = 0; j<NPTS; j++)
for(k = 0; k<NPTS; k++)
y[j][k] = y[j][k]/NUM_SIMS;
lsd = lsd_initialize((NPTS)*(NPTS) /* nobs */, 2 /* domain dimension */,
1 /* range dimension */, &inform);
for(j = 0; j<NPTS; j++)
for(k = 0; k<NPTS; k++)
lsd->xdata[j + k*(NPTS)][0] = (j+.5)/(NPTS),
lsd->xdata[j + k*(NPTS)][1] = (k+.5)/(NPTS),
lsd->ydata[j + k*(NPTS)][0] = y[j][k]*(NPTS*NPTS);
return lsd;
}
void print_graph_data(struct RBF * rbf )
{
int i, j, NUM_GRAPH = 15;
double v[1], x[2];
FILE *fp;
fp = fopen("rbf_graph_3.txt" , "w");
for(i=0; i<NUM_GRAPH; i++){
for(j=0; j<NUM_GRAPH; j++){
x[0] = i/(NUM_GRAPH-1.), x[1] = j/(NUM_GRAPH-1.);
rbf_value(rbf, x, 0, 0, v);
fprintf(fp, " %f \t", v[0]);
}
fprintf(fp, "\n");
}
fprintf(fp, "\n");
for(i=0; i<NPTS; i++){
for(j=0; j<NPTS; j++){
rbf_value(rbf, glob_lsd_1->xdata[j+NPTS*i], 0, 0, v);
fprintf(fp, " %f \t", v[0]);
}
fprintf(fp, "\n");
}
fprintf(fp, "\n");
for(i=0; i<NPTS; i++){
for(j=0; j<NPTS; j++){
v[0] = glob_lsd_1->ydata[j+NPTS*i][0];
fprintf(fp, " %f \t", v[0]);
}
fprintf(fp, "\n");
}
fclose(fp);
}
TYPE_export(void) get_mat (int nfuns, int nvars, double x[], double g[],
long int igrad[], long int jgrad[],
double grad[])
{ int i, j, nnz;
double h = 1.e-7, *gx, *gxph, *jacob_col;
gx = (double *) CALLOC(nfuns+1, sizeof(double));
gxph = (double *) CALLOC(nfuns+1, sizeof(double));
jacob_col = (double *) CALLOC(nfuns+1, sizeof(double));
lsq_error_const(gx, x);
nnz = 0;
jgrad[0] = 0;
for(i=0; i<nvars; i++){
x[i] = x[i] +h;
lsq_error_const(gxph, x);
for(j=0; j<nfuns; j++)
jacob_col[j] = (gxph[j] - gx[j])/h;
for(j=0; j<nfuns; j++)
if(jacob_col[j] != 0.){
grad[nnz] = jacob_col[j];
igrad[nnz++] = j;}
jgrad[i+1] = nnz;
x[i] = x[i] -h;
}
free(gx), free(gxph), free(jacob_col);
}
TYPE_export(void) const_P (double x[], double g[],
long int igrad[], long int jgrad[],
double grad[], long int nzg)
{ int i, nfuns = 124, nvars = 63;
double h = 1.e-7, temp;
nfuns = NFUNS, nvars = NVARS;
if(nzg >0) get_mat(nfuns, nvars, x, g, igrad, jgrad, grad);
else{
for(i=1; i<nvars; i++){
x[i] += h;
temp = (lsq_nlin(x)-g[1])/h;
x[i] -= h;
grad[jgrad[i]] = temp;
}
}
}
double lsq_nlin(double *x)
{ /* code to compute the value of the objective function only.
This function is called by the Jacobian generating function
const_P */
double temp;
int info;
x_to_rbf(x, glob_rbf_1, centers_1, rad_coef_1,
center_parameter_1, domain_scale_1, const_coef_1);
temp = sqrt(rbf_error(glob_rbf_1, glob_lsd_1, &info));
return temp;
}
The routines listed in this section are sufficient to allow the user to compute least squares fits to large and complicated data sets. At this level, we have attempted to hide all the distracting details of radial basis functions. We have identified the routines which will allow the user to compute sophisticated surface fits and interpolants without having to dive into the minutiae of radial basis functions. To this end we provide an initialization routine so that the user can initialize the area where he will store the data he is trying to fit. We also provide the all-purpose rbf_fit routine which will produce a least squares fit to the data from (a somewhat intelligent) choice of RBF parameters. Likewise, we provide the interpolation routine rbf_interpolate. Having computed a least-squares fit or an interpolant we must be able to evaluate the fit and asses the error. The routines rbf_error and rbf_value satisfy this need.
TYPE_export(struct LSD *) lsd_initialize(int nobs, int domain_dim,
int range_dim, int *inform)
.......................................................................
Purpose: This function allocates and initializes all the pointers
in the LSD data structure.
Return Value: A pointer to the newly initialized LSD structure
Arguments:
nobs Number of observations (input)
domain_dim Domain dimension (input)
range_dim Range dimension (input)
inform Error code (output)
Fatal Errors (which return a NULL pointer):
inform Error Condition
1 nobs <= 0
2 domain_dim <= 0
3 range_dim <= 0
4 allocation failure for some member of LSD
.......................................................................
TYPE_export(struct RBF *) rbf_fit(struct LSD * lsd, int max_centers,
int * inform)
.......................................................................
Purpose: This function returns the result
of applying local and global optimization
techniques to the standard regression
problem. The data for the regression is in lsd.
This function is a good starting point to explore
the rbf fit properties. The default parameters
are chosen to balance the computational burden against
the attempt to find a good fit.
See rbf_fit_opt_int and rbf_fit_radfcn for functions
which can change the defaults set up in this function.
Return Value: A pointer to the solution of the RBF regression
problem. The default radial function is the gaussian.
Arguments:
lsd The regression data (input)
max_centers An upperbound on the number of centers
to be used in the search for a good RBF fit (input)
inform Error code (output)
Fatal Errors:
inform Error Condition
1 problem is mis-specified (ie integer
data in lsd is less than or equal 0.
2 max_centers = 0.
3 allocation error
4 allocation error from rbf_initialize
5 allocation error from rbf_set_values
.......................................................................
TYPE_export(struct RBF *) rbf_interpolate(struct LSD * lsd_in,
double rcond, int *rank, int *inform)
.......................................................................
Purpose: This function computes an interpolant
to the data in lsd_in. The user can set the radial
function by invoking rbf_fit_radfcn.
Return Value: The rbf interpolant.
Arguments:
lsd_in Pointer to LSD structure (input)
rcond This parameter is used to determine
the effective rank of the
(nobs x 1+nobs) matrix for
the interpolation problem.
Suggested value for rcond is 1.e-6;
This variable must be positive. (input)
rank The numerical effective rank of the
interpolation problem which depends
on rcond (output)
inform Error code (output)
Fatal Errors:
inform Error Condition
1 rcond is not > 0
2 Error from linear lsq solver dgelsx
(check integer members of the
structures)
3 Allocation failure (check integer
members of the structures)
.......................................................................
TYPE_export(double) rbf_error(struct RBF * rbf_in, struct LSD * lsd_in,
int *inform)
.......................................................................
Purpose: This function computes the average weighted
squared error
(sum_{i=1}^nobs w_i ||F(x_i) - y_i||^2)/(nobs*range_dim)
between the values in the rbf F specified in rbf_in
and the Least Squares Data set specified in lsd_in
Return Value: Weighted average squared error
Arguments:
rbf_in Pointer to RBF structure (input)
lsd_in Pointer to LSD structure (input)
inform Error code (output)
Fatal Errors:
inform Error Condition
1 The two structures are incompatible
due to different specification of
domain_dim
range_dim
members
2 Allocation error
.......................................................................
TYPE_export(void) rbf_value(struct RBF * rbf_in, double *x,
int ider1, int ider2, double *value)
.......................................................................
Purpose: This function computes the value
D_{(ider1, ider2)}F(x)
(of the function specified in rbf_in)
and returns it in the vector value.
Return Value: void
Arguments:
rbf_in Pointer to RBF structure (input)
x Vector of length rbf_in->domain_dim (input)
ider1 Partial derivative specification (input)
ider2 Partial derivative specification (input)
value Value of D_{(ider1, ider2)}F(x) (output)
Fatal Errors:
No Fatal errors are defined since this function is the heart
of many other computations and error handling would adversely
affect performance.
.......................................................................
The second level of RBFpack requires the user to become more involved with the RBF parameters. In particular to user rbf_lsq_fit (which computes a least-squares fit given the nonlinear RBF parameters) requires the user to specify some of the RBF parameters (such as the center parameters, the scale parameters, the centers, and of course the number of the centers). The user is assisted in this effort by rbf_initialize (which sets up the pointers for the RBF) and rbf_set_values which generates trial values for the RBF parameters.
TYPE_export(void) rbf_lsq_fit(struct RBF * rbf_in, struct LSD * lsd_in,
double rcond, int *rank, int *inform)
.......................................................................
Purpose: This function computes a least squares fit
to the data in lsd_in. It assumes that all
the values in rbf_in have been determined with
the exception of the members {rad_coef const_coef}.
This function then computes the appropriate coefficient
vaues and places them in rbf_in.
Return Value: void
Arguments:
rbf_in Pointer to RBF structure (input/output)
lsd_in Pointer to LSD structure (input)
rcond This parameter is used to determine
the effective rank of the
(nobs x 1+num_centers) matrix for
the least squares problem.
Suggested value for rcond is 1.e-6;
This variable must be positive. (input)
rank The numerical effective rank of the
least squares problem which depends
on rcond (output)
inform Error code (output)
Fatal Errors:
inform Error Condition
1 The two structures are incompatible
due to different specification of
domain_dim
range_dim
members
2 rcond is not > 0
3 Error from linear lsq solver dgelsx
(check integer members of the
structures)
4 Allocation failure (check integer
members of the structures)
.......................................................................
TYPE_export(struct RBF *) rbf_initialize(int num_centers, int domain_dim,
int range_dim,
double (STDCALL *radial_fcn)(double, double, int),
int level, double seed, int * inform)
.......................................................................
Purpose: This function allocates and initializes all the pointers
in the RBF data structure. Depending on the specified
level some of the numerical data in the structure may
also be initialized.
Return Value: A pointer to the newly initialized RBF structure
Arguments:
num_centers Number of centers (input)
domain_dim Domain dimension (input)
range_dim Range dimension (input)
radial_fcn Name of a function (input)
level Level of numerical values placed in (input)
the rbf structure (input)
{
level <= 0 All double data set to zero in rbf structure
level = 1 Center parameters and domain scale are
set to 1.1 and 1.2 respectively
level >= 2 (level 1) and centers are filled out
randomly in [0,1]^(domain_dim)
}
seed A seed for the random number generator
inform Error code (output)
Fatal Errors (which return a NULL pointer):
inform Error Condition
1 num_centers <= 0
2 domain_dim <= 0
3 range_dim <= 0
4 allocation failure for some member of RBF
.......................................................................
TYPE_export(void) rbf_set_values(struct RBF * rbf_in, struct LSD * lsd_in,
double seed, int *inform)
.......................................................................
Purpose: This function selects centers and domain_scale values
based on the xdata in lsd_in. The centers are randomly
chosen in the smallest hyper-cube containing the xdata
values. The domain_scale values are selected to equalize
the apparent width of each dimension.
Return Value: void
Arguments:
rbf_in Pointer to RBF structure (input)
lsd_in Pointer to LSD structure (input)
seed Seed for the random number generator (input)
inform Error code (output)
Fatal Errors:
inform Error Condition
1 The two structures are incompatible
due to different specification of
domain_dim
range_dim
members
2 Allocation failure (check value of domain_dim)
.......................................................................
There are four built-in radial functions. They are the gaussian, hardy_multiquadric, inverse_hardy_multiquadric, and the radial_poly. The user can add his own radial function. RBFpack supports user supplied radial functions by allowing the user to code and specify a new radial function. All the user must do is provide code based on the template of the gaussian listed below. The value r2 is the (scaled) distance squared, p is a parameter used to change the shape of the (univariate) radial function, and k is the derivative (with respect to r2). Note, there is a routine which can be used to check the correct programming of the partials (rbf_fd_der) which uses finite differences to compute the partials.
TYPE_export(double) gaussian(double r2, double p, int k)
{
/* Purpose: This function computes the Gaussian
radial function defined by
f(r2) = exp(-p*p*r2)
*/
double value, f_value=0.;
value = exp(-p*p*r2);
if (k == 0) f_value = value;
if (k == 1) f_value = -p*p*value;
if (k == 2) f_value = p*p*p*p*value;
return f_value;
}
TYPE_export(double) gaussian(double r2, double p, int k)
.......................................................................
Purpose: This function computes the Gaussian
radial function defined by
f(r2) = exp(-p*p*r2)
as well as its first and second derivatives in r2.
This function is used as a pointer to a function in
rbf->radial_fcn.
Return Value: D^k f(r2).
Arguments:
r2 r2 (usually the squared scaled distance from
scnorm) (input)
p A scalling parameter (input)
This parameter is supplied most often
by center_parameter members.
k The derivative value (input)
Only values 0, 1, 2 are expected,
all others return 0.
Fatal Errors:
No Fatal errors are defined since this function is the heart
of many other computations and error handling would adversely
affect performance.
.......................................................................
TYPE_export(double) radial_poly(double r2, double p, int k)
.......................................................................
Purpose: This function computes the (compactly supported)
radial function defined by
f(r2) =
(1-p*sqrt(r2))^4*(4*p*sqrt(r2) + 1) if p*sqrt(r2) < 1
0 if p*sqrt(r2) >= 1
as well as its first and second derivatives in r2.
This function is used as a pointer to a function in
rbf->radial_fcn.
Return Value: D^k f(r2).
Arguments:
r2 r2 (usually the squared scaled distance from
scnorm) (input)
p A scalling parameter (input)
This parameter is supplied most often
by center_parameter members.
k The derivative value (input)
Only values 0, 1, 2 are expected all
others return 0.
Fatal Errors:
No Fatal errors are defined since this function is the heart
of many other computations and error handling would adversely
affect performance.
.......................................................................
TYPE_export(double) hardy_multiquadric(double r2, double p, int k)
.......................................................................
Purpose: This function computes the hardy multiquadric
radial function defined by
f(r2) = sqrt(p**2 + r2)
as well as its first and second derivatives in r2.
This function is used as a pointer to a function in
rbf->radial_fcn.
Return Value: D^k f(r2).
Arguments:
r2 r2 (usually the squared scaled distance from
scnorm) (input)
p A scalling parameter (input)
This parameter is supplied most often
by center_parameter members.
k The derivative value (input)
Only values 0, 1, 2 are expected,
all others return 0.
Fatal Errors:
No Fatal errors are defined since this function is the heart
of many other computations and error handling would adversely
affect performance.
.......................................................................
TYPE_export(double) inverse_hardy_multiquadric(double r2, double p, int k)
.......................................................................
Purpose: This function computes the inverse hardy multiquadric
radial function defined by
f(r2) = 1/sqrt(p**2 + r2)
as well as its first and second derivatives in r2.
This function is used as a pointer to a function in
rbf->radial_fcn.
Return Value: D^k f(r2).
Arguments:
r2 r2 (usually the squared scaled distance from
scnorm) (input)
p A scaling parameter (input)
This parameter is supplied most often
by center_parameter members.
k The derivative value (input)
Only values 0, 1, 2 are expected,
all others return 0.
Fatal Errors:
No Fatal errors are defined since this function is the heart
of many other computations and error handling would adversely
affect performance.
.......................................................................
These routines allow the user to change the various parameters in the rbf_fit routine.
TYPE_export(void) rbf_fit_opt_int(char * par_name, int value,
int * inform)
.......................................................................
Purpose: This function sets integer parameters for rbf_fit.
Setting Optimize = 0 forces no optimization. In this
case the function rbf_fit relies solely on global search.
Setting Optimize = 1 forces nonlinear optimization in
conjunction with global search. Otherwise the
function attempts to choose the most efficient method.
Note: The default values are reset after each call to
rbf_fit. Thus they must be explicitly overridden before
each invocation of rbf_fit.
Return Value: Void.
Arguments:
par_name The name of the parameter to be set (input)
value The value of the parameter (input)
par_name value default
"Optimize" 0, 1 Not set
"Num_tries" 1, 2, ... 150 This parameter controls
the global search. Larger
values will produce smaller
sums of squares and take
more time.
"Max_it" 1, 2, ... 30 This parameter controls the
nonlinear optimizer. Larger
values will produce smaller
sums of squares and take
more time.
Fatal Errors:
inform Error
1 Optimize value not 0 or 1
2 Num_tries value not positive
3 Max_it value not positive
4 Unrecognized option.
.......................................................................
TYPE_export(void) rbf_fit_opt_radfcn(double (STDCALL *radial_fcn)
(double, double, int))
.......................................................................
Purpose: To change the default radial function used in rbf_fit
or rbf_interpolate.
Note: The default values are reset after each call to
rbf_fit or rbf_interpolate. Thus they must be explicitly
overridden before each invocation of rbf_fit.
Return Value: Void.
Arguments:
radial_fcn The pointer to the new radial function. (input)
The default radial function is the gaussian.
This package internally supports three other
radial functions:
radial_poly
hardy_multiquadric
inverse_hardy_multiquadric
The user can supply his own radial function
provided
it is properly prototyped and implemented as any of
the above 3 functions.
Fatal Errors: None
.......................................................................
These routines are designed to read or write data to files.
TYPE_export(int) set_io_file(char * file_name, char * s, int *inform)
.......................................................................
Purpose: This function opens or closes input or output
files.
Return Value: 0 for correct function
1 for unrecognized character string (s) input
-1 for file open failure
Arguments:
file_name Name of file where data is to be sent (input)
s Option character string (input)
{The only strings recognized by the function
are:
"open output file"
"open output file append"
"open input file"
"close output file"
"close input file"
}
inform Error code (output)
Fatal Errors:
inform Error Condition
1 Output file could not be opened
2 Input file could not be opened
3 Unrecognizable second argument
.......................................................................
TYPE_export(void) rbf_view(struct RBF * rbf_in, int level, int * inform)
.......................................................................
Purpose: This function prints out various levels of information
contained in the RBF structure.
Return Value: void
Arguments:
rbf_in Pointer to RBF structure (input)
level Level of printed information (nested) (input)
{
level = 0 All integer information is printed from
rbf_in
level > 0 All coefficients are printed out
level > 1 All centers are printed out
level > 2 All center parameters are printed out
level > 3 All "domain scale" info is printed out
}
inform Error code (output)
Fatal Errors:
inform Error Condition
1 Output file not set (see set_io_file)
.......................................................................
TYPE_export(void) rbf_error_stat(struct RBF * rbf_in, struct LSD * lsd_in,
int level, double * stat_out, int *inform)
.......................................................................
Purpose: This function computes statistics on the Least
Squares Fit and if requested it prints out these
statistics.
Return Value: void
Arguments:
rbf_in Pointer to RBF structure (input)
lsd_in Pointer to LSD structure (input)
level level of printed output required (input)
{
level < 0 no printed output
level <= 0 stat_out is printed out
}
stat_out A vector of length 5 which contains statistics
of the fit
{
stat_out[0] = rsq (R-squre)
stat_out[1] = sse (normalized sum of
squares error)
stat_out[2] = sst (sum of squares total)
stat_out[3] = mae (maximum absolute error)
stat_out[4] = aae (average absolute error)
Note rsq = 1 - (sse/sst) and under
usual conditions is between 0 and 1.
}
inform Error code (output)
Fatal Errors:
inform Error Condition
1 The two structures are incompatible
due to different specification of
domain_dim
range_dim
members
2 Allocation failure (check value of range_dim)
3 Output file not set
4 Allocation error in rbf_error
.......................................................................
TYPE_export(struct LSD *) lsd_read(char *file_name, int option,
int *inform)
.......................................................................
Purpose: This function reads data from the data file
specified in file_name. The first line of the file
must contain only the three integers
nobs domain_dim range_dim
Return Value: (LSD *) A pointer to a Least Squares Data set
Arguments:
file_name The file_name of the data set (input)
option Option describing the format of the
file to be read(input)
{
option = 0
Data formated with data numbers and weights ie
0 .5 .3 .2 <-- xdata
0 .1 <-- ydata
0 2. <-- weight
option = 1
data
number xdata ydata weight
0 .5 .3 .2 .1 2
1 .5 .2 .2 .2 7
option = 2
data
number xdata ydata
0 .5 .3 .2 .1
The weights are set to 1 in option 2.
}
inform Error code (output)
Fatal Errors:
inform Error Condition
1 set_io_file failed to open the designated
input file
2 lsd_initialize failed to allocate the LSD
data structure.
3 set_io_file failed to close the input file
.......................................................................
TYPE_export(void) lsd_view(struct LSD * lsd_in, int level,int * inform)
.......................................................................
Purpose: This function prints out various levels of information
contained in the LSD structure.
Return Value: void
Arguments:
lsd_in Pointer to RBF structure (input)
level Level of printed information (nested) (input)
{
level = 0 All integer information is printed from
lsd_in
level > 0 The first level data points are printed out
with their weights.
}
inform Error code (output)
Fatal Errors:
inform Error Condition
1 Output file not set (see set_io_file)
.......................................................................
These routines allow one to optimize the RBF parameters in the context of a least-squares problem.
TYPE_export(struct RBF *) rbf_optimize(struct RBF * rbf_in,
struct LSD * lsd,
int itmax, int * inform)
.......................................................................
Purpose: This function returns the result
of applying optimization to the standard regression
problem. The starting point for the nonlinear
regression is contained in rbf_in and the data for
the regression is in lsd.
Return Value: A pointer to the solution of the NL regression
problem.
Arguments:
rbf_in The initial guess for the Nonlinear
regression proble. (input)
lsd The regression data (input)
itmax An upperbound on the number of line
searches (input)
in the optimization algorithm (input)
inform Error code (output)
Fatal Errors:
inform Error Condition
1 rbf_in has some int members <= 0
2 rbf_copy failed
3 allocation error
J<0 error from grg2 whose value is J
.......................................................................
TYPE_export(struct RBF *) step_up_optimize(struct RBF *rbf_in,
struct LSD * lsd,
double * lbound, double * ubound, int num_steps,
int num_fcn_evals_per_step, int itmax, int *inform)
.......................................................................
Purpose: This function returns the result
of applying a step_up algorithm followed by
optimization to the standard regression
problem. The starting point for the
regression is contained in rbf_in and the data for
the regression is in lsd.
Return Value: A pointer to the solution of the step up
Nonlinear regression problem.
Arguments:
rbf_in The initial guess for the Nonlinear
regression problem. (input)
lsd The regression data (input)
lbound A pointer to the lower bounds for the
centers and the center_parameter. (input)
ubound A pointer to the lower bounds for the
centers and the center_parameter. (input)
num_fcn_evals_per_step The number of random centers tried
per step in the step_up algorithm (input)
itmax An upper bound on the number of line searches
in the optimization algorithm. (input)
inform Error code (output)
Fatal Errors:
inform Error Condition
1 rbf_in has some int members <= 0
2 rbf_copy failed
.......................................................................
TYPE_export(struct RBF *) step_up(struct RBF *rbf_in, struct LSD * lsd,
double * lbound, double * ubound, int num_steps,
int num_fcn_evals_per_step, int *inform)
.......................................................................
Purpose: This function returns the result
of applying a step_up algorithm followed by
optimization to the standard regression
problem. The starting point for the
regression is contained in rbf_in and the data for
the regression is in lsd.
Return Value: A pointer to the solution of the step up
Nonlinear regression problem.
Arguments:
rbf_in The initial guess for the Nonlinear
regression problem. (input)
lsd The regression data (input)
lbound A pointer to the lower bounds for the
centers and the center_parameter. (input)
ubound A pointer to the lower bounds for the
centers and the center_parameter. (input)
num_fcn_evals_per_step The number of random centers tried
per step in the step_up algorithm (input)
itmax An upper bound on the number of line searches
in the optimization algorithm. (input)
inform Error code (output)
Fatal Errors:
inform Error Condition
1 rbf_in has some int members <= 0
2 rbf_copy failed
.......................................................................
TYPE_export(void) rbf_to_x(struct RBF * rbf, double *x, int centers,
int rad_coef, int center_parameter,
int domain_scale, int const_coef)
.......................................................................
Purpose: This function moves the parameters in rbf
into the appropriate positions in x.
Return Value: void
Arguments:
rbf The radial basis function which will
receive the parameters in x (output)
x The parameters to be copied into
rbf. (input)
centers Variable indicating whether the centers
are to be varied in the NL regression
problem (input)
rad_coef Variable indicating whether the radial
coefficients are to be varied in the NL
regression problem (input)
center_parameter Variable indicating whether the center
parameters are to be varied in the NL
regression problem (input)
domain_scale Variable indicating whether the domain
scale are to be varied in the NL
regression problem (input)
const_coef Variable indicating whether the constant
coefficients are to be varied in the NL
regression problem (input)
.......................................................................
TYPE_export(void) x_to_rbf(double *x, struct RBF * rbf, int centers,
int rad_coef, int center_parameter,
int domain_scale, int const_coef)
.......................................................................
Purpose: This function moves the parameters pointed to
by x into the appropriate positions into
the rbf.
Return Value: void
Arguments:
x The parameters to be copied into
rbf. (input)
rbf The radial basis function which will
receive the parameters in x (output)
centers Variable indicating whether the centers
are to be varied in the NL regression
problem (input)
rad_coef Variable indicating whether the radial
coefficients are to be varied in the NL
regression problem (input)
center_parameter Variable indicating whether the center
parameters are to be varied in the NL
regression problem (input)
domain_scale Variable indicating whether the domain
scale are to be varied in the NL
regression problem (input)
const_coef Variable indicating whether the constant
coefficients are to be varied in the NL
regression problem (input)
.......................................................................
TYPE_export(void) lsq_error(double g[], double x[])
.......................................................................
Purpose: This function computes the square root of the
least squares error given the values x and
the global variables centers, rad_coef, center_parameter,
domain_scale, const_coef. This version assumes that
rad_coef = 0 and it computes the linear variables by calling
rbf_lsq_fit. This function is designed to be called
by grg2.
Return Value: void
Arguments:
g The objective value for the NL regression
is stored in g[1] (output)
x The input independent variables for the
grg2 function (input)
Global Variables Used:
glob_lsd Pointer to global Least Squares Data set
glob_rbf Pointer to global rbf
centers = 1 Centers enter the optimization
rad_coef = 0 coefficients do not explicitly enter
center_parameter = 1 center_parameter enter the optimization
domain_scale = 0 domain_scale is not optimized
const_coef = 0 const_coef does not explicitly enter the
optimization problem.
.......................................................................
TYPE_export(struct LSD *) test_data(struct RBF * rbf, int nobs,
double seed, int choice,
int * inform)
.......................................................................
Purpose: This function returns an LSD * which is generated by known
data so that we can test the least squares functions.
Return Value: A pointer to the LSD test data.
Arguments:
rbf Radial basis function
nobs Number of observations (input)
domain_dim Domain dimension (input)
range_dim Range dimension (input)
seed Seed for the random number generator (input)
choice Choice of data set generated. Unless otherwise
stated the weights are 1 and the xdata values are
randomly generated in the unit cube.
Choice Data
0 1
1 ydata[i] = rbf(lsd->xdata[i])
2 ydata[i] = rbf(lsd->xdata[i])
with weights alternately 0 and 1
3 ydata[i] = rbf(lsd->xdata[i])
with weights alternately 1+i
4 same as 3 but also writes the data
to the file lsq_data.txt
inform Error code (output)
Fatal Errors (which return a NULL pointer):
inform Error Condition
1 nobs <= 0
2 rbf->domain_dim <= 0
3 rbf->range_dim <= 0
4 allocation failure for some member of LSD
5 error from set_io_file
.......................................................................
TYPE_export(double) rbf_compare(struct RBF * rbf1, struct RBF * rbf2)
.......................................................................
Purpose: This function compares two rbfs. If they are
different then the function returns a positive
number indicating the difference.
Return Value: -1. if one of the two rbfs are NULL
0. if the rbfs point to identical data
Otherwise a positive value is returned
indicating that the two rbfs represent
different functions
Arguments:
rbf1 The first rbf (input)
rbf2 The second rbf (input)
.......................................................................
TYPE_export(struct RBF *) rbf_copy(struct RBF * rbf_in, int * inform)
.......................................................................
Purpose: This function produces a copy of rbf_in.
Return Value: Pointer to a copy of rbf_in.
Arguments:
rbf_in The RBF to be copied (input)
inform Error code (output)
Fatal Errors:
inform Error Condition
1 rbf_in is inconsistently specified
2 rbf_initialize failed
.......................................................................
TYPE_export(void) copy_rbf(struct RBF * rbf_in,
struct RBF * rbf_out)
.......................................................................
Purpose: This function produces a copy of rbf_in in rbf_out.
it is assumed that rbf_out has been appropriately
initialized. This function is useful if one wants
to allocate an rbf once and reuse the allocation by
storing "smaller" rbfs in the space.
Return Value: Void.
Arguments:
rbf_in The rbf to be copied (input)
rbf_out The copied version of rbf_in (out_put)
rbf_out should be initialized so that the
following equalities and inequalities hold.
rbf_in-> num_centers <= rbf_out->num_centers
rbf_in-> domain_dim = rbf_out->domain_dim
rbf_in-> range_dim = rbf_out-> range_dim
Fatal Errors:
No fatal errors.
.......................................................................
TYPE_export(void) rbf_fd_der(struct RBF *rbf_in, double *x, int ider1,
int ider2, double *value)
.......................................................................
Purpose: This function approximates the value
D_{(ider1, ider2)}F(x)
(of the function specified in rbf_in
using finite differences) and returns
it in the vector value.
Return Value: void
Arguments:
rbf_in Pointer to RBF structure (input)
x Vector of length rbf_in->domain_dim (input)
ider1 Partial derivative specification (input)
ider2 Partial derivative specification (input)
value Value of D_{(ider1, ider2)}F(x) (output)
Fatal Errors:
No Fatal errors are defined since this function is the heart
of many other computations and error handling would adversely
affect performance.
.......................................................................
TYPE_export(void) free_LSD(struct LSD * lsd_in)
.......................................................................
Purpose: This function frees all the space allocated in
the LSD structure.
Return Value: void
Arguments:
rbf_in The pointer to the LSD structure (input)
No Fatal Errors
.......................................................................
TYPE_export(void) free_RBF(struct RBF * rbf_in)
.......................................................................
Purpose: This function frees all the space allocated in
the RBF structure.
Return Value: void
Arguments:
rbf_in The pointer to the RBF structure (input)
No Fatal Errors
.......................................................................
TYPE_export(void) get_box(struct LSD * lsd_in, double * min_box, double *max_box)
.......................................................................
Purpose: This function computes the lower left
(in min_box) and upper right (in max_box)
coordinates of the smallest box containing
the xdata values in lsd_in.
Return Value: void
Arguments:
min_box A vector of length lsd_in->domain_data (output)
The coordinates in min_box satisfy
min_box[j] = min{lsd_in->xdata[i][j]:
1<= i <= lsd_in->nobs}
max_box A vector of length lsd_in->domain_data (output)
The coordinates in min_box satisfy
max_box[j] = max{lsd_in->xdata[i][j]: 1<= i <= lsd_in->nobs}
.......................................................................
struct RBF /* Radial Basis Funcion */ {
int num_centers;
int domain_dim;
int range_dim;
double (STDCALL *radial_fcn)(double, double, int);
double ** rad_coef;
double * const_coef;
double ** centers;
double * center_parameter;
double * domain_scale;
};
struct LSD /*Least Squares Data */ {
int nobs;
int domain_dim;
int range_dim;
double ** xdata;
double ** ydata;
double * weight;
};
Using RBFpack is easy provided the user understands the Surface Fitting Paradigm. This paradigm is:
|
Data (LSD *) ® |
Surface (RBF *) ® |
Evaluate / Operate |
The easiest way to envision an appropriate data set for RBFpack is to think of a table or spreadsheet of data arranged in a rectangle. The first M columns represent independent values, while the last N columns contain the dependent values. We illustrate this with surface values drawn from the function
(w, x, y, z) ® (w + x + y + z, 2w + x + y + z)
which has a 4 dimensional domain and a 2 dimensional range. Some data from this surface can be found in the table below:
|
Observation # |
Ind Var 1 |
Ind Var 2 |
Ind Var 3 |
Ind Var 4 |
Dep Var 1 |
Dep Var 2 |
|
1 |
.3 |
.2 |
.4 |
.5 |
1.4 |
1.7 |
|
2 |
.5 |
.5 |
.4 |
.1 |
1.5 |
2.0 |
|
3 |
.1 |
.1 |
.2 |
.4 |
.7 |
0.8 |
RBFpack has a data structure (LSD) to easily store and retrieve these values. Thus, the user must first initialize an LSD data structure of the correct size using lsd_struct and then stuff the problem data into the structure. In particular, lsd for the above data set would be initialized by the call lsd_struct(3, 4, 2) so that lsd.nobs = 3, lsd.domain_dim = 4, lsd.range_dim = 2, lsd.xdata is a (lsd.nobs by lsd.domain_dim) matrix, lsd.ydata is a (lsd.nobs by lsd.range_dim) matrix, and lsd.weight is a (lsd.nobs by 1) matrix.
Once having stored the data in an LSD data structure, the next operation is to produce an RBF approximation to this data. There are three main methods for doing this. The function, rbf_fit, is recommended if a least squares fit is desired. This function attempts to find a global, nonlinear least-squares fit to the data by varying the centers and the center parameters. A sophisticated local minimizer (LSGRG2) is used to complement the global techniques. It is possible that rbf_fit could take too much time to arrive at an answer. In this case, we recommend rbf_lsq_fit. This will produce a least squares fit using linear techniques. The last option is to use the interpolation routine rbf_interpolate. This function produces an interpolating RBF to the data and hence is not appropriate for noisy data.
At this point, it is assumed that an RBF has been computed and stored in an RBF data structure (see rbf_struct). Having a mathematical representation of the surface is nice, but it is not very helpful if there are no evaluation or differentiation routines. Thus, one uses rbf_value to get function values and derivatives. A function is available to compute the least squares error, rbf_error, and finally there is a prototype m-file for graphing certain RBFs, rbf_graph_m.
RBFpack is designed to produce a surface from discrete data. Once produced, we can evaluate, differentiate, or view the error between the surface and data. The general methodology is to minimize the least squares error between the data and a candidate surface (using the RBF representations mentioned above). Thus the general sequence of events is:
|
Data (lsd_struct) ® |
Surface (rbf_struct) ® |
Evaluate / Operate |
Consider the problem of approximating the surface (x, y, sin(x+y^2)) from the data
|
X values |
Y values |
Sin(x+y^2) |
|
0 |
0 |
0 |
|
0.5 |
0 |
0.479426 |
|
1 |
0 |
0.841471 |
|
0 |
0.5 |
0.247404 |
|
0.5 |
0.5 |
0.681639 |
|
1 |
0.5 |
0.948985 |
|
0 |
1 |
0.841471 |
|
0.5 |
1 |
0.997495 |
|
1 |
1 |
0.909297 |
This data has nine observations. The domain dimension is two and the range dimension is one. We present this data to the RBFpack functions through a Least Squares Data structure called LSD.
The following m-file follows the paradigm described above. First, the LSD data structure is allocated and filled out. Second, the RBF approximation is computed and third, this approximation is graphed, the error is computed, and the RBF is evaluated at (0,0). Notice that at most 3 centers are required for the fit. The routine rbf_fit can be very time-consuming. For that reason we have given the user the option of guiding the routine by furnishing a second argument which places an upper bound on the number of centers allowed. This m-file is provided with the RBFpack product under the name example_1.m.
function example_1()
num_points = 15;
% allocate the lsd data structure
lsd = lsd_struct(num_points*num_points, 2, 1);
x = [0:num_points-1]/(num_points-1);
for i = 1:num_points
for j = 1:num_points
lsd.xdata(j + (i-1)*num_points, :) = [x(i);x(j)];
end
end
lsd.ydata = sin(lsd.xdata(:,1) + lsd.xdata(:,2).^2);
% Get the approximation
rbf = rbf_fit(lsd, 5);
% Get the error
'the error is: '
rbf_error(rbf,lsd)
% Graph the result
rbf_graph_m(rbf, lsd, [0,0], [1,1], 'Surface Plot of RBF computed using rbf\_fit');
% Get the value of the RBF at the origin (0,0)
'The value at (0,0) is:'
rbf_value(rbf, [0,0])
This example uses the same data set as above and computes the least squares fit based on 10 centers. The location of the centers is randomly generated and not optimized. Comparing the error of the fit in this example with that of Example 1, we see that the error is smaller in Example 1 even though the number of centers is more than tripled. This illustrates the power of optimizing the parameters in the RBF fit. This m-file is provided with the RBFpack product under the name example_2.m.
functionexample_2() num_points = 15; % allocate the lsd data structure lsd = lsd_struct(num_points*num_points, 2, 1); x = [0:num_points-1]/(num_points-1); for i = 1:num_points for j = 1:num_points lsd.xdata(j + (i-1)*num_points, :) = [x(i);x(j)]; end end lsd.ydata = sin(lsd.xdata(:,1) + lsd.xdata(:,2).^2); % Get the approximation rbf = rbf_lsq_fit(lsd, 10) %Examine the rbf values rbf.rad_coef rbf.centers rbf.center_parameter rbf.domain_scale rbf % Get the error 'the error is: ' rbf_error(rbf,lsd) % Graph the result rbf_graph_m(rbf, lsd, [0,0], [1,1], 'Surface Plot of RBF computed using rbf\_lsq\_fit');
Interpolation is appropriate when the data is not noisy. The interpolation routine returns an RBF with as many centers as data points. Thus for very large data sets, least squares may be appropriate even if the data is not noisy. We also demonstrate the use of rbf_fit_opt_radfcn to change the radial function from 'gaussian' to 'radial_poly'. This m-file is provided with the RBFpack product under the name example_3.m.
functionexample_3() num_points = 5; % allocate the lsd data structure lsd = lsd_struct(num_points*num_points, 2, 1); x = [0:num_points-1]/(num_points-1); for i = 1:num_points for j = 1:num_points lsd.xdata(j + (i-1)*num_points, :) = [x(i);x(j)]; end end lsd.ydata = sin(lsd.xdata(:,1) + lsd.xdata(:,2).^2); % Get the approximation rbf_fit_opt_radfcn('radial_poly'); rbf = rbf_interpolate(lsd); % Get the error 'the error is: ' rbf_error(rbf,lsd) % Graph the result rbf_graph_m(rbf, lsd, [0,0], [1,1], 'Surface Plot of RBF computed by rbf\_interpolate'); % Get the value of the RBF at the origin (0,0) 'The value at (0,0) is:' rbf_value(rbf, [0,0])
This example resets the default radial function from 'gaussian' to 'radial_poly'. We then try to fit data with noise first by using rbf_lsq_fit and then by rbf_fit. When executed, graphs of both functions are displayed. This m-file is provided with the RBFpack product under the name example_4.m.
functionexample_4() num_points = 30; % allocate the lsd data structure lsd = lsd_struct(num_points, 1, 1); x = [0:num_points-1]/(num_points-1); for j = 1:num_points lsd.xdata(j , 1) = x(j); end lsd.ydata = sin(6*lsd.xdata(:,1).^2)+ (rand(num_points, 1)-.5); % Get the approximation rbf_fit_opt_radfcn('radial_poly'); rbf = rbf_lsq_fit(lsd,9); % Get the error 'the error is: ' rbf_error(rbf,lsd) % Graph the result rbf_graph_m(rbf, lsd, 0, 1, 'Plot of data and RBF computed by rbf\_lsq\_fit'); figure rbf_fit_opt_radfcn('radial_poly'); rbf = rbf_fit(lsd, 3) % Graph the result rbf_graph_m(rbf, lsd, 0, 1, 'Plot of data and RBF computed by rbf\_fit'); % Get the error 'the error is: ' rbf_error(rbf,lsd)
This example illustrates the