Source Edit

Interpolation

This module implements various interpolation routines. See also:

rbf module for RBF interpolation of scattered data in arbitrary dimensions.

1D interpolation

Hermite spline (recommended): cubic spline that works with many types of values. Accepts derivatives if available.
Cubic spline: cubic spline that only works with floats.
Linear spline: Linear spline that works with many types of values.

Extrapolation is supported for all 1D interpolators by passing the type of extrapolation as an argument of eval. The default is to use the interpolator's native method to extrapolate. This means that Linear does linear extrapolation, while Hermite and Cubic performs cubic extrapolation. Other options are using a constant value, using the values of the edges, linear extrapolation and raising an error if the x-value is outside the domain.

Example:

import src/numericalnim/interpolate
import numericalnim, std/[math, sequtils]

let x = linspace(0.0, 1.0, 10)
let y = x.mapIt(sin(it))

let interp = newHermiteSpline(x, y)

let val = interp.eval(0.314)

let valExtrapolate = interp.eval(1.1, Edge)

2D interpolation

Bicubic: Works on gridded data.
Bilinear: Works on gridded data.
Nearest neighbour: Works on gridded data.
Barycentric: Works on scattered data.

Example:

import src/numericalnim/interpolate
import arraymancer, numericalnim

var z = newTensor[float](10, 10)
for i, x in linspace(0.0, 1.0, 10):
  for j, y in linspace(0.0, 1.0, 10):
    z[i, j] = x * y

let xlim = (0.0, 1.0)
let ylim = (0.0, 1.0)
let interp = newBicubicSpline(z, xlim, ylim)

let val = interp.eval(0.314, 0.628)

3D interpolation

Trilinear: works on gridded data

Example:

import src/numericalnim/interpolate
import arraymancer, numericalnim

var f = newTensor[float](10, 10, 10)
for i, x in linspace(0.0, 1.0, 10):
  for j, y in linspace(0.0, 1.0, 10):
    for k, z in linspace(0.0, 1.0, 10):
      f[i, j, k] = x * y * z

let xlim = (0.0, 1.0)
let ylim = (0.0, 1.0)
let zlim = (0.0, 1.0)
let interp = newTrilinearSpline(f, xlim, ylim, zlim)

let val = interp.eval(0.314, 0.628, 0.414)

Imports

utils, common/commonTypes, private/arraymancerOverloads, rbf

Types

Eval2DHandler[T] = proc (self: Interpolator2DType[T]; x, y: float): T {.nimcall.}: Source Edit
Eval3DHandler[T] = proc (self: Interpolator3DType[T]; x, y, z: float): T {. nimcall.}: Source Edit
EvalHandler[T] = proc (self: InterpolatorType[T]; x: float): T {.nimcall.}: Source Edit
EvalUnstructured2DHandler[T; U] = proc (self: InterpolatorUnstructured2DType[T, U]; x, y: float): U {.nimcall.}: Source Edit
ExtrapolateKind = enum Constant, Edge, Linear, Native, Error: Source Edit
Interpolator2DType[T] = ref object z*, xGrad*, yGrad*, xyGrad*: Tensor[T] alphaCache*: Table[(int, int), Tensor[T]] dx*, dy*: float xLim*, yLim*: tuple[lower: float, upper: float] eval_handler*: Eval2DHandler[T]: Source Edit
Interpolator3DType[T] = ref object f*: Tensor[T] dx*, dy*, dz*: float xLim*, yLim*, zLim*: tuple[lower: float, upper: float] eval_handler*: Eval3DHandler[T]: Source Edit
InterpolatorType[T] = ref object X*: seq[float] Y*: seq[T] coeffs_f*: Tensor[float] coeffs_T*: Tensor[T] high*: int len*: int eval_handler*: EvalHandler[T] deriveval_handler*: EvalHandler[T]: Source Edit
InterpolatorUnstructured2DType[T; U] = ref object values*: Tensor[T] points*: Tensor[U] dt*: DelaunayTriangulation z*, gradX*, gradY*: Table[(T, T), U] boundPoints*: array[4, Vector2] eval_handler*: EvalUnstructured2DHandler[T, U]: Source Edit

Procs

proc derivEval[T; U](interpolator: InterpolatorType[T]; x: float; extrap: ExtrapolateKind = Native; extrapValue: U = missing()): T: Evaluates the derivative of an interpolator.
x: The point to evaluate the interpolator at.

extrap: The extrapolation method to use. Available options are:

Constant: Set all points outside the range of the interpolator to extrapValue.

Edge: Use the value of the left/right edge.

Linear: Uses linear extrapolation using the two points closest to the edge.

Native (default): Uses the native method of the interpolator to extrapolate. For Linear1D it will be a linear extrapolation, and for Cubic and Hermite splines it will be cubic extrapolation.

Error: Raises an ValueError if x is outside the range.

extrapValue: The extrapolation value to use when extrap = Constant.

Beware: Native extrapolation for the cubic splines can very quickly diverge if the extrapolation is too far away from the interpolation points.
Source Edit
proc derivEval[T; U](spline: InterpolatorType[T]; x: openArray[float]; extrap: ExtrapolateKind = Native; extrapValue: U = missing()): seq[T]: Evaluates the derivative of an interpolator at all points in x. Source Edit
proc derivEval_cubicspline[T](spline: InterpolatorType[T]; x: float): T: Source Edit
proc derivEval_hermitespline[T](spline: InterpolatorType[T]; x: float): T: Source Edit
proc derivEval_linear1d[T](spline: InterpolatorType[T]; x: float): T: Source Edit
proc eval[T; U](interpolator: InterpolatorType[T]; x: float; extrap: ExtrapolateKind = Native; extrapValue: U = missing()): T: Evaluates an interpolator.
x: The point to evaluate the interpolator at.

extrap: The extrapolation method to use. Available options are:

Constant: Set all points outside the range of the interpolator to extrapValue.

Edge: Use the value of the left/right edge.

Linear: Uses linear extrapolation using the two points closest to the edge.

Native (default): Uses the native method of the interpolator to extrapolate. For Linear1D it will be a linear extrapolation, and for Cubic and Hermite splines it will be cubic extrapolation.

Error: Raises an ValueError if x is outside the range.

extrapValue: The extrapolation value to use when extrap = Constant.

Beware: Native extrapolation for the cubic splines can very quickly diverge if the extrapolation is too far away from the interpolation points.
Source Edit
proc eval[T; U](spline: InterpolatorType[T]; x: openArray[float]; extrap: ExtrapolateKind = Native; extrapValue: U = missing()): seq[ T]: Evaluates an interpolator at all points in x. Source Edit
proc eval_barycentric2d[T, U](self: InterpolatorUnstructured2DType[T, U]; x, y: float): U: Source Edit
proc eval_bicubic[T](self: Interpolator2DType[T]; x, y: float): T {.nimcall.}: Source Edit
proc eval_bilinear[T](self: Interpolator2DType[T]; x, y: float): T {.nimcall.}: Source Edit
proc eval_cubicspline[T](spline: InterpolatorType[T]; x: float): T: Source Edit
proc eval_hermitespline[T](spline: InterpolatorType[T]; x: float): T: Source Edit
proc eval_linear1d[T](spline: InterpolatorType[T]; x: float): T: Source Edit
proc eval_nearestneigh[T](self: Interpolator2DType[T]; x, y: float): T {.nimcall.}: Source Edit
proc eval_trilinear[T](self: Interpolator3DType[T]; x, y, z: float): T {.nimcall.}: Source Edit
proc findInterval(list: openArray[float]; x: float): int {.inline, ...raises: [], tags: [], forbids: [].}: Source Edit
proc newBarycentric2D[T: SomeFloat; U](points: Tensor[T]; values: Tensor[U]): InterpolatorUnstructured2DType[ T, U]: Barycentric interpolation of scattered points in 2D.

Inputs:

points: Tensor of shape (nPoints, 2) with the coordinates of all points.

values: Tensor of shape (nPoints) with the function values.

Returns:

Interpolator object that can be evaluated using interp.eval(x, y.

Source Edit
proc newBicubicSpline[T](z: Tensor[T]; xlim, ylim: (float, float)): Interpolator2DType[ T]: Returns a bicubic spline for regularly gridded data.
z: Tensor with the function values. x corrensponds to the rows and y to the columns. Must be sorted so ascendingly in both variables.

xlim: the lowest and highest x-value

ylim: the lowest and highest y-value

Source Edit
proc newBilinearSpline[T](z: Tensor[T]; xlim, ylim: (float, float)): Interpolator2DType[ T]: Returns a bilinear spline for regularly gridded data.
z: Tensor with the function values. x corrensponds to the rows and y to the columns. Must be sorted so ascendingly in both variables.

xlim: the lowest and highest x-value

ylim: the lowest and highest y-value

Source Edit
proc newCubicSpline[T: SomeFloat](X: openArray[float]; Y: openArray[T]): InterpolatorType[ T]: Returns a cubic spline. Source Edit
proc newHermiteSpline[T](X: openArray[float]; Y, dY: openArray[T]): InterpolatorType[ T]: Constructs a cubic Hermite spline using x, y and derivatives of y. Source Edit
proc newHermiteSpline[T](X: openArray[float]; Y: openArray[T]): InterpolatorType[ T]: Constructs a cubic Hermite spline by approximating the derivatives. Source Edit
proc newLinear1D[T](X: openArray[float]; Y: openArray[T]): InterpolatorType[T]: Constructs a linear interpolator. Source Edit
proc newNearestNeighbour2D[T](z: Tensor[T]; xlim, ylim: (float, float)): Interpolator2DType[ T]: Returns a nearest neighbour interpolator for regularly gridded data.
z: Tensor with the function values. x corrensponds to the rows and y to the columns. Must be sorted so ascendingly in both variables.

xlim: the lowest and highest x-value

ylim: the lowest and highest y-value

Source Edit
proc newTrilinearSpline[T](f: Tensor[T]; xlim, ylim, zlim: (float, float)): Interpolator3DType[ T]: Returns a trilinear spline for regularly gridded data.
z: Tensor with the function values. x corrensponds to the first dimension, y to the second and z to the third. Must be sorted so ascendingly in both variables.

xlim: the lowest and highest x-value

ylim: the lowest and highest y-value

zlim: the lowest and highest z-value

Source Edit
proc toDerivNumContextProc[T](spline: InterpolatorType[T]): NumContextProc[T, float]: Convert interpolator's derivative to NumContextProc. Source Edit
proc toDerivProc[T](spline: InterpolatorType[T]): InterpolatorProc[T]: Returns a proc to evaluate the derivative of the interpolator. Source Edit
proc toProc[T](spline: InterpolatorType[T]): InterpolatorProc[T]: Returns a proc to evaluate the interpolator. Source Edit

Converters

converter toNumContextProc[T](spline: InterpolatorType[T]): NumContextProc[T, float]: Convert interpolator to NumContextProc. Source Edit

Templates

template eval[T, U](interpolator: InterpolatorUnstructured2DType[T, U]; x, y: T): untyped: Evaluate interpolator. Source Edit
template eval[T](interpolator: Interpolator2DType[T]; x, y: float): untyped: Evaluate interpolator. Source Edit
template eval[T](interpolator: Interpolator3DType[T]; x, y, z: float): untyped: Evaluate interpolator. Source Edit