#!/usr/bin/env python3 """ Kelly-optimal allocation between exercise (throughput therapy) and fasting/GNG (hepatic therapy) using Ben's CGM + meal + exercise data. Model: max E[log(improvement)] = a·log(1 + t_e) + b·log(1 + t_g) subject to t_e + t_g + t_eat = 16h (waking hours) We estimate a and b empirically from the data, then solve for optimal allocation. """ import csv import re from datetime import datetime, timedelta from collections import defaultdict import math # ============================================================ # 1. Parse CGM data # ============================================================ cgm = {} # datetime -> mg/dL with open('/home/ben/downloads/Blood Glucose-2026-02-19-2026-03-17.csv') as f: reader = csv.DictReader(f) for row in reader: dt = datetime.strptime(row['Date/Time'], '%Y-%m-%d %H:%M:%S') bg = float(row['Blood Glucose (mg/dL)']) cgm[dt] = bg print(f"CGM readings: {len(cgm)}") # ============================================================ # 2. Parse meal log # ============================================================ meals = [] # (date, time, description) with open('/home/ben/ava/ben/health/meal-log.md') as f: content = f.read() current_date = None for line in content.strip().split('\n'): line = line.strip() if not line: continue # Match date headers like "# March 14" or "# Feb 28" date_match = re.match(r'^#\s+(March|Feb|February)\s+(\d+)', line) if date_match: month_str, day = date_match.groups() month = 3 if month_str == 'March' else 2 current_date = datetime(2026, month, int(day)).date() continue # Match meal entries like "14:25 description" meal_match = re.match(r'^(\d{1,2}):(\d{2})\s+(.*)', line) if meal_match and current_date: h, m, desc = meal_match.groups() meal_time = datetime.combine(current_date, datetime.min.time().replace(hour=int(h), minute=int(m))) meals.append((current_date, meal_time, desc.strip())) meals.sort(key=lambda x: x[1]) print(f"Meals: {len(meals)}") # ============================================================ # 3. Exercise data (from intervals.icu, previously extracted) # ============================================================ # Manually entered from our previous analysis exercise_sessions = [ # (date, type, duration_min, calories, start_time_approx) ('2026-02-20', 'VirtualRide', 35, 197, '06:00'), ('2026-02-22', 'VirtualRide', 35, 210, '06:00'), ('2026-03-01', 'VirtualRide', 29, 149, '06:00'), ('2026-03-07', 'VirtualRide', 42, 245, '06:00'), ('2026-03-08', 'VirtualRide', 33, 186, '06:00'), ('2026-03-09', 'VirtualRide', 40, 240, '06:00'), ('2026-03-09', 'WeightTraining', 21, 241, '06:45'), ('2026-03-14', 'VirtualRide', 60, 449, '06:00'), ('2026-03-15', 'Walk', 20, 50, '10:00'), ('2026-03-16', 'VirtualRide', 35, 200, '06:00'), ] exercise_by_date = defaultdict(lambda: {'duration_min': 0, 'calories': 0}) for date_str, etype, dur, cal, start in exercise_sessions: d = datetime.strptime(date_str, '%Y-%m-%d').date() exercise_by_date[d]['duration_min'] += dur exercise_by_date[d]['calories'] += cal # ============================================================ # 4. Build daily summaries # ============================================================ # For each day, compute: # - Exercise hours (t_e) # - Eating window (t_eat) # - GNG hours estimate (t_g) = 16 - t_e - t_eat (approx) # - TTGNG: ~11-12h after last meal (from our analysis) # - Overnight glucose (4-6am next day) as outcome metric # - Daily mean glucose as outcome metric # - Time in range start_date = datetime(2026, 2, 19).date() end_date = datetime(2026, 3, 17).date() daily_data = [] for day_offset in range((end_date - start_date).days + 1): d = start_date + timedelta(days=day_offset) # Get CGM for this day day_readings = [(dt, bg) for dt, bg in cgm.items() if dt.date() == d] if len(day_readings) < 100: # need decent coverage continue day_readings.sort() bgs = [bg for _, bg in day_readings] mean_bg = sum(bgs) / len(bgs) tir = sum(1 for bg in bgs if 70 <= bg <= 140) / len(bgs) * 100 time_above_140 = sum(1 for bg in bgs if bg > 140) / len(bgs) * 100 # Overnight glucose (4-6am) overnight_readings = [bg for dt, bg in day_readings if 4 <= dt.hour < 6] overnight_bg = sum(overnight_readings) / len(overnight_readings) if overnight_readings else None # Next day overnight (the overnight AFTER this day's meals) next_d = d + timedelta(days=1) next_overnight = [(dt, bg) for dt, bg in cgm.items() if dt.date() == next_d and 4 <= dt.hour < 6] next_overnight_bg = sum(bg for _, bg in next_overnight) / len(next_overnight) if next_overnight else None # Exercise ex = exercise_by_date.get(d, {'duration_min': 0, 'calories': 0}) t_e = ex['duration_min'] / 60 # hours # Meals today day_meals = [(mt, desc) for dm, mt, desc in meals if dm == d] if day_meals: first_meal = min(mt for mt, _ in day_meals) last_meal = max(mt for mt, _ in day_meals) eating_window = (last_meal - first_meal).total_seconds() / 3600 t_eat = eating_window + 1 # add ~1h for digestion after last meal # Estimate GNG hours: from ~12h after last meal to next first meal # Or for a given day: hours since midnight that were in GNG state # Better: hours in the 24h day spent in GNG # GNG starts ~11-12h after last meal of PREVIOUS day # Let's compute fasting duration from last meal previous day prev_d = d - timedelta(days=1) prev_meals = [(mt, desc) for dm, mt, desc in meals if dm == prev_d] if prev_meals: prev_last_meal = max(mt for mt, _ in prev_meals) fast_duration = (first_meal - prev_last_meal).total_seconds() / 3600 ttgng = 11.5 # our estimate gng_hours = max(0, fast_duration - ttgng) else: # No previous day meal data - might be a multi-day fast fast_duration = None gng_hours = None n_meals = len(day_meals) else: # No meals = fasting day t_eat = 0 eating_window = 0 fast_duration = 24 # at least gng_hours = max(0, 24 - 11.5) # ~12.5h of GNG n_meals = 0 first_meal = None last_meal = None # Categorize meals as protein-dominant or carb-containing carb_keywords = ['rice', 'bread', 'toast', 'pasta', 'spaghetti', 'pancake', 'burrito', 'cookie', 'cake', 'ice cream', 'croissant', 'noodle', 'corn', 'potato', 'sweet potato', 'mamaliga', 'fig', 'fruit', 'banana', 'apple', 'mango', 'honey', 'sugar', 'syrup', 'jam', 'pita', 'roll', 'cornbread', 'quinoa', 'oat'] protein_keywords = ['steak', 'chicken', 'beef', 'salmon', 'shrimp', 'fish', 'egg', 'yogurt', 'protein', 'sardine', 'pork', 'ham', 'burger', 'gyro', 'tilapia', 'meatball'] carb_meals = 0 protein_meals = 0 for mt, desc in day_meals: desc_lower = desc.lower() has_carb = any(kw in desc_lower for kw in carb_keywords) has_protein = any(kw in desc_lower for kw in protein_keywords) if has_carb: carb_meals += 1 if has_protein and not has_carb: protein_meals += 1 daily_data.append({ 'date': d, 'mean_bg': mean_bg, 'tir': tir, 'pct_above_140': time_above_140, 'overnight_bg': overnight_bg, 'next_overnight_bg': next_overnight_bg, 't_exercise_h': t_e, 'exercise_cal': ex['calories'], 't_eat_h': t_eat if day_meals else 0, 'eating_window_h': eating_window if day_meals else 0, 'fast_duration_h': fast_duration, 'gng_hours': gng_hours, 'n_meals': n_meals, 'carb_meals': carb_meals, 'protein_meals': protein_meals, }) print(f"\nDays with full data: {len(daily_data)}") # ============================================================ # 5. Estimate potency coefficients # ============================================================ # We want to find a and b in: improvement ∝ a·log(1 + t_e) + b·log(1 + t_g) # # "Improvement" proxy: lower mean BG, lower overnight BG, higher TIR # We'll use daily mean BG as the outcome (lower = better) # # Actually, let's think about this differently. # The compounding model says each day's "therapy dose" contributes to # long-term improvement. We need to estimate: # - How much does 1h of exercise improve next-day metrics? # - How much does 1h of GNG improve next-day metrics? print("\n" + "="*70) print("DAILY DATA SUMMARY") print("="*70) print(f"{'Date':>10} {'Mean':>5} {'ON_BG':>5} {'TIR':>5} {'Ex_h':>5} {'Eat_h':>5} {'Fast_h':>6} {'GNG_h':>5} {'Meals':>5} {'Carb':>4} {'Prot':>4}") for dd in daily_data: gng = f"{dd['gng_hours']:.1f}" if dd['gng_hours'] is not None else "?" fast = f"{dd['fast_duration_h']:.1f}" if dd['fast_duration_h'] is not None else "?" on_bg = f"{dd['overnight_bg']:.0f}" if dd['overnight_bg'] is not None else "?" print(f"{str(dd['date']):>10} {dd['mean_bg']:>5.1f} {on_bg:>5} {dd['tir']:>5.1f} {dd['t_exercise_h']:>5.1f} {dd['t_eat_h']:>5.1f} {fast:>6} {gng:>5} {dd['n_meals']:>5} {dd['carb_meals']:>4} {dd['protein_meals']:>4}") # ============================================================ # 6. Simple regression: next-day overnight BG ~ exercise + GNG + carbs # ============================================================ print("\n" + "="*70) print("REGRESSION: What predicts better next-day outcomes?") print("="*70) # Filter to days where we have all the data complete = [d for d in daily_data if d['gng_hours'] is not None and d['next_overnight_bg'] is not None and d['fast_duration_h'] is not None] print(f"\nComplete days for regression: {len(complete)}") if len(complete) >= 5: # Manual OLS for: next_overnight_bg = c0 + c1*t_exercise + c2*gng_hours + c3*carb_meals + c4*time_trend # Using numpy-free approach n = len(complete) # Add time trend (day index) to control for overall improvement base_date = complete[0]['date'] # Prepare variables Y = [d['next_overnight_bg'] for d in complete] X_exercise = [d['t_exercise_h'] for d in complete] X_gng = [d['gng_hours'] for d in complete] X_carb = [d['carb_meals'] for d in complete] X_protein = [d['protein_meals'] for d in complete] X_trend = [(d['date'] - base_date).days for d in complete] # Correlations def corr(a, b): n = len(a) mean_a = sum(a)/n mean_b = sum(b)/n cov = sum((a[i]-mean_a)*(b[i]-mean_b) for i in range(n)) std_a = (sum((x-mean_a)**2 for x in a))**0.5 std_b = (sum((x-mean_b)**2 for x in b))**0.5 if std_a == 0 or std_b == 0: return 0 return cov / (std_a * std_b) print(f"\nCorrelations with NEXT-DAY overnight BG (lower = better):") print(f" Exercise hours: r = {corr(X_exercise, Y):+.3f}") print(f" GNG hours: r = {corr(X_gng, Y):+.3f}") print(f" Carb meals: r = {corr(X_carb, Y):+.3f}") print(f" Protein meals: r = {corr(X_protein, Y):+.3f}") print(f" Time trend: r = {corr(X_trend, Y):+.3f}") # Same for same-day mean BG Y_mean = [d['mean_bg'] for d in complete] print(f"\nCorrelations with SAME-DAY mean BG:") print(f" Exercise hours: r = {corr(X_exercise, Y_mean):+.3f}") print(f" GNG hours: r = {corr(X_gng, Y_mean):+.3f}") print(f" Carb meals: r = {corr(X_carb, Y_mean):+.3f}") print(f" Protein meals: r = {corr(X_protein, Y_mean):+.3f}") print(f" Time trend: r = {corr(X_trend, Y_mean):+.3f}") # TIR Y_tir = [d['tir'] for d in complete] print(f"\nCorrelations with SAME-DAY TIR (higher = better):") print(f" Exercise hours: r = {corr(X_exercise, Y_tir):+.3f}") print(f" GNG hours: r = {corr(X_gng, Y_tir):+.3f}") print(f" Carb meals: r = {corr(X_carb, Y_tir):+.3f}") print(f" Protein meals: r = {corr(X_protein, Y_tir):+.3f}") print(f" Time trend: r = {corr(X_trend, Y_tir):+.3f}") # ============================================================ # 7. Kelly optimal allocation # ============================================================ print("\n" + "="*70) print("KELLY-OPTIMAL ALLOCATION MODEL") print("="*70) # Group days by therapy mode exercise_days = [d for d in complete if d['t_exercise_h'] > 0] rest_days = [d for d in complete if d['t_exercise_h'] == 0] high_gng = [d for d in complete if d['gng_hours'] and d['gng_hours'] > 5] low_gng = [d for d in complete if d['gng_hours'] is not None and d['gng_hours'] <= 5] # Dose-response for exercise print(f"\nExercise dose-response (next-day overnight BG):") if exercise_days: ex_next_on = [d['next_overnight_bg'] for d in exercise_days] print(f" Exercise days (n={len(exercise_days)}): next overnight = {sum(ex_next_on)/len(ex_next_on):.1f}") if rest_days: rest_next_on = [d['next_overnight_bg'] for d in rest_days] print(f" Rest days (n={len(rest_days)}): next overnight = {sum(rest_next_on)/len(rest_next_on):.1f}") print(f"\nGNG dose-response (next-day overnight BG):") if high_gng: high_next_on = [d['next_overnight_bg'] for d in high_gng] print(f" High GNG >5h (n={len(high_gng)}): next overnight = {sum(high_next_on)/len(high_next_on):.1f}") if low_gng: low_next_on = [d['next_overnight_bg'] for d in low_gng] print(f" Low GNG ≤5h (n={len(low_gng)}): next overnight = {sum(low_next_on)/len(low_next_on):.1f}") # ============================================================ # 8. Estimate a and b using marginal effects # ============================================================ # For each day, compute the "improvement" = base_overnight - actual_overnight # where base = overall average overnight # Then regress improvement on log(1 + t_e) and log(1 + t_g) print("\n" + "="*70) print("POTENCY ESTIMATION") print("="*70) # Use a rolling approach: for each day, what's the improvement in # the WEEKLY metrics (not just next day)? # Better: compute 3-day rolling average improvement # Group into ~4-day windows windows = [] for i in range(0, len(complete) - 2, 3): chunk = complete[i:i+3] avg_exercise = sum(d['t_exercise_h'] for d in chunk) / len(chunk) avg_gng = sum(d['gng_hours'] for d in chunk) / len(chunk) avg_next_on = sum(d['next_overnight_bg'] for d in chunk) / len(chunk) avg_mean = sum(d['mean_bg'] for d in chunk) / len(chunk) avg_carb = sum(d['carb_meals'] for d in chunk) / len(chunk) mid_date = chunk[len(chunk)//2]['date'] windows.append({ 'date': mid_date, 'avg_exercise': avg_exercise, 'avg_gng': avg_gng, 'avg_next_overnight': avg_next_on, 'avg_mean_bg': avg_mean, 'avg_carb': avg_carb, }) print(f"\n3-day windows: {len(windows)}") for w in windows: print(f" {w['date']}: ex={w['avg_exercise']:.2f}h gng={w['avg_gng']:.1f}h next_ON={w['avg_next_overnight']:.1f} mean={w['avg_mean_bg']:.1f} carbs={w['avg_carb']:.1f}") # ============================================================ # 9. Direct estimation via comparison # ============================================================ print("\n" + "="*70) print("DIRECT POTENCY COMPARISON") print("="*70) # Compare days with similar carb intake but different exercise/GNG print("\nControlling for carb meals (0-1 carb meals):") low_carb_days = [d for d in complete if d['carb_meals'] <= 1] if low_carb_days: ex_lc = [d for d in low_carb_days if d['t_exercise_h'] > 0] rest_lc = [d for d in low_carb_days if d['t_exercise_h'] == 0] if ex_lc: print(f" + Exercise (n={len(ex_lc)}): mean {sum(d['mean_bg'] for d in ex_lc)/len(ex_lc):.1f}, next ON {sum(d['next_overnight_bg'] for d in ex_lc)/len(ex_lc):.1f}") if rest_lc: print(f" - Exercise (n={len(rest_lc)}): mean {sum(d['mean_bg'] for d in rest_lc)/len(rest_lc):.1f}, next ON {sum(d['next_overnight_bg'] for d in rest_lc)/len(rest_lc):.1f}") print("\nControlling for exercise (rest days only):") rest_only = [d for d in complete if d['t_exercise_h'] == 0] if rest_only: high_gng_rest = [d for d in rest_only if d['gng_hours'] and d['gng_hours'] > 5] low_gng_rest = [d for d in rest_only if d['gng_hours'] is not None and d['gng_hours'] <= 5] if high_gng_rest: print(f" High GNG (n={len(high_gng_rest)}): mean {sum(d['mean_bg'] for d in high_gng_rest)/len(high_gng_rest):.1f}, next ON {sum(d['next_overnight_bg'] for d in high_gng_rest)/len(high_gng_rest):.1f}") if low_gng_rest: print(f" Low GNG (n={len(low_gng_rest)}): mean {sum(d['mean_bg'] for d in low_gng_rest)/len(low_gng_rest):.1f}, next ON {sum(d['next_overnight_bg'] for d in low_gng_rest)/len(low_gng_rest):.1f}") # ============================================================ # 10. Estimate Kelly fractions # ============================================================ print("\n" + "="*70) print("KELLY FRACTION ESTIMATION") print("="*70) # Method: treat each therapy as a "bet" that produces a "return" # Exercise: 1h of exercise produces X mg/dL improvement # GNG: 1h of GNG produces Y mg/dL improvement # Kelly fraction = edge/odds = marginal_return / variance # Compute marginal returns per hour if complete: # Exercise marginal effect (per hour) ex_days_data = [(d['t_exercise_h'], d['mean_bg'], d['next_overnight_bg']) for d in complete if d['t_exercise_h'] > 0] no_ex_data = [(d['mean_bg'], d['next_overnight_bg']) for d in complete if d['t_exercise_h'] == 0] if ex_days_data and no_ex_data: avg_ex_mean = sum(d[1] for d in ex_days_data) / len(ex_days_data) avg_noex_mean = sum(d[0] for d in no_ex_data) / len(no_ex_data) avg_ex_hours = sum(d[0] for d in ex_days_data) / len(ex_days_data) ex_effect_per_hour = (avg_noex_mean - avg_ex_mean) / avg_ex_hours print(f"\nExercise effect on same-day mean BG: {ex_effect_per_hour:+.1f} mg/dL per hour") avg_ex_on = sum(d[2] for d in ex_days_data) / len(ex_days_data) avg_noex_on = sum(d[1] for d in no_ex_data) / len(no_ex_data) ex_on_effect = (avg_noex_on - avg_ex_on) / avg_ex_hours print(f"Exercise effect on next-day overnight: {ex_on_effect:+.1f} mg/dL per hour") # GNG marginal effect (per hour) gng_data = [(d['gng_hours'], d['mean_bg'], d['next_overnight_bg']) for d in complete if d['gng_hours'] is not None and d['gng_hours'] > 0] if gng_data: # Split at median GNG hours gng_sorted = sorted(gng_data, key=lambda x: x[0]) mid = len(gng_sorted) // 2 low_half = gng_sorted[:mid] high_half = gng_sorted[mid:] avg_low_gng = sum(d[0] for d in low_half) / len(low_half) avg_high_gng = sum(d[0] for d in high_half) / len(high_half) avg_low_mean = sum(d[1] for d in low_half) / len(low_half) avg_high_mean = sum(d[1] for d in high_half) / len(high_half) avg_low_on = sum(d[2] for d in low_half) / len(low_half) avg_high_on = sum(d[2] for d in high_half) / len(high_half) gng_delta_h = avg_high_gng - avg_low_gng gng_effect_mean = (avg_low_mean - avg_high_mean) / gng_delta_h if gng_delta_h > 0 else 0 gng_effect_on = (avg_low_on - avg_high_on) / gng_delta_h if gng_delta_h > 0 else 0 print(f"\nGNG effect on same-day mean BG: {gng_effect_mean:+.1f} mg/dL per hour") print(f"GNG effect on next-day overnight: {gng_effect_on:+.1f} mg/dL per hour") # Potency ratio a/b print(f"\n--- POTENCY COEFFICIENTS ---") # Using next-day overnight as the compounding metric a = abs(ex_on_effect) if ex_on_effect else 0 b = abs(gng_effect_on) if gng_effect_on else 0 if a + b > 0: print(f" a (exercise potency): {a:.2f} mg/dL improvement per hour") print(f" b (GNG potency): {b:.2f} mg/dL improvement per hour") print(f" Ratio a/b: {a/b:.2f}" if b > 0 else " b = 0") # Kelly optimal allocation # max: a·log(1+t_e) + b·log(1+t_g) s.t. t_e + t_g = T # FOC: a/(1+t_e) = b/(1+t_g) # => t_e = (a·T + a - b) / (a + b) # => t_g = (b·T + b - a) / (a + b) T = 16 # waking hours available # But we need to subtract eating time # Average eating window = ~5-8h on eating days # Let's model for different eating windows for t_eat in [0, 2, 4, 6, 8]: T_available = T - t_eat t_e_opt = (a * T_available + a - b) / (a + b) t_g_opt = (b * T_available + b - a) / (a + b) # Clamp to [0, T_available] t_e_opt = max(0, min(T_available, t_e_opt)) t_g_opt = T_available - t_e_opt print(f"\n If eating window = {t_eat}h ({T_available}h available):") print(f" Optimal exercise: {t_e_opt:.1f}h") print(f" Optimal GNG: {t_g_opt:.1f}h") print(f" Exercise fraction: {t_e_opt/T_available*100:.0f}%") # ============================================================ # 11. Week-over-week improvement model # ============================================================ print("\n" + "="*70) print("WEEK-OVER-WEEK COMPOUNDING MODEL") print("="*70) # Group by week and see if exercise-heavy vs GNG-heavy weeks improve more weekly = defaultdict(lambda: {'exercise_h': 0, 'gng_h': 0, 'mean_bgs': [], 'overnight_bgs': [], 'carb_meals': 0, 'days': 0}) for d in daily_data: # ISO week week = d['date'].isocalendar()[1] weekly[week]['exercise_h'] += d['t_exercise_h'] weekly[week]['mean_bgs'].append(d['mean_bg']) if d['overnight_bg']: weekly[week]['overnight_bgs'].append(d['overnight_bg']) if d['gng_hours'] is not None: weekly[week]['gng_h'] += d['gng_hours'] weekly[week]['carb_meals'] += d['carb_meals'] weekly[week]['days'] += 1 weeks_sorted = sorted(weekly.items()) print(f"\nWeek Ex_h GNG_h Mean Overnight Carbs Days") prev_overnight = None for wk, wd in weeks_sorted: mean = sum(wd['mean_bgs'])/len(wd['mean_bgs']) if wd['mean_bgs'] else 0 onight = sum(wd['overnight_bgs'])/len(wd['overnight_bgs']) if wd['overnight_bgs'] else 0 delta = f" (Δ{onight - prev_overnight:+.1f})" if prev_overnight else "" prev_overnight = onight print(f" {wk} {wd['exercise_h']:.1f} {wd['gng_h']:.1f} {mean:.1f} {onight:.1f}{delta} {wd['carb_meals']} {wd['days']}") # ============================================================ # 12. Simulation: project outcomes under different allocations # ============================================================ print("\n" + "="*70) print("SIMULATION: 4-week projections under different strategies") print("="*70) # Current state: overnight BG ~ 114, mean BG ~ 115 # Use estimated per-hour effects to project current_overnight = 114 current_mean = 115 # Assumptions based on data: # - Exercise: a mg/dL improvement per hour, diminishing returns (log) # - GNG: b mg/dL improvement per hour, diminishing returns (log) # - Carb penalty: ~2 mg/dL per carb-heavy meal to overnight # - Each week's improvement compounds print(f"\nCurrent baseline: overnight {current_overnight}, mean {current_mean}") print(f"Exercise potency (a): {a:.2f} mg/dL/h") print(f"GNG potency (b): {b:.2f} mg/dL/h") strategies = [ ("Current protocol (1h ex, ~6h GNG/day, 4 eat days)", 4, 1.0, 6, 3, 2), # (name, eat_days/wk, ex_h/day, gng_h/eat_day, fast_days, carb_meals/eat_day) ("Max exercise (1.5h ex, ~4h GNG, 5 eat days)", 5, 1.5, 4, 2, 2), ("Max fasting (0.5h ex, ~8h GNG, 3 eat days)", 3, 0.5, 8, 4, 1), ("Balanced (1h ex, ~8h GNG, 4 eat days, low carb)", 4, 1.0, 8, 3, 1), ("Your Mar 14 template daily (1h ex, protein-dominant)", 5, 1.0, 6, 2, 0), ] for name, eat_days, ex_h, gng_h_eat, fast_days, carb_per_day in strategies: overnight = current_overnight mean_bg = current_mean for week in range(4): # Weekly therapy dose weekly_ex = eat_days * ex_h weekly_gng = eat_days * gng_h_eat + fast_days * 16 # fast days = ~16h GNG weekly_carb = eat_days * carb_per_day # Improvement (diminishing returns via log) ex_improvement = a * math.log(1 + weekly_ex) gng_improvement = b * math.log(1 + weekly_gng) carb_penalty = 0.3 * weekly_carb # estimated overnight -= (ex_improvement + gng_improvement - carb_penalty) * 0.5 # weekly scaling mean_bg -= (ex_improvement + gng_improvement - carb_penalty) * 0.6 # Floor effects (can't go below healthy baseline) overnight = max(85, overnight) mean_bg = max(90, mean_bg) print(f"\n {name}:") print(f" 4-week projected overnight: {overnight:.0f} (Δ{overnight - current_overnight:+.0f})") print(f" 4-week projected mean BG: {mean_bg:.0f} (Δ{mean_bg - current_mean:+.0f})") print("\n" + "="*70) print("CAVEATS") print("="*70) print(""" 1. n is small (~15 days with complete data). Confidence intervals are wide. 2. Time trend confounds everything — you're improving steadily, so later days (with more exercise) look better partly due to cumulative effect. 3. Food composition is the dominant variable, not exercise or fasting duration. The Kelly model treats exercise and GNG as the two levers, but in reality there's a third lever (food quality) that dominates both. 4. The simulation assumes linear compounding, which almost certainly breaks down as you approach normal glucose levels. 5. Single-day exercise effects are small; the real exercise benefit is cumulative adaptation over weeks/months (GLUT4 density, mitochondrial biogenesis, muscle mass) — which our 4-week dataset can't capture. """)