From d585d157927eb19bd34daf0d47e8e5aae6f670b2 Mon Sep 17 00:00:00 2001 From: Pete Schwamb Date: Thu, 15 Oct 2026 12:35:33 -0500 Subject: [PATCH 1/2] Optional algorithm experiments, off by default Adds AlgorithmExperiments, an options struct carried on AlgorithmInput (with a protocol-extension default of .none, so existing conformers are unchanged) and passed to generatePrediction. With .none every experiment code path is skipped and the forecast and dose are standard LoopAlgorithm. The experiments are the leading candidates from LoopEval's closed-loop counterfactual evaluation: - retrospectiveCorrectionDropScale / RiseScale: asymmetric RC gains, standard and integral (C01, asymmetric IRC, with integral RC on) - integralRetrospectiveCorrectionMaxVelocity: optional cap on the IRC correction rate (part of C01 as scored) - postLowRiseCut: within 12 h of a reading below 70 and while glucose is below 180, scale the rising side of RC (C20) - volatilityBand: lower the near-term forecast by recent glucose volatility, tapering to zero by 2 h, so the low guard sees it (C29) - calmHighLicence: a larger automatic-bolus application factor on a calm, carb-free high (C23/C25) - negativeInsulinDamper: the Loop and Learn customization (C37) AlgorithmExperiments.stack(medianVolatility:) builds the LoopEval stack (C30) and .asymmetricIntegralRC the C01 configuration. GlucoseVolatility.sigma5 is the causal volatility estimate the band and the licence key on. Tests: existing fixtures decode with no experiments; explicit .none, unity gains and experiments whose gates are closed reproduce the standard output; one unit test per mechanism. --- Sources/LoopAlgorithm/AlgorithmInput.swift | 7 + .../LoopAlgorithm/AlgorithmInputFixture.swift | 6 + .../Experiments/AlgorithmExperiments.swift | 173 +++++++++++++++ .../ExperimentalForecastAdjustments.swift | 139 ++++++++++++ .../Experiments/GlucoseVolatility.swift | 43 ++++ Sources/LoopAlgorithm/LoopAlgorithm.swift | 46 +++- .../IntegralRetrospectiveCorrection.swift | 22 +- .../StandardRetrospectiveCorrection.swift | 14 +- .../AlgorithmExperimentsTests.swift | 205 ++++++++++++++++++ 9 files changed, 646 insertions(+), 9 deletions(-) create mode 100644 Sources/LoopAlgorithm/Experiments/AlgorithmExperiments.swift create mode 100644 Sources/LoopAlgorithm/Experiments/ExperimentalForecastAdjustments.swift create mode 100644 Sources/LoopAlgorithm/Experiments/GlucoseVolatility.swift create mode 100644 Tests/LoopAlgorithmTests/AlgorithmExperimentsTests.swift diff --git a/Sources/LoopAlgorithm/AlgorithmInput.swift b/Sources/LoopAlgorithm/AlgorithmInput.swift index e2139ab..7e7e884 100644 --- a/Sources/LoopAlgorithm/AlgorithmInput.swift +++ b/Sources/LoopAlgorithm/AlgorithmInput.swift @@ -34,4 +34,11 @@ public protocol AlgorithmInput { var recommendationType: DoseRecommendationType { get } var automaticBolusApplicationFactor: Double? { get } // Defaults to 0.4 var gradualTransitionsThreshold: Double? { get } + /// Optional algorithm experiments. Defaults to `.none` (standard LoopAlgorithm) for any conforming + /// type that does not provide it. + var experiments: AlgorithmExperiments { get } +} + +extension AlgorithmInput { + public var experiments: AlgorithmExperiments { .none } } diff --git a/Sources/LoopAlgorithm/AlgorithmInputFixture.swift b/Sources/LoopAlgorithm/AlgorithmInputFixture.swift index 4355778..98d9471 100644 --- a/Sources/LoopAlgorithm/AlgorithmInputFixture.swift +++ b/Sources/LoopAlgorithm/AlgorithmInputFixture.swift @@ -35,6 +35,12 @@ public struct AlgorithmInputFixture: AlgorithmInput { public var recommendationType: DoseRecommendationType = .automaticBolus public var automaticBolusApplicationFactor: Double? public var gradualTransitionsThreshold: Double? + /// Optional experiments; absent in existing fixtures, which therefore run standard LoopAlgorithm. + public var algorithmExperiments: AlgorithmExperiments? + + public var experiments: AlgorithmExperiments { + algorithmExperiments ?? .none + } public var recommendationInsulinModel: InsulinModel { recommendationInsulinType.insulinModel diff --git a/Sources/LoopAlgorithm/Experiments/AlgorithmExperiments.swift b/Sources/LoopAlgorithm/Experiments/AlgorithmExperiments.swift new file mode 100644 index 0000000..7005f27 --- /dev/null +++ b/Sources/LoopAlgorithm/Experiments/AlgorithmExperiments.swift @@ -0,0 +1,173 @@ +// +// AlgorithmExperiments.swift +// LoopAlgorithm +// +// Copyright © 2026 LoopKit Authors. All rights reserved. +// + +import Foundation + +/// Optional, off-by-default algorithm experiments. +/// +/// Each experiment is a candidate change that was evaluated in LoopEval (closed-loop counterfactual +/// replay of real CGM and insulin data). They exist so the candidates can be run for testing in real +/// apps. **None of them changes behavior unless it is explicitly enabled**: with `AlgorithmExperiments.none` +/// (the default everywhere) every code path below is skipped and LoopAlgorithm produces exactly the +/// standard forecast and dose. +/// +/// The experiments, and the LoopEval candidate each one implements: +/// +/// | Property | LoopEval | What it does | +/// |---|---|---| +/// | `retrospectiveCorrectionDropScale` / `...RiseScale` | C01 (with IRC on) | Weight the falling (sensitivity) side of retrospective correction differently from the rising side | +/// | `integralRetrospectiveCorrectionMaxVelocity` | part of C01 as scored | Cap on how fast integral RC may bend the forecast | +/// | `postLowRiseCut` | C20 | After a low, treat an unexplained rise as transient until glucose proves it is a real high | +/// | `volatilityBand` | C29 | Lower the near-term forecast by the recent volatility of glucose, so the low guard sees it | +/// | `calmHighLicence` | C23 / C25 | A high with low volatility and no carbs on board takes a larger automatic-bolus application factor | +/// | `negativeInsulinDamper` | C37 | Damp forecast rises while the recent insulin effect is pushing glucose up (negative insulin) | +/// +/// The LoopEval "deployable stack" (C30) is `volatilityBand` + `calmHighLicence` + `postLowRiseCut` +/// together, with per-person constants; see `AlgorithmExperiments.stack(...)`. +public struct AlgorithmExperiments: Codable, Equatable, Sendable { + + /// Multiplier on the retrospective-correction effect when the current discrepancy is negative + /// (glucose falling faster than predicted). Applies to standard and integral RC. 1.0 = standard. + public var retrospectiveCorrectionDropScale: Double + + /// Multiplier on the retrospective-correction effect when the current discrepancy is positive + /// (glucose rising faster than predicted). Applies to standard and integral RC. 1.0 = standard. + public var retrospectiveCorrectionRiseScale: Double + + /// Ceiling on the integral-RC correction rate, in mg/dL per minute. nil = no ceiling (standard). + /// Only used when integral retrospective correction is enabled. + public var integralRetrospectiveCorrectionMaxVelocity: Double? + + public var postLowRiseCut: PostLowRiseCut? + public var volatilityBand: VolatilityBand? + public var calmHighLicence: CalmHighLicence? + public var negativeInsulinDamper: NegativeInsulinDamper? + + public init( + retrospectiveCorrectionDropScale: Double = 1.0, + retrospectiveCorrectionRiseScale: Double = 1.0, + integralRetrospectiveCorrectionMaxVelocity: Double? = nil, + postLowRiseCut: PostLowRiseCut? = nil, + volatilityBand: VolatilityBand? = nil, + calmHighLicence: CalmHighLicence? = nil, + negativeInsulinDamper: NegativeInsulinDamper? = nil + ) { + self.retrospectiveCorrectionDropScale = retrospectiveCorrectionDropScale + self.retrospectiveCorrectionRiseScale = retrospectiveCorrectionRiseScale + self.integralRetrospectiveCorrectionMaxVelocity = integralRetrospectiveCorrectionMaxVelocity + self.postLowRiseCut = postLowRiseCut + self.volatilityBand = volatilityBand + self.calmHighLicence = calmHighLicence + self.negativeInsulinDamper = negativeInsulinDamper + } + + /// No experiments: standard LoopAlgorithm. + public static let none = AlgorithmExperiments() + + /// True when no experiment is enabled. + public var isNone: Bool { self == .none } + + var needsVolatility: Bool { volatilityBand != nil || calmHighLicence != nil } + + /// C01, asymmetric integral RC, as scored in LoopEval. Use together with + /// `useIntegralRetrospectiveCorrection = true`. + public static let asymmetricIntegralRC = AlgorithmExperiments( + retrospectiveCorrectionDropScale: 1.5, + retrospectiveCorrectionRiseScale: 0.5, + integralRetrospectiveCorrectionMaxVelocity: 4.0 + ) + + /// C30, the LoopEval deployable stack. `medianVolatility` is the person's own median of + /// `GlucoseVolatility.sigma5` (mg/dL per 5 min) over the last few weeks; LoopEval found 4 weeks + /// of history sufficient (9 % drift). + public static func stack(medianVolatility: Double) -> AlgorithmExperiments { + AlgorithmExperiments( + postLowRiseCut: PostLowRiseCut(), + volatilityBand: VolatilityBand(gain: Swift.min(1.0, 6.7 / Swift.max(medianVolatility, 1e-6))), + calmHighLicence: CalmHighLicence(maxVolatility: medianVolatility) + ) + } +} + +extension AlgorithmExperiments { + + /// C20: post-low-gated RC rise-cut. Within `window` of a glucose reading below `lowThreshold`, + /// and while current glucose is below `releaseGlucose`, the rising side of retrospective + /// correction is scaled by `riseScale` (replacing `retrospectiveCorrectionRiseScale`). The + /// rebound after a low is mostly rescue carbs and counter-regulation, not a lasting need for + /// insulin; above `releaseGlucose` the high is treated as real and RC is normal. + public struct PostLowRiseCut: Codable, Equatable, Sendable { + public var riseScale: Double + public var window: TimeInterval + public var lowThreshold: Double // mg/dL + public var releaseGlucose: Double // mg/dL + + public init(riseScale: Double = 0.3, window: TimeInterval = 720 * 60, + lowThreshold: Double = 70, releaseGlucose: Double = 180) { + self.riseScale = riseScale + self.window = window + self.lowThreshold = lowThreshold + self.releaseGlucose = releaseGlucose + } + } + + /// C29: volatility-widened lower forecast. Each predicted point τ minutes ahead is lowered by + /// `gain × σ5 × (τ/5)^exponent` up to `horizon`, tapering to zero by `taperEnd`. Eventual glucose is + /// untouched; the predicted minimum, which drives the low guard and suspends, sees the volatility. + /// With `carbGate`, the band is off while carbs are on board. + public struct VolatilityBand: Codable, Equatable, Sendable { + public var gain: Double + public var horizon: TimeInterval + public var taperEnd: TimeInterval + public var exponent: Double + public var carbGate: Bool + + public init(gain: Double, horizon: TimeInterval = 60 * 60, taperEnd: TimeInterval = 120 * 60, + exponent: Double = 0.71, carbGate: Bool = true) { + self.gain = gain + self.horizon = horizon + self.taperEnd = taperEnd + self.exponent = exponent + self.carbGate = carbGate + } + } + + /// C23 / C25: calm-high licence. When current glucose is at or above `minGlucose`, recent + /// volatility σ5 is at or below `maxVolatility` (the person's own median), and (with `carbGate`) + /// no carbs are on board, the automatic-bolus application factor is multiplied by + /// `applicationFactorScale`, capped at 1. Only affects the automatic-bolus dosing strategy. + public struct CalmHighLicence: Codable, Equatable, Sendable { + public var applicationFactorScale: Double + public var minGlucose: Double // mg/dL + public var maxVolatility: Double // mg/dL per 5 min + public var carbGate: Bool + + public init(applicationFactorScale: Double = 2.0, minGlucose: Double = 180, + maxVolatility: Double, carbGate: Bool = true) { + self.applicationFactorScale = applicationFactorScale + self.minGlucose = minGlucose + self.maxVolatility = maxVolatility + self.carbGate = carbGate + } + } + + /// C37: Negative Insulin Damper (Loop and Learn customization). The glucose effect of insulin + /// delivered up to `lag` ago is summed over its upward steps (the push from negative insulin, + /// i.e. delivery below scheduled basal); the larger that push relative to an anchor of + /// 0.8 × peak activity time × basal × ISF, the more forecast rises are damped. + public struct NegativeInsulinDamper: Codable, Equatable, Sendable { + public var anchorAlpha: Double + public var marginalSlope: Double + public var lag: TimeInterval + + public init(anchorAlpha: Double = 0.75, marginalSlope: Double = 0.05, lag: TimeInterval = 15 * 60) { + self.anchorAlpha = anchorAlpha + self.marginalSlope = marginalSlope + self.lag = lag + } + } +} diff --git a/Sources/LoopAlgorithm/Experiments/ExperimentalForecastAdjustments.swift b/Sources/LoopAlgorithm/Experiments/ExperimentalForecastAdjustments.swift new file mode 100644 index 0000000..5d284bc --- /dev/null +++ b/Sources/LoopAlgorithm/Experiments/ExperimentalForecastAdjustments.swift @@ -0,0 +1,139 @@ +// +// ExperimentalForecastAdjustments.swift +// LoopAlgorithm +// +// Copyright © 2026 LoopKit Authors. All rights reserved. +// + +import Foundation + +extension AlgorithmExperiments.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 AlgorithmExperiments.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 { + + /// Applies the forecast-side experiments, in order: negative insulin damper, then volatility band. + /// Returns `prediction.glucose` unchanged when neither is enabled. + public static func applyExperimentalForecastAdjustments( + to prediction: LoopPrediction, + experiments: AlgorithmExperiments, + volatility: Double?, + basal: [AbsoluteScheduleValue], + sensitivity: [AbsoluteScheduleValue], + insulinModel: InsulinModel + ) -> [PredictedGlucoseValue] { + var glucose = prediction.glucose + if let damper = experiments.negativeInsulinDamper { + glucose = negativeInsulinDamped(glucose, damper: damper, dosesRelativeToBasal: prediction.dosesRelativeToBasal, + basal: basal, sensitivity: sensitivity, insulinModel: insulinModel) + } + if let band = experiments.volatilityBand, let volatility, + !(band.carbGate && (prediction.activeCarbs ?? 0) > 0) { + glucose = volatilityBanded(glucose, band: band, volatility: volatility) + } + return glucose + } + + static func volatilityBanded(_ glucose: [PredictedGlucoseValue], band: AlgorithmExperiments.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)) + } + } + + /// Port of the Loop and Learn Negative Insulin Damper (LoopDataManager.computeNegativeInsulinDamper). + static func negativeInsulinDamped(_ glucose: [PredictedGlucoseValue], damper: AlgorithmExperiments.NegativeInsulinDamper, + dosesRelativeToBasal: [BasalRelativeDose], basal: [AbsoluteScheduleValue], + sensitivity: [AbsoluteScheduleValue], insulinModel: InsulinModel) -> [PredictedGlucoseValue] { + let unit = LoopUnit.milligramsPerDeciliter + guard let anchorDate = glucose.first?.startDate else { return glucose } + let cutoff = anchorDate.addingTimeInterval(-damper.lag) + + // Doses as of `lag` ago: later doses dropped, basal deliveries clipped at the cutoff (net units pro-rated), + // boluses kept whole. + var lagged: [BasalRelativeDose] = [] + for dose in dosesRelativeToBasal where dose.startDate < cutoff { + if case .bolus = dose.type { lagged.append(dose); continue } + if dose.endDate <= cutoff { lagged.append(dose); continue } + var clipped = dose + let full = dose.endDate.timeIntervalSince(dose.startDate) + clipped.endDate = cutoff + clipped.volume = full > 0 ? dose.volume * cutoff.timeIntervalSince(dose.startDate) / full : dose.volume + lagged.append(clipped) + } + let effects = lagged.glucoseEffects(insulinSensitivityHistory: sensitivity, from: anchorDate.addingTimeInterval(-.minutes(5))) + var positiveSum = 0.0 + if effects.count > 1 { + for i in 1.. 0 else { return glucose } // a 0 U/h basal segment disables the damper + + let a0 = damper.anchorAlpha, ms = damper.marginalSlope + let lin = (1.0 - a0) / anchorPoint + let tp = (1 - ms) / (2 * lin) + let alpha: Double + if positiveSum < tp { + alpha = 1 - lin * positiveSum + } else { + let tv = (1 - lin * tp) * tp + alpha = (tv + ms * (positiveSum - tp)) / positiveSum + } + let damping = Swift.max(0, 1 - Swift.max(ms, alpha)) + guard damping > 0 else { return glucose } + + let keep = 1 - damping + var damped: [PredictedGlucoseValue] = [] + damped.reserveCapacity(glucose.count) + var value = 0.0 + for (i, point) in glucose.enumerated() { + let x = point.quantity.doubleValue(for: unit) + if i == 0 { value = x; damped.append(point); continue } + let delta = x - glucose[i - 1].quantity.doubleValue(for: unit) + value += delta > 0 ? keep * delta : delta + damped.append(PredictedGlucoseValue(startDate: point.startDate, quantity: LoopQuantity(unit: unit, doubleValue: value))) + } + return damped + } + + static func peakActivity(of model: InsulinModel) -> TimeInterval { + if let preset = model as? ExponentialInsulinModelPreset { return preset.peakActivity } + if let exponential = model as? ExponentialInsulinModel { return exponential.peakActivityTime } + return .minutes(75) + } +} diff --git a/Sources/LoopAlgorithm/Experiments/GlucoseVolatility.swift b/Sources/LoopAlgorithm/Experiments/GlucoseVolatility.swift new file mode 100644 index 0000000..4a6edfc --- /dev/null +++ b/Sources/LoopAlgorithm/Experiments/GlucoseVolatility.swift @@ -0,0 +1,43 @@ +// +// GlucoseVolatility.swift +// LoopAlgorithm +// +// Copyright © 2026 LoopKit Authors. All rights reserved. +// + +import Foundation + +/// Causal short-term glucose volatility, used by the `volatilityBand` and `calmHighLicence` experiments. +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: the exponentially weighted RMS of the 5-minute glucose increment, in mg/dL per 5 min, + /// over the last ~5 h of readings. Readings 2–7 min apart contribute (scaled to 5 min); a gap of + /// more than 7 min restarts the estimate. With no usable pair the floor sqrt(2 × noise²) is returned. + public static func sigma5(_ glucoseHistory: [GlucoseType]) -> Double { + let unit = LoopUnit.milligramsPerDeciliter + let floorVariance = 2.0 * noise * noise + var variance = Double.nan + var previous: (date: Date, value: Double)? = nil + for sample in glucoseHistory.suffix(60) { + let value = sample.quantity.doubleValue(for: unit) + if let previous { + let minutes = sample.startDate.timeIntervalSince(previous.date) / 60.0 + if minutes > 2.0 && minutes <= 7.0 { + let d5 = (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.startDate, value) + } + return variance.isNaN ? floorVariance.squareRoot() : variance.squareRoot() + } +} diff --git a/Sources/LoopAlgorithm/LoopAlgorithm.swift b/Sources/LoopAlgorithm/LoopAlgorithm.swift index 2e711a1..6790ef4 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, + experiments: AlgorithmExperiments = .none ) -> LoopPrediction where CarbType: CarbEntry, GlucoseType: GlucoseSampleValue, InsulinDoseType: InsulinDose { var prediction: [PredictedGlucoseValue] = [] @@ -254,10 +255,25 @@ public struct LoopAlgorithm { let rc: RetrospectiveCorrection + // Experiments (all no-ops at their defaults): asymmetric RC gains, and the post-low rise-cut, + // which replaces the rise gain while its gate is open. + let rcDropScale = experiments.retrospectiveCorrectionDropScale + var rcRiseScale = experiments.retrospectiveCorrectionRiseScale + if let cut = experiments.postLowRiseCut, cut.isActive(glucoseHistory: glucoseHistory) { + rcRiseScale = cut.riseScale + } + if useIntegralRetrospectiveCorrection { - rc = IntegralRetrospectiveCorrection(effectDuration: LoopMath.retrospectiveCorrectionEffectDuration) + rc = IntegralRetrospectiveCorrection( + effectDuration: LoopMath.retrospectiveCorrectionEffectDuration, + dropScale: rcDropScale, + riseScale: rcRiseScale, + maxVelocity: experiments.integralRetrospectiveCorrectionMaxVelocity) } else { - rc = StandardRetrospectiveCorrection(effectDuration: LoopMath.retrospectiveCorrectionEffectDuration) + rc = StandardRetrospectiveCorrection( + effectDuration: LoopMath.retrospectiveCorrectionEffectDuration, + dropScale: rcDropScale, + riseScale: rcRiseScale) } if let latestGlucose = glucoseHistory.last { @@ -550,9 +566,22 @@ public struct LoopAlgorithm { includingPositiveVelocityAndRC: input.includePositiveVelocityAndRC, useMidAbsorptionISF: input.useMidAbsorptionISF, carbAbsorptionModel: input.carbAbsorptionModel.model, - gradualTransitionsThreshold: input.gradualTransitionsThreshold + gradualTransitionsThreshold: input.gradualTransitionsThreshold, + experiments: input.experiments ) + // Experiments that adjust the forecast before dosing (no-ops unless enabled). + let volatility: Double? = input.experiments.needsVolatility ? GlucoseVolatility.sigma5(input.glucoseHistory) : nil + if !input.experiments.isNone { + prediction.glucose = applyExperimentalForecastAdjustments( + to: prediction, + experiments: input.experiments, + volatility: volatility, + basal: input.basal, + sensitivity: input.sensitivity, + insulinModel: input.recommendationInsulinModel) + } + let sensitivityForDosing: [AbsoluteScheduleValue] if input.useMidAbsorptionISF { sensitivityForDosing = input.sensitivity @@ -587,9 +616,16 @@ public struct LoopAlgorithm { target: input.target) result = .success(.init(manual: recommendation)) case .automaticBolus: + var applicationFactor = input.automaticBolusApplicationFactor ?? defaultBolusPartialApplicationFactor + if let licence = input.experiments.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..ad2d88e --- /dev/null +++ b/Tests/LoopAlgorithmTests/AlgorithmExperimentsTests.swift @@ -0,0 +1,205 @@ +// +// 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() { + XCTAssertTrue(loadScenario("suspend").experiments.isNone) + XCTAssertTrue(loadScenario("carbs_with_isf_change").experiments.isNone) + } + + func testExplicitNoneMatchesStandard() { + for name in ["suspend", "carbs_with_isf_change"] { + let standard = loadScenario(name) + var explicit = standard + explicit.algorithmExperiments = AlgorithmExperiments() + assertSameOutput(LoopAlgorithm.run(input: standard), LoopAlgorithm.run(input: explicit)) + } + } + + func testUnityScalesMatchStandard() { + var input = loadScenario("carbs_with_isf_change") + let standard = LoopAlgorithm.run(input: input) + input.algorithmExperiments = AlgorithmExperiments(retrospectiveCorrectionDropScale: 1.0, retrospectiveCorrectionRiseScale: 1.0) + assertSameOutput(standard, LoopAlgorithm.run(input: input)) + } + + func testGatedExperimentsAreInertWhenTheirGateIsClosed() { + // Mock: glucose rising 100 -> 160, no low, below 180, no carbs. The post-low cut has no low to + // key on and the calm-high licence needs >= 180, so both must leave the output unchanged. + let now = ISO8601DateFormatter().date(from: "2024-01-03T12:00:00+0000")! + var input = AlgorithmInputFixture.mock(for: now) + input.recommendationType = .automaticBolus + let standard = LoopAlgorithm.run(input: input) + input.algorithmExperiments = AlgorithmExperiments( + postLowRiseCut: .init(), + calmHighLicence: .init(maxVolatility: 100)) + assertSameOutput(standard, LoopAlgorithm.run(input: input)) + } + + // 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 = AlgorithmExperiments.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.algorithmExperiments = AlgorithmExperiments(postLowRiseCut: .init()) + 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 testVolatilityBandShape() { + let now = Date() + let glucose = (0...48).map { PredictedGlucoseValue(startDate: now.addingTimeInterval(.minutes(5 * Double($0))), quantity: .glucose(120)) } + let band = AlgorithmExperiments.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.algorithmExperiments = AlgorithmExperiments(calmHighLicence: .init(maxVolatility: 5)) + let licensed = LoopAlgorithm.run(input: input).recommendation!.automatic!.bolusUnits! + XCTAssertGreaterThan(standard, 0) + XCTAssertEqual(licensed, min(2 * standard, standard / 0.4), accuracy: 0.051) // AF 0.4 -> 0.8, pump-rounded + // Volatile high: no licence. + input.algorithmExperiments = AlgorithmExperiments(calmHighLicence: .init(maxVolatility: 0.5)) + XCTAssertEqual(LoopAlgorithm.run(input: input).recommendation!.automatic!.bolusUnits!, standard, accuracy: 1e-9) + } + + // MARK: - Negative insulin damper + + func testNegativeInsulinDamperDampsRisesOnly() { + let now = Date() + // Three hours of suspended basal (scheduled 1 U/h) ending 20 min ago: negative insulin pushing glucose up. + let suspend = BasalRelativeDose(type: .basal(scheduledRate: 1.0), startDate: now.addingTimeInterval(.hours(-3)), + endDate: now.addingTimeInterval(.minutes(-20)), volume: 0) + let values: [Double] = [120, 130, 140, 135, 145] + let glucose = values.enumerated().map { PredictedGlucoseValue(startDate: now.addingTimeInterval(.minutes(5 * Double($0.offset))), quantity: .glucose($0.element)) } + let basal = [AbsoluteScheduleValue(startDate: now.addingTimeInterval(.hours(-10)), endDate: now.addingTimeInterval(.hours(6)), value: 1.0)] + let isf = [AbsoluteScheduleValue(startDate: now.addingTimeInterval(.hours(-10)), endDate: now.addingTimeInterval(.hours(6)), value: LoopQuantity.glucose(50))] + let damped = LoopAlgorithm.negativeInsulinDamped(glucose, damper: .init(), dosesRelativeToBasal: [suspend], basal: basal, + sensitivity: isf, insulinModel: ExponentialInsulinModelPreset.rapidActingAdult) + let d = damped.map { $0.quantity.doubleValue(for: unit) } + XCTAssertEqual(d[0], 120, accuracy: 1e-9) + XCTAssertLessThan(d[1] - d[0], 10) // rise damped + XCTAssertEqual(d[3] - d[2], -5, accuracy: 1e-9) // fall kept whole + // No negative insulin: unchanged. + let none = LoopAlgorithm.negativeInsulinDamped(glucose, damper: .init(), dosesRelativeToBasal: [], basal: basal, + sensitivity: isf, insulinModel: ExponentialInsulinModelPreset.rapidActingAdult) + XCTAssertEqual(none.map { $0.quantity.doubleValue(for: unit) }, values) + } + + // MARK: - Codable + + func testExperimentsRoundTrip() throws { + let stack = AlgorithmExperiments.stack(medianVolatility: 5.0) + let data = try JSONEncoder().encode(stack) + XCTAssertEqual(try JSONDecoder().decode(AlgorithmExperiments.self, from: data), stack) + XCTAssertEqual(stack.volatilityBand?.gain, 1.0) // min(1, 6.7 / 5) + XCTAssertEqual(AlgorithmExperiments.stack(medianVolatility: 10).volatilityBand!.gain, 0.67, accuracy: 1e-12) + } +} From 8b5957c3ea59fc5d52cad4cfd1c27d12d463c8fa Mon Sep 17 00:00:00 2001 From: Pete Schwamb Date: Thu, 15 Oct 2026 13:08:29 -0500 Subject: [PATCH 2/2] Reduce the experiments to two fixed choices: asymmetric integral RC and the stack The experiments will be offered to people who choose to run them and report results; fewer choices pool more people per experiment. The public API is now AlgorithmExperiment, an enum with two cases and no tunable parameters: - asymmetricIntegralRC (LoopEval C01): integral RC switched on by the experiment, falling side x1.5, rising side x0.5, rate capped at 4 mg/dL/min - stack (LoopEval C30): post-low RC rise-cut + volatility band + calm-high licence, keyed to the person's volatility baseline AlgorithmInput gains experiment and volatilityBaseline (both default nil). The component parameters are internal so they cannot be mixed. The negative insulin damper is removed. GlucoseVolatility.baseline(_:) computes the person's median sigma5 from a glucose history; when no baseline is supplied the stack computes it from the input history (needs a day of readings) or runs its post-low cut alone. Glucose is thinned to 5-minute spacing before sigma5: with 1-minute CGM the estimator otherwise sees only increments it ignores and never leaves its floor. Identical to the evaluated estimator on 5-minute data. --- Sources/LoopAlgorithm/AlgorithmInput.swift | 14 +- .../LoopAlgorithm/AlgorithmInputFixture.swift | 15 +- .../Experiments/AlgorithmExperiments.swift | 219 ++++++------------ .../ExperimentalForecastAdjustments.swift | 95 +------- .../Experiments/GlucoseVolatility.swift | 50 +++- Sources/LoopAlgorithm/LoopAlgorithm.swift | 41 ++-- .../AlgorithmExperimentsTests.swift | 136 ++++++----- 7 files changed, 234 insertions(+), 336 deletions(-) diff --git a/Sources/LoopAlgorithm/AlgorithmInput.swift b/Sources/LoopAlgorithm/AlgorithmInput.swift index 7e7e884..2bec7ec 100644 --- a/Sources/LoopAlgorithm/AlgorithmInput.swift +++ b/Sources/LoopAlgorithm/AlgorithmInput.swift @@ -34,11 +34,17 @@ public protocol AlgorithmInput { var recommendationType: DoseRecommendationType { get } var automaticBolusApplicationFactor: Double? { get } // Defaults to 0.4 var gradualTransitionsThreshold: Double? { get } - /// Optional algorithm experiments. Defaults to `.none` (standard LoopAlgorithm) for any conforming - /// type that does not provide it. - var experiments: AlgorithmExperiments { 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 experiments: AlgorithmExperiments { .none } + public var experiment: AlgorithmExperiment? { nil } + public var volatilityBaseline: Double? { nil } } diff --git a/Sources/LoopAlgorithm/AlgorithmInputFixture.swift b/Sources/LoopAlgorithm/AlgorithmInputFixture.swift index 98d9471..4af45be 100644 --- a/Sources/LoopAlgorithm/AlgorithmInputFixture.swift +++ b/Sources/LoopAlgorithm/AlgorithmInputFixture.swift @@ -35,12 +35,9 @@ public struct AlgorithmInputFixture: AlgorithmInput { public var recommendationType: DoseRecommendationType = .automaticBolus public var automaticBolusApplicationFactor: Double? public var gradualTransitionsThreshold: Double? - /// Optional experiments; absent in existing fixtures, which therefore run standard LoopAlgorithm. - public var algorithmExperiments: AlgorithmExperiments? - - public var experiments: AlgorithmExperiments { - algorithmExperiments ?? .none - } + /// 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 @@ -154,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) } @@ -190,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 { @@ -212,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 index 7005f27..1f4f714 100644 --- a/Sources/LoopAlgorithm/Experiments/AlgorithmExperiments.swift +++ b/Sources/LoopAlgorithm/Experiments/AlgorithmExperiments.swift @@ -7,167 +7,90 @@ import Foundation -/// Optional, off-by-default algorithm experiments. +/// An optional algorithm experiment. `nil` (the default everywhere) is standard LoopAlgorithm. /// -/// Each experiment is a candidate change that was evaluated in LoopEval (closed-loop counterfactual -/// replay of real CGM and insulin data). They exist so the candidates can be run for testing in real -/// apps. **None of them changes behavior unless it is explicitly enabled**: with `AlgorithmExperiments.none` -/// (the default everywhere) every code path below is skipped and LoopAlgorithm produces exactly the -/// standard forecast and dose. -/// -/// The experiments, and the LoopEval candidate each one implements: -/// -/// | Property | LoopEval | What it does | -/// |---|---|---| -/// | `retrospectiveCorrectionDropScale` / `...RiseScale` | C01 (with IRC on) | Weight the falling (sensitivity) side of retrospective correction differently from the rising side | -/// | `integralRetrospectiveCorrectionMaxVelocity` | part of C01 as scored | Cap on how fast integral RC may bend the forecast | -/// | `postLowRiseCut` | C20 | After a low, treat an unexplained rise as transient until glucose proves it is a real high | -/// | `volatilityBand` | C29 | Lower the near-term forecast by the recent volatility of glucose, so the low guard sees it | -/// | `calmHighLicence` | C23 / C25 | A high with low volatility and no carbs on board takes a larger automatic-bolus application factor | -/// | `negativeInsulinDamper` | C37 | Damp forecast rises while the recent insulin effect is pushing glucose up (negative insulin) | -/// -/// The LoopEval "deployable stack" (C30) is `volatilityBand` + `calmHighLicence` + `postLowRiseCut` -/// together, with per-person constants; see `AlgorithmExperiments.stack(...)`. -public struct AlgorithmExperiments: Codable, Equatable, Sendable { - - /// Multiplier on the retrospective-correction effect when the current discrepancy is negative - /// (glucose falling faster than predicted). Applies to standard and integral RC. 1.0 = standard. - public var retrospectiveCorrectionDropScale: Double - - /// Multiplier on the retrospective-correction effect when the current discrepancy is positive - /// (glucose rising faster than predicted). Applies to standard and integral RC. 1.0 = standard. - public var retrospectiveCorrectionRiseScale: Double - - /// Ceiling on the integral-RC correction rate, in mg/dL per minute. nil = no ceiling (standard). - /// Only used when integral retrospective correction is enabled. - public var integralRetrospectiveCorrectionMaxVelocity: Double? - - public var postLowRiseCut: PostLowRiseCut? - public var volatilityBand: VolatilityBand? - public var calmHighLicence: CalmHighLicence? - public var negativeInsulinDamper: NegativeInsulinDamper? - - public init( - retrospectiveCorrectionDropScale: Double = 1.0, - retrospectiveCorrectionRiseScale: Double = 1.0, - integralRetrospectiveCorrectionMaxVelocity: Double? = nil, - postLowRiseCut: PostLowRiseCut? = nil, - volatilityBand: VolatilityBand? = nil, - calmHighLicence: CalmHighLicence? = nil, - negativeInsulinDamper: NegativeInsulinDamper? = nil - ) { - self.retrospectiveCorrectionDropScale = retrospectiveCorrectionDropScale - self.retrospectiveCorrectionRiseScale = retrospectiveCorrectionRiseScale - self.integralRetrospectiveCorrectionMaxVelocity = integralRetrospectiveCorrectionMaxVelocity - self.postLowRiseCut = postLowRiseCut - self.volatilityBand = volatilityBand - self.calmHighLicence = calmHighLicence - self.negativeInsulinDamper = negativeInsulinDamper - } +/// 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 +} - /// No experiments: standard LoopAlgorithm. - public static let none = AlgorithmExperiments() +/// 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 - /// True when no experiment is enabled. - public var isNone: Bool { self == .none } + static let none = ExperimentParameters() var needsVolatility: Bool { volatilityBand != nil || calmHighLicence != nil } - /// C01, asymmetric integral RC, as scored in LoopEval. Use together with - /// `useIntegralRetrospectiveCorrection = true`. - public static let asymmetricIntegralRC = AlgorithmExperiments( - retrospectiveCorrectionDropScale: 1.5, - retrospectiveCorrectionRiseScale: 0.5, - integralRetrospectiveCorrectionMaxVelocity: 4.0 - ) - - /// C30, the LoopEval deployable stack. `medianVolatility` is the person's own median of - /// `GlucoseVolatility.sigma5` (mg/dL per 5 min) over the last few weeks; LoopEval found 4 weeks - /// of history sufficient (9 % drift). - public static func stack(medianVolatility: Double) -> AlgorithmExperiments { - AlgorithmExperiments( - postLowRiseCut: PostLowRiseCut(), - volatilityBand: VolatilityBand(gain: Swift.min(1.0, 6.7 / Swift.max(medianVolatility, 1e-6))), - calmHighLicence: CalmHighLicence(maxVolatility: medianVolatility) - ) - } -} - -extension AlgorithmExperiments { - - /// C20: post-low-gated RC rise-cut. Within `window` of a glucose reading below `lowThreshold`, - /// and while current glucose is below `releaseGlucose`, the rising side of retrospective - /// correction is scaled by `riseScale` (replacing `retrospectiveCorrectionRiseScale`). The - /// rebound after a low is mostly rescue carbs and counter-regulation, not a lasting need for - /// insulin; above `releaseGlucose` the high is treated as real and RC is normal. - public struct PostLowRiseCut: Codable, Equatable, Sendable { - public var riseScale: Double - public var window: TimeInterval - public var lowThreshold: Double // mg/dL - public var releaseGlucose: Double // mg/dL - - public init(riseScale: Double = 0.3, window: TimeInterval = 720 * 60, - lowThreshold: Double = 70, releaseGlucose: Double = 180) { - self.riseScale = riseScale - self.window = window - self.lowThreshold = lowThreshold - self.releaseGlucose = releaseGlucose + /// `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) + } } } - /// C29: volatility-widened lower forecast. Each predicted point τ minutes ahead is lowered by - /// `gain × σ5 × (τ/5)^exponent` up to `horizon`, tapering to zero by `taperEnd`. Eventual glucose is - /// untouched; the predicted minimum, which drives the low guard and suspends, sees the volatility. - /// With `carbGate`, the band is off while carbs are on board. - public struct VolatilityBand: Codable, Equatable, Sendable { - public var gain: Double - public var horizon: TimeInterval - public var taperEnd: TimeInterval - public var exponent: Double - public var carbGate: Bool + init() {} - public init(gain: Double, horizon: TimeInterval = 60 * 60, taperEnd: TimeInterval = 120 * 60, - exponent: Double = 0.71, carbGate: Bool = true) { - self.gain = gain - self.horizon = horizon - self.taperEnd = taperEnd - self.exponent = exponent - self.carbGate = carbGate - } + struct PostLowRiseCut: Equatable { + var riseScale = 0.3 + var window: TimeInterval = .hours(12) + var lowThreshold = 70.0 // mg/dL + var releaseGlucose = 180.0 // mg/dL } - /// C23 / C25: calm-high licence. When current glucose is at or above `minGlucose`, recent - /// volatility σ5 is at or below `maxVolatility` (the person's own median), and (with `carbGate`) - /// no carbs are on board, the automatic-bolus application factor is multiplied by - /// `applicationFactorScale`, capped at 1. Only affects the automatic-bolus dosing strategy. - public struct CalmHighLicence: Codable, Equatable, Sendable { - public var applicationFactorScale: Double - public var minGlucose: Double // mg/dL - public var maxVolatility: Double // mg/dL per 5 min - public var carbGate: Bool - - public init(applicationFactorScale: Double = 2.0, minGlucose: Double = 180, - maxVolatility: Double, carbGate: Bool = true) { - self.applicationFactorScale = applicationFactorScale - self.minGlucose = minGlucose - self.maxVolatility = maxVolatility - self.carbGate = carbGate - } + struct VolatilityBand: Equatable { + var gain: Double + var horizon: TimeInterval = .minutes(60) + var taperEnd: TimeInterval = .minutes(120) + var exponent = 0.71 + var carbGate = true } - /// C37: Negative Insulin Damper (Loop and Learn customization). The glucose effect of insulin - /// delivered up to `lag` ago is summed over its upward steps (the push from negative insulin, - /// i.e. delivery below scheduled basal); the larger that push relative to an anchor of - /// 0.8 × peak activity time × basal × ISF, the more forecast rises are damped. - public struct NegativeInsulinDamper: Codable, Equatable, Sendable { - public var anchorAlpha: Double - public var marginalSlope: Double - public var lag: TimeInterval - - public init(anchorAlpha: Double = 0.75, marginalSlope: Double = 0.05, lag: TimeInterval = 15 * 60) { - self.anchorAlpha = anchorAlpha - self.marginalSlope = marginalSlope - self.lag = lag - } + 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 index 5d284bc..020b2a2 100644 --- a/Sources/LoopAlgorithm/Experiments/ExperimentalForecastAdjustments.swift +++ b/Sources/LoopAlgorithm/Experiments/ExperimentalForecastAdjustments.swift @@ -7,7 +7,7 @@ import Foundation -extension AlgorithmExperiments.PostLowRiseCut { +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 { @@ -22,7 +22,7 @@ extension AlgorithmExperiments.PostLowRiseCut { } } -extension AlgorithmExperiments.CalmHighLicence { +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 @@ -31,29 +31,9 @@ extension AlgorithmExperiments.CalmHighLicence { extension LoopAlgorithm { - /// Applies the forecast-side experiments, in order: negative insulin damper, then volatility band. - /// Returns `prediction.glucose` unchanged when neither is enabled. - public static func applyExperimentalForecastAdjustments( - to prediction: LoopPrediction, - experiments: AlgorithmExperiments, - volatility: Double?, - basal: [AbsoluteScheduleValue], - sensitivity: [AbsoluteScheduleValue], - insulinModel: InsulinModel - ) -> [PredictedGlucoseValue] { - var glucose = prediction.glucose - if let damper = experiments.negativeInsulinDamper { - glucose = negativeInsulinDamped(glucose, damper: damper, dosesRelativeToBasal: prediction.dosesRelativeToBasal, - basal: basal, sensitivity: sensitivity, insulinModel: insulinModel) - } - if let band = experiments.volatilityBand, let volatility, - !(band.carbGate && (prediction.activeCarbs ?? 0) > 0) { - glucose = volatilityBanded(glucose, band: band, volatility: volatility) - } - return glucose - } - - static func volatilityBanded(_ glucose: [PredictedGlucoseValue], band: AlgorithmExperiments.VolatilityBand, + /// 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 @@ -71,69 +51,4 @@ extension LoopAlgorithm { quantity: LoopQuantity(unit: unit, doubleValue: point.quantity.doubleValue(for: unit) + offset)) } } - - /// Port of the Loop and Learn Negative Insulin Damper (LoopDataManager.computeNegativeInsulinDamper). - static func negativeInsulinDamped(_ glucose: [PredictedGlucoseValue], damper: AlgorithmExperiments.NegativeInsulinDamper, - dosesRelativeToBasal: [BasalRelativeDose], basal: [AbsoluteScheduleValue], - sensitivity: [AbsoluteScheduleValue], insulinModel: InsulinModel) -> [PredictedGlucoseValue] { - let unit = LoopUnit.milligramsPerDeciliter - guard let anchorDate = glucose.first?.startDate else { return glucose } - let cutoff = anchorDate.addingTimeInterval(-damper.lag) - - // Doses as of `lag` ago: later doses dropped, basal deliveries clipped at the cutoff (net units pro-rated), - // boluses kept whole. - var lagged: [BasalRelativeDose] = [] - for dose in dosesRelativeToBasal where dose.startDate < cutoff { - if case .bolus = dose.type { lagged.append(dose); continue } - if dose.endDate <= cutoff { lagged.append(dose); continue } - var clipped = dose - let full = dose.endDate.timeIntervalSince(dose.startDate) - clipped.endDate = cutoff - clipped.volume = full > 0 ? dose.volume * cutoff.timeIntervalSince(dose.startDate) / full : dose.volume - lagged.append(clipped) - } - let effects = lagged.glucoseEffects(insulinSensitivityHistory: sensitivity, from: anchorDate.addingTimeInterval(-.minutes(5))) - var positiveSum = 0.0 - if effects.count > 1 { - for i in 1.. 0 else { return glucose } // a 0 U/h basal segment disables the damper - - let a0 = damper.anchorAlpha, ms = damper.marginalSlope - let lin = (1.0 - a0) / anchorPoint - let tp = (1 - ms) / (2 * lin) - let alpha: Double - if positiveSum < tp { - alpha = 1 - lin * positiveSum - } else { - let tv = (1 - lin * tp) * tp - alpha = (tv + ms * (positiveSum - tp)) / positiveSum - } - let damping = Swift.max(0, 1 - Swift.max(ms, alpha)) - guard damping > 0 else { return glucose } - - let keep = 1 - damping - var damped: [PredictedGlucoseValue] = [] - damped.reserveCapacity(glucose.count) - var value = 0.0 - for (i, point) in glucose.enumerated() { - let x = point.quantity.doubleValue(for: unit) - if i == 0 { value = x; damped.append(point); continue } - let delta = x - glucose[i - 1].quantity.doubleValue(for: unit) - value += delta > 0 ? keep * delta : delta - damped.append(PredictedGlucoseValue(startDate: point.startDate, quantity: LoopQuantity(unit: unit, doubleValue: value))) - } - return damped - } - - static func peakActivity(of model: InsulinModel) -> TimeInterval { - if let preset = model as? ExponentialInsulinModelPreset { return preset.peakActivity } - if let exponential = model as? ExponentialInsulinModel { return exponential.peakActivityTime } - return .minutes(75) - } } diff --git a/Sources/LoopAlgorithm/Experiments/GlucoseVolatility.swift b/Sources/LoopAlgorithm/Experiments/GlucoseVolatility.swift index 4a6edfc..46de30e 100644 --- a/Sources/LoopAlgorithm/Experiments/GlucoseVolatility.swift +++ b/Sources/LoopAlgorithm/Experiments/GlucoseVolatility.swift @@ -7,7 +7,7 @@ import Foundation -/// Causal short-term glucose volatility, used by the `volatilityBand` and `calmHighLicence` experiments. +/// Causal short-term glucose volatility, used by the `stack` experiment. public enum GlucoseVolatility { /// EWMA weight on the previous variance. @@ -15,20 +15,49 @@ public enum GlucoseVolatility { /// CGM noise level; the variance is floored at 2 × noise². public static let noise = 1.4 - /// σ5: the exponentially weighted RMS of the 5-minute glucose increment, in mg/dL per 5 min, - /// over the last ~5 h of readings. Readings 2–7 min apart contribute (scaled to 5 min); a gap of - /// more than 7 min restarts the estimate. With no usable pair the floor sqrt(2 × noise²) is returned. + /// σ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 - for sample in glucoseHistory.suffix(60) { - let value = sample.quantity.doubleValue(for: unit) + var out: [Double] = [] + for sample in samples { if let previous { - let minutes = sample.startDate.timeIntervalSince(previous.date) / 60.0 + let minutes = sample.date.timeIntervalSince(previous.date) / 60.0 if minutes > 2.0 && minutes <= 7.0 { - let d5 = (value - previous.value) * (5.0 / minutes) + 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) @@ -36,8 +65,9 @@ public enum GlucoseVolatility { variance = .nan } } - previous = (sample.startDate, value) + previous = sample + out.append(variance.isNaN ? floorVariance.squareRoot() : variance.squareRoot()) } - return variance.isNaN ? floorVariance.squareRoot() : variance.squareRoot() + return out } } diff --git a/Sources/LoopAlgorithm/LoopAlgorithm.swift b/Sources/LoopAlgorithm/LoopAlgorithm.swift index 6790ef4..3713995 100644 --- a/Sources/LoopAlgorithm/LoopAlgorithm.swift +++ b/Sources/LoopAlgorithm/LoopAlgorithm.swift @@ -180,7 +180,7 @@ public struct LoopAlgorithm { useMidAbsorptionISF: Bool = false, carbAbsorptionModel: CarbAbsorptionComputable = PiecewiseLinearAbsorption(), gradualTransitionsThreshold: Double? = 40.0, - experiments: AlgorithmExperiments = .none + experiment: AlgorithmExperiment? = nil ) -> LoopPrediction where CarbType: CarbEntry, GlucoseType: GlucoseSampleValue, InsulinDoseType: InsulinDose { var prediction: [PredictedGlucoseValue] = [] @@ -255,20 +255,21 @@ public struct LoopAlgorithm { let rc: RetrospectiveCorrection - // Experiments (all no-ops at their defaults): asymmetric RC gains, and the post-low rise-cut, - // which replaces the rise gain while its gate is open. - let rcDropScale = experiments.retrospectiveCorrectionDropScale - var rcRiseScale = experiments.retrospectiveCorrectionRiseScale - if let cut = experiments.postLowRiseCut, cut.isActive(glucoseHistory: glucoseHistory) { + // 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 { + if useIntegralRetrospectiveCorrection || experimentParameters.forceIntegralRetrospectiveCorrection { rc = IntegralRetrospectiveCorrection( effectDuration: LoopMath.retrospectiveCorrectionEffectDuration, dropScale: rcDropScale, riseScale: rcRiseScale, - maxVelocity: experiments.integralRetrospectiveCorrectionMaxVelocity) + maxVelocity: experimentParameters.integralRetrospectiveCorrectionMaxVelocity) } else { rc = StandardRetrospectiveCorrection( effectDuration: LoopMath.retrospectiveCorrectionEffectDuration, @@ -567,19 +568,19 @@ public struct LoopAlgorithm { useMidAbsorptionISF: input.useMidAbsorptionISF, carbAbsorptionModel: input.carbAbsorptionModel.model, gradualTransitionsThreshold: input.gradualTransitionsThreshold, - experiments: input.experiments + experiment: input.experiment ) - // Experiments that adjust the forecast before dosing (no-ops unless enabled). - let volatility: Double? = input.experiments.needsVolatility ? GlucoseVolatility.sigma5(input.glucoseHistory) : nil - if !input.experiments.isNone { - prediction.glucose = applyExperimentalForecastAdjustments( - to: prediction, - experiments: input.experiments, - volatility: volatility, - basal: input.basal, - sensitivity: input.sensitivity, - insulinModel: input.recommendationInsulinModel) + // 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] @@ -617,7 +618,7 @@ public struct LoopAlgorithm { result = .success(.init(manual: recommendation)) case .automaticBolus: var applicationFactor = input.automaticBolusApplicationFactor ?? defaultBolusPartialApplicationFactor - if let licence = input.experiments.calmHighLicence, let volatility, + if let licence = experimentParameters.calmHighLicence, let volatility, licence.isActive(currentGlucose: latestGlucose.quantity.doubleValue(for: .milligramsPerDeciliter), volatility: volatility, activeCarbs: prediction.activeCarbs ?? 0) { diff --git a/Tests/LoopAlgorithmTests/AlgorithmExperimentsTests.swift b/Tests/LoopAlgorithmTests/AlgorithmExperimentsTests.swift index ad2d88e..cb9ebed 100644 --- a/Tests/LoopAlgorithmTests/AlgorithmExperimentsTests.swift +++ b/Tests/LoopAlgorithmTests/AlgorithmExperimentsTests.swift @@ -35,37 +35,55 @@ final class AlgorithmExperimentsTests: XCTestCase { // MARK: - Default is standard LoopAlgorithm func testExistingFixturesDecodeWithoutExperiments() { - XCTAssertTrue(loadScenario("suspend").experiments.isNone) - XCTAssertTrue(loadScenario("carbs_with_isf_change").experiments.isNone) + XCTAssertNil(loadScenario("suspend").experiment) + XCTAssertNil(loadScenario("carbs_with_isf_change").experiment) } - func testExplicitNoneMatchesStandard() { - for name in ["suspend", "carbs_with_isf_change"] { - let standard = loadScenario(name) - var explicit = standard - explicit.algorithmExperiments = AlgorithmExperiments() - assertSameOutput(LoopAlgorithm.run(input: standard), LoopAlgorithm.run(input: explicit)) - } - } - - func testUnityScalesMatchStandard() { - var input = loadScenario("carbs_with_isf_change") + 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.algorithmExperiments = AlgorithmExperiments(retrospectiveCorrectionDropScale: 1.0, retrospectiveCorrectionRiseScale: 1.0) + input.experiment = .stack + input.volatilityBaseline = 5 assertSameOutput(standard, LoopAlgorithm.run(input: input)) } - func testGatedExperimentsAreInertWhenTheirGateIsClosed() { - // Mock: glucose rising 100 -> 160, no low, below 180, no carbs. The post-low cut has no low to - // key on and the calm-high licence needs >= 180, so both must leave the output unchanged. + 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) - input.recommendationType = .automaticBolus let standard = LoopAlgorithm.run(input: input) - input.algorithmExperiments = AlgorithmExperiments( - postLowRiseCut: .init(), - calmHighLicence: .init(maxVolatility: 100)) + 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 @@ -106,7 +124,7 @@ final class AlgorithmExperimentsTests: XCTestCase { func testPostLowRiseCutGate() { let now = Date() - let cut = AlgorithmExperiments.PostLowRiseCut() + 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") @@ -121,7 +139,7 @@ final class AlgorithmExperimentsTests: XCTestCase { 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.algorithmExperiments = AlgorithmExperiments(postLowRiseCut: .init()) + input.experiment = .stack let cut = LoopAlgorithm.run(input: input) XCTAssertLessThan(cut.predictedGlucose.last!.quantity.doubleValue(for: unit), standard.predictedGlucose.last!.quantity.doubleValue(for: unit)) @@ -140,10 +158,28 @@ final class AlgorithmExperimentsTests: XCTestCase { 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 = AlgorithmExperiments.VolatilityBand(gain: 1.0) + 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) @@ -161,45 +197,29 @@ final class AlgorithmExperimentsTests: XCTestCase { 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.algorithmExperiments = AlgorithmExperiments(calmHighLicence: .init(maxVolatility: 5)) + input.experiment = .stack + input.volatilityBaseline = 5 let licensed = LoopAlgorithm.run(input: input).recommendation!.automatic!.bolusUnits! XCTAssertGreaterThan(standard, 0) - XCTAssertEqual(licensed, min(2 * standard, standard / 0.4), accuracy: 0.051) // AF 0.4 -> 0.8, pump-rounded - // Volatile high: no licence. - input.algorithmExperiments = AlgorithmExperiments(calmHighLicence: .init(maxVolatility: 0.5)) - XCTAssertEqual(LoopAlgorithm.run(input: input).recommendation!.automatic!.bolusUnits!, standard, accuracy: 1e-9) - } - - // MARK: - Negative insulin damper - - func testNegativeInsulinDamperDampsRisesOnly() { - let now = Date() - // Three hours of suspended basal (scheduled 1 U/h) ending 20 min ago: negative insulin pushing glucose up. - let suspend = BasalRelativeDose(type: .basal(scheduledRate: 1.0), startDate: now.addingTimeInterval(.hours(-3)), - endDate: now.addingTimeInterval(.minutes(-20)), volume: 0) - let values: [Double] = [120, 130, 140, 135, 145] - let glucose = values.enumerated().map { PredictedGlucoseValue(startDate: now.addingTimeInterval(.minutes(5 * Double($0.offset))), quantity: .glucose($0.element)) } - let basal = [AbsoluteScheduleValue(startDate: now.addingTimeInterval(.hours(-10)), endDate: now.addingTimeInterval(.hours(6)), value: 1.0)] - let isf = [AbsoluteScheduleValue(startDate: now.addingTimeInterval(.hours(-10)), endDate: now.addingTimeInterval(.hours(6)), value: LoopQuantity.glucose(50))] - let damped = LoopAlgorithm.negativeInsulinDamped(glucose, damper: .init(), dosesRelativeToBasal: [suspend], basal: basal, - sensitivity: isf, insulinModel: ExponentialInsulinModelPreset.rapidActingAdult) - let d = damped.map { $0.quantity.doubleValue(for: unit) } - XCTAssertEqual(d[0], 120, accuracy: 1e-9) - XCTAssertLessThan(d[1] - d[0], 10) // rise damped - XCTAssertEqual(d[3] - d[2], -5, accuracy: 1e-9) // fall kept whole - // No negative insulin: unchanged. - let none = LoopAlgorithm.negativeInsulinDamped(glucose, damper: .init(), dosesRelativeToBasal: [], basal: basal, - sensitivity: isf, insulinModel: ExponentialInsulinModelPreset.rapidActingAdult) - XCTAssertEqual(none.map { $0.quantity.doubleValue(for: unit) }, values) + // 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 testExperimentsRoundTrip() throws { - let stack = AlgorithmExperiments.stack(medianVolatility: 5.0) - let data = try JSONEncoder().encode(stack) - XCTAssertEqual(try JSONDecoder().decode(AlgorithmExperiments.self, from: data), stack) - XCTAssertEqual(stack.volatilityBand?.gain, 1.0) // min(1, 6.7 / 5) - XCTAssertEqual(AlgorithmExperiments.stack(medianVolatility: 10).volatilityBand!.gain, 0.67, accuracy: 1e-12) + 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]) } }