RBF interpolation with example from wikipedia https://en.wikipedia.org/wiki/Radial_basis_function_interpolation

rbf_interpolation_1d.nim

import ggplotnim, math
import arraymancer

const ε = 3
proc φ(r: float): float =
  result = exp(-pow((ε.float * r), 2.0))

proc toMatrix(n: int, start, stop: float): Tensor[float] =
  result = zeros[float]([n, n])
  let xs = linspace(start, stop, n)
  for i in 0 ..< n:
    for j in 0 ..< n:
      result[i, j] = φ(abs(xs[i] - xs[j]))

proc interp(ws, xs, xFine: Tensor[float]): Tensor[float] =
  result = zeros[float](xFine.size.int)
  for i in 0 ..< xFine.size:
    var res = 0.0
    for j in 0 ..< ws.size:
      let wk = ws[j]
      res += wk * φ(abs(xFine[i] - xs[j]))
    result[i] = res

proc rbf(xs, ys: Tensor[float]): Tensor[float] =
  doAssert xs.rank == 1
  let A = toMatrix(xs.size, xs.min, xs.max)
  result = solve(A, ys)

proc testFn(x: float): float =
  result = exp(x * cos(3 * PI * x))

# generate a grid to interpolate from
let xs = linspace(0.0, 1.0, 15)
let ys = xs.map_inline(testFn(x))
# and a fine mesh for comparison
let xFine = linspace(0.0, 1.0, 1000)
let yFine = xFine.map_inline(testFn(x))

let ws = rbf(xs, ys)
let yInterp = interp(ws, xs, xFine)

let df = toDf(xs, ys)
let dfFine = toDf({ "x" : xFine,
                    "yReal" : yFine,
                    "yInterp" : yInterp })
  .gather(["yReal", "yInterp"], "type", "y")
ggplot(dfFine, aes("x", "y", color = "type")) +
  geom_line() +
  geom_point(data = df, aes = aes("xs", "ys"), color = "black") +
  ggsave("/t/rbf_interp.png")