Bound IntegralRC correction rate to a settings-free physiological ceiling

Sources/LoopAlgorithm/RetrospectiveCorrection/IntegralRetrospectiveCorrection.swift

@@ -45,20 +45,36 @@ public class IntegralRetrospectiveCorrection: RetrospectiveCorrection {
     static let integralGain: Double = ((1 - integralForget) / integralForget) *
         (persistentDiscrepancyGain - currentDiscrepancyGain)
     static let proportionalGain: Double = currentDiscrepancyGain - integralGain
-
+
+    /// Default ceiling on the RC correction rate. Bounds how fast the integral RC may bend
+    /// the forecast, so a large or spurious sustained discrepancy can't wind the correction
+    /// up into over-/under-dosing. Chosen at the top of physiologically-plausible *sustained*
+    /// unmodeled glucose velocity (~4 mg/dL/min); unlike a bound scaled by ISF/basal/target,
+    /// it does not depend on the user's dosing settings — the plausible velocity of unmodeled
+    /// physiology is the same regardless of a person's insulin needs or target range.
+    public static let defaultMaxCorrectionVelocity = LoopQuantity(
+        unit: .milligramsPerDeciliterPerSecond, doubleValue: 4.0 / 60.0)
+
+    /// Ceiling applied to the correction rate (see `defaultMaxCorrectionVelocity`); nil disables it.
+    public let maxCorrectionVelocity: LoopQuantity?
+
     /// All math is performed with glucose expressed in mg/dL
     private let unit = LoopUnit.milligramsPerDeciliter
-    
+
     /// State variables reported in diagnostic issue report
     var recentDiscrepancyValues: [Double] = []
     var integralCorrectionEffectDuration: TimeInterval?
     var proportionalCorrection: Double = 0.0
     var integralCorrection: Double = 0.0
     var differentialCorrection: Double = 0.0
+    /// Correction rate actually used (after the `maxCorrectionVelocity` clamp), for diagnostics.
+    var correctionVelocity: LoopQuantity?
     var currentDate: Date = Date()
 
-    public init(effectDuration: TimeInterval) {
+    public init(effectDuration: TimeInterval,
+                maxCorrectionVelocity: LoopQuantity? = IntegralRetrospectiveCorrection.defaultMaxCorrectionVelocity) {
         self.effectDuration = effectDuration
+        self.maxCorrectionVelocity = maxCorrectionVelocity
     }
 
     /**
@@ -158,8 +174,18 @@ public class IntegralRetrospectiveCorrection: RetrospectiveCorrection {
 
         let retrospectionTimeInterval = currentDiscrepancy.endDate.timeIntervalSince(currentDiscrepancy.startDate)
         let discrepancyTime = max(retrospectionTimeInterval, retrospectiveCorrectionGroupingInterval)
-            let velocity = LoopQuantity(unit: .milligramsPerDeciliterPerSecond, doubleValue: scaledCorrection / discrepancyTime)
-
+        var velocityValue = scaledCorrection / discrepancyTime
+
+        // Bound the correction rate to a settings-free physiological ceiling (see
+        // `defaultMaxCorrectionVelocity`), so a large or spurious sustained discrepancy
+        // can't wind the integral up into over-/under-dosing.
+        if let maxCorrectionVelocity {
+            let cap = abs(maxCorrectionVelocity.doubleValue(for: .milligramsPerDeciliterPerSecond))
+            velocityValue = min(max(velocityValue, -cap), cap)
+        }
+        let velocity = LoopQuantity(unit: .milligramsPerDeciliterPerSecond, doubleValue: velocityValue)
+        correctionVelocity = velocity
+
         // Update array of glucose correction effects
         glucoseCorrectionEffect = startingGlucose.decayEffect(atRate: velocity, for: integralCorrectionEffectDuration!)
 

Tests/LoopAlgorithmTests/Mocks/IntegralRetrospectiveCorrectionTests.swift

@@ -45,4 +45,81 @@ final class IntegralRetrospectiveCorrectionTests: XCTestCase {
         XCTAssertEqual(effect.last?.quantity.doubleValue(for: .milligramsPerDeciliter) ?? 0, 110, accuracy: 0.05)
         XCTAssertEqual(effect.last?.startDate, dateFormatter.date(from: "2015-07-13T13:00:00")!)
     }
+
+    // MARK: - Correction-rate clamp (settings-free)
+    //
+    // A physiological ceiling on the RC correction rate replaces the deployed-LoopKit
+    // integral clamp, which scaled the bound by ISF×basal and the target range. See
+    // `defaultMaxCorrectionVelocity`.
+
+    private let capMgdlPerSec = 4.0 / 60.0
+
+    /// `count` contiguous 5-min discrepancies of `valuePerStep` mg/dL, ending at `endingAt`.
+    private func windup(endingAt: Date, count: Int, valuePerStep: Double) -> [GlucoseChange] {
+        (0..<count).reversed().map { i in
+            let end = endingAt.addingTimeInterval(.minutes(-Double(i) * 5))
+            return GlucoseChange(startDate: end.addingTimeInterval(.minutes(-5)),
+                                 endDate: end, quantity: .glucose(valuePerStep))
+        }
+    }
+
+    private func compute(_ irc: IntegralRetrospectiveCorrection, from: SimpleGlucoseValue,
+                         _ discrepancies: [GlucoseChange]) -> [GlucoseEffect] {
+        irc.computeEffect(
+            startingAt: from,
+            retrospectiveGlucoseDiscrepanciesSummed: discrepancies,
+            recencyInterval: TimeInterval(minutes: 15),
+            retrospectiveCorrectionGroupingInterval: LoopMath.retrospectiveCorrectionGroupingInterval)
+    }
+
+    func testCorrectionRateClampedToCeiling() {
+        let start = dateFormatter.date(from: "2015-07-13T12:00:00")!
+        let g = SimpleGlucoseValue(startDate: start, quantity: .glucose(150))
+        // Huge sustained under-prediction so the correction rate exceeds the ceiling.
+        let discrepancies = windup(endingAt: start, count: 18, valuePerStep: 250)
+
+        let unclamped = IntegralRetrospectiveCorrection(
+            effectDuration: LoopMath.retrospectiveCorrectionEffectDuration, maxCorrectionVelocity: nil)
+        let unclampedEffect = compute(unclamped, from: g, discrepancies)
+        XCTAssertGreaterThan(unclamped.correctionVelocity!.doubleValue(for: .milligramsPerDeciliterPerSecond),
+                             capMgdlPerSec, "windup must exceed the ceiling so the clamp is exercised")
+
+        // Default ceiling (4 mg/dL/min).
+        let clamped = IntegralRetrospectiveCorrection(effectDuration: LoopMath.retrospectiveCorrectionEffectDuration)
+        let clampedEffect = compute(clamped, from: g, discrepancies)
+        XCTAssertEqual(clamped.correctionVelocity!.doubleValue(for: .milligramsPerDeciliterPerSecond),
+                       capMgdlPerSec, accuracy: 1e-9)
+        // Clamping the rate reduces the forecast excursion.
+        XCTAssertLessThan(clampedEffect.last!.quantity.doubleValue(for: .milligramsPerDeciliter),
+                          unclampedEffect.last!.quantity.doubleValue(for: .milligramsPerDeciliter))
+    }
+
+    func testCorrectionRateClampIsSymmetricForNegativeDiscrepancies() {
+        let start = dateFormatter.date(from: "2015-07-13T12:00:00")!
+        let g = SimpleGlucoseValue(startDate: start, quantity: .glucose(90))
+        let discrepancies = windup(endingAt: start, count: 18, valuePerStep: -250)
+
+        let clamped = IntegralRetrospectiveCorrection(effectDuration: LoopMath.retrospectiveCorrectionEffectDuration)
+        _ = compute(clamped, from: g, discrepancies)
+        XCTAssertEqual(clamped.correctionVelocity!.doubleValue(for: .milligramsPerDeciliterPerSecond),
+                       -capMgdlPerSec, accuracy: 1e-9)
+    }
+
+    func testCorrectionRateUnclampedWhenBelowCeiling() {
+        let start = dateFormatter.date(from: "2015-07-13T12:00:00")!
+        let g = SimpleGlucoseValue(startDate: start, quantity: .glucose(120))
+        // Small windup: rate stays under the ceiling, so the clamp must not alter anything.
+        let discrepancies = windup(endingAt: start, count: 3, valuePerStep: 5)
+
+        let clamped = IntegralRetrospectiveCorrection(effectDuration: LoopMath.retrospectiveCorrectionEffectDuration)
+        let e1 = compute(clamped, from: g, discrepancies)
+        XCTAssertLessThan(abs(clamped.correctionVelocity!.doubleValue(for: .milligramsPerDeciliterPerSecond)),
+                          capMgdlPerSec)
+
+        let noClamp = IntegralRetrospectiveCorrection(
+            effectDuration: LoopMath.retrospectiveCorrectionEffectDuration, maxCorrectionVelocity: nil)
+        let e2 = compute(noClamp, from: g, discrepancies)
+        XCTAssertEqual(e1.last!.quantity.doubleValue(for: .milligramsPerDeciliter),
+                       e2.last!.quantity.doubleValue(for: .milligramsPerDeciliter), accuracy: 1e-9)
+    }
 }