This library interpolates cubic splines for one-dimensional sets of data.
It computes "natural splines", which means the second derivative at the endpoints is set to zero.
For more information on the mathematics used, check out this paper.
Compute a spline, given the minimum and maximum values of x
and a
list of data for f(x)
. The data should be evenly spaced and in order of
increasing x
.
For example, if we have the data
f(2) = 1
f(3) = 5.2
f(4) = 3.2
f(5) = 0.8
Then we would generate a spline by calling
fSpline = withRange 2 6 [ 1, 5.2, 3.2, 0.8 ]
If there is only one data point, then the spline will be a horizontal line passing through that point. If the data is empty, the spline will be zero everywhere.
Same as withRange
, except instead of passing the endpoint as the
second argument, you pass the x-distance between data points.
fSpline = withDelta 2 1 [ 1, 5.2, 3.2, 0.8 ]
-- equivalent to withRange 2 6 [ 1, 5.2, 3.2, 0.8 ]
Given x1
and a spline f(x)
, return f(x1)
Return the first derivative of the curve at the given point
Return the second derivative of the curve at the given point
Return a record containing the value, slope, and concavity of the curve at the given point
module Interpolate.Cubic
exposing
( Spline
, LocalCurve
, withRange
, withDelta
, valueAt
, slopeAt
, concavityAt
, curveAt
)
{-| This library interpolates cubic splines for one-dimensional sets of data.
It computes "natural splines", which means the second derivative at the endpoints
is set to zero.
For more information on the mathematics used, check out
[this paper.](http://web.archive.org/web/20090408054627/http://online.redwoods.cc.ca.us/instruct/darnold/laproj/Fall98/SkyMeg/Proj.PDF)
# Creating splines
@docs withRange, withDelta, Spline
# Interpolating
@docs valueAt, slopeAt, concavityAt, curveAt, LocalCurve
-}
import Array exposing (Array)
{-| Given `x1` and a spline `f(x)`, return `f(x1)`
-}
valueAt : Float -> Spline -> Float
valueAt =
evaluate cubic
{-| Return the first derivative of the curve at the given point
-}
slopeAt : Float -> Spline -> Float
slopeAt =
evaluate firstDerivative
{-| Return the second derivative of the curve at the given point
-}
concavityAt : Float -> Spline -> Float
concavityAt =
evaluate secondDerivative
{-| Return a record containing the value, slope, and concavity of the curve
at the given point
-}
curveAt : Float -> Spline -> LocalCurve
curveAt =
let
curve x coeff =
{ value = cubic x coeff
, slope = firstDerivative x coeff
, concavity = secondDerivative x coeff
}
in
evaluate curve
cubic : Float -> Coefficients -> Float
cubic x { a, b, c, d } =
a * x ^ 3 + b * x ^ 2 + c * x + d
firstDerivative : Float -> Coefficients -> Float
firstDerivative x { a, b, c, d } =
3 * a * x ^ 2 + 2 * b * x + c
secondDerivative : Float -> Coefficients -> Float
secondDerivative x { a, b, c, d } =
6 * a * x + 2 * b
evaluate : (Float -> Coefficients -> a) -> Float -> Spline -> a
evaluate f x (Spline spline) =
let
maxIndex =
Array.length spline.coefficients - 1
index =
(x - spline.start)
/ spline.dx
|> floor
|> clamp 0 maxIndex
offset =
x - spline.dx * (toFloat index) - spline.start
coeff =
Array.get index spline.coefficients
|> Maybe.withDefault { a = 0, b = 0, c = 0, d = 0 }
in
f offset coeff
{-| -}
type alias LocalCurve =
{ value : Float
, slope : Float
, concavity : Float
}
{-| Compute a spline, given the minimum and maximum values of `x` and a
list of data for `f(x)`. The data should be evenly spaced and in order of
increasing `x`.
For example, if we have the data
f(2) = 1
f(3) = 5.2
f(4) = 3.2
f(5) = 0.8
Then we would generate a spline by calling
fSpline = withRange 2 6 [ 1, 5.2, 3.2, 0.8 ]
If there is only one data point, then the spline will be a horizontal line
passing through that point. If the data is empty, the spline will be zero
everywhere.
-}
withRange : Float -> Float -> List Float -> Spline
withRange start end heights =
let
n =
List.length heights |> toFloat
dx =
(end - start) / (n - 1)
in
if 1 < n then
withDelta start dx heights
else
{ coefficients = degenerateCoefficients heights
, start = 0
, dx = 0
}
|> Spline
{-| Same as `withRange`, except instead of passing the endpoint as the
second argument, you pass the x-distance between data points.
fSpline = withDelta 2 1 [ 1, 5.2, 3.2, 0.8 ]
-- equivalent to withRange 2 6 [ 1, 5.2, 3.2, 0.8 ]
-}
withDelta : Float -> Float -> List Float -> Spline
withDelta start dx heights =
{ coefficients = findCoefficients dx heights
, start = start
, dx = dx
}
|> Spline
degenerateCoefficients : List Float -> Array Coefficients
degenerateCoefficients heights =
{ a = 0
, b = 0
, c = 0
, d = List.head heights |> Maybe.withDefault 0
}
|> Array.repeat 1
findCoefficients : Float -> List Float -> Array Coefficients
findCoefficients dx heights =
let
concavity y0 y1 y2 =
6 * (y0 + y2 - 2 * y1) / dx ^ 2
sweep y ( uPrev, vPrev ) =
( 1 / (4 - uPrev), (y - vPrev) / (4 - uPrev) )
backSub ( u, v ) mNext =
v - u * mNext
piece ( y0, m0 ) ( y1, m1 ) =
{ a = (m1 - m0) / 6 / dx
, b = m0 / 2
, c = (y1 - y0) / dx - (m1 + 2 * m0) * dx / 6
, d = y0
}
in
mapTriple concavity heights
|> List.scanl sweep ( 0, 0 )
|> scanr backSub 0
|> List.map2 (,) heights
|> mapPair piece
|> Array.fromList
mapTriple : (a -> a -> a -> b) -> List a -> List b
mapTriple f x0 =
let
tail =
List.tail x0
doubleTail =
tail |> Maybe.andThen List.tail
in
case ( tail, doubleTail ) of
( Just x1, Just x2 ) ->
List.map3 f x0 x1 x2
otherwise ->
[]
mapPair : (a -> a -> b) -> List a -> List b
mapPair f x0 =
case (List.tail x0) of
Just x1 ->
List.map2 f x0 x1
otherwise ->
[]
scanr : (a -> b -> b) -> b -> List a -> List b
scanr f init =
List.reverse >> List.scanl f init >> List.reverse
{-| -}
type Spline
= Spline
{ coefficients : Array Coefficients
, start : Float
, dx : Float
}
type alias Coefficients =
{ a : Float
, b : Float
, c : Float
, d : Float
}