radialScatteredFit¶

Computes an approximation to scattered data in \(ℜ^n\) for \(n \geq 1\) using radial-basis functions.

Synopsis¶

radialScatteredFit (abscissae, fdata, numCenters)

Required Arguments¶

float abscissae[[]] (Input): Array of size dimension × numPoints containing the abscissae of the data points. The argument abscissae[i][j] is the abscissa value of the (i+1)-th data point in the (j+1)-th dimension.
float fdata[] (Input): Array with numPoints components containing the ordinates for the problem.
int numCenters (Input): The number of centers to be used when computing the radial-basis fit. The argument numCenters should be less than or equal to numPoints.

Return Value¶

The structure that represents the radial-basis fit. If a fit cannot be computed, then None is returned.

Optional Arguments¶

centers (Input)

User-supplied centers. See the Description section of this function for details.

centersRatio, float (Input)

The desired ratio of centers placed on an evenly spaced grid to the total number of centers. The condition that the same number of centers placed on a grid for each dimension must be equal. Thus, the actual number of centers placed on a grid is usually less than centersRatio × numCenters, but will never be more than centersRatio × numCenters. The remaining centers are randomly chosen from the set of abscissae given in abscissae.

Default: centersRatio = 0.5

randomSeed, int (Input)

The value of the random seed used when determining the random subset of abscissae to use as centers. By changing the value of seed on different calls to radialScatteredFit, with the same data set, a different set of random centers will be chosen. Setting randomSeed to zero forces the random number seed to be based on the system clock, so a possibly different set of centers will be chosen each time the program is executed.

Default: randomSeed = 234579

supplyBasis, float radialFunction (distance) (Input)

User-supplied function to compute the values of the radial functions.

Default: Hardy multiquadric

supplyDelta, float (Input)

The delta used in the default basis function.

\[\phi(r) = \sqrt{r^2 + \delta^2}\]

Default: supplyDelta = 1

weights, float[]

This option requires the user to provide the weights.

Default: all weights equal one

noSvd

This option forces the use of a QR decomposition instead of a singular value decomposition. This may result in space savings for large problems.

Description¶

The function radialScatteredFit computed a least-squares fit to scattered data in \(ℜ^d\) where d = dimension. More precisely, let n = ndata, x = abscissae, f = fdata, and d = dimension. Then we have

\[x^0, \ldots, x^{n-1} \subset \Re^d f_0, \ldots, f_{n-1} \subset \Re^1\]

This function computes a function F which approximates the above data in the sense that it minimizes the sum-of-squares error

\[\sum_{i=0}^{n-1} w_i \left(F\left(x^i\right)-f_i\right)^2\]

where w = weights. Of course, we must restrict the functional form of F. This is done as follows:

\[F(x) := \sum_{j=0}^{k-1} \alpha_j \phi \left(\|x-c\|^2 + \delta^2\right)^{1/2}\]

The function ɸ is called the radial function. It maps \(ℜ^1\) into \(ℜ^1\), only defined for the nonnegative reals. For the purpose of this routine, the user-supplied function

\[\phi(r) = \left(r^2 + \delta^2\right)^{1/2}\]

Note that the value of delta is defaulted to 1. It can be set by the user by using the keyword supplyDelta. The parameter δ is used to scale the problem. Generally choose δ to be near the minimum spacing of the centers.

The default basis function is called the Hardy multiquadric, and it is defined as

\[\phi(r) = \left(r^2 + \delta^2\right)^{1/2}\]

A key feature of this routine is the user’s control over the selection of the basis function.

To obtain the default selection of centers, we first compute the number of centers that will be on a grid and how many are on a random subset of the abscissae. Next, we compute those centers on a grid. Finally, a random subset of abscissa are obtained determining where the centers are placed. Let us examine the selection of centers in more detail.

First, we restrict the computed grid to have the same number of grid values in each of the dimension directions. Then, the number of centers placed on a grid, numGridded, is computed as follows:

\[\alpha = (\texttt{centersRatio}) (\texttt{numCenters})\]

\[\beta = ⌊\alpha^{1/dimension}⌋\]

\[\texttt{numGridded} = \beta^{dimension}\]

Note that there are β grid values in each of the dimension directions. Then we have

\[\texttt{numRandom} = (\texttt{numCenters}) - (\texttt{numGridded})\]

Now we know how many centers will be placed on a grid and how many will be placed on a random subset of the abscissae. The gridded centers are computed such that they are equally spaced in each of the dimension directions. The last problem is to compute a random subset, without replacement, of the abscissa. The selection is based on a random seed. The default seed is 234579. The user can change this using the optional argument randomSeed. Once the subset is computed, we use the abscissae as centers.

Since the selection of good centers for a specific problem is an unsolved problem at this time, we have given the ultimate flexibility to the user. That is, you can select your own centers using the keyword centers. As a rule of thumb, the centers should be interspersed with the abscissae.

The return value for this function is the structure, which contains all the information necessary to evaluate the fit. This is then passed to the function radialEvaluate to produce values of the fitted function.

Examples¶

Example 1¶

This example, generates data from a function and contaminates it with noise on a grid of 10 equally spaced points.The fit is evaluated on a finer grid and compared with the actual function values.

from __future__ import print_function
from numpy import *
from pyimsl.math.constant import constant
from pyimsl.math.radialEvaluate import radialEvaluate
from pyimsl.math.radialScatteredFit import radialScatteredFit
from pyimsl.math.randomSeedSet import randomSeedSet
from pyimsl.math.randomUniform import randomUniform


def F(x):
    return sin(2.0 * pi * x)


ndata = 10
num_centers = 5
noise_size = 0.25
fdata = empty(ndata)
xdata = empty(ndata)
xdata2 = empty((ndata * 2, 1), dtype=double)
pi = constant("pi")
randomSeedSet(234579)
noise = randomUniform(ndata)

# Set up the sampled data points with noise.
for i in range(0, ndata):
    xdata[i] = float(i) / (ndata - 1)
    fdata[i] = F(xdata[i]) + noise_size * (1.0 - 2.0 * noise[i])

# Compute the radial fit.
radial_fit = radialScatteredFit(xdata, fdata, num_centers)

# Compare result to the original function at twice as many values as
# there were original data points.
for i in range(0, ndata * 2):
    xdata2[i, 0] = float(i) / (2 * (ndata - 1))

# Evaluate the fit at these new points.
fdata2 = radialEvaluate(xdata2, radial_fit)
print("    I     TRUE       APPROX     ERROR")
for i in range(0, 2 * ndata):
    print("%5d %10.5f %10.5f %10.5f"
          % (i + 1, F(xdata2[i, 0]), fdata2[i], F(xdata2[i, 0]) - fdata2[i]))

Output¶

    I     TRUE       APPROX     ERROR
  0.00000   -0.12821    0.12821
  0.34202    0.40184   -0.05982
  0.64279    0.79627   -0.15348
  0.86603    1.04735   -0.18132
  0.98481    1.15482   -0.17001
  0.98481    1.12662   -0.14181
  0.86603    0.97878   -0.11275
  0.64279    0.73465   -0.09187
  0.34202    0.42350   -0.08148
  0.00000    0.07850   -0.07850
 -0.34202   -0.26554   -0.07648
 -0.64279   -0.57454   -0.06825
 -0.86603   -0.81744   -0.04858
 -0.98481   -0.96821   -0.01660
 -0.98481   -1.00729    0.02249
 -0.86603   -0.92258    0.05655
 -0.64279   -0.70962    0.06684
 -0.34202   -0.37134    0.02932
 -0.00000    0.08297   -0.08297
  0.34202    0.63895   -0.29693

Example 2¶

This example generates data from a function and contaminates it with noise.We fit this data successively on grids of size 10, 20, …, 100. Now interpolate and print the 2-norm of the difference between the interpolated result and actual function values. Note that double precision is used for higher accuracy.

from __future__ import print_function
from numpy import *
from pyimsl.math.constant import constant
from pyimsl.math.radialEvaluate import radialEvaluate
from pyimsl.math.radialScatteredFit import radialScatteredFit
from pyimsl.math.randomSeedSet import randomSeedSet
from pyimsl.math.randomUniform import randomUniform
from pyimsl.math.vectorNorm import vectorNorm


def G(x, y):
    return exp(y / 2.0) * sin(x) - cos(y / 2.0)


def radial_function(r):
    return log(1.0 + r)


nrandom = 200
noise_size = 1.0
pi = constant("pi")

# Get the random numbers used for the noise.
randomSeedSet(234579)
noise = randomUniform(nrandom + 1)
for i in range(0, nrandom):
    noise[i] = 1.0 - 2.0 * noise[i]
print("     NDATA         || Error ||_2")

for ndata in range(10, 101, 10):
    num_centers = ndata
    fdata = empty((ndata), dtype='double')
    xydata = empty((ndata, 2), dtype='double')
    # Set up the sampled data points with noise.
    for i in range(0, 2 * ndata, 2):
        ii = int(i / 2)
        xydata[ii, 0] = 3.0 * noise[i]
        xydata[ii, 1] = 3.0 * noise[i + 1]
        fdata[ii] = G(xydata[ii, 0], xydata[ii, 1]) + noise_size * noise[i]
    # Compute the radial fit.
    ratio = 0.5
    radial_struct = radialScatteredFit(xydata, fdata,
                                       num_centers, centersRatio=ratio,
                                       supplyBasis=radial_function)
    fit = radialEvaluate(xydata, radial_struct)
    for i in range(0, ndata):
        fit[i] = fit[i] - fdata[i]
    print("%8d %17.8f" % (ndata, vectorNorm(fit)))

Output¶

     NDATA         || Error ||_2
      0.00000000
      0.00000000
      0.00000000
      0.00000000
      0.00000000
      0.00000000
      0.00000000
      0.00000000
      0.00000000
      0.00000000

Fatal Errors¶

IMSL_STOP_USER_FCN

Request from user supplied function to stop algorithm.

User flag = “#”.