diff --git a/src/FSharp.Stats/Interpolation/Akima.fs b/src/FSharp.Stats/Interpolation/Akima.fs index d9805459..6bbc0434 100644 --- a/src/FSharp.Stats/Interpolation/Akima.fs +++ b/src/FSharp.Stats/Interpolation/Akima.fs @@ -3,9 +3,16 @@ open System open FSharp.Stats +/// +/// 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. +/// module Akima = - type SplineCoef = { + /// + /// Subsplines differ from regular splines because they are discontinuous in the second derivative. + /// + type SubSplineCoef = { /// sample points (N+1), sorted ascending XValues : float [] /// Zero order spline coefficients (N) @@ -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 + /// + /// Computes coefficients for piecewise interpolating subsplines. + /// + /// x values that are sorted ascending + /// function value at x values + /// first derivatives at x values + /// Coefficients that define the interpolating function. + /// + /// + /// // 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 + /// + /// + /// + /// 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. let interpolateHermiteSorted (xValues:float []) (yValues:float []) (firstDerivatives:float []) = if xValues.Length <> yValues.Length || xValues.Length <> firstDerivatives.Length then failwith "input arrays differ in length" @@ -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 + /// + /// Computes coefficients for piecewise interpolating subsplines. + /// + /// function value at x values + /// Coefficients that define the interpolating function. + /// + /// + /// // 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 /// + /// + /// + /// 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. let interpolate (xValues:float []) (yValues:float []) = if xValues.Length <> yValues.Length then failwith "input arrays differ in length" @@ -82,29 +125,102 @@ module Akima = tmp interpolateHermiteSorted xValues yValues dd + /// + /// Returns function that takes x value and predicts the corresponding interpolating y value. + /// + /// Interpolation functions coefficients. + /// X value of which the y value should be predicted. + /// Function that takes an x value and returns function value. + /// + /// + /// // 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 + /// + /// + /// 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. + 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])) + /// + /// Returns function that takes x value and predicts the corresponding slope. + /// + /// Interpolation functions coefficients. + /// X value of which the slope should be predicted. + /// Function that takes an x value and returns slope. + /// + /// + /// // 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 + /// + /// + /// 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. + 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]) + /// + /// Returns function that takes x value and predicts the corresponding interpolating curvature. + /// + /// Interpolation functions coefficients. + /// X value of which the curvature should be predicted. + /// Function that takes an x value and returns curvature. + /// + /// + /// // 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 + /// + /// + /// 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. + 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 @@ -113,13 +229,24 @@ module Akima = tmp integral - /// - let integrate (splineCoeffs:SplineCoef) xVal = + /// + /// Returns integral from interpolating function from x=0 to x=xVal. + /// + /// Interpolation functions coefficients. + /// X value up to which the integral should be calculated. + /// Integral (area under the curve) from x=0 to x=xVal + 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.))) - - /// + + /// + /// Returns integral from interpolating function from x=xVal1 to x=xVal2. + /// + /// Interpolation functions coefficients. + /// X value from where the integral should be calculated. + /// X value up to which the integral should be calculated. + /// Integral (area under the curve) from x=xVal1 to x=xVal2 let definiteIntegral (integrateF: float -> float) xVal1 xVal2 = integrateF xVal2 - integrateF xVal1 diff --git a/src/FSharp.Stats/Interpolation/CubicSpline.fs b/src/FSharp.Stats/Interpolation/CubicSpline.fs index cb7beec8..3312c69d 100644 --- a/src/FSharp.Stats/Interpolation/CubicSpline.fs +++ b/src/FSharp.Stats/Interpolation/CubicSpline.fs @@ -6,7 +6,7 @@ module CubicSpline = open FSharp.Stats type BoundaryCondition = - ///most used spline variant: f'' at borders is set to 0 + ///most used spline variant: f'' at borders is set to 0 | Natural ///f' at first point is the same as f' at the last point | Periodic @@ -18,14 +18,37 @@ module CubicSpline = | Quadratic ///f' at first and last knot are set by user | Clamped - + + /// + /// Contains x data, y data (c0), slopes (c1), curvatures (c2), and the third derivative at each knot + /// type CubicSplineCoef = { XData : vector + /// + /// vector of [a0;b0;c0;d0;a1;b1;...;d(n-2)] where f_n(x) = (an)x^3 + (bn)x^2 + (cn)x + (dn) + /// C0_3 : vector} with static member Create x c = {XData=x;C0_3=c} - /// computes all coefficients for piecewise interpolating splines. In the form of [a0;b0;c0;d0;a1;b1;...;d(n-2)]. - /// where: fn(x) = (an)x^3+(bn)x^2+(cn)x+(dn) + /// + /// Computes coefficients for piecewise interpolating splines. In the form of [a0;b0;c0;d0;a1;b1;...;d(n-2)]. + /// where: fn(x) = (an)x^3+(bn)x^2+(cn)x+(dn) + /// + /// Condition that defines how slopes and curvatures are defined at the left- and rightmost knot + /// function value at x values + /// Coefficients that define the interpolating function. + /// + /// + /// // 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 = + /// CubicSpline.interpolate CubicSpline.BoundaryConditions.Natural xData yData + /// + /// let interpolate (boundaryCondition: BoundaryCondition) (xValues: Vector) (yValues: Vector) = //f(x) = ax³+bx²+cx+d //f'(x) = 3ax²+2bx+c @@ -196,8 +219,33 @@ module CubicSpline = let coeffs = Algebra.LinearAlgebra.SolveLinearSystem A b CubicSplineCoef.Create xValues coeffs - - ///Fits a cubic spline with given coefficients. Only defined within the range of the given xValues + + /// + /// Returns function that takes x value (that lies within the range of input x values) and predicts the corresponding interpolating y value. + /// + /// Interpolation functions coefficients. + /// X value of which the y value should be predicted. + /// Function that takes an x value and returns function value. + /// + /// + /// // 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 = + /// CubicSpline.interpolate CubicSpline.BoundaryConditions.Natural xData yData + /// + /// // get interpolating value + /// let func = CubicSpline.predictWithinRange(coefLinSpl) + /// + /// // get function value at x=3.4 + /// func 3.4 + /// + /// + /// + /// Only defined within the range of the given xValues! let predictWithinRange (coefficients: CubicSplineCoef) x = let sortedX = coefficients.XData |> Seq.sort let intervalNumber = @@ -220,7 +268,32 @@ module CubicSpline = yValue - ///fits a cubic spline fit even outside of the interval defined in xValues by linear interpolation of slope given by border knots + /// + /// Returns function that takes x value and predicts the corresponding interpolating y value. + /// + /// Interpolation functions coefficients. + /// X value of which the y value should be predicted. + /// Function that takes an x value and returns function value. + /// + /// + /// // 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 = + /// CubicSpline.interpolate CubicSpline.BoundaryConditions.Natural xData yData + /// + /// // get interpolating function + /// let func = CubicSpline.predict(coefLinSpl) + /// + /// // get function value at x=3.4 + /// func 3.4 + /// + /// + /// + /// x values outside of the xValue range are predicted by straight lines defined by the nearest knot! let predict (coefficients: CubicSplineCoef) x = let sortedX = coefficients.XData |> Seq.sort let xHead = coefficients.XData |> Seq.head @@ -279,15 +352,96 @@ module CubicSpline = thirdDerivative | _ -> failwithf "for cubic splines no derivative > 3 is defined" + /// + /// Returns function that takes x value and predicts the corresponding slope. + /// + /// Interpolation functions coefficients. + /// X value of which the slope should be predicted. + /// Function that takes an x value and returns the function slope. + /// + /// + /// // 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 = + /// CubicSpline.interpolate CubicSpline.BoundaryConditions.Natural xData yData + /// + /// // get slope function + /// let func = CubicSpline.getFirstDerivative(coefLinSpl) + /// + /// // get slope at x=3.4 + /// func 3.4 + /// + /// + /// + /// x values outside of the xValue range are predicted by straight lines defined by the nearest knot! let getFirstDerivative (coefficients: CubicSplineCoef) x = getDerivative 1 coefficients x + /// + /// Returns function that takes x value and predicts the corresponding curvature. + /// + /// Interpolation functions coefficients. + /// X value of which the curvature should be predicted. + /// Function that takes an x value and returns the function curvature. + /// + /// + /// // 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 = + /// CubicSpline.interpolate CubicSpline.BoundaryConditions.Natural xData yData + /// + /// // get curvature function + /// let func = CubicSpline.getSecondDerivative(coefLinSpl) + /// + /// // get curvature at x=3.4 + /// func 3.4 + /// + /// + /// + /// x values outside of the xValue range are predicted by straight lines defined by the nearest knot! let getSecondDerivative (coefficients: CubicSplineCoef) x = getDerivative 2 coefficients x + /// + /// Returns function that takes x value and predicts the corresponding third derivative. + /// + /// Interpolation functions coefficients. + /// X value of which the y value should be predicted. + /// Function that takes an x value and returns the function third derivative. + /// + /// + /// // 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 = + /// CubicSpline.interpolate CubicSpline.BoundaryConditions.Natural xData yData + /// + /// // get third derivative function + /// let func = CubicSpline.getThirdDerivative(coefLinSpl) + /// + /// // get third derivative at x=3.4 + /// func 3.4 + /// + /// + /// + /// x values outside of the xValue range are predicted by straight lines defined by the nearest knot! let getThirdDerivative (coefficients: CubicSplineCoef) x = getDerivative 3 coefficients x + /// + /// Hermite cubic splines are defined by the function values and their slopes (first derivatives). If the slopws are unknown, they must be estimated. + /// module Hermite = type HermiteCoef = { @@ -302,6 +456,26 @@ module CubicSpline = XValues=xValues;YValues=yValues;Slopes=slopes } + /// + /// Computes coefficients for piecewise interpolating splines. + /// The x data has to be sorted ascending. + /// + /// function value at x values + /// Coefficients that define the interpolating function. + /// + /// + /// // 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 = + /// CubicSpline.Hermite.interpolate xData yData + /// + /// + /// + /// Second derivative (curvature) is NOT continuous at knots to allow higher flexibility to reduce oscillations! let interpolate (xData: Vector) (yData: Vector) = let slopes = Vector.init xData.Length (fun i -> @@ -316,7 +490,27 @@ module CubicSpline = ) HermiteCoef.Create xData yData slopes - ///if the knots are monotone in/decreasing, the spline also is monotone (http://www.korf.co.uk/spline.pdf) + /// + /// Computes coefficients for piecewise interpolating splines. + /// If the knots are monotone in/decreasing, the spline also is monotone (http://www.korf.co.uk/spline.pdf) + /// The x data has to be sorted ascending + /// + /// function value at x values + /// Coefficients that define the interpolating function. + /// + /// + /// // e.g. days since a certain event + /// let xData = vector [|0.;1.;5.;4.;3.;|] + /// // some measured feature + /// let yData = vector [|1.;5.;6.;13.;13.1|] + /// + /// // get coefficients for piecewise interpolating cubic polynomials + /// let coefficients = + /// CubicSpline.Hermite.interpolateTryMonotonicity xData yData + /// + /// + /// + /// Second derivative (curvature) is NOT continuous at knots to allow higher flexibility to reduce oscillations! let interpolateTryMonotonicity (xData: Vector) (yData: Vector) = let calcSlope i = let s1 = (xData.[i+1] - xData.[i]) / (yData.[i+1] - yData.[i]) @@ -345,8 +539,32 @@ module CubicSpline = let slopes = loop 1 [slopeAtFstKnot] HermiteCoef.Create xData yData (slopes |> vector) - ///calculates a function to interpolate between the datapoints with given slopes (yData'). - ///the data has to be sorted ascending + + /// + /// Returns function that takes x value and predicts the corresponding interpolating y value. + /// + /// Interpolation functions coefficients. + /// X value of which the y value should be predicted. + /// Function that takes an x value and returns function value. + /// + /// + /// // 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 = + /// CubicSpline.Hermite.interpolate xData yData + /// + /// // get interpolating function + /// let func = CubicSpline.Hermite.predict coefLinSpl + /// + /// // get function value at x=3.4 + /// func 3.4 + /// + /// + /// x values outside of the xValue range are predicted by straight lines defined by the nearest knot! let predict (coef: HermiteCoef) x = let n = coef.XValues.Length @@ -469,23 +687,23 @@ module CubicSpline = Akima.interpolate xValues yValues [] - let interpolateAtX (splineCoeffs: Akima.SplineCoef) xVal = + let interpolateAtX (splineCoeffs: Akima.SubSplineCoef) xVal = Akima.predict splineCoeffs xVal [] - let firstDerivative (splineCoeffs: Akima.SplineCoef) xVal = + let firstDerivative (splineCoeffs: Akima.SubSplineCoef) xVal = Akima.getFirstDerivative splineCoeffs xVal [] - let secondDerivative (splineCoeffs: Akima.SplineCoef) xVal = + let secondDerivative (splineCoeffs: Akima.SubSplineCoef) xVal = Akima.getSecondDerivative splineCoeffs xVal [] - let computeIndefiniteIntegral (splineCoeffs: Akima.SplineCoef) = + let computeIndefiniteIntegral (splineCoeffs: Akima.SubSplineCoef) = Akima.computeIndefiniteIntegral splineCoeffs [] - let integrate (splineCoeffs: Akima.SplineCoef) xVal = + let integrate (splineCoeffs: Akima.SubSplineCoef) xVal = Akima.integrate splineCoeffs xVal []