Skills › Scientific Agent Skills › neurokit2 › eda_pipeline.py
eda_pipeline.py
301 lines · skills/neurokit2/scripts/eda_pipeline.py at 5f2bc03
1#!/usr/bin/env python3
2"""Run an explicit, bounded EDA cleaning/decomposition/SCR workflow."""
4from __future__ import annotations
6import argparse
7import statistics
8from typing import Any
10from _common import (
11 MAX_CSV_BYTES,
12 MAX_ROWS,
13 NEUROKIT2_VERSION,
14 PINNED_INSTALL,
15 CliError,
16 checked_input_file,
17 dataframe_rows,
18 emit_json,
19 finite_float,
20 json_number,
21 read_numeric_columns,
22 require_deidentified,
23 run_cli,
24 sha256_bytes,
25 write_csv,
26)
28CLEAN_METHODS = ("neurokit", "biosppy", "none")
29PHASIC_METHODS = ("highpass", "smoothmedian", "cvxeda", "sparseda")
30PEAK_METHODS = ("neurokit", "gamboa2008", "kim2004", "vanhalem2020", "nabian2018")
33def _finite_values(values: Any) -> list[float]:
34 result: list[float] = []
35 for value in values:
36 number = json_number(value)
37 if number is not None:
38 result.append(float(number))
39 return result
42def _stats(values: Any) -> dict[str, float | int | None]:
43 finite = _finite_values(values)
44 return {
45 "count": len(finite),
46 "maximum": max(finite) if finite else None,
47 "mean": statistics.fmean(finite) if finite else None,
48 "median": statistics.median(finite) if finite else None,
49 "minimum": min(finite) if finite else None,
50 }
53def build_parser() -> argparse.ArgumentParser:
54 parser = argparse.ArgumentParser(
55 description=(
56 "Clean, explicitly decompose, and detect SCRs in bounded "
57 "deidentified EDA CSV data or a reproducible synthetic signal using "
58 "NeuroKit2 0.2.13."
59 )
60 )
61 source = parser.add_mutually_exclusive_group(required=True)
62 source.add_argument("--input", help="local .csv input")
63 source.add_argument("--synthetic", action="store_true")
64 parser.add_argument("--column", default="EDA")
65 parser.add_argument("--sampling-rate", type=float, required=True, help="Hz")
66 parser.add_argument("--clean-method", choices=CLEAN_METHODS, default="neurokit")
67 parser.add_argument("--phasic-method", choices=PHASIC_METHODS, default="highpass")
68 parser.add_argument("--peak-method", choices=PEAK_METHODS, default="neurokit")
69 parser.add_argument(
70 "--amplitude-min",
71 type=float,
72 default=0.1,
73 help="relative-to-largest SCR threshold for neurokit/kim2004",
74 )
75 parser.add_argument(
76 "--sympathetic-method",
77 choices=("none", "posada", "ghiasi"),
78 default="none",
79 )
80 parser.add_argument("--unit", default="unspecified")
81 parser.add_argument(
82 "--duration", type=float, default=60.0, help="synthetic seconds"
83 )
84 parser.add_argument("--scr-number", type=int, default=6)
85 parser.add_argument("--seed", type=int, default=42)
86 parser.add_argument("--root", default=".")
87 parser.add_argument("--max-rows", type=int, default=MAX_ROWS)
88 parser.add_argument(
89 "--deidentified",
90 action="store_true",
91 help="required with --input",
92 )
93 parser.add_argument("--signals-output", help="optional processed local .csv")
94 parser.add_argument("--output", help="optional local .json report")
95 parser.add_argument("--force", action="store_true")
96 return parser
99def main() -> None:
100 args = build_parser().parse_args()
101 sampling_rate = finite_float(
102 args.sampling_rate,
103 name="--sampling-rate",
104 minimum=1.0,
105 maximum=5000.0,
106 )
107 amplitude_min = finite_float(
108 args.amplitude_min,
109 name="--amplitude-min",
110 minimum=0.0,
111 maximum=1.0,
112 )
113 if not 1 <= args.max_rows <= MAX_ROWS:
114 raise CliError(f"--max-rows must be between 1 and {MAX_ROWS}")
115 if not 0 <= args.scr_number <= 10_000:
116 raise CliError("--scr-number must be between 0 and 10000")
117 if not 0 <= args.seed < 2**32:
118 raise CliError("--seed must be an unsigned 32-bit integer")
119 if not args.unit or len(args.unit) > 32:
120 raise CliError("--unit must be a nonempty string of at most 32 characters")
122 import neurokit2 as nk
123 import numpy as np
124 import pandas as pd
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 raw = 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 if not sampling_rate.is_integer():
156 raise CliError(
157 "--sampling-rate must be an integer for NeuroKit2 synthetic EDA"
158 )
159 row_count = round(duration * sampling_rate)
160 if not 1 <= row_count <= args.max_rows:
161 raise CliError(f"synthetic row count must be between 1 and {args.max_rows}")
162 try:
163 raw = nk.eda_simulate(
164 duration=row_count / sampling_rate,
165 length=row_count,
166 sampling_rate=int(sampling_rate),
167 scr_number=args.scr_number,
168 random_state=args.seed,
169 )
170 except (TypeError, ValueError, RuntimeError) as exc:
171 raise CliError(f"synthetic EDA generation failed: {exc}") from exc
172 row_count = len(raw)
173 source = "neurokit2_synthetic"
175 duration_s = len(raw) / sampling_rate
176 try:
177 clean = nk.eda_clean(
178 raw,
179 sampling_rate=sampling_rate,
180 method=args.clean_method,
181 )
182 components = nk.eda_phasic(
183 clean,
184 sampling_rate=sampling_rate,
185 method=args.phasic_method,
186 )
187 peak_signals, peak_info = nk.eda_peaks(
188 components["EDA_Phasic"],
189 sampling_rate=sampling_rate,
190 method=args.peak_method,
191 amplitude_min=amplitude_min,
192 )
193 except (TypeError, ValueError, RuntimeError, IndexError) as exc:
194 raise CliError(f"EDA processing failed: {exc}") from exc
196 signals = pd.concat(
197 [
198 pd.DataFrame({"EDA_Raw": raw, "EDA_Clean": clean}),
199 components.reset_index(drop=True),
200 peak_signals.reset_index(drop=True),
201 ],
202 axis=1,
203 )
204 amplitudes = _finite_values(peak_info.get("SCR_Amplitude", []))
205 scr_count = len(peak_info.get("SCR_Peaks", []))
206 warnings: list[str] = []
207 threshold_applied = args.peak_method in {"neurokit", "kim2004"}
208 if not threshold_applied:
209 warnings.append(
210 "--amplitude-min is ignored by the selected NeuroKit2 peak detector"
211 )
212 if args.phasic_method == "cvxeda":
213 warnings.append(
214 "cvxEDA requires the optional cvxopt dependency; report its version "
215 "and validate decomposition parameters"
216 )
217 if source == "deidentified_local_csv" and args.unit == "unspecified":
218 warnings.append(
219 "EDA units were not declared; do not interpret amplitudes as "
220 "microsiemens without acquisition metadata"
221 )
222 if scr_count == 0:
223 warnings.append(
224 "no SCR peaks were detected; inspect signal orientation and quality"
225 )
226 sympathetic: dict[str, float | int | None] | None = None
227 if args.sympathetic_method != "none":
228 if duration_s <= 64:
229 warnings.append(
230 "EDA sympathetic index skipped: this CLI requires more than 64 seconds"
231 )
232 else:
233 try:
234 result = nk.eda_sympathetic(
235 clean,
236 sampling_rate=sampling_rate,
237 method=args.sympathetic_method,
238 show=False,
239 )
240 sympathetic = {
241 str(key): json_number(value) for key, value in result.items()
242 }
243 except (TypeError, ValueError, RuntimeError, ZeroDivisionError) as exc:
244 warnings.append(f"EDA sympathetic index failed: {exc}")
246 report: dict[str, Any] = {
247 "duration_s": duration_s,
248 "input_unit": args.unit,
249 "methods": {
250 "clean": args.clean_method,
251 "decomposition": args.phasic_method,
252 "peak_detection": args.peak_method,
253 "peak_threshold_relative_to_largest": (
254 amplitude_min if threshold_applied else None
255 ),
256 "peak_threshold_applied": threshold_applied,
257 "peak_threshold_reference": {
258 "neurokit": "maximum_candidate_prominence",
259 "kim2004": "maximum_candidate_phasic_peak_height",
260 }.get(args.peak_method),
261 "sympathetic": args.sympathetic_method,
262 },
263 "neurokit2_version": installed,
264 "output_schema_observed": {
265 "info_keys": sorted(str(key) for key in peak_info),
266 "signal_columns": [str(column) for column in signals.columns],
267 },
268 "path_redacted": True,
269 "research_use_only": True,
270 "row_count": row_count,
271 "sampling_rate_hz": sampling_rate,
272 "scr": {
273 "amplitude": _stats(amplitudes),
274 "count": scr_count,
275 },
276 "source": source,
277 "sympathetic": sympathetic,
278 "tonic": _stats(signals["EDA_Tonic"]),
279 "warning": (
280 "EDA features are method-, sensor-, site-, unit-, and population-dependent. "
281 "This output is not diagnosis, monitoring, or medical-device validation."
282 ),
283 "warnings": warnings,
284 }
285 if args.signals_output:
286 payload = write_csv(
287 args.signals_output,
288 [str(column) for column in signals.columns],
289 dataframe_rows(signals),
290 root=args.root,
291 force=args.force,
292 )
293 report["signals_csv"] = {
294 "path_redacted": True,
295 "sha256": sha256_bytes(payload),
296 }
297 emit_json(report, output=args.output, root=args.root, force=args.force)
300if __name__ == "__main__":
301 raise SystemExit(run_cli(main))