#!/usr/bin/env python3 """Safety-first multimodal correlation engine for HealthManager. The engine computes exploratory, phase-stratified Spearman associations from complete cases only. Missing symptom days are never interpreted as symptom-free. Only aggregate statistics are persisted; daily values and source names remain in the protected source tables. """ from __future__ import annotations import argparse import hashlib import itertools import json import math import random import re import sqlite3 import sys from collections import defaultdict from dataclasses import dataclass, replace from datetime import date, datetime, timedelta from pathlib import Path from typing import Any, Callable, Iterable from scipy.stats import spearmanr BASE = Path.home() / ".hermes" / "assets" / "Gesundheit" DB = BASE / "health_data.db" MIN_PAIRS = 14 MIN_TARGET_COVERAGE = 0.6 LAGS = range(0, 4) METHOD = "spearman_7d_block_permutation_complete_case_phase_stratified_bh_v1" MAX_BLOCK_PERMUTATIONS = 999 SYMPTOM_DIMENSIONS = { "aphthen_mundulzera", "gi_darm", "muedigkeit_fatigue", "haut", "augen", "gelenke", "vaskulaer_thrombose_warnzeichen", } LAB_FEATURE_ALLOWLIST = { ("c_reaktives_protein_crp", "mg_l"): "lab:crp_mg_l", ("crp", "mg_l"): "lab:crp_mg_l", ("d_dimer", "ug_l"): "lab:d_dimer_ug_l", ("fibrinogen", "g_l"): "lab:fibrinogen_g_l", ("faktor_viii", "percent"): "lab:factor_viii_percent", ("thrombozyten", "tsd_ul"): "lab:platelets_tsd_ul", ("leukozyten", "tsd_ul"): "lab:leukocytes_tsd_ul", ("ferritin", "ng_ml"): "lab:ferritin_ng_ml", } APPLE_METRICS = ( "step_count", "active_energy", "sleep_analysis", "resting_heart_rate", "heart_rate_variability", "blood_oxygen_saturation", "respiratory_rate", "physical_effort", ) ALLOWED_PREDICTORS = ( {"nutrition:histamine_score"} | {f"apple:{metric}" for metric in APPLE_METRICS} | set(LAB_FEATURE_ALLOWLIST.values()) ) ALLOWED_TARGETS = {"symptom_total"} | {f"symptom_{dimension}" for dimension in SYMPTOM_DIMENSIONS} ALLOWED_PHASES = {"unknown", "baseline_pre_treatment", "early_treatment", "stable_treatment"} ALLOWED_STATUSES = {"computed", "insufficient_n", "insufficient_coverage"} @dataclass(frozen=True) class CorrelationResult: predictor: str target: str lag_days: int medication_phase: str n: int eligible_target_days: int expected_target_days: int target_coverage: float missing_pairs: int rho: float | None p_value: float | None q_value: float | None status: str method: str quality_flags: tuple[str, ...] interpretation: str def conn(db_path: Path) -> sqlite3.Connection: connection = sqlite3.connect(db_path) connection.row_factory = sqlite3.Row return connection def _ensure_row_factory(connection: sqlite3.Connection) -> None: connection.row_factory = sqlite3.Row def _float(value: Any) -> float | None: try: if value is None or str(value).strip() == "": return None number = float(str(value).replace(",", ".").strip()) return number if math.isfinite(number) else None except (TypeError, ValueError): return None def _canonical_day(value: Any) -> str | None: text = str(value or "").strip() if re.fullmatch(r"\d{4}-\d{2}-\d{2}", text): try: return date.fromisoformat(text).isoformat() except ValueError: return None if not re.fullmatch( r"\d{4}-\d{2}-\d{2}T\d{2}:\d{2}:\d{2}(?:\.\d{1,6})?(?:Z|[+-]\d{2}:\d{2})?", text, ): return None try: parsed = datetime.fromisoformat(text[:-1] + "+00:00" if text.endswith("Z") else text) except ValueError: return None return parsed.date().isoformat() def _day_add(day: str, lag: int) -> str: return (date.fromisoformat(day) + timedelta(days=lag)).isoformat() def _slug(text: str) -> str: normalized = str(text or "").casefold().translate(str.maketrans({"ä": "ae", "ö": "oe", "ü": "ue", "ß": "ss"})) return re.sub(r"[^a-z0-9]+", "_", normalized).strip("_") def _unit_slug(text: str) -> str: return _slug(str(text or "").replace("µ", "u").replace("μ", "u").replace("%", " percent ")) def _severity(value: Any) -> float | None: text = str(value or "").casefold().strip() exact = { "0": 0.0, "1": 1.0, "2": 2.0, "3": 3.0, "keine (0)": 0.0, "none (0)": 0.0, "leicht (1)": 1.0, "mild (1)": 1.0, "mittel (2)": 2.0, "moderat (2)": 2.0, "schwer (3)": 3.0, "stark (3)": 3.0, "hoch (3)": 3.0, } return exact.get(text) def load_symptom_targets(connection: sqlite3.Connection) -> dict[str, dict[str, float]]: """Load only complete structured quick logs; absence stays missing.""" _ensure_row_factory(connection) by_day: dict[str, dict[str, list[float | None]]] = defaultdict(lambda: defaultdict(list)) for row in connection.execute( """SELECT id,datum,symptom,schwergrad FROM symptom_log WHERE kontext='daily_quick_score' AND datum IS NOT NULL ORDER BY id""" ): dimension = _slug(str(row["symptom"])) if dimension not in SYMPTOM_DIMENSIONS: continue day = str(row["datum"])[:10] # Count every known-dimension source row before parsing. One invalid, # identical or conflicting duplicate invalidates the complete day. by_day[day][dimension].append(_severity(row["schwergrad"])) result: dict[str, dict[str, float]] = { "symptom_total": {}, **{f"symptom_{dimension}": {} for dimension in SYMPTOM_DIMENSIONS}, } for day, dimensions in by_day.items(): if set(dimensions) != SYMPTOM_DIMENSIONS: continue if any(len(values) != 1 or values[0] is None for values in dimensions.values()): continue complete_values = { dimension: float(values[0]) for dimension, values in dimensions.items() if values[0] is not None } result["symptom_total"][day] = sum(complete_values.values()) for dimension, value in complete_values.items(): result[f"symptom_{dimension}"][day] = value return result def load_nutrition_features(connection: sqlite3.Connection) -> dict[str, dict[str, float]]: """Use only the rule-based score whose unit semantics are explicit. Macro and nutrient columns stay excluded until the separate YAZIO nutrient-unit audit is complete. """ _ensure_row_factory(connection) series: dict[str, float] = {} for row in connection.execute( """SELECT datum,histamine_score,item_count,histamine_unknown_count FROM nutrition_daily_summary_v2 WHERE datum IS NOT NULL""" ): if int(row["item_count"] or 0) <= 0 or int(row["histamine_unknown_count"] or 0) != 0: continue value = _float(row["histamine_score"]) if value is not None: series[str(row["datum"])[:10]] = value return {"nutrition:histamine_score": series} def load_apple_features( connection: sqlite3.Connection, points_loader: Callable[[str], dict[str, Any]] | None = None, ) -> dict[str, dict[str, float]]: loader = points_loader if loader is None: _ensure_row_factory(connection) sys.path.insert(0, str(Path(__file__).resolve().parent)) from apple_health_analytics import daily_points placeholders = ",".join("?" for _ in APPLE_METRICS) rows_by_metric: dict[str, list[sqlite3.Row]] = defaultdict(list) rows = connection.execute( f"""SELECT id,metric,value,unit,start_date,end_date,source_name,file_name,file_hash FROM apple_health_records WHERE value IS NOT NULL AND start_date IS NOT NULL AND metric IN ({placeholders}) ORDER BY metric,start_date,end_date,source_name,id""", APPLE_METRICS, ) for row in rows: rows_by_metric[str(row["metric"])].append(row) def from_rows(metric: str) -> dict[str, Any]: return daily_points(metric, _rows=rows_by_metric.get(metric, [])) loader = from_rows features: dict[str, dict[str, float]] = {} for metric in APPLE_METRICS: series: dict[str, float] = {} for day, point in loader(metric).items(): if getattr(point, "quality", None) != "direct": continue value = _float(getattr(point, "value", None)) if value is not None: series[day] = value features[f"apple:{metric}"] = series return features def load_lab_features(connection: sqlite3.Connection) -> dict[str, dict[str, float]]: """Load only exact, strictly provenance-verified same-day lab observations.""" _ensure_row_factory(connection) features: dict[str, dict[str, float]] = defaultdict(dict) rows = connection.execute( """SELECT parameter_name,wert,einheit,abnahme_datum,befund_datum FROM laborwerte WHERE lower(trim(COALESCE(validierungsstatus,'')))='validiert' AND verified_against_original=1 AND reference_range_source='scanned_original' AND trim(COALESCE(einheit,''))<>'' ORDER BY id""" ) candidates: dict[tuple[str, str], list[str]] = defaultdict(list) exact_number = re.compile(r"^[+-]?\d+(?:[.,]\d+)?$") for row in rows: safe_feature = LAB_FEATURE_ALLOWLIST.get( (_slug(str(row["parameter_name"])), _unit_slug(str(row["einheit"]))) ) if not safe_feature: continue # Prefer a valid specimen date; if absent or invalid, fall back to a # valid report date. Rows without either valid ISO day are excluded. day = _canonical_day(row["abnahme_datum"]) or _canonical_day(row["befund_datum"]) if day is None: continue # Count every provenance-valid source row before parsing/censoring. candidates[(safe_feature, day)].append(str(row["wert"] or "").strip()) for (safe_feature, day), raw_values in candidates.items(): if len(raw_values) != 1 or not exact_number.fullmatch(raw_values[0]): continue value = _float(raw_values[0]) if value is not None: features[safe_feature][day] = value return dict(features) def medication_phase_function(connection: sqlite3.Connection) -> Callable[[str], str]: _ensure_row_factory(connection) disruptive = connection.execute( """SELECT COUNT(*) FROM medication_administrations WHERE (lower(COALESCE(medication_name,'')) LIKE '%hyrimoz%' OR lower(COALESCE(medication_name,'')) LIKE '%adalimumab%') AND lower(trim(COALESCE(event_type,''))) IN ('stopped','paused','restart','dose_change','abgesetzt','pausiert','neustart')""" ).fetchone()[0] administration_days = [ day for row in connection.execute( """SELECT datum FROM medication_administrations WHERE (lower(COALESCE(medication_name,'')) LIKE '%hyrimoz%' OR lower(COALESCE(medication_name,'')) LIKE '%adalimumab%') AND lower(trim(COALESCE(event_type,''))) IN ('administered','verabreicht')""" ) if (day := _canonical_day(row["datum"])) is not None ] start = date.fromisoformat(min(administration_days)) if administration_days and not disruptive else None def phase(day: str) -> str: if start is None: return "unknown" current = date.fromisoformat(day[:10]) if current < start: return "baseline_pre_treatment" if current < start + timedelta(days=56): return "early_treatment" return "stable_treatment" return phase def benjamini_hochberg(p_values: Iterable[float]) -> list[float]: values = list(p_values) count = len(values) if not values: return [] order = sorted(range(count), key=values.__getitem__) adjusted = [1.0] * count running = 1.0 for rank_index in range(count - 1, -1, -1): original_index = order[rank_index] rank = rank_index + 1 running = min(running, values[original_index] * count / rank) adjusted[original_index] = min(max(running, values[original_index]), 1.0) return adjusted def _spearman_rho(xs: list[float], ys: list[float]) -> float: return float(spearmanr(xs, ys).statistic) def _calendar_week_blocks(values: list[float], days: list[str]) -> list[list[float]]: if len(values) != len(days): raise ValueError("values and days must have equal length") by_week: dict[str, list[tuple[str, float]]] = defaultdict(list) for day_text, value in zip(days, values, strict=True): day = date.fromisoformat(day_text) week_start = (day - timedelta(days=day.weekday())).isoformat() by_week[week_start].append((day_text, value)) return [ [value for _day, value in sorted(by_week[week_start])] for week_start in sorted(by_week) ] def _block_permutation_p( xs: list[float], ys: list[float], pair_days: list[str], seed_text: str ) -> float: """Two-sided calendar-week block permutation preserving within-week order.""" observed = abs(_spearman_rho(xs, ys)) blocks = _calendar_week_blocks(ys, pair_days) indexes = tuple(range(len(blocks))) if len(blocks) < 2: return 1.0 if len(blocks) <= 7: permutations = list(itertools.permutations(indexes)) else: seed = int(hashlib.sha256(seed_text.encode()).hexdigest()[:16], 16) rng = random.Random(seed) permutations = [indexes] seen = {indexes} while len(permutations) < MAX_BLOCK_PERMUTATIONS: candidate = list(indexes) rng.shuffle(candidate) item = tuple(candidate) if item not in seen: seen.add(item) permutations.append(item) extreme = 0 for permutation in permutations: permuted = [value for block_index in permutation for value in blocks[block_index]] if abs(_spearman_rho(xs, permuted)) >= observed - 1e-12: extreme += 1 return max(extreme, 1) / len(permutations) def _result_for_pairs( predictor: str, target: str, lag: int, phase_name: str, eligible: int, expected: int, coverage: float, xs: list[float], ys: list[float], pair_days: list[str], ) -> CorrelationResult: n = len(xs) flags = ["complete_case", "phase_stratified", "7d_block_permutation", "no_causality"] common = dict( predictor=predictor, target=target, lag_days=lag, medication_phase=phase_name, n=n, eligible_target_days=eligible, expected_target_days=expected, target_coverage=round(coverage, 6), missing_pairs=max(eligible - n, 0), q_value=None, method=METHOD, ) if coverage < MIN_TARGET_COVERAGE: flags.append("low_target_coverage") return CorrelationResult( **common, rho=None, p_value=None, status="insufficient_coverage", quality_flags=tuple(flags), interpretation="Zielabdeckung unter 60 Prozent; keine Kausalität.", ) if n < MIN_PAIRS or len(set(xs)) < 3 or len(set(ys)) < 3: flags.append("low_n_or_variance") return CorrelationResult( **common, rho=None, p_value=None, status="insufficient_n", quality_flags=tuple(flags), interpretation="Unzureichende vollständige Datenpaare oder Varianz; keine Kausalität.", ) rho = _spearman_rho(xs, ys) p_value = _block_permutation_p( xs, ys, pair_days, f"{predictor}|{target}|{lag}|{phase_name}" ) return CorrelationResult( **common, rho=round(rho, 6), p_value=round(p_value, 6), status="computed", quality_flags=tuple(flags), interpretation="Explorative phasenstratifizierte Rangassoziation mit 7-Tage-Blockpermutation; keine Kausalität.", ) def compute_results( connection: sqlite3.Connection, apple_points_loader: Callable[[str], dict[str, Any]] | None = None, ) -> list[CorrelationResult]: targets = load_symptom_targets(connection) features: dict[str, dict[str, float]] = {} features.update(load_nutrition_features(connection)) features.update(load_apple_features(connection, apple_points_loader)) features.update(load_lab_features(connection)) phase_for = medication_phase_function(connection) results: list[CorrelationResult] = [] for target_name, target_series in targets.items(): if not target_series: continue phases = sorted({phase_for(day) for day in target_series}) or ["unknown"] for predictor, predictor_series in sorted(features.items()): if not predictor_series: continue for lag in LAGS: for phase_name in phases: shifted_predictor_days = [_day_add(day, lag) for day in predictor_series] start = max(min(target_series), min(shifted_predictor_days)) end = min(max(target_series), max(shifted_predictor_days)) if start > end: continue expected_days: list[str] = [] current = date.fromisoformat(start) end_date = date.fromisoformat(end) while current <= end_date: outcome_day = current.isoformat() predictor_day = _day_add(outcome_day, -lag) if phase_for(outcome_day) == phase_name and phase_for(predictor_day) == phase_name: expected_days.append(outcome_day) current += timedelta(days=1) if not expected_days: continue eligible_days = [day for day in expected_days if day in target_series] coverage = len(eligible_days) / len(expected_days) xs: list[float] = [] ys: list[float] = [] pair_days: list[str] = [] for outcome_day in eligible_days: predictor_day = _day_add(outcome_day, -lag) if predictor_day not in predictor_series: continue xs.append(predictor_series[predictor_day]) ys.append(target_series[outcome_day]) pair_days.append(outcome_day) results.append( _result_for_pairs( predictor, target_name, lag, phase_name, len(eligible_days), len(expected_days), coverage, xs, ys, pair_days, ) ) grouped: dict[tuple[str, str], list[int]] = defaultdict(list) for index, result in enumerate(results): if result.status == "computed" and result.p_value is not None: grouped[(result.target, result.medication_phase)].append(index) for indexes in grouped.values(): adjusted = benjamini_hochberg([ results[index].p_value if results[index].p_value is not None else 1.0 for index in indexes ]) for index, q_value in zip(indexes, adjusted): signal = q_value <= 0.1 and results[index].n >= 30 interpretation = ( "Exploratives Signal nach Mehrfachtestkorrektur; keine Kausalität." if signal else "Kein belastbares korrigiertes Signal; explorativ und keine Kausalität." ) results[index] = replace(results[index], q_value=round(q_value, 6), interpretation=interpretation) return results def ensure_schema(connection: sqlite3.Connection) -> None: connection.execute( """CREATE TABLE IF NOT EXISTS multimodal_correlation_results ( id INTEGER PRIMARY KEY AUTOINCREMENT, predictor TEXT NOT NULL, target TEXT NOT NULL, lag_days INTEGER NOT NULL, medication_phase TEXT NOT NULL, n INTEGER NOT NULL, eligible_target_days INTEGER NOT NULL, expected_target_days INTEGER NOT NULL DEFAULT 0, target_coverage REAL NOT NULL DEFAULT 0, missing_pairs INTEGER NOT NULL, rho REAL, p_value REAL, q_value REAL, status TEXT NOT NULL, method TEXT NOT NULL, quality_flags TEXT NOT NULL, interpretation TEXT NOT NULL, computed_at TEXT NOT NULL DEFAULT CURRENT_TIMESTAMP, UNIQUE(predictor,target,lag_days,medication_phase) )""" ) columns = {row[1] for row in connection.execute("PRAGMA table_info(multimodal_correlation_results)")} if "expected_target_days" not in columns: connection.execute( "ALTER TABLE multimodal_correlation_results ADD COLUMN expected_target_days INTEGER NOT NULL DEFAULT 0" ) if "target_coverage" not in columns: connection.execute( "ALTER TABLE multimodal_correlation_results ADD COLUMN target_coverage REAL NOT NULL DEFAULT 0" ) connection.execute( "CREATE INDEX IF NOT EXISTS idx_multimodal_correlations_status ON multimodal_correlation_results(status,target,medication_phase)" ) def _validate_result_for_storage(result: CorrelationResult) -> None: if result.predictor not in ALLOWED_PREDICTORS: raise ValueError("predictor is not allowlisted") if result.target not in ALLOWED_TARGETS: raise ValueError("target is not allowlisted") if result.medication_phase not in ALLOWED_PHASES or result.status not in ALLOWED_STATUSES: raise ValueError("phase or status is not allowlisted") if result.lag_days not in LAGS or result.method != METHOD: raise ValueError("lag or method is not allowlisted") counts = ( result.n, result.eligible_target_days, result.expected_target_days, result.missing_pairs, ) if any(isinstance(value, bool) or not isinstance(value, int) or value < 0 for value in counts): raise ValueError("counts must be non-negative integers") if result.expected_target_days < 1 or not ( result.n <= result.eligible_target_days <= result.expected_target_days ): raise ValueError("count invariants are invalid") if result.missing_pairs != result.eligible_target_days - result.n: raise ValueError("missing-pair count is inconsistent") if isinstance(result.target_coverage, bool) or not isinstance(result.target_coverage, (int, float)): raise ValueError("target coverage must be numeric") expected_coverage = result.eligible_target_days / result.expected_target_days if not math.isfinite(float(result.target_coverage)) or not math.isclose( float(result.target_coverage), expected_coverage, rel_tol=0.0, abs_tol=5e-7 ): raise ValueError("target coverage is invalid or inconsistent") def bounded(value: float | None, lower: float, upper: float, name: str) -> None: if value is None: return if isinstance(value, bool) or not isinstance(value, (int, float)): raise ValueError(f"{name} must be numeric") if not math.isfinite(float(value)) or not lower <= float(value) <= upper: raise ValueError(f"{name} is outside its valid range") bounded(result.rho, -1.0, 1.0, "rho") bounded(result.p_value, 0.0, 1.0, "p_value") bounded(result.q_value, 0.0, 1.0, "q_value") base_flags = ( "complete_case", "phase_stratified", "7d_block_permutation", "no_causality" ) if result.status == "insufficient_coverage": expected_flags = base_flags + ("low_target_coverage",) expected_interpretation = "Zielabdeckung unter 60 Prozent; keine Kausalität." if result.target_coverage >= MIN_TARGET_COVERAGE: raise ValueError("coverage status is inconsistent") if any(value is not None for value in (result.rho, result.p_value, result.q_value)): raise ValueError("insufficient-coverage result contains inferential values") elif result.status == "insufficient_n": expected_flags = base_flags + ("low_n_or_variance",) expected_interpretation = ( "Unzureichende vollständige Datenpaare oder Varianz; keine Kausalität." ) if result.target_coverage < MIN_TARGET_COVERAGE: raise ValueError("insufficient-n status is inconsistent with coverage") if any(value is not None for value in (result.rho, result.p_value, result.q_value)): raise ValueError("insufficient-n result contains inferential values") else: expected_flags = base_flags if result.target_coverage < MIN_TARGET_COVERAGE or result.n < MIN_PAIRS: raise ValueError("computed status is inconsistent with evidence thresholds") if any(value is None for value in (result.rho, result.p_value, result.q_value)): raise ValueError("computed result lacks inferential values") signal = result.q_value is not None and result.q_value <= 0.1 and result.n >= 30 expected_interpretation = ( "Exploratives Signal nach Mehrfachtestkorrektur; keine Kausalität." if signal else "Kein belastbares korrigiertes Signal; explorativ und keine Kausalität." ) if result.quality_flags != expected_flags: raise ValueError("quality flags are not generated by the safe engine") if result.interpretation != expected_interpretation: raise ValueError("interpretation is not generated by the safe engine") def store_results(connection: sqlite3.Connection, results: Iterable[CorrelationResult]) -> None: materialized = list(results) for result in materialized: _validate_result_for_storage(result) ensure_schema(connection) connection.execute("DELETE FROM multimodal_correlation_results") connection.executemany( """INSERT INTO multimodal_correlation_results (predictor,target,lag_days,medication_phase,n,eligible_target_days, expected_target_days,target_coverage,missing_pairs,rho,p_value,q_value,status, method,quality_flags,interpretation,computed_at) VALUES(?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,CURRENT_TIMESTAMP)""", [ ( result.predictor, result.target, result.lag_days, result.medication_phase, result.n, result.eligible_target_days, result.expected_target_days, result.target_coverage, result.missing_pairs, result.rho, result.p_value, result.q_value, result.status, result.method, json.dumps(result.quality_flags), result.interpretation, ) for result in materialized ], ) def main(argv: list[str] | None = None) -> int: parser = argparse.ArgumentParser(description=__doc__) parser.add_argument("--db", type=Path, required=True) args = parser.parse_args(argv) if not args.db.exists(): raise SystemExit(f"Database not found: {args.db}") connection = conn(args.db) try: results = compute_results(connection) store_results(connection, results) connection.commit() computed = sum(result.status == "computed" for result in results) insufficient = sum(result.status == "insufficient_n" for result in results) print(json.dumps({"status": "ok", "computed": computed, "insufficient": insufficient})) finally: connection.close() return 0 if __name__ == "__main__": raise SystemExit(main())