Skills › Scientific Agent Skills › neurokit2 › ecg_hrv_pipeline.py
ecg_hrv_pipeline.py
322 lines · skills/neurokit2/scripts/ecg_hrv_pipeline.py at 5f2bc03
1#!/usr/bin/env python3
2"""Run a bounded ECG processing and duration-aware HRV workflow."""
4from __future__ import annotations
6import argparse
7from typing import Any
9from _common import (
10 MAX_CSV_BYTES,
11 MAX_ROWS,
12 NEUROKIT2_VERSION,
13 PINNED_INSTALL,
14 CliError,
15 checked_input_file,
16 dataframe_rows,
17 emit_json,
18 finite_float,
19 json_number,
20 parse_name_list,
21 read_numeric_columns,
22 require_deidentified,
23 run_cli,
24 sha256_bytes,
25 write_csv,
26)
28ECG_METHODS = (
29 "neurokit",
30 "pantompkins1985",
31 "hamilton2002",
32 "elgendi2010",
33 "engzeemod2012",
34)
35HRV_DOMAINS = {"time", "frequency", "nonlinear"}
38def _row_to_json(frame: Any) -> dict[str, float | int | None]:
39 if len(frame) != 1:
40 raise CliError("expected one-row HRV output")
41 return {str(column): json_number(frame.iloc[0][column]) for column in frame.columns}
44def _summary(values: Any) -> dict[str, float | int | None]:
45 finite = values[values.notna()] if hasattr(values, "notna") else values
46 if len(finite) == 0:
47 return {"count": 0, "maximum": None, "mean": None, "minimum": None}
48 return {
49 "count": len(finite),
50 "maximum": json_number(finite.max()),
51 "mean": json_number(finite.mean()),
52 "minimum": json_number(finite.min()),
53 }
56def _artifact_counts(info: dict[str, Any]) -> dict[str, int]:
57 result: dict[str, int] = {}
58 for name in ("ectopic", "missed", "extra", "longshort"):
59 value = info.get(f"ECG_fixpeaks_{name}", [])
60 try:
61 result[name] = len(value)
62 except TypeError:
63 result[name] = 0
64 return result
67def build_parser() -> argparse.ArgumentParser:
68 parser = argparse.ArgumentParser(
69 description=(
70 "Process a bounded deidentified ECG CSV or a reproducible synthetic "
71 "signal with NeuroKit2 0.2.13, report peak correction and quality, "
72 "and gate HRV domains by duration and beat count."
73 )
74 )
75 source = parser.add_mutually_exclusive_group(required=True)
76 source.add_argument("--input", help="local .csv input")
77 source.add_argument("--synthetic", action="store_true")
78 parser.add_argument("--column", default="ECG")
79 parser.add_argument("--sampling-rate", type=float, required=True, help="Hz")
80 parser.add_argument("--method", choices=ECG_METHODS, default="neurokit")
81 parser.add_argument(
82 "--domains",
83 default="time",
84 help="comma-separated subset of time,frequency,nonlinear",
85 )
86 parser.add_argument(
87 "--duration", type=float, default=60.0, help="synthetic seconds"
88 )
89 parser.add_argument("--heart-rate", type=float, default=70.0, help="synthetic BPM")
90 parser.add_argument("--seed", type=int, default=42)
91 parser.add_argument("--root", default=".")
92 parser.add_argument("--max-rows", type=int, default=MAX_ROWS)
93 parser.add_argument(
94 "--deidentified",
95 action="store_true",
96 help="required with --input",
97 )
98 parser.add_argument("--signals-output", help="optional processed local .csv")
99 parser.add_argument("--output", help="optional local .json report")
100 parser.add_argument("--force", action="store_true")
101 return parser
104def main() -> None:
105 args = build_parser().parse_args()
106 sampling_rate = finite_float(
107 args.sampling_rate,
108 name="--sampling-rate",
109 minimum=20.0,
110 maximum=20_000.0,
111 )
112 domains = set(parse_name_list(args.domains, name="--domains"))
113 unknown = sorted(domains - HRV_DOMAINS)
114 if unknown:
115 raise CliError(f"unknown HRV domains: {', '.join(unknown)}")
116 if not domains:
117 raise CliError("select at least one HRV domain")
118 if not 1 <= args.max_rows <= MAX_ROWS:
119 raise CliError(f"--max-rows must be between 1 and {MAX_ROWS}")
120 if not 0 <= args.seed < 2**32:
121 raise CliError("--seed must be an unsigned 32-bit integer")
123 import neurokit2 as nk
124 import numpy as np
126 installed = str(nk.__version__)
127 if installed != NEUROKIT2_VERSION:
128 raise CliError(
129 f"this reproducible workflow requires NeuroKit2 {NEUROKIT2_VERSION}; "
130 f"found {installed}. Install with: {PINNED_INSTALL}"
131 )
133 if args.input:
134 require_deidentified(args.deidentified)
135 path = checked_input_file(
136 args.input,
137 root=args.root,
138 suffixes={".csv"},
139 max_bytes=MAX_CSV_BYTES,
140 )
141 selected, _, row_count = read_numeric_columns(
142 path,
143 [args.column],
144 max_rows=args.max_rows,
145 )
146 ecg = np.asarray(selected[args.column], dtype=float)
147 source = "deidentified_local_csv"
148 else:
149 duration = finite_float(
150 args.duration,
151 name="--duration",
152 minimum=5.0,
153 maximum=3600.0,
154 )
155 heart_rate = finite_float(
156 args.heart_rate,
157 name="--heart-rate",
158 minimum=20.0,
159 maximum=240.0,
160 )
161 if not sampling_rate.is_integer():
162 raise CliError(
163 "--sampling-rate must be an integer for NeuroKit2 synthetic ECG"
164 )
165 row_count = round(duration * sampling_rate)
166 if not 1 <= row_count <= args.max_rows:
167 raise CliError(f"synthetic row count must be between 1 and {args.max_rows}")
168 try:
169 ecg = nk.ecg_simulate(
170 duration=row_count / sampling_rate,
171 length=row_count,
172 sampling_rate=int(sampling_rate),
173 heart_rate=heart_rate,
174 random_state=args.seed,
175 )
176 except (TypeError, ValueError, RuntimeError) as exc:
177 raise CliError(f"synthetic ECG generation failed: {exc}") from exc
178 row_count = len(ecg)
179 source = "neurokit2_synthetic"
181 duration_s = len(ecg) / sampling_rate
182 if duration_s < 5:
183 raise CliError("ECG must contain at least 5 seconds for this pipeline")
184 try:
185 signals, info = nk.ecg_process(
186 ecg,
187 sampling_rate=sampling_rate,
188 method=args.method,
189 )
190 except (TypeError, ValueError, RuntimeError, IndexError) as exc:
191 raise CliError(f"ECG processing failed: {exc}") from exc
192 peaks = info.get("ECG_R_Peaks")
193 if peaks is None or len(peaks) < 3:
194 raise CliError("fewer than three R-peaks were detected")
195 beat_count = len(peaks)
196 peak_span_s = float(peaks[-1] - peaks[0]) / sampling_rate
198 warnings: list[str] = []
199 hrv: dict[str, dict[str, float | int | None]] = {}
200 if "time" in domains:
201 if beat_count < 20:
202 warnings.append("time-domain HRV skipped: fewer than 20 detected beats")
203 else:
204 try:
205 hrv["time"] = _row_to_json(
206 nk.hrv_time(info, sampling_rate=sampling_rate)
207 )
208 except (TypeError, ValueError, RuntimeError, ZeroDivisionError) as exc:
209 warnings.append(f"time-domain HRV failed: {exc}")
210 if peak_span_s < 300:
211 warnings.append(
212 "time-domain HRV is shorter than the conventional 5-minute "
213 "short-term recording; metric-specific validation is required"
214 )
215 if "frequency" in domains:
216 if peak_span_s < 120 or beat_count < 50:
217 warnings.append(
218 "frequency-domain HRV skipped: require at least 120 seconds between "
219 "first and last corrected peaks and "
220 "50 detected beats in this conservative CLI"
221 )
222 else:
223 try:
224 hrv["frequency"] = _row_to_json(
225 nk.hrv_frequency(
226 info, sampling_rate=sampling_rate,
227 psd_method="welch", normalize=False, interpolation_rate=4,
228 )
229 )
230 except (TypeError, ValueError, RuntimeError, ZeroDivisionError) as exc:
231 warnings.append(f"frequency-domain HRV failed: {exc}")
232 warnings.append(
233 "ULF is not validated by this bounded pipeline; VLF interpretation "
234 "from short records is uncertain. Frequency band powers use "
235 "normalize=False (ms^2), not maximum-PSD normalization."
236 )
237 if peak_span_s < 300:
238 warnings.append(
239 "frequency-domain HRV is below the conventional 5-minute "
240 "short-term window; do not interpret VLF/ULF and justify bands"
241 )
242 if "nonlinear" in domains:
243 if beat_count < 100:
244 warnings.append(
245 "nonlinear HRV skipped: fewer than 100 detected beats; many "
246 "entropy/fractal metrics need substantially more"
247 )
248 else:
249 try:
250 hrv["nonlinear"] = _row_to_json(
251 nk.hrv_nonlinear(info, sampling_rate=sampling_rate)
252 )
253 except (TypeError, ValueError, RuntimeError, ZeroDivisionError) as exc:
254 warnings.append(f"nonlinear HRV failed: {exc}")
256 artifact_counts = _artifact_counts(info)
257 corrected_total = sum(artifact_counts.values())
258 if corrected_total:
259 warnings.append(
260 "R-peak corrections were applied by ecg_process(); inspect the raw "
261 "trace, uncorrected peaks, and correction categories"
262 )
263 if source == "deidentified_local_csv":
264 warnings.append(
265 "input amplitude units are not inferred; record sensor units, lead, "
266 "gain, filters, and clock metadata separately"
267 )
269 report: dict[str, Any] = {
270 "artifact_correction": {
271 "algorithm": "Lipponen-Tarvainen via ecg_process(correct_artifacts=True)",
272 "category_counts": artifact_counts,
273 "corrected_event_count_sum": corrected_total,
274 "uncorrected_peak_count": len(info.get("ECG_R_Peaks_Uncorrected", peaks)),
275 },
276 "detected_r_peaks": beat_count,
277 "duration_s": duration_s,
278 "hrv_peak_span_s": peak_span_s,
279 "hrv_interval_validity": "not established by automatic peak correction",
280 "hrv": hrv,
281 "hrv_frequency_settings": {
282 "psd_method": "welch",
283 "normalize": False,
284 "interpolation_rate_hz": 4,
285 } if "frequency" in domains else None,
286 "hrv_requested_domains": sorted(domains),
287 "ecg_method": args.method,
288 "neurokit2_version": installed,
289 "output_schema_observed": {
290 "info_keys": sorted(str(key) for key in info),
291 "signal_columns": [str(column) for column in signals.columns],
292 },
293 "path_redacted": True,
294 "quality": _summary(signals["ECG_Quality"]),
295 "research_use_only": True,
296 "row_count": row_count,
297 "sampling_rate_hz": sampling_rate,
298 "source": source,
299 "warning": (
300 "This output is for research and education, not diagnosis, patient "
301 "monitoring, or medical-device validation. Validate the detector, "
302 "sensor, population, protocol, and artifact policy for the intended use."
303 ),
304 "warnings": warnings,
305 }
306 if args.signals_output:
307 payload = write_csv(
308 args.signals_output,
309 [str(column) for column in signals.columns],
310 dataframe_rows(signals),
311 root=args.root,
312 force=args.force,
313 )
314 report["signals_csv"] = {
315 "path_redacted": True,
316 "sha256": sha256_bytes(payload),
317 }
318 emit_json(report, output=args.output, root=args.root, force=args.force)
321if __name__ == "__main__":
322 raise SystemExit(run_cli(main))