From 6d91d6cab805c68f8b46388e8e3a99cef068e0f8 Mon Sep 17 00:00:00 2001
From: "github-actions[bot]"
<41898282+github-actions[bot]@users.noreply.github.com>
Date: Mon, 13 Apr 2026 13:06:18 +0000
Subject: [PATCH 1/2] feat: add integrate function to LinearSpline, Step, and
CubicSpline
Implements definite integral computation for three piecewise interpolation
types, addressing part of issue #291.
- LinearSpline.integrate: exact integral of piecewise linear function
- Step.integrate: exact integral of step function
- CubicSpline.integrate: exact integral of piecewise cubic polynomial
All three use closed-form antiderivatives for numerical precision.
Adds 5 new tests covering known analytical results (y=2x, y=x^2, step).
1198/1198 tests pass.
Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
---
src/FSharp.Stats/Interpolation.fs | 101 ++++++++++++++++++++++
tests/FSharp.Stats.Tests/Interpolation.fs | 52 +++++++++++
2 files changed, 153 insertions(+)
diff --git a/src/FSharp.Stats/Interpolation.fs b/src/FSharp.Stats/Interpolation.fs
index 29f58812..9e258ca6 100644
--- a/src/FSharp.Stats/Interpolation.fs
+++ b/src/FSharp.Stats/Interpolation.fs
@@ -567,6 +567,34 @@ module Interpolation =
let k = leftSegmentIdx lsc.XValues x
lsc.C1.[k]
+ ///
+ /// Returns the definite integral of the linear spline from xVal1 to xVal2.
+ ///
+ /// Linear spline coefficients given as input x values, intersects, and slopes.
+ /// Lower bound of integration.
+ /// Upper bound of integration.
+ /// Definite integral (signed area under the curve) from xVal1 to xVal2.
+ /// xVal1 and xVal2 should lie within the range of the input x values; values outside the range are extrapolated using the nearest segment.
+ let integrate (lsc: LinearSplineCoef) xVal1 xVal2 =
+ if xVal1 = xVal2 then 0.
+ else
+ // Integral of segment k from x1 to x2:
+ // C0[k]*(x2-x1) + C1[k]*((x2-xk)^2 - (x1-xk)^2)/2
+ let segmentIntegral k x1 x2 =
+ let xk = lsc.XValues.[k]
+ lsc.C0.[k] * (x2 - x1) + lsc.C1.[k] * ((x2 - xk) * (x2 - xk) - (x1 - xk) * (x1 - xk)) / 2.
+ let k1 = leftSegmentIdx lsc.XValues xVal1
+ let k2 = leftSegmentIdx lsc.XValues xVal2
+ if k1 = k2 then
+ segmentIntegral k1 xVal1 xVal2
+ else
+ let firstPart = segmentIntegral k1 xVal1 lsc.XValues.[k1 + 1]
+ let lastPart = segmentIntegral k2 lsc.XValues.[k2] xVal2
+ let middleParts =
+ [ k1 + 1 .. k2 - 1 ]
+ |> List.sumBy (fun k -> segmentIntegral k lsc.XValues.[k] lsc.XValues.[k + 1])
+ firstPart + middleParts + lastPart
+
///
/// Module to create linear splines from x,y coordinates. x,y coordinates are interpolated by straight lines between two knots.
///
@@ -741,6 +769,36 @@ module Interpolation =
let differentiate (lsc: StepCoef) x =
0.
+ ///
+ /// Returns the definite integral of the step function from xVal1 to xVal2.
+ ///
+ /// Step function coefficients given as input x values and intersects.
+ /// Lower bound of integration.
+ /// Upper bound of integration.
+ /// Definite integral (signed area under the step function) from xVal1 to xVal2.
+ /// xVal1 and xVal2 should lie within the range of the input x values.
+ let integrate (lsc: StepCoef) xVal1 xVal2 =
+ if xVal1 = xVal2 then 0.
+ else
+ let n = lsc.XValues.Length
+ // Find interval index k such that XValues[k] <= x < XValues[k+1], clamped to [0, n-2]
+ let getInterval x =
+ if x >= lsc.XValues.[n - 1] then n - 2
+ elif x <= lsc.XValues.[0] then 0
+ else
+ lsc.XValues |> Array.findIndex (fun xk -> xk > x) |> fun idx -> idx - 1
+ let k1 = getInterval xVal1
+ let k2 = getInterval xVal2
+ if k1 = k2 then
+ lsc.C0.[k1] * (xVal2 - xVal1)
+ else
+ let firstPart = lsc.C0.[k1] * (lsc.XValues.[k1 + 1] - xVal1)
+ let lastPart = lsc.C0.[k2] * (xVal2 - lsc.XValues.[k2])
+ let middleParts =
+ [ k1 + 1 .. k2 - 1 ]
+ |> List.sumBy (fun k -> lsc.C0.[k] * (lsc.XValues.[k + 1] - lsc.XValues.[k]))
+ firstPart + middleParts + lastPart
+
///
/// Module to create piecewise cubic polynomials (cubic subsplines) from x,y coordinates.
@@ -1741,6 +1799,49 @@ module Interpolation =
let getThirdDerivative (coefficients: CubicSplineCoef) x =
getDerivative 3 coefficients x
+ ///
+ /// Returns the definite integral of the cubic spline from xVal1 to xVal2.
+ ///
+ /// Interpolation functions coefficients.
+ /// Lower bound of integration.
+ /// Upper bound of integration.
+ /// Definite integral (signed area under the curve) from xVal1 to xVal2.
+ /// xVal1 and xVal2 should lie within the range of the input x values; values outside are handled by extrapolating the nearest segment's polynomial.
+ let integrate (coefficients: CubicSplineCoef) xVal1 xVal2 =
+ if xVal1 = xVal2 then 0.
+ else
+ let sortedX = coefficients.XData |> Seq.sort |> Array.ofSeq
+ let n = sortedX.Length - 1 // number of intervals
+
+ // Find interval index k such that sortedX[k] <= x < sortedX[k+1], clamped to [0, n-1]
+ let getInterval x =
+ if x >= sortedX.[n] then n - 1
+ elif x < sortedX.[0] then 0
+ else
+ sortedX |> Array.findIndex (fun xk -> xk > x) |> fun idx -> idx - 1
+
+ // Antiderivative of the polynomial for interval k, evaluated at x:
+ // F_k(x) = a*x^4/4 + b*x^3/3 + c*x^2/2 + d*x
+ let antideriv k x =
+ let a = coefficients.C0_3.[4 * k + 0]
+ let b = coefficients.C0_3.[4 * k + 1]
+ let c = coefficients.C0_3.[4 * k + 2]
+ let d = coefficients.C0_3.[4 * k + 3]
+ a * x * x * x * x / 4. + b * x * x * x / 3. + c * x * x / 2. + d * x
+
+ let i1 = getInterval xVal1
+ let i2 = getInterval xVal2
+
+ if i1 = i2 then
+ antideriv i1 xVal2 - antideriv i1 xVal1
+ else
+ let firstPart = antideriv i1 sortedX.[i1 + 1] - antideriv i1 xVal1
+ let lastPart = antideriv i2 xVal2 - antideriv i2 sortedX.[i2]
+ let middleParts =
+ [ i1 + 1 .. i2 - 1 ]
+ |> List.sumBy (fun k -> antideriv k sortedX.[k + 1] - antideriv k sortedX.[k])
+ firstPart + middleParts + lastPart
+
///
/// Hermite cubic splines are defined by the function values and their slopes (first derivatives). If the slopws are unknown, they must be estimated.
///
diff --git a/tests/FSharp.Stats.Tests/Interpolation.fs b/tests/FSharp.Stats.Tests/Interpolation.fs
index 6a805481..dd9c3140 100644
--- a/tests/FSharp.Stats.Tests/Interpolation.fs
+++ b/tests/FSharp.Stats.Tests/Interpolation.fs
@@ -142,6 +142,58 @@ let BezierInterpolationTests =
+[]
+let integrationTests =
+ testList "Interpolation.integrate" [
+
+ testCase "LinearSpline.integrate linear function" <| fun () ->
+ // y = 2x at {0,2,4} => integral [0,4] = [x^2]_0^4 = 16
+ let coefs = LinearSpline.interpolate [|0.;2.;4.|] [|0.;4.;8.|]
+ Expect.floatClose Accuracy.high (LinearSpline.integrate coefs 0. 4.) 16.0
+ "integral of y=2x from 0 to 4 should be 16"
+ // [1,3]: segments [0,2] and [2,4], so partial cross-segment
+ // ∫[1,3] 2x dx = [x^2]_1^3 = 9 - 1 = 8
+ Expect.floatClose Accuracy.high (LinearSpline.integrate coefs 1. 3.) 8.0
+ "integral of y=2x from 1 to 3 should be 8"
+
+ testCase "LinearSpline.integrate returns zero for equal bounds" <| fun () ->
+ let coefs = LinearSpline.interpolate [|0.;1.;2.|] [|1.;2.;3.|]
+ Expect.floatClose Accuracy.high (LinearSpline.integrate coefs 1. 1.) 0.0
+ "integral with equal bounds should be 0"
+
+ testCase "Step.integrate constant segments" <| fun () ->
+ // y = 2 on [0,1), y = 3 on [1,2), y = 4 on [2,3)
+ let coefs = Step.interpolate [|0.;1.;2.;3.|] [|2.;3.;4.;5.|]
+ // ∫[0,3] = 2*1 + 3*1 + 4*1 = 9
+ Expect.floatClose Accuracy.high (Step.integrate coefs 0. 3.) 9.0
+ "integral of step function from 0 to 3 should be 9"
+ // ∫[0.5,2.5] = 2*0.5 + 3*1 + 4*0.5 = 1 + 3 + 2 = 6
+ Expect.floatClose Accuracy.high (Step.integrate coefs 0.5 2.5) 6.0
+ "partial integral of step function from 0.5 to 2.5 should be 6"
+
+ testCase "CubicSpline.integrate quadratic function" <| fun () ->
+ // Quadratic boundary condition reproduces y=x^2 exactly
+ let t = vector [| 1.; 2.; 3.; 4. |]
+ let u = vector [| 1.; 4.; 9.; 16. |] // y = x^2
+ let coefs = CubicSpline.interpolate CubicSpline.Quadratic t u
+ // ∫[1,4] x^2 dx = [x^3/3]_1^4 = 64/3 - 1/3 = 21
+ Expect.floatClose Accuracy.high (CubicSpline.integrate coefs 1. 4.) 21.0
+ "integral of y=x^2 from 1 to 4 should be 21"
+ // ∫[1,2] x^2 dx = 8/3 - 1/3 = 7/3
+ Expect.floatClose Accuracy.high (CubicSpline.integrate coefs 1. 2.) (7. / 3.)
+ "integral of y=x^2 from 1 to 2 should be 7/3"
+
+ testCase "CubicSpline.integrate returns zero for equal bounds" <| fun () ->
+ let t = vector [| 0.; 1.; 2. |]
+ let u = vector [| 0.; 1.; 4. |]
+ let coefs = CubicSpline.interpolate CubicSpline.Natural t u
+ Expect.floatClose Accuracy.high (CubicSpline.integrate coefs 1. 1.) 0.0
+ "integral with equal bounds should be 0"
+ ]
+
+
+
+
From d7fd78d687c6006bd6a44ff0a510e90164191cee Mon Sep 17 00:00:00 2001
From: "github-actions[bot]"
Date: Mon, 13 Apr 2026 13:06:22 +0000
Subject: [PATCH 2/2] ci: trigger checks