# Interpolation

Interpolation is the process of constructing a function that takes on specified values at specified points. In engineering and science, one often has a number of data points, obtained by sampling or experimentation, which represent the values of a function for a limited number of values of the independent variable. It is often required to interpolate (i.e. estimate) the value of that function for an intermediate value of the independent variable.

Smile package `smile.interpolation`

provides a variety of
algorithms on 1d and 2d data. These algorithms implements the interface
`Interpolation`

(or `Interpolation2D`

for 2d data),
of which the method `interpolate`

takes a value and return
an interploated value.

## Piecewise Linear Interpolation

`LinearInterpolation`

is quick and easy, but it is not very precise. Another disadvantage
is that the interpolant is not differentiable at the control points.

```
val x = Array(0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0)
val y = Array(0.0, 0.8415, 0.9093, 0.1411, -0.7568, -0.9589, -0.2794)
val linear = new LinearInterpolation(x, y)
val data = (0 to 60).map { i =>
val x = i * 0.1
val y = linear.interpolate(x)
Array(x, y)
}.toArray
val canvas = line(data)
show(canvas)
```

```
double[] x = {0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0}
double[] y = {0.0, 0.8415, 0.9093, 0.1411, -0.7568, -0.9589, -0.2794}
var linear = new LinearInterpolation(x, y)
double[][] data = new double[61][2]
for (int i = 0; i < data.length; i++) {
data[i][0] = i * 0.1;
data[i][1] = linear.interpolate(data[i][0]);
}
var canvas = LinePlot.plot(data, Line.Style.SOLID)
canvas.window()
```

`BilinearInterpolation`

is an extension of linear interpolation for
interpolating functions of two variables on a
regular grid. The key idea is to perform linear interpolation first in one
direction, and then again in the other direction.

## Cubic Spline Interpolation

Spline interpolation uses low-degree polynomials in each of the intervals, and chooses the polynomial pieces such that they fit smoothly together. The resulting function is called a spline. The natural cubic spline is piecewise cubic and twice continuously differentiable. Furthermore, its second derivative is zero at the end points.

Like polynomial interpolation, spline interpolation incurs a smaller error than linear interpolation and the interpolant is smoother. However, the interpolant is easier to evaluate than the high-degree polynomials used in polynomial interpolation. It also does not suffer from Runge's phenomenon.

```
val cubic = new CubicSplineInterpolation1D(x, y)
val data = (0 to 60).map { i =>
val x = i * 0.1
val y = cubic.interpolate(x)
Array(x, y)
}.toArray
canvas.line("Cubic", data, RED)
show(canvas)
```

```
var cubic = new CubicSplineInterpolation1D(x, y)
double[][] data = new double[61][2]
for (int i = 0; i < data.length; i++) {
data[i][0] = i * 0.1;
data[i][1] = cubic.interpolate(data[i][0]);
}
canvas.line("Cubic", data, Color.RED)
canvas.window()
```

For 2d grid, we provide two implementations: `CubicSplineInterpolation2D`

and
`BicubicInterpolation`

. `CubicSplineInterpolation2D`

is
similar to one-dimensional splines as it guarantees the
continuity of the first and second function derivatives. In contrast,
`BicubicInterpolation`

guarantees continuity of only
gradient and cross-derivative. Second derivatives
could be discontinuous.

In image processing, bicubic interpolation is often chosen over bilinear interpolation or nearest neighbor in image resampling, when speed is not an issue. Images resampled with bicubic interpolation are smoother and have fewer interpolation artifacts.

```
val x1 = Array(1950.0, 1960, 1970, 1980, 1990)
val x2 = Array(10.0, 20, 30)
val y = Array(
Array(150.697, 199.592, 187.625),
Array(179.323, 195.072, 250.287),
Array(203.212, 179.092, 322.767),
Array(226.505, 153.706, 426.730),
Array(249.633, 120.281, 598.243)
)
val canvas = heatmap(y, Palette.jet(256))
canvas.setTitle("Original")
show(canvas)
val cubic = new CubicSplineInterpolation2D(x1, x2, y)
val data = Array.ofDim[Double](101, 101)
for (i <- 0 to 100; j <- 0 to 100)
data(i)(j) = cubic.interpolate(1950 + i*0.4, 10 + j*0.2)
val cubicPlot = heatmap(data, Palette.jet(256))
cubicPlot.setTitle("Cubic")
show(cubicPlot)
val bicubic = new BicubicInterpolation(x1, x2, y)
for (i <- 0 to 100; j <- 0 to 100)
data(i)(j) = bicubic.interpolate(1950 + i*0.4, 10 + j*0.2)
val bicubicPlot = heatmap(data, Palette.jet(256))
bicubicPlot.setTitle("Bicubic")
show(bicubicPlot)
```

```
double[] x1 = {1950.0, 1960, 1970, 1980, 1990}
double[] x2 = {10.0, 20, 30}
double[][] y = {
{150.697, 199.592, 187.625},
{179.323, 195.072, 250.287},
{203.212, 179.092, 322.767},
{226.505, 153.706, 426.730},
{249.633, 120.281, 598.243}
}
var canvas = Heatmap.plot(y, Palette.jet(256))
canvas.setTitle("Original")
canvas.window()
var cubic = new CubicSplineInterpolation2D(x1, x2, y)
double[][] data = new double[101][101]
for (int i = 0; i <= 100; i++) {
for (int j = 0; j <= 100; j++) {
data[i][j] = cubic.interpolate(1950 + i*0.4, 10 + j*0.2);
}
}
var cubicPlot = Heatmap.plot(data, Palette.jet(256))
cubicPlot.setTitle("Cubic")
cubicPlot.window()
val bicubic = new BicubicInterpolation(x1, x2, y)
for (int i = 0; i <= 100; i++) {
for (int j = 0; j <= 100; j++) {
data[i][j] = cubic.interpolate(1950 + i*0.4, 10 + j*0.2);
}
}
var bicubicPlot = Heatmap.plot(data, Palette.jet(256))
bicubicPlot.setTitle("Bicubic")
bicubicPlot.window()
```

## Kriging Interpolation

Kriging interpolation can be used for the data points irregularly distributed in space. Kriging belongs to the family of linear least squares estimation algorithms, also known as Gauss-Markov estimation or Gaussian process regression. We implement ordinary kriging for interpolation with power variogram.

```
val kriging = new KrigingInterpolation1D(x, y)
val data = Array.ofDim[Double](61, 2)
for (i <- 0 to 60) {
data(i)(0) = i * 0.1
data(i)(1) = kriging.interpolate(data(i)(0))
}
canvas.line("Kriging", data, CYAN)
show(canvas)
```

```
var kriging = new KrigingInterpolation1D(x, y)
double[][] data = new double[61][2]
for (int i = 0; i < data.length; i++) {
data[i][0] = i * 0.1;
data[i][1] = kriging.interpolate(data[i][0]);
}
canvas.line("Kriging", data, Color.CYAN)
canvas.window()
```

## RBF Interpolation

Radial basis function interpolation is a popular method for the data points are
irregularly distributed in space. In its basic form, radial basis function
interpolation is in the form
`y(x) = Σ w`

,
where the approximating function y(x) is represented as a sum of N radial
basis functions φ, each associated with a different center c_{i} φ(||x-c_{i}||)_{i}, and weighted
by an appropriate coefficient w_{i}. For distance, one usually chooses
euclidean distance. The weights w_{i} can
be estimated using the matrix methods of linear least squares, because
the approximating function is linear in the weights.

The points c_{i} often called the centers or collocation points
of the RBF interpolant. Note also that the centers c_{i} can be
located at arbitrary points in the domain, and do not require a grid.
For certain RBF exponential convergence has been shown. Radial basis
functions were successfully applied to problems as diverse as computer
graphics, neural networks, for the solution of differential equations
via collocation methods and many other problems.

Other popular choices for φ comprise the Gaussian function and the so
called thin plate splines. Thin plate splines result from the solution of
a variational problem. The advantage of the thin plate splines is that
their conditioning is invariant under scalings. Gaussians, multi-quadrics
and inverse multi-quadrics are infinitely smooth and and involve a scale
or shape parameter, r_{0} > 0.
Decreasing r_{0} tends to
flatten the basis function. For a given function, the quality of
approximation may strongly depend on this parameter. In particular,
increasing r_{0} has the effect of better conditioning
(the separation distance of the scaled points increases).

A variant on RBF interpolation is normalized radial basis function (NRBF) interpolation, in which we require the sum of the basis functions to be unity. NRBF arises more naturally from a Bayesian statistical perspective. However, there is no evidence that either the NRBF method is consistently superior to the RBF method, or vice versa.

```
val rbf = new RBFInterpolation1D(x, y, new smile.math.rbf.GaussianRadialBasis())
val data = Array.ofDim[Double](61, 2)
for (i <- 0 to 60) {
data(i)(0) = i * 0.1
data(i)(1) = rbf.interpolate(data(i)(0))
}
canvas.line("RBF", data, GREEN)
show(canvas)
```

```
var rbf = new RBFInterpolation1D(x, y, new smile.math.rbf.GaussianRadialBasis())
double[][] data = new double[61][2]
for (int i = 0; i < data.length; i++) {
data[i][0] = i * 0.1;
data[i][1] = rbf.interpolate(data[i][0]);
}
canvas.line("RBF", data, Color.GREEN)
canvas.window()
```

## Shepard Interpolation

Shepard interpolation is a special case of normalized radial basis function
interpolation if the function φ(r) goes to infinity as r → 0, and is
finite for r > 0. In this case, the weights w_{i} are just equal to
the respective function values y_{i}. So we need not solve linear
equations and thus it works for very large N. An example of such φ
is φ(r) = r^{-p} with (typically) 1 < p ≤ 3.

```
val shepard = new ShepardInterpolation1D(x, y, 3)
val data = Array.ofDim[Double](61, 2)
for (i <- 0 to 60) {
data(i)(0) = i * 0.1
data(i)(1) = shepard.interpolate(data(i)(0))
}
canvas.line("Shepard", data, PINK)
show(canvas)
```

```
var shepard = new ShepardInterpolation1D(x, y, 3)
double[][] data = new double[61][2]
for (int i = 0; i < data.length; i++) {
data[i][0] = i * 0.1;
data[i][1] = shepard.interpolate(data[i][0]);
}
canvas.line("Shepard", data, Color.PINK)
canvas.window()
```

Shepard interpolation is rarely as accurate as the well-tuned application of other radial basis functions. However, it is simple, fast, and often sufficient for quick and dirty applications.

## Laplace Interpolation

Laplace interpolation can restore missing or unmeasured values on a 2-dimensional evenly spaced regular grid. In some sense, Laplace interpolation produces the smoothest possible interpolant, which are obtained by solving a very sparse linear equations with biconjugate gradient method.

```
val zz = Array.ofDim[Double](101, 101)
val ww = Array.ofDim[Double](101, 101)
for (i <- 0 to 100; j <- 0 to 100) {
zz(i)(j) = if (java.lang.Math.random() < 0.2) Double.NaN else data(i)(j)
ww(i)(j) = zz(i)(j)
}
val canvas = heatmap(ww, Palette.jet(256))
canvas.setTitle("Missing Values")
show(canvas)
LaplaceInterpolation.interpolate(zz)
val canvas = heatmap(zz, Palette.jet(256))
canvas.setTitle("Laplace")
show(canvas)
```

```
double[][] zz = new double[101][101]
double[][] ww = new double[101][101]
for (int i = 0; i <= 100; i++) {
for (int j = 0; j <= 100; j++) {
zz[i][j] = Math.random() < 0.2 ? Double.NaN : data[i][j];
ww[i][j] = zz[i][j];
}
}
var canvas = Heatmap.plot(ww, Palette.jet(256))
canvas.setTitle("Missing Values")
canvas.window()
LaplaceInterpolation.interpolate(zz)
var canvas = Heatmap.plot(zz, Palette.jet(256))
canvas.setTitle("Laplace")
canvas.window()
```