Merge pull request #375 from fslaborg/repo-assist/improve-interpolation-integrate-20260413-a4be8a75855ed951

[Repo Assist] feat: add `integrate` to LinearSpline, Step, and CubicSpline

Author: annacwe

src/FSharp.Stats/Interpolation.fs

@@ -567,6 +567,36 @@ module Interpolation =
             let k = leftSegmentIdx lsc.XValues x
             lsc.C1.[k]
 
+        /// <summary>
+        ///   Returns the definite integral of the linear spline from xVal1 to xVal2.
+        /// </summary>
+        /// <param name="lsc">Linear spline coefficients given as input x values, intersects, and slopes.</param>
+        /// <param name="xVal1">Lower bound of integration.</param>
+        /// <param name="xVal2">Upper bound of integration.</param>
+        /// <returns>Definite integral (signed area under the curve) from xVal1 to xVal2.</returns>
+        /// <remarks>xVal1 and xVal2 should lie within the range of the input x values; values outside the range are extrapolated using the nearest segment.</remarks>
+        let rec getIntegralBetween (lsc: LinearSplineCoef) xVal1 xVal2 =
+            if xVal1 > xVal2 then
+                - getIntegralBetween lsc xVal2 xVal1
+            elif 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
+
     /// <summary>
     ///   Module to create linear splines from x,y coordinates. x,y coordinates are interpolated by straight lines between two knots.
     /// </summary>
@@ -741,6 +771,38 @@ module Interpolation =
         let differentiate (lsc: StepCoef) x =
             0.
 
+        /// <summary>
+        ///   Returns the definite integral of the step function from xVal1 to xVal2.
+        /// </summary>
+        /// <param name="lsc">Step function coefficients given as input x values and intersects.</param>
+        /// <param name="xVal1">Lower bound of integration.</param>
+        /// <param name="xVal2">Upper bound of integration.</param>
+        /// <returns>Definite integral (signed area under the step function) from xVal1 to xVal2.</returns>
+        /// <remarks>xVal1 and xVal2 should lie within the range of the input x values.</remarks>
+        let rec getIntegralBetween (lsc: StepCoef) xVal1 xVal2 =
+            if xVal1 > xVal2 then
+                - getIntegralBetween lsc xVal2 xVal1
+            elif 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
+
 
     /// <summary>
     ///   Module to create piecewise cubic polynomials (cubic subsplines) from x,y coordinates. 
@@ -1741,6 +1803,51 @@ module Interpolation =
         let getThirdDerivative (coefficients: CubicSplineCoef) x =
             getDerivative 3 coefficients x
 
+        /// <summary>
+        ///   Returns the definite integral of the cubic spline from xVal1 to xVal2.
+        /// </summary>
+        /// <param name="coefficients">Interpolation functions coefficients.</param>
+        /// <param name="xVal1">Lower bound of integration.</param>
+        /// <param name="xVal2">Upper bound of integration.</param>
+        /// <returns>Definite integral (signed area under the curve) from xVal1 to xVal2.</returns>
+        /// <remarks>xVal1 and xVal2 should lie within the range of the input x values; values outside are handled by extrapolating the nearest segment's polynomial.</remarks>
+        let rec getIntegralBetween (coefficients: CubicSplineCoef) xVal1 xVal2 =
+            if xVal1 > xVal2 then
+                - getIntegralBetween coefficients xVal2 xVal1
+            elif 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
+
         /// <summary>
         /// Hermite cubic splines are defined by the function values and their slopes (first derivatives). If the slopws are unknown, they must be estimated.
         /// </summary>

tests/FSharp.Stats.Tests/Interpolation.fs

@@ -142,6 +142,166 @@ let BezierInterpolationTests =
 
 
 
+[<Tests>]
+let integrationTests =
+    testList "Interpolation.integrate" [
+
+        testCase "LinearSpline.getIntegralBetween positive 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.getIntegralBetween 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.getIntegralBetween coefs 1. 3.) 8.0
+                "integral of y=2x from 1 to 3 should be 8"
+            // ∫[0.5,1.5] 2x dx = 1.5^2 - 0.5^2 = 2.25 - 0.25 = 2 
+            Expect.floatClose Accuracy.high (LinearSpline.getIntegralBetween coefs 0.5 1.5) 2.0
+                "integral of y=2x from 0.5 to 1.5 should be 2"
+        
+        testCase "LinearSpline.getIntegralBetween reversed limits" <| fun () ->
+            let coefs = LinearSpline.interpolate [|0.;2.;4.|] [|0.;4.;8.|]
+            let forward = LinearSpline.getIntegralBetween coefs 1. 3.
+            let reverse = LinearSpline.getIntegralBetween coefs 3. 1.
+
+            Expect.floatClose Accuracy.high forward (-reverse)
+                "integrating with reversed limits should negate the result"
+
+        testCase "LinearSpline.getIntegralBetween is additive" <| fun () ->
+            let coefs = LinearSpline.interpolate [|0.;2.;4.|] [|0.;4.;8.|]
+            let whole = LinearSpline.getIntegralBetween coefs 1. 3.
+            let partials = LinearSpline.getIntegralBetween coefs 1. 2. + LinearSpline.getIntegralBetween coefs 2. 3.
+            Expect.floatClose Accuracy.high whole partials
+                "integral should be additive across knot boundaries"
+
+        testCase "LinearSpline.getIntegralBetween returns zero for equal bounds" <| fun () ->
+            let coefs = LinearSpline.interpolate [|0.;1.;2.|] [|1.;2.;3.|]
+            Expect.floatClose Accuracy.high (LinearSpline.getIntegralBetween coefs 1. 1.) 0.0
+                "integral with equal bounds should be 0"
+        
+        testCase "LinearSpline.getIntegralBetween negative linear function" <| fun () ->
+            let coefs = LinearSpline.interpolate [|0.;1.|] [|-1.;0.|]
+            // ∫[0,1] (x-1) dx = -1/2
+            Expect.floatClose Accuracy.high (LinearSpline.getIntegralBetween coefs 0. 1.) -0.5
+                "integral of negative linear function should be negative"
+        
+        testCase "LinearSpline.getIntegralBetween across zero crossing" <| fun () ->
+            let coefs = LinearSpline.interpolate [|0.;1.;2.|] [|-1.;0.;1.|]
+            // ∫[0,2] (x-1) dx = 0
+            Expect.floatClose Accuracy.high (LinearSpline.getIntegralBetween coefs 0. 2.) 0.0
+                "equivalent positive and negative areas should cancel"
+
+        testCase "Step.getIntegralBetween constant positive 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.getIntegralBetween 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.getIntegralBetween coefs 0.5 2.5) 6.0
+                "partial integral of step function from 0.5 to 2.5 should be 6"
+            // ∫[1,2] = 3 * 1 = 3
+            Expect.floatClose Accuracy.high (Step.getIntegralBetween coefs 1. 2.) 3.0
+                "integral within exact segment boundaries should be 3"
+            Expect.floatClose Accuracy.high (Step.getIntegralBetween coefs 1.5 1.5) 0.0
+                "integral over zero-width interval should be zero"
+            // Uses first interval value for values below x = 0
+            // ∫[-1,0] = 2 * 1 = 2
+            Expect.floatClose Accuracy.high (Step.getIntegralBetween coefs -1 0) 2.0
+                "integral before the first breakpoint should use first segment value"
+        
+        testCase "Step.getIntegralBetween is additive" <| fun () ->
+            let coefs = Step.interpolate [|0.;1.;2.;3.|] [|2.;3.;4.;5.|]
+            let whole = Step.getIntegralBetween coefs 0.5 2.5
+            let partials = (Step.getIntegralBetween coefs 0.5 1.5) + (Step.getIntegralBetween coefs 1.5 2.5)
+            Expect.floatClose Accuracy.high whole partials
+                "integral should be additive acreoss sub-intervals"
+        
+        testCase "Step.getIntegralBetween reversed limits" <| fun () ->
+            let coefs = Step.interpolate [|0.;1.;2.;3.|] [|2.;3.;4.;5.|]
+
+            let forward = Step.getIntegralBetween coefs 0.5 2.5
+            let reverse = Step.getIntegralBetween coefs 2.5 0.5
+
+            Expect.floatClose Accuracy.high reverse (-forward)
+                "integrating with reversed limits should negate the result"
+        
+        testCase "Step.getIntegralBetween constant negative segments" <| fun () ->
+            // y = -2 on [0,1), y = -3 on [1,2)
+            let coefs = Step.interpolate [|0.;1.;2.|] [|-2.;-3.;-3.|]
+            //∫[0,2] = (-2)*1 + (-3)*1 = -5
+            Expect.floatClose Accuracy.high (Step.getIntegralBetween coefs 0. 2.) -5.0
+                "Integral of negative step function should be negative"
+
+        testCase "Step.getIntegralBetween across zero crossing" <| fun () ->
+            // y= -2 on [0,1), y = 2 on [1,2)
+            let coefs = Step.interpolate [|0.;1.;2.|] [|-2.; 2.; 2.|]
+            // ∫[0,2] = (-2)*1 + 2*1 = 0
+            Expect.floatClose Accuracy.high (Step.getIntegralBetween coefs 0. 2.) 0.0
+                "equivalent positive and negative areas should cancel"
+        
+        testCase "CubicSpline.getIntegralBetween positive 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.getIntegralBetween 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.getIntegralBetween coefs 1. 2.) (7. / 3.)
+                "integral of y=x^2 from 1 to 2 should be 7/3"
+            // ∫[1.5,3.5] x² dx = (3.5³ - 1.5³)/3
+            Expect.floatClose Accuracy.high (CubicSpline.getIntegralBetween coefs 1.5 3.5) ((3.5 ** 3. - 1.5 ** 3.) / 3.)
+                "integral across several spline intervals should be correct"
+
+        testCase "CubicSpline.getIntegralBetween is additive" <| fun () ->
+            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
+            let whole= CubicSpline.getIntegralBetween coefs 1. 4.
+            let partials = CubicSpline.getIntegralBetween coefs 1. 2 + CubicSpline.getIntegralBetween coefs 2. 3. + CubicSpline.getIntegralBetween coefs 3. 4.
+            Expect.floatClose Accuracy.high whole partials
+                "integral should be additive across sub-intervals"
+
+        testCase "CubicSpline.getIntegralBetween 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.getIntegralBetween coefs 1. 1.) 0.0
+                "integral with equal bounds should be 0"
+        
+        testCase "CubicSpline.getIntegralBetween reversed limits" <| fun () ->
+            let t = vector [| 1.; 2.; 3.; 4. |]
+            let u = vector [| 1.; 4.; 9.; 16. |]
+            let coefs = CubicSpline.interpolate CubicSpline.Quadratic t u
+
+            let forward = CubicSpline.getIntegralBetween coefs 1. 4.
+            let reverse = CubicSpline.getIntegralBetween coefs 4. 1.
+
+            Expect.floatClose Accuracy.high forward (-reverse)
+                "integrating with reversed limits should negate the result"
+
+        testCase "CubicSpline.getIntegralBetween negative values" <| fun () ->
+            let t = vector [|1.;2.;3.;4.|]
+            let u = vector [|-1.;-4.;-9.;-16.|]
+            let coefs = CubicSpline.interpolate CubicSpline.Quadratic t u
+            // ∫[1,4] -x² dx = -21
+            Expect.floatClose Accuracy.high (CubicSpline.getIntegralBetween coefs 1. 4.) -21.0
+                "integral of negative quadratic should be negative"
+        
+        testCase "CubicSpline.getIntegralBetween positive and negative values" <| fun () ->
+            let t = vector [|-2.;-1.;0.;1.;2.|]
+            let u = vector [|-2.;-1.;0.;1.;2.|]
+            let coefs = CubicSpline.interpolate CubicSpline.Quadratic t u
+            // ∫[-2,2] x dx = 0
+            Expect.floatClose Accuracy.high (CubicSpline.getIntegralBetween coefs -2. 2.) 0.0
+                "positive and negative contributions should cancel"
+    ]
+
+
+
+