diff --git a/Sources/LoopAlgorithm/AlgorithmInput.swift b/Sources/LoopAlgorithm/AlgorithmInput.swift index e2139ab..2bec7ec 100644 --- a/Sources/LoopAlgorithm/AlgorithmInput.swift +++ b/Sources/LoopAlgorithm/AlgorithmInput.swift @@ -34,4 +34,17 @@ public protocol AlgorithmInput { var recommendationType: DoseRecommendationType { get } var automaticBolusApplicationFactor: Double? { get } // Defaults to 0.4 var gradualTransitionsThreshold: Double? { get } + /// Optional algorithm experiment; nil is standard LoopAlgorithm. Conforming types that do not + /// provide it get nil. + var experiment: AlgorithmExperiment? { get } + /// The person's volatility baseline for the `stack` experiment: the median σ5 (mg/dL per 5 min) + /// over their last four weeks of glucose, from `GlucoseVolatility.baseline(_:)`. When nil it is + /// computed from `glucoseHistory` if that holds at least a day of readings; otherwise the stack runs + /// its post-low rise-cut only. + var volatilityBaseline: Double? { get } +} + +extension AlgorithmInput { + public var experiment: AlgorithmExperiment? { nil } + public var volatilityBaseline: Double? { nil } } diff --git a/Sources/LoopAlgorithm/AlgorithmInputFixture.swift b/Sources/LoopAlgorithm/AlgorithmInputFixture.swift index 4355778..4af45be 100644 --- a/Sources/LoopAlgorithm/AlgorithmInputFixture.swift +++ b/Sources/LoopAlgorithm/AlgorithmInputFixture.swift @@ -35,6 +35,9 @@ public struct AlgorithmInputFixture: AlgorithmInput { public var recommendationType: DoseRecommendationType = .automaticBolus public var automaticBolusApplicationFactor: Double? public var gradualTransitionsThreshold: Double? + /// Optional experiment; absent in existing fixtures, which therefore run standard LoopAlgorithm. + public var experiment: AlgorithmExperiment? + public var volatilityBaseline: Double? public var recommendationInsulinModel: InsulinModel { recommendationInsulinType.insulinModel @@ -148,6 +151,8 @@ extension AlgorithmInputFixture: Codable { self.automaticBolusApplicationFactor = try container.decodeIfPresent(Double.self, forKey: .automaticBolusApplicationFactor) self.gradualTransitionsThreshold = try container.decodeIfPresent(Double.self, forKey: .gradualTransitionsThreshold) ?? 40.0 + self.experiment = try container.decodeIfPresent(AlgorithmExperiment.self, forKey: .experiment) + self.volatilityBaseline = try container.decodeIfPresent(Double.self, forKey: .volatilityBaseline) } @@ -184,6 +189,8 @@ extension AlgorithmInputFixture: Codable { try container.encode(recommendationType.rawValue, forKey: .recommendationType) try container.encode(automaticBolusApplicationFactor, forKey: .automaticBolusApplicationFactor) try container.encode(gradualTransitionsThreshold, forKey: .gradualTransitionsThreshold) + try container.encodeIfPresent(experiment, forKey: .experiment) + try container.encodeIfPresent(volatilityBaseline, forKey: .volatilityBaseline) } private enum CodingKeys: String, CodingKey { @@ -206,6 +213,8 @@ extension AlgorithmInputFixture: Codable { case recommendationType case automaticBolusApplicationFactor case gradualTransitionsThreshold + case experiment + case volatilityBaseline } } diff --git a/Sources/LoopAlgorithm/Experiments/AlgorithmExperiments.swift b/Sources/LoopAlgorithm/Experiments/AlgorithmExperiments.swift new file mode 100644 index 0000000..1f4f714 --- /dev/null +++ b/Sources/LoopAlgorithm/Experiments/AlgorithmExperiments.swift @@ -0,0 +1,96 @@ +// +// AlgorithmExperiments.swift +// LoopAlgorithm +// +// Copyright © 2026 LoopKit Authors. All rights reserved. +// + +import Foundation + +/// An optional algorithm experiment. `nil` (the default everywhere) is standard LoopAlgorithm. +/// +/// There are deliberately only two, so that people who choose to try one pool into large groups whose +/// results can be compared. Each is a fixed configuration that was evaluated in LoopEval (closed-loop +/// counterfactual replay of real CGM and insulin data, scored against the standard algorithm swept over +/// insulin needs); neither has tunable parameters. +public enum AlgorithmExperiment: String, Codable, CaseIterable, Sendable { + + /// LoopEval C01. Integral retrospective correction (switched on by the experiment, whatever the + /// integral-RC setting), with the falling side of the correction weighted 1.5× and the rising side + /// 0.5×, and the correction rate capped at 4 mg/dL per minute. It reacts faster to glucose + /// running below the forecast and more slowly to glucose running above it. + case asymmetricIntegralRC + + /// LoopEval C30, three parts that act in different situations: + /// - **Post-low rise-cut** (C20). For 12 h after a reading below 70 mg/dL, and while glucose is + /// below 180, the rising side of retrospective correction is weighted 0.3×: the rebound after a + /// low (rescue carbs, counter-regulation) is not treated as a lasting need for insulin. + /// - **Volatility band** (C29). While no carbs are on board, the first two hours of the forecast are + /// lowered in proportion to recent glucose volatility, so the low guard sees it. Eventual glucose + /// is unchanged. + /// - **Calm-high licence** (C23/C25). When glucose is 180 or above, volatility is at or below the + /// person's own median and no carbs are on board, the automatic-bolus application factor is + /// doubled (capped at 1). Automatic-bolus dosing only. + /// + /// The band depth and the licence threshold are keyed to the person's median volatility; supply it + /// as `AlgorithmInput.volatilityBaseline` (see `GlucoseVolatility.baseline(_:)`). + case stack +} + +/// The internal configuration each experiment expands to. Not public: the experiments are fixed. +struct ExperimentParameters: Equatable { + var retrospectiveCorrectionDropScale = 1.0 + var retrospectiveCorrectionRiseScale = 1.0 + var integralRetrospectiveCorrectionMaxVelocity: Double? = nil // mg/dL per minute + var forceIntegralRetrospectiveCorrection = false + var postLowRiseCut: PostLowRiseCut? = nil + var volatilityBand: VolatilityBand? = nil + var calmHighLicence: CalmHighLicence? = nil + + static let none = ExperimentParameters() + + var needsVolatility: Bool { volatilityBand != nil || calmHighLicence != nil } + + /// `volatilityBaseline` is the person's median σ5 (mg/dL per 5 min); only the stack uses it. + init(experiment: AlgorithmExperiment?, volatilityBaseline: Double?) { + switch experiment { + case nil: + break + case .asymmetricIntegralRC: + retrospectiveCorrectionDropScale = 1.5 + retrospectiveCorrectionRiseScale = 0.5 + integralRetrospectiveCorrectionMaxVelocity = 4.0 + forceIntegralRetrospectiveCorrection = true + case .stack: + postLowRiseCut = PostLowRiseCut() + if let baseline = volatilityBaseline, baseline.isFinite, baseline > 0 { + volatilityBand = VolatilityBand(gain: Swift.min(1.0, 6.7 / baseline)) + calmHighLicence = CalmHighLicence(maxVolatility: baseline) + } + } + } + + init() {} + + struct PostLowRiseCut: Equatable { + var riseScale = 0.3 + var window: TimeInterval = .hours(12) + var lowThreshold = 70.0 // mg/dL + var releaseGlucose = 180.0 // mg/dL + } + + struct VolatilityBand: Equatable { + var gain: Double + var horizon: TimeInterval = .minutes(60) + var taperEnd: TimeInterval = .minutes(120) + var exponent = 0.71 + var carbGate = true + } + + struct CalmHighLicence: Equatable { + var applicationFactorScale = 2.0 + var minGlucose = 180.0 // mg/dL + var maxVolatility: Double // mg/dL per 5 min + var carbGate = true + } +} diff --git a/Sources/LoopAlgorithm/Experiments/ExperimentalForecastAdjustments.swift b/Sources/LoopAlgorithm/Experiments/ExperimentalForecastAdjustments.swift new file mode 100644 index 0000000..020b2a2 --- /dev/null +++ b/Sources/LoopAlgorithm/Experiments/ExperimentalForecastAdjustments.swift @@ -0,0 +1,54 @@ +// +// ExperimentalForecastAdjustments.swift +// LoopAlgorithm +// +// Copyright © 2026 LoopKit Authors. All rights reserved. +// + +import Foundation + +extension ExperimentParameters.PostLowRiseCut { + /// True when current glucose is below `releaseGlucose` and some reading within `window` before it + /// was below `lowThreshold`. + func isActive(glucoseHistory: [GlucoseType]) -> Bool { + let unit = LoopUnit.milligramsPerDeciliter + guard let latest = glucoseHistory.last, + latest.quantity.doubleValue(for: unit) < releaseGlucose else { return false } + for sample in glucoseHistory.reversed() { + if latest.startDate.timeIntervalSince(sample.startDate) > window { break } + if sample.quantity.doubleValue(for: unit) < lowThreshold { return true } + } + return false + } +} + +extension ExperimentParameters.CalmHighLicence { + func isActive(currentGlucose: Double, volatility: Double, activeCarbs: Double) -> Bool { + if carbGate && activeCarbs > 0 { return false } + return volatility.isFinite && currentGlucose >= minGlucose && volatility <= maxVolatility + } +} + +extension LoopAlgorithm { + + /// The volatility band: each predicted point τ minutes ahead is lowered by + /// `gain × σ5 × (τ/5)^exponent` up to `horizon`, tapering to zero by `taperEnd`. + static func volatilityBanded(_ glucose: [PredictedGlucoseValue], band: ExperimentParameters.VolatilityBand, + volatility: Double) -> [PredictedGlucoseValue] { + guard band.gain > 0, let t0 = glucose.first?.startDate else { return glucose } + let unit = LoopUnit.milligramsPerDeciliter + let horizon = band.horizon.minutes + let taperEnd = Swift.max(band.taperEnd.minutes, horizon + 1) + let atHorizon = pow(Swift.max(horizon, 5.0) / 5.0, band.exponent) + return glucose.map { point in + let tau = point.startDate.timeIntervalSince(t0).minutes + var shape = 0.0 + if tau <= 0 { shape = 0 } + else if tau <= horizon { shape = pow(Swift.max(tau, 5.0) / 5.0, band.exponent) } + else if tau < taperEnd { shape = atHorizon * (taperEnd - tau) / (taperEnd - horizon) } + let offset = -band.gain * volatility * shape + return PredictedGlucoseValue(startDate: point.startDate, + quantity: LoopQuantity(unit: unit, doubleValue: point.quantity.doubleValue(for: unit) + offset)) + } + } +} diff --git a/Sources/LoopAlgorithm/Experiments/GlucoseVolatility.swift b/Sources/LoopAlgorithm/Experiments/GlucoseVolatility.swift new file mode 100644 index 0000000..46de30e --- /dev/null +++ b/Sources/LoopAlgorithm/Experiments/GlucoseVolatility.swift @@ -0,0 +1,73 @@ +// +// GlucoseVolatility.swift +// LoopAlgorithm +// +// Copyright © 2026 LoopKit Authors. All rights reserved. +// + +import Foundation + +/// Causal short-term glucose volatility, used by the `stack` experiment. +public enum GlucoseVolatility { + + /// EWMA weight on the previous variance. + public static let lambda = 0.875 + /// CGM noise level; the variance is floored at 2 × noise². + public static let noise = 1.4 + + /// σ5 now: the exponentially weighted RMS of the 5-minute glucose increment, in mg/dL per 5 min, + /// over the last ~5 h of readings. + public static func sigma5(_ glucoseHistory: [GlucoseType]) -> Double { + let series = sigma5Series(fiveMinute(glucoseHistory).suffix(60)) + return series.last ?? (2.0 * noise * noise).squareRoot() + } + + /// The person's volatility baseline for the `stack` experiment: the median of σ5 over a glucose + /// history. LoopEval used each person's own median; four weeks of history gave a stable value + /// (about 9 % drift between successive four-week windows). Returns nil with less than a day of readings. + public static func baseline(_ glucoseHistory: [GlucoseType]) -> Double? { + let series = sigma5Series(fiveMinute(glucoseHistory)) + guard series.count >= 288 else { return nil } + let sorted = series.sorted() + let mid = sorted.count / 2 + return sorted.count % 2 == 0 ? (sorted[mid - 1] + sorted[mid]) / 2 : sorted[mid] + } + + /// Thin to 5-minute spacing: keep a reading only when it is at least 4 minutes after the last one + /// kept. With 5-minute CGM every reading is kept; with 1-minute CGM the estimator would otherwise + /// see only 1-minute increments, which it ignores, and never leave its floor. + static func fiveMinute(_ glucoseHistory: [GlucoseType]) -> [(date: Date, value: Double)] { + let unit = LoopUnit.milligramsPerDeciliter + var kept: [(date: Date, value: Double)] = [] + for sample in glucoseHistory { + if let last = kept.last, sample.startDate.timeIntervalSince(last.date) < .minutes(4) { continue } + kept.append((sample.startDate, sample.quantity.doubleValue(for: unit))) + } + return kept + } + + /// σ5 after each reading. Readings 2–7 min apart contribute (scaled to 5 min); a gap of more than + /// 7 min restarts the estimate, which reads as the floor until the next usable pair. + static func sigma5Series(_ samples: S) -> [Double] where S.Element == (date: Date, value: Double) { + let floorVariance = 2.0 * noise * noise + var variance = Double.nan + var previous: (date: Date, value: Double)? = nil + var out: [Double] = [] + for sample in samples { + if let previous { + let minutes = sample.date.timeIntervalSince(previous.date) / 60.0 + if minutes > 2.0 && minutes <= 7.0 { + let d5 = (sample.value - previous.value) * (5.0 / minutes) + variance = variance.isNaN + ? Swift.max(d5 * d5, floorVariance) + : Swift.max(lambda * variance + (1 - lambda) * d5 * d5, floorVariance) + } else if minutes > 7.0 { + variance = .nan + } + } + previous = sample + out.append(variance.isNaN ? floorVariance.squareRoot() : variance.squareRoot()) + } + return out + } +} diff --git a/Sources/LoopAlgorithm/LoopAlgorithm.swift b/Sources/LoopAlgorithm/LoopAlgorithm.swift index 2e711a1..3713995 100644 --- a/Sources/LoopAlgorithm/LoopAlgorithm.swift +++ b/Sources/LoopAlgorithm/LoopAlgorithm.swift @@ -179,7 +179,8 @@ public struct LoopAlgorithm { includingPositiveVelocityAndRC: Bool = true, useMidAbsorptionISF: Bool = false, carbAbsorptionModel: CarbAbsorptionComputable = PiecewiseLinearAbsorption(), - gradualTransitionsThreshold: Double? = 40.0 + gradualTransitionsThreshold: Double? = 40.0, + experiment: AlgorithmExperiment? = nil ) -> LoopPrediction where CarbType: CarbEntry, GlucoseType: GlucoseSampleValue, InsulinDoseType: InsulinDose { var prediction: [PredictedGlucoseValue] = [] @@ -254,10 +255,26 @@ public struct LoopAlgorithm { let rc: RetrospectiveCorrection - if useIntegralRetrospectiveCorrection { - rc = IntegralRetrospectiveCorrection(effectDuration: LoopMath.retrospectiveCorrectionEffectDuration) + // Experiments (no-ops when `experiment` is nil): asymmetric RC gains, and the post-low + // rise-cut, which replaces the rise gain while its gate is open. + let experimentParameters = ExperimentParameters(experiment: experiment, volatilityBaseline: nil) + let rcDropScale = experimentParameters.retrospectiveCorrectionDropScale + var rcRiseScale = experimentParameters.retrospectiveCorrectionRiseScale + if let cut = experimentParameters.postLowRiseCut, cut.isActive(glucoseHistory: glucoseHistory) { + rcRiseScale = cut.riseScale + } + + if useIntegralRetrospectiveCorrection || experimentParameters.forceIntegralRetrospectiveCorrection { + rc = IntegralRetrospectiveCorrection( + effectDuration: LoopMath.retrospectiveCorrectionEffectDuration, + dropScale: rcDropScale, + riseScale: rcRiseScale, + maxVelocity: experimentParameters.integralRetrospectiveCorrectionMaxVelocity) } else { - rc = StandardRetrospectiveCorrection(effectDuration: LoopMath.retrospectiveCorrectionEffectDuration) + rc = StandardRetrospectiveCorrection( + effectDuration: LoopMath.retrospectiveCorrectionEffectDuration, + dropScale: rcDropScale, + riseScale: rcRiseScale) } if let latestGlucose = glucoseHistory.last { @@ -550,9 +567,22 @@ public struct LoopAlgorithm { includingPositiveVelocityAndRC: input.includePositiveVelocityAndRC, useMidAbsorptionISF: input.useMidAbsorptionISF, carbAbsorptionModel: input.carbAbsorptionModel.model, - gradualTransitionsThreshold: input.gradualTransitionsThreshold + gradualTransitionsThreshold: input.gradualTransitionsThreshold, + experiment: input.experiment ) + // Experiments that adjust the forecast before dosing (no-ops when `experiment` is nil). + let experimentParameters = ExperimentParameters( + experiment: input.experiment, + volatilityBaseline: input.experiment == .stack + ? (input.volatilityBaseline ?? GlucoseVolatility.baseline(input.glucoseHistory)) : nil) + let volatility: Double? = experimentParameters.needsVolatility ? GlucoseVolatility.sigma5(input.glucoseHistory) : nil + if let band = experimentParameters.volatilityBand, let volatility, + !(band.carbGate && (prediction.activeCarbs ?? 0) > 0) + { + prediction.glucose = volatilityBanded(prediction.glucose, band: band, volatility: volatility) + } + let sensitivityForDosing: [AbsoluteScheduleValue] if input.useMidAbsorptionISF { sensitivityForDosing = input.sensitivity @@ -587,9 +617,16 @@ public struct LoopAlgorithm { target: input.target) result = .success(.init(manual: recommendation)) case .automaticBolus: + var applicationFactor = input.automaticBolusApplicationFactor ?? defaultBolusPartialApplicationFactor + if let licence = experimentParameters.calmHighLicence, let volatility, + licence.isActive(currentGlucose: latestGlucose.quantity.doubleValue(for: .milligramsPerDeciliter), + volatility: volatility, activeCarbs: prediction.activeCarbs ?? 0) + { + applicationFactor = Swift.min(1.0, applicationFactor * licence.applicationFactorScale) + } let recommendation = recommendAutomaticDose( for: correction, - applicationFactor: input.automaticBolusApplicationFactor ?? defaultBolusPartialApplicationFactor, + applicationFactor: applicationFactor, neutralBasalRate: scheduledBasalRate, activeInsulin: prediction.activeInsulin!, maxBolus: input.maxBolus, diff --git a/Sources/LoopAlgorithm/RetrospectiveCorrection/IntegralRetrospectiveCorrection.swift b/Sources/LoopAlgorithm/RetrospectiveCorrection/IntegralRetrospectiveCorrection.swift index d850b9d..d8b3a79 100644 --- a/Sources/LoopAlgorithm/RetrospectiveCorrection/IntegralRetrospectiveCorrection.swift +++ b/Sources/LoopAlgorithm/RetrospectiveCorrection/IntegralRetrospectiveCorrection.swift @@ -57,8 +57,18 @@ public class IntegralRetrospectiveCorrection: RetrospectiveCorrection { var differentialCorrection: Double = 0.0 var currentDate: Date = Date() - public init(effectDuration: TimeInterval) { + /// Experimental asymmetric gains (see `AlgorithmExperiments`): multiply the correction when the + /// current discrepancy run is negative (drop) or positive (rise). 1.0 / 1.0 is standard IRC. + let dropScale: Double + let riseScale: Double + /// Experimental ceiling on the correction rate, mg/dL per minute. nil = none (standard IRC). + let maxVelocity: Double? + + public init(effectDuration: TimeInterval, dropScale: Double = 1.0, riseScale: Double = 1.0, maxVelocity: Double? = nil) { self.effectDuration = effectDuration + self.dropScale = dropScale + self.riseScale = riseScale + self.maxVelocity = maxVelocity } /** @@ -158,7 +168,15 @@ public class IntegralRetrospectiveCorrection: RetrospectiveCorrection { let retrospectionTimeInterval = currentDiscrepancy.endDate.timeIntervalSince(currentDiscrepancy.startDate) let discrepancyTime = max(retrospectionTimeInterval, retrospectiveCorrectionGroupingInterval) - let velocity = LoopQuantity(unit: .milligramsPerDeciliterPerSecond, doubleValue: scaledCorrection / discrepancyTime) + if dropScale != 1.0 || riseScale != 1.0 { + scaledCorrection *= currentDiscrepancyValue < 0 ? dropScale : riseScale + } + var velocityValue = scaledCorrection / discrepancyTime + if let maxVelocity { + let cap = abs(maxVelocity) / 60.0 + velocityValue = Swift.min(Swift.max(velocityValue, -cap), cap) + } + let velocity = LoopQuantity(unit: .milligramsPerDeciliterPerSecond, doubleValue: velocityValue) // Update array of glucose correction effects glucoseCorrectionEffect = startingGlucose.decayEffect(atRate: velocity, for: integralCorrectionEffectDuration!) diff --git a/Sources/LoopAlgorithm/RetrospectiveCorrection/StandardRetrospectiveCorrection.swift b/Sources/LoopAlgorithm/RetrospectiveCorrection/StandardRetrospectiveCorrection.swift index 0ee2aff..f5aa5fd 100644 --- a/Sources/LoopAlgorithm/RetrospectiveCorrection/StandardRetrospectiveCorrection.swift +++ b/Sources/LoopAlgorithm/RetrospectiveCorrection/StandardRetrospectiveCorrection.swift @@ -23,8 +23,15 @@ public class StandardRetrospectiveCorrection: RetrospectiveCorrection { /// All math is performed with glucose expressed in mg/dL private let unit = LoopUnit.milligramsPerDeciliter - public init(effectDuration: TimeInterval) { + /// Experimental asymmetric gains (see `AlgorithmExperiments`): multiply the correction when the + /// current discrepancy is negative (drop) or positive (rise). 1.0 / 1.0 is standard RC. + let dropScale: Double + let riseScale: Double + + public init(effectDuration: TimeInterval, dropScale: Double = 1.0, riseScale: Double = 1.0) { self.effectDuration = effectDuration + self.dropScale = dropScale + self.riseScale = riseScale } public func computeEffect( @@ -43,7 +50,10 @@ public class StandardRetrospectiveCorrection: RetrospectiveCorrection { } // Standard retrospective correction math - let currentDiscrepancyValue = currentDiscrepancy.quantity.doubleValue(for: unit) + var currentDiscrepancyValue = currentDiscrepancy.quantity.doubleValue(for: unit) + if dropScale != 1.0 || riseScale != 1.0 { + currentDiscrepancyValue *= currentDiscrepancyValue < 0 ? dropScale : riseScale + } totalGlucoseCorrectionEffect = LoopQuantity(unit: unit, doubleValue: currentDiscrepancyValue) let retrospectionTimeInterval = currentDiscrepancy.endDate.timeIntervalSince(currentDiscrepancy.startDate) diff --git a/Tests/LoopAlgorithmTests/AlgorithmExperimentsTests.swift b/Tests/LoopAlgorithmTests/AlgorithmExperimentsTests.swift new file mode 100644 index 0000000..cb9ebed --- /dev/null +++ b/Tests/LoopAlgorithmTests/AlgorithmExperimentsTests.swift @@ -0,0 +1,225 @@ +// +// AlgorithmExperimentsTests.swift +// LoopAlgorithmTests +// +// Copyright © 2026 LoopKit Authors. All rights reserved. +// + +import XCTest +@testable import LoopAlgorithm + +final class AlgorithmExperimentsTests: XCTestCase { + + private let unit = LoopUnit.milligramsPerDeciliter + + private func loadScenario(_ name: String) -> AlgorithmInputFixture { + let decoder = JSONDecoder() + decoder.dateDecodingStrategy = .iso8601 + let url = Bundle.module.url(forResource: name + "_scenario", withExtension: "json", subdirectory: "Fixtures")! + return try! decoder.decode(AlgorithmInputFixture.self, from: try! Data(contentsOf: url)) + } + + private func glucoseSeries(_ values: [Double], endingAt end: Date, spacing: TimeInterval = .minutes(5)) -> [FixtureGlucoseSample] { + values.enumerated().map { i, v in + FixtureGlucoseSample(startDate: end.addingTimeInterval(-Double(values.count - 1 - i) * spacing), quantity: .glucose(v)) + } + } + + private func assertSameOutput(_ a: AlgorithmOutput, _ b: AlgorithmOutput, + file: StaticString = #filePath, line: UInt = #line) { + XCTAssertEqual(a.recommendation, b.recommendation, file: file, line: line) + XCTAssertEqual(a.predictedGlucose.map { $0.quantity.doubleValue(for: unit) }, + b.predictedGlucose.map { $0.quantity.doubleValue(for: unit) }, file: file, line: line) + } + + // MARK: - Default is standard LoopAlgorithm + + func testExistingFixturesDecodeWithoutExperiments() { + XCTAssertNil(loadScenario("suspend").experiment) + XCTAssertNil(loadScenario("carbs_with_isf_change").experiment) + } + + func testGatedExperimentsAreInertWhenTheirGateIsClosed() { + // Mock: glucose rising 100 -> 160, no low, below 180. With carbs on board the band is gated off, + // there is no low for the post-low cut, and the licence needs >= 180: the stack must change nothing. + let now = ISO8601DateFormatter().date(from: "2024-01-03T12:00:00+0000")! + var input = AlgorithmInputFixture.mock(for: now) + input.recommendationType = .automaticBolus + input.carbEntries = [FixtureCarbEntry(absorptionTime: .hours(3), startDate: now.addingTimeInterval(.minutes(-30)), quantity: .carbs(30))] + let standard = LoopAlgorithm.run(input: input) + input.experiment = .stack + input.volatilityBaseline = 5 + assertSameOutput(standard, LoopAlgorithm.run(input: input)) + } + + func testStackWithoutABaselineRunsOnlyThePostLowCut() { + // No baseline and under a day of glucose: band and licence are off, so with no low nothing changes. + let now = ISO8601DateFormatter().date(from: "2024-01-03T12:00:00+0000")! + var input = AlgorithmInputFixture.mock(for: now) + let standard = LoopAlgorithm.run(input: input) + input.experiment = .stack + assertSameOutput(standard, LoopAlgorithm.run(input: input)) + XCTAssertNil(ExperimentParameters(experiment: .stack, volatilityBaseline: nil).volatilityBand) + XCTAssertNotNil(ExperimentParameters(experiment: .stack, volatilityBaseline: nil).postLowRiseCut) + } + + func testAsymmetricIntegralRCTurnsIntegralRCOn() { + var input = loadScenario("carbs_with_isf_change") + input.experiment = .asymmetricIntegralRC + input.useIntegralRetrospectiveCorrection = false + let forced = LoopAlgorithm.run(input: input) + input.useIntegralRetrospectiveCorrection = true + assertSameOutput(forced, LoopAlgorithm.run(input: input)) + let p = ExperimentParameters(experiment: .asymmetricIntegralRC, volatilityBaseline: nil) + XCTAssertEqual(p.retrospectiveCorrectionDropScale, 1.5) + XCTAssertEqual(p.retrospectiveCorrectionRiseScale, 0.5) + XCTAssertEqual(p.integralRetrospectiveCorrectionMaxVelocity, 4.0) + XCTAssertNil(p.postLowRiseCut) + } + + func testStackParameters() { + let p = ExperimentParameters(experiment: .stack, volatilityBaseline: 5) + XCTAssertEqual(p.volatilityBand?.gain, 1.0) // min(1, 6.7 / 5) + XCTAssertEqual(ExperimentParameters(experiment: .stack, volatilityBaseline: 10).volatilityBand!.gain, 0.67, accuracy: 1e-12) + XCTAssertEqual(p.calmHighLicence?.maxVolatility, 5) + XCTAssertEqual(p.postLowRiseCut, ExperimentParameters.PostLowRiseCut()) + XCTAssertEqual(ExperimentParameters(experiment: nil, volatilityBaseline: 5), ExperimentParameters.none) + } + + // MARK: - Retrospective correction gains + + func testStandardRCRiseAndDropScales() { + let now = Date() + let latest = SimpleGlucoseValue(startDate: now, quantity: .glucose(150)) + func total(_ discrepancy: Double, drop: Double, rise: Double) -> Double { + let rc = StandardRetrospectiveCorrection(effectDuration: .minutes(60), dropScale: drop, riseScale: rise) + let change = GlucoseChange(startDate: now.addingTimeInterval(.minutes(-30)), endDate: now, quantity: .glucose(discrepancy)) + _ = rc.computeEffect(startingAt: latest, retrospectiveGlucoseDiscrepanciesSummed: [change], + recencyInterval: .minutes(15), retrospectiveCorrectionGroupingInterval: .minutes(30)) + return rc.totalGlucoseCorrectionEffect!.doubleValue(for: unit) + } + XCTAssertEqual(total(10, drop: 1.5, rise: 0.5), 5, accuracy: 1e-9) + XCTAssertEqual(total(-10, drop: 1.5, rise: 0.5), -15, accuracy: 1e-9) + XCTAssertEqual(total(10, drop: 1, rise: 1), 10, accuracy: 1e-9) + } + + func testIntegralRCVelocityCap() { + let now = Date() + let latest = SimpleGlucoseValue(startDate: now, quantity: .glucose(250)) + // A 30-min discrepancy of +200 mg/dL is a correction rate far above 4 mg/dL/min. + let change = GlucoseChange(startDate: now.addingTimeInterval(.minutes(-30)), endDate: now, quantity: .glucose(200)) + let capped = IntegralRetrospectiveCorrection(effectDuration: .minutes(60), maxVelocity: 4.0) + let effect = capped.computeEffect(startingAt: latest, retrospectiveGlucoseDiscrepanciesSummed: [change], + recencyInterval: .minutes(15), retrospectiveCorrectionGroupingInterval: .minutes(30)) + let uncapped = IntegralRetrospectiveCorrection(effectDuration: .minutes(60)) + let effectU = uncapped.computeEffect(startingAt: latest, retrospectiveGlucoseDiscrepanciesSummed: [change], + recencyInterval: .minutes(15), retrospectiveCorrectionGroupingInterval: .minutes(30)) + let rise: ([GlucoseEffect]) -> Double = { $0.last!.quantity.doubleValue(for: self.unit) - $0.first!.quantity.doubleValue(for: self.unit) } + XCTAssertLessThan(rise(effect), rise(effectU)) + // First 5-min step can rise by at most 4 mg/dL/min * 5 min. + XCTAssertLessThanOrEqual(effect[1].quantity.doubleValue(for: unit) - effect[0].quantity.doubleValue(for: unit), 20.0 + 1e-6) + } + + // MARK: - Post-low rise-cut gate + + func testPostLowRiseCutGate() { + let now = Date() + let cut = ExperimentParameters.PostLowRiseCut() + XCTAssertTrue(cut.isActive(glucoseHistory: glucoseSeries([90, 65, 80, 120], endingAt: now))) + XCTAssertFalse(cut.isActive(glucoseHistory: glucoseSeries([90, 65, 120, 190], endingAt: now)), "released above 180") + XCTAssertFalse(cut.isActive(glucoseHistory: glucoseSeries([90, 100, 110, 120], endingAt: now)), "no low") + let old = [FixtureGlucoseSample(startDate: now.addingTimeInterval(.hours(-13)), quantity: .glucose(60))] + + glucoseSeries([100, 110, 120], endingAt: now) + XCTAssertFalse(cut.isActive(glucoseHistory: old), "low outside the 12 h window") + } + + func testPostLowRiseCutLowersTheForecastAfterALow() { + // Rebound after a low, with insulin on board: the unexplained rise is cut, so the forecast is lower. + let now = ISO8601DateFormatter().date(from: "2024-01-03T12:00:00+0000")! + var input = AlgorithmInputFixture.mock(for: now) + input.glucoseHistory = glucoseSeries([62, 66, 75, 90, 110, 130], endingAt: now.addingTimeInterval(.minutes(-1))) + let standard = LoopAlgorithm.run(input: input) + input.experiment = .stack + let cut = LoopAlgorithm.run(input: input) + XCTAssertLessThan(cut.predictedGlucose.last!.quantity.doubleValue(for: unit), + standard.predictedGlucose.last!.quantity.doubleValue(for: unit)) + } + + // MARK: - Volatility + + func testSigma5() { + let now = Date() + XCTAssertEqual(GlucoseVolatility.sigma5(glucoseSeries([120, 120, 120, 120], endingAt: now)), (2 * 1.4 * 1.4).squareRoot(), accuracy: 1e-9) + // Constant 10 mg/dL steps: the EWMA of 100 is 100. + XCTAssertEqual(GlucoseVolatility.sigma5(glucoseSeries([100, 110, 120, 130, 140], endingAt: now)), 10, accuracy: 1e-9) + // A gap > 7 min restarts the estimate: only the steps after the gap count. + let gapped = glucoseSeries([100, 140, 180], endingAt: now.addingTimeInterval(.minutes(-30))) + + glucoseSeries([100, 102], endingAt: now) + XCTAssertEqual(GlucoseVolatility.sigma5(gapped), 2.0, accuracy: 1e-9) // one 2 mg/dL step; 4 > floor 3.92 + } + + func testSigma5OnOneMinuteCGM() { + // The same 10 mg/dL per 5 min trend sampled every minute must read like 5-minute data, not the floor. + let now = Date() + let oneMinute = (0...40).map { FixtureGlucoseSample(startDate: now.addingTimeInterval(.minutes(Double($0 - 40))), quantity: .glucose(100 + 2 * Double($0))) } + XCTAssertEqual(GlucoseVolatility.sigma5(oneMinute), 10, accuracy: 1e-9) + } + + func testBaselineIsTheMedianOfSigma5() { + let now = Date() + XCTAssertNil(GlucoseVolatility.baseline(glucoseSeries(Array(repeating: 120, count: 100), endingAt: now)), "under a day") + // A calm stretch then a longer +/-10 stretch: the median is volatile, above the floor and at most 10. + let flat = Array(repeating: 120.0, count: 100) + let zig = (0..<200).map { $0 % 2 == 0 ? 120.0 : 130.0 } + let b = GlucoseVolatility.baseline(glucoseSeries(flat + zig, endingAt: now))! + XCTAssertGreaterThan(b, (2 * 1.4 * 1.4).squareRoot()) + XCTAssertLessThanOrEqual(b, 10) + } + + func testVolatilityBandShape() { + let now = Date() + let glucose = (0...48).map { PredictedGlucoseValue(startDate: now.addingTimeInterval(.minutes(5 * Double($0))), quantity: .glucose(120)) } + let band = ExperimentParameters.VolatilityBand(gain: 1.0) + let banded = LoopAlgorithm.volatilityBanded(glucose, band: band, volatility: 4.0) + let offset: (Int) -> Double = { banded[$0].quantity.doubleValue(for: self.unit) - 120 } + XCTAssertEqual(offset(0), 0, accuracy: 1e-9) + XCTAssertEqual(offset(12), -4.0 * pow(12.0, 0.71), accuracy: 1e-9) // 60 min: full depth + XCTAssertEqual(offset(18), -4.0 * pow(12.0, 0.71) / 2, accuracy: 1e-9) // 90 min: half way down the taper + XCTAssertEqual(offset(24), 0, accuracy: 1e-9) // 120 min: gone + XCTAssertEqual(offset(48), 0, accuracy: 1e-9) // eventual untouched + } + + // MARK: - Calm-high licence + + func testCalmHighLicenceDoublesTheApplicationFactorOnACalmHigh() { + let now = ISO8601DateFormatter().date(from: "2024-01-03T12:00:00+0000")! + var input = AlgorithmInputFixture.mock(for: now) + input.recommendationType = .automaticBolus + input.glucoseHistory = glucoseSeries([220, 220, 221, 220], endingAt: now.addingTimeInterval(.minutes(-1))) + let standard = LoopAlgorithm.run(input: input).recommendation!.automatic!.bolusUnits! + input.experiment = .stack + input.volatilityBaseline = 5 + let licensed = LoopAlgorithm.run(input: input).recommendation!.automatic!.bolusUnits! + XCTAssertGreaterThan(standard, 0) + // AF 0.4 -> 0.8; the band (on, no carbs) lowers only the near-term forecast, so the bolus roughly doubles. + XCTAssertGreaterThan(licensed, 1.5 * standard) + XCTAssertLessThanOrEqual(licensed, 2 * standard + 0.051) + // Volatile relative to the person's baseline: no licence. (The band is still on, and only lowers the dose.) + input.volatilityBaseline = 0.5 + XCTAssertLessThanOrEqual(LoopAlgorithm.run(input: input).recommendation!.automatic!.bolusUnits!, standard + 1e-9) + } + + // MARK: - Codable + + func testFixtureRoundTripCarriesTheExperiment() throws { + var input = loadScenario("suspend") + input.experiment = .stack + input.volatilityBaseline = 4.2 + let encoder = JSONEncoder(); encoder.dateEncodingStrategy = .iso8601 + let decoder = JSONDecoder(); decoder.dateDecodingStrategy = .iso8601 + let back = try decoder.decode(AlgorithmInputFixture.self, from: encoder.encode(input)) + XCTAssertEqual(back.experiment, .stack) + XCTAssertEqual(back.volatilityBaseline, 4.2) + XCTAssertEqual(AlgorithmExperiment.allCases, [.asymmetricIntegralRC, .stack]) + } +}