Skip to content

Commit

Permalink
add XML documentation Interpolation
Browse files Browse the repository at this point in the history
  • Loading branch information
bvenn committed Jul 14, 2023
1 parent b7e9877 commit 7f0104e
Show file tree
Hide file tree
Showing 2 changed files with 374 additions and 29 deletions.
155 changes: 141 additions & 14 deletions src/FSharp.Stats/Interpolation/Akima.fs
Original file line number Diff line number Diff line change
Expand Up @@ -3,9 +3,16 @@
open System
open FSharp.Stats

/// <summary>
/// Module to create piecewise cubic polynomials (cubic subsplines) from x,y coordinates.
/// Akima subsplines are more flexible than standard cubic splines because the are NOT continuous in the function curvature, thereby diminishing oscillating behaviour.
/// </summary>
module Akima =

type SplineCoef = {
/// <summary>
/// Subsplines differ from regular splines because they are discontinuous in the second derivative.
/// </summary>
type SubSplineCoef = {
/// sample points (N+1), sorted ascending
XValues : float []
/// Zero order spline coefficients (N)
Expand All @@ -21,8 +28,27 @@ module Akima =
XValues=xValues;C0=c0;C1=c1;C2=c2;C3=c3
}

// For reference see: http://www.dorn.org/uni/sls/kap06/f08_01.htm
/// <summary>
/// Computes coefficients for piecewise interpolating subsplines.
/// </summary>
/// <param name="xValues">x values that are sorted ascending</param>
/// <param name="yValues">function value at x values</param>
/// <param name="firstDerivatives">first derivatives at x values</param>
/// <returns>Coefficients that define the interpolating function.</returns>
/// <example>
/// <code>
/// // e.g. days since a certain event
/// let xData = vector [|0.;1.;5.;4.;3.;|]
/// // some measured feature
/// let yData = vector [|1.;5.;4.;13.;17.|]
///
/// // get coefficients for piecewise interpolating cubic polynomials
/// let coefficients =
/// Akima.interpolate xData yData
///
/// </code>
/// </example>
/// <remarks>Second derivative (curvature) is NOT continuous at knots to allow higher flexibility to reduce oscillations! For reference see: http://www.dorn.org/uni/sls/kap06/f08_0204.htm.</remarks>
let interpolateHermiteSorted (xValues:float []) (yValues:float []) (firstDerivatives:float []) =
if xValues.Length <> yValues.Length || xValues.Length <> firstDerivatives.Length then
failwith "input arrays differ in length"
Expand All @@ -41,13 +67,30 @@ module Akima =
c1.[i] <- firstDerivatives.[i]
c2.[i] <- (3.*(yValues.[i+1] - yValues.[i])/w - 2. * firstDerivatives.[i] - firstDerivatives.[i+1])/w
c3.[i] <- (2.*(yValues.[i] - yValues.[i+1])/w + firstDerivatives.[i] + firstDerivatives.[i+1])/ww
SplineCoef.Create xValues c0 c1 c2 c3
SubSplineCoef.Create xValues c0 c1 c2 c3

let internal leftSegmentIdx arr value =
LinearSpline.leftSegmentIdx arr value

//For reference see: http://www.dorn.org/uni/sls/kap06/f08_0204.htm
/// <summary>
/// Computes coefficients for piecewise interpolating subsplines.
/// </summary>
/// <param name="yValues">function value at x values</param>
/// <returns>Coefficients that define the interpolating function.</returns>
/// <example>
/// <code>
/// // e.g. days since a certain event
/// let xData = vector [|0.;1.;5.;4.;3.;|]
/// // some measured feature
/// let yData = vector [|1.;5.;4.;13.;17.|]
///
/// // get coefficients for piecewise interpolating cubic polynomials
/// let coefficients =
/// Akima.interpolate xData yData
///
/// </code>
/// </example>
/// <remarks>Second derivative (curvature) is NOT continuous at knots to allow higher flexibility to reduce oscillations! For reference see: http://www.dorn.org/uni/sls/kap06/f08_0204.htm.</remarks>
let interpolate (xValues:float []) (yValues:float []) =
if xValues.Length <> yValues.Length then
failwith "input arrays differ in length"
Expand Down Expand Up @@ -82,29 +125,102 @@ module Akima =
tmp
interpolateHermiteSorted xValues yValues dd

/// <summary>
/// Returns function that takes x value and predicts the corresponding interpolating y value.
/// </summary>
/// <param name="splineCoeffs">Interpolation functions coefficients.</param>
/// <param name="xVal">X value of which the y value should be predicted.</param>
/// <returns>Function that takes an x value and returns function value.</returns>
/// <example>
/// <code>
/// // e.g. days since a certain event
/// let xData = vector [|0.;1.;5.;4.;3.;|]
/// // some measured feature
/// let yData = vector [|1.;5.;4.;13.;17.|]
///
/// // get coefficients for piecewise interpolating cubic polynomials
/// let coefficients =
/// Akima.interpolate xData yData
///
/// // get function to predict y value
/// let interpolFunc =
/// Akima.predict coefficients
///
let predict (splineCoeffs: SplineCoef) xVal =
/// // get function value at x=3.4
/// interpolFunc 3.4
/// </code>
/// </example>
/// <remarks>Second derivative (curvature) is NOT continuous at knots to allow higher flexibility to reduce oscillations! For reference see: http://www.dorn.org/uni/sls/kap06/f08_0204.htm.</remarks>
let predict (splineCoeffs: SubSplineCoef) xVal =
let k = leftSegmentIdx splineCoeffs.XValues xVal
let x = xVal - splineCoeffs.XValues.[k]
splineCoeffs.C0.[k] + x*(splineCoeffs.C1.[k] + x*(splineCoeffs.C2.[k] + x*splineCoeffs.C3.[k]))


/// <summary>
/// Returns function that takes x value and predicts the corresponding slope.
/// </summary>
/// <param name="splineCoeffs">Interpolation functions coefficients.</param>
/// <param name="xVal">X value of which the slope should be predicted.</param>
/// <returns>Function that takes an x value and returns slope.</returns>
/// <example>
/// <code>
/// // e.g. days since a certain event
/// let xData = vector [|0.;1.;5.;4.;3.;|]
/// // some measured feature
/// let yData = vector [|1.;5.;4.;13.;17.|]
///
let getFirstDerivative (splineCoeffs: SplineCoef) xVal =
/// // get coefficients for piecewise interpolating cubic polynomials
/// let coefficients =
/// Akima.interpolate xData yData
///
/// // get function to predict slope
/// let interpolFunc =
/// Akima.getFirstDerivative coefficients
///
/// // get slope at x=3.4
/// interpolFunc 3.4
/// </code>
/// </example>
/// <remarks>Second derivative (curvature) is NOT continuous at knots to allow higher flexibility to reduce oscillations! For reference see: http://www.dorn.org/uni/sls/kap06/f08_0204.htm.</remarks>
let getFirstDerivative (splineCoeffs: SubSplineCoef) xVal =
let k = leftSegmentIdx splineCoeffs.XValues xVal
let x = xVal - splineCoeffs.XValues.[k]
splineCoeffs.C1.[k] + x*(2.*splineCoeffs.C2.[k] + x*3.*splineCoeffs.C3.[k])


/// <summary>
/// Returns function that takes x value and predicts the corresponding interpolating curvature.
/// </summary>
/// <param name="splineCoeffs">Interpolation functions coefficients.</param>
/// <param name="xVal">X value of which the curvature should be predicted.</param>
/// <returns>Function that takes an x value and returns curvature.</returns>
/// <example>
/// <code>
/// // e.g. days since a certain event
/// let xData = vector [|0.;1.;5.;4.;3.;|]
/// // some measured feature
/// let yData = vector [|1.;5.;4.;13.;17.|]
///
/// // get coefficients for piecewise interpolating cubic polynomials
/// let coefficients =
/// Akima.interpolate xData yData
///
let getSecondDerivative (splineCoeffs:SplineCoef) xVal =
/// // get function to predict curvature
/// let interpolFunc =
/// Akima.getSecondDerivative coefficients
///
/// // get curvature at x=3.4
/// interpolFunc 3.4
/// </code>
/// </example>
/// <remarks>Second derivative (curvature) is NOT continuous at knots to allow higher flexibility to reduce oscillations! For reference see: http://www.dorn.org/uni/sls/kap06/f08_0204.htm.</remarks>
let getSecondDerivative (splineCoeffs: SubSplineCoef) xVal =
let k = leftSegmentIdx splineCoeffs.XValues xVal
let x = xVal - splineCoeffs.XValues.[k]
2.*splineCoeffs.C2.[k] + x*6.*splineCoeffs.C3.[k]


///
let computeIndefiniteIntegral (splineCoeffs:SplineCoef) =
let internal computeIndefiniteIntegral (splineCoeffs: SubSplineCoef) =
let integral =
let tmp = Array.zeroCreate splineCoeffs.C0.Length
for i = 0 to tmp.Length-2 do
Expand All @@ -113,13 +229,24 @@ module Akima =
tmp
integral

///
let integrate (splineCoeffs:SplineCoef) xVal =
/// <summary>
/// Returns integral from interpolating function from x=0 to x=xVal.
/// </summary>
/// <param name="splineCoeffs">Interpolation functions coefficients.</param>
/// <param name="xVal">X value up to which the integral should be calculated.</param>
/// <returns>Integral (area under the curve) from x=0 to x=xVal</returns>
let integrate (splineCoeffs: SubSplineCoef) xVal =
let integral = computeIndefiniteIntegral splineCoeffs
let k = leftSegmentIdx splineCoeffs.XValues xVal
let x = xVal - splineCoeffs.XValues.[k]
integral.[k] + x*(splineCoeffs.C0.[k] + x*(splineCoeffs.C1.[k]/2. + x*(splineCoeffs.C2.[k]/3. + x*splineCoeffs.C3.[k]/4.)))

///

/// <summary>
/// Returns integral from interpolating function from x=xVal1 to x=xVal2.
/// </summary>
/// <param name="splineCoeffs">Interpolation functions coefficients.</param>
/// <param name="xVal1">X value from where the integral should be calculated.</param>
/// <param name="xVal2">X value up to which the integral should be calculated.</param>
/// <returns>Integral (area under the curve) from x=xVal1 to x=xVal2</returns>
let definiteIntegral (integrateF: float -> float) xVal1 xVal2 =
integrateF xVal2 - integrateF xVal1
Loading

0 comments on commit 7f0104e

Please sign in to comment.