add clamped spline interpolation

fslaborg · Aug 2, 2023 · 8ab0975 · 8ab0975
1 parent 7de8b82
commit 8ab0975
Show file tree

Hide file tree

Showing 3 changed files with 26 additions and 18 deletions.
diff --git a/README.md b/README.md
@@ -22,13 +22,18 @@ FSharp.Stats is a multipurpose project for statistical testing, linear algebra,
 
 #### Fitting
   - <a href="https://fslab.org/FSharp.Stats/Fitting.html#Linear-Regression">Linear regression</a>
+    - <a href="https://fslab.org/FSharp.Stats/Fitting.html#Simple-Linear-Regression">Simple linear regression (weighted and constrained)</a>
+    - <a href="https://fslab.org/FSharp.Stats/Fitting.html#Polynomial-Regression">Polynomial regression (weighted and constrained)</a>
   - <a href="https://fslab.org/FSharp.Stats/Fitting.html#Nonlinear-Regression">Nonlinear regression</a>
   - <a href="https://fslab.org/FSharp.Stats/Fitting.html#Smoothing-spline">Spline regression</a>
   - <a href="https://fslab.org/FSharp.Stats/GoodnessOfFit.html">Goodness of fit</a>
 
 #### Interpolation
+  - <a href="https://fslab.org/FSharp.Stats/Interpolation.html#Polynomial-Interpolation">Linear spline interpolation</a>
   - <a href="https://fslab.org/FSharp.Stats/Interpolation.html#Polynomial-Interpolation">Polynomial interpolation</a>
-  - <a href="https://fslab.org/FSharp.Stats/Interpolation.html#Cubic-interpolating-Spline">Spline interpolation</a>
+  - <a href="https://fslab.org/FSharp.Stats/Interpolation.html#Cubic-interpolating-Spline">Cubic spline interpolation</a>
+  - <a href="https://fslab.org/FSharp.Stats/Interpolation.html">Akima subspline interpolation</a>
+  - <a href="https://fslab.org/FSharp.Stats/Interpolation.html">Hermite subspline interpolation</a>
 
 
 #### Signal processing

diff --git a/docs/Interpolation.fsx b/docs/Interpolation.fsx
@@ -67,13 +67,13 @@ let testDataX = [|1. .. 10.|]
 
 let testDataY = [|0.;-1.;0.;0.;0.;0.;1.;1.;3.;3.5|]
 
-
-let coefLinear     = Interpolation.interpolate(testDataX,testDataY,InterpolationMethod.LinearSpline)
-let coefAkima      = Interpolation.interpolate(testDataX,testDataY,InterpolationMethod.AkimaSubSpline)
-let coefCubicNa    = Interpolation.interpolate(testDataX,testDataY,InterpolationMethod.CubicSpline Interpolation.CubicSpline.BoundaryCondition.Natural)
-let coefCubicPe    = Interpolation.interpolate(testDataX,testDataY,InterpolationMethod.CubicSpline Interpolation.CubicSpline.BoundaryCondition.Periodic)
-let coefCubicNo    = Interpolation.interpolate(testDataX,testDataY,InterpolationMethod.CubicSpline Interpolation.CubicSpline.BoundaryCondition.NotAKnot)
-let coefCubicPa    = Interpolation.interpolate(testDataX,testDataY,InterpolationMethod.CubicSpline Interpolation.CubicSpline.BoundaryCondition.Parabolic)
+let coefLinear     = Interpolation.interpolate(testDataX,testDataY,InterpolationMethod.LinearSpline) // Straight lines passing all points
+let coefAkima      = Interpolation.interpolate(testDataX,testDataY,InterpolationMethod.AkimaSubSpline) // Akima cubic subspline
+let coefCubicNa    = Interpolation.interpolate(testDataX,testDataY,InterpolationMethod.CubicSpline Interpolation.CubicSpline.BoundaryCondition.Natural) // cubic spline with f'' at borders is set to 0
+let coefCubicPe    = Interpolation.interpolate(testDataX,testDataY,InterpolationMethod.CubicSpline Interpolation.CubicSpline.BoundaryCondition.Periodic) // cubic spline with equal f' at borders
+let coefCubicNo    = Interpolation.interpolate(testDataX,testDataY,InterpolationMethod.CubicSpline Interpolation.CubicSpline.BoundaryCondition.NotAKnot) // cubic spline with continous f''' at second and penultimate knot
+let coefCubicPa    = Interpolation.interpolate(testDataX,testDataY,InterpolationMethod.CubicSpline Interpolation.CubicSpline.BoundaryCondition.Parabolic) // cubic spline with quadratic polynomial at borders
+let coefCubicCl    = Interpolation.interpolate(testDataX,testDataY,InterpolationMethod.CubicSpline (Interpolation.CubicSpline.BoundaryCondition.Clamped (0,-1))) // cubic spline with border f' set to 0 and -1
 let coefHermite    = Interpolation.interpolate(testDataX,testDataY,InterpolationMethod.HermiteSpline)
 let coefPolynomial = Interpolation.interpolate(testDataX,testDataY,InterpolationMethod.Polynomial)
 
@@ -86,6 +86,7 @@ let interpolationComparison =
         [1. .. 0.01 .. 10.] |> List.map (fun x -> x,Interpolation.predict(coefCubicPe)    x) |> Chart.Line |> Chart.withTraceInfo "Cubic_periodic"
         [1. .. 0.01 .. 10.] |> List.map (fun x -> x,Interpolation.predict(coefCubicNo)    x) |> Chart.Line |> Chart.withTraceInfo "Cubic_notaknot"
         [1. .. 0.01 .. 10.] |> List.map (fun x -> x,Interpolation.predict(coefCubicPa)    x) |> Chart.Line |> Chart.withTraceInfo "Cubic_parabolic"
+        [1. .. 0.01 .. 10.] |> List.map (fun x -> x,Interpolation.predict(coefCubicCl)    x) |> Chart.Line |> Chart.withTraceInfo "Cubic_clamped"
         [1. .. 0.01 .. 10.] |> List.map (fun x -> x,Interpolation.predict(coefHermite)    x) |> Chart.Line |> Chart.withTraceInfo "Hermite"
         [1. .. 0.01 .. 10.] |> List.map (fun x -> x,Interpolation.predict(coefPolynomial) x) |> Chart.Line |> Chart.withTraceInfo "Polynomial"
     ]

diff --git a/src/FSharp.Stats/Interpolation.fs b/src/FSharp.Stats/Interpolation.fs
@@ -650,7 +650,7 @@ module Interpolation =
             ///first and last polynomial are quadratic, not cubic
             | Quadratic
             ///f' at first and last knot are set by user
-            | Clamped
+            | Clamped of (float*float)
 
         /// <summary>
         /// Contains x data, y data (c0), slopes (c1), curvatures (c2), and the third derivative at each knot
@@ -912,15 +912,17 @@ module Interpolation =
 
                     [|(tmp1,0.);(tmp2,0.)|]
 
-                | Clamped -> //user defined border f''
-                    failwith "Not implemented yet. Slopes m1 and m2 have to be set by user"
-                    ////first condition: f''0(x0) = m1
-                    //tmp1.[0] <- 6. * firstX
-                    //tmp1.[1] <- 2.
-                    ////second condition: f''n-1(xn) = m2
-                    //tmp2.[4 * (intervalNumber - 1) + 0] <- 6. * lastX
-                    //tmp2.[4 * (intervalNumber - 1) + 1] <- 2.
-                    //[|(tmp1,m1);(tmp2,m2)|]
+                | Clamped (m1,m2)-> //user defined border f'
+                    //first condition: f'0(x0) = m1
+                    tmp1.[0] <- 3. * (pown firstX 2)
+                    tmp1.[1] <- 2. * firstX
+                    tmp1.[2] <- 1.
+
+                    //second condition: f'n-1(xn) = m2
+                    tmp2.[4 * (intervalNumber - 1) + 0] <- 3. * (pown lastX 2)
+                    tmp2.[4 * (intervalNumber - 1) + 1] <- 2. * lastX
+                    tmp2.[4 * (intervalNumber - 1) + 2] <- 1.
+                    [|(tmp1,m1);(tmp2,m2)|]
 
             let (equations,solutions) =
                 let interPolConstraints =