Skills › Scientific Agent Skills › neurokit2 › generate_synthetic.py
generate_synthetic.py
221 lines · skills/neurokit2/scripts/generate_synthetic.py at 5f2bc03
1#!/usr/bin/env python3
2"""Generate deterministic, dependency-free synthetic biosignal CSV fixtures."""
4from __future__ import annotations
6import argparse
7import math
8import random
9from collections.abc import Sequence
10from typing import Any
12from _common import (
13 MAX_ROWS,
14 CliError,
15 emit_json,
16 finite_float,
17 parse_name_list,
18 run_cli,
19 sha256_bytes,
20 write_csv,
21)
23MODALITIES = ("ecg", "ppg", "rsp", "eda", "emg", "trigger")
24UNITS = {
25 "time_s": "s",
26 "ecg": "arbitrary_unit",
27 "ppg": "arbitrary_unit",
28 "rsp": "arbitrary_unit",
29 "eda": "arbitrary_unit",
30 "emg": "arbitrary_unit",
31 "trigger": "binary",
32}
35def _gaussian(phase: float, center: float, width: float) -> float:
36 distance = min(abs(phase - center), 1.0 - abs(phase - center))
37 return math.exp(-0.5 * (distance / width) ** 2)
40def _event_times(duration: float, interval: float) -> list[float]:
41 first = interval
42 return [
43 first + index * interval
44 for index in range(max(0, int((duration - first) // interval) + 1))
45 if first + index * interval < duration
46 ]
49def generate_rows(
50 *,
51 duration: float,
52 sampling_rate: float,
53 heart_rate: float,
54 respiratory_rate: float,
55 modalities: Sequence[str],
56 seed: int,
57 event_interval: float,
58) -> tuple[list[str], list[list[float]], list[float]]:
59 """Return deterministic analytic fixtures, not validated physiology."""
61 samples = round(duration * sampling_rate)
62 if not 1 <= samples <= MAX_ROWS:
63 raise CliError(f"synthetic row count must be between 1 and {MAX_ROWS}")
64 rng = random.Random(seed)
65 events = _event_times(duration, event_interval)
66 cardiac_period = 60.0 / heart_rate
67 respiratory_period = 60.0 / respiratory_rate
68 columns = ["time_s", *modalities]
69 rows: list[list[float]] = []
71 for sample in range(samples):
72 time_s = sample / sampling_rate
73 cardiac_phase = (time_s % cardiac_period) / cardiac_period
74 respiratory_phase = (time_s % respiratory_period) / respiratory_period
75 values: dict[str, float] = {}
77 if "ecg" in modalities:
78 values["ecg"] = (
79 0.12 * _gaussian(cardiac_phase, 0.18, 0.025)
80 - 0.15 * _gaussian(cardiac_phase, 0.36, 0.012)
81 + 1.00 * _gaussian(cardiac_phase, 0.40, 0.010)
82 - 0.25 * _gaussian(cardiac_phase, 0.43, 0.014)
83 + 0.30 * _gaussian(cardiac_phase, 0.68, 0.055)
84 + rng.gauss(0.0, 0.008)
85 )
86 if "ppg" in modalities:
87 pulse = max(0.0, math.sin(math.pi * cardiac_phase)) ** 2.5
88 notch = 0.12 * _gaussian(cardiac_phase, 0.72, 0.035)
89 values["ppg"] = pulse - notch + rng.gauss(0.0, 0.006)
90 if "rsp" in modalities:
91 values["rsp"] = math.sin(2 * math.pi * respiratory_phase) + rng.gauss(
92 0.0, 0.01
93 )
94 if "eda" in modalities:
95 phasic = 0.0
96 for event_time in events:
97 elapsed = time_s - event_time
98 if elapsed >= 0:
99 phasic += (1.0 - math.exp(-elapsed / 0.45)) * math.exp(
100 -elapsed / 2.2
101 )
102 values["eda"] = (
103 1.0 + 0.002 * time_s + 0.15 * phasic + rng.gauss(0.0, 0.0015)
104 )
105 if "emg" in modalities:
106 active = any(0.5 <= time_s - event_time < 1.2 for event_time in events)
107 scale = 0.16 if active else 0.01
108 values["emg"] = rng.gauss(0.0, scale)
109 if "trigger" in modalities:
110 values["trigger"] = (
111 1.0
112 if any(
113 0 <= time_s - event_time < max(0.02, 2.0 / sampling_rate)
114 for event_time in events
115 )
116 else 0.0
117 )
118 rows.append([time_s, *(values[name] for name in modalities)])
119 return columns, rows, events
122def build_parser() -> argparse.ArgumentParser:
123 parser = argparse.ArgumentParser(
124 description=(
125 "Generate a deterministic local CSV of analytic biosignal-like fixtures. "
126 "The output is for software tests and education, not physiological validation."
127 )
128 )
129 parser.add_argument("--output", required=True, help="local .csv output path")
130 parser.add_argument("--root", default=".", help="existing local I/O boundary")
131 parser.add_argument("--duration", type=float, default=30.0, help="seconds")
132 parser.add_argument("--sampling-rate", type=float, default=250.0, help="Hz")
133 parser.add_argument("--heart-rate", type=float, default=70.0, help="beats/min")
134 parser.add_argument(
135 "--respiratory-rate", type=float, default=15.0, help="breaths/min"
136 )
137 parser.add_argument(
138 "--modalities",
139 default="ecg,rsp,eda,trigger",
140 help=f"comma-separated subset of: {', '.join(MODALITIES)}",
141 )
142 parser.add_argument("--event-interval", type=float, default=5.0, help="seconds")
143 parser.add_argument("--seed", type=int, default=42)
144 parser.add_argument("--report", help="optional local .json report path")
145 parser.add_argument("--force", action="store_true")
146 return parser
149def main() -> None:
150 args = build_parser().parse_args()
151 duration = finite_float(
152 args.duration, name="--duration", minimum=0.1, maximum=3600.0
153 )
154 sampling_rate = finite_float(
155 args.sampling_rate,
156 name="--sampling-rate",
157 minimum=1.0,
158 maximum=5000.0,
159 )
160 heart_rate = finite_float(
161 args.heart_rate, name="--heart-rate", minimum=20.0, maximum=240.0
162 )
163 respiratory_rate = finite_float(
164 args.respiratory_rate,
165 name="--respiratory-rate",
166 minimum=2.0,
167 maximum=80.0,
168 )
169 event_interval = finite_float(
170 args.event_interval,
171 name="--event-interval",
172 minimum=0.5,
173 maximum=3600.0,
174 )
175 modalities = parse_name_list(args.modalities, name="--modalities")
176 unknown = sorted(set(modalities) - set(MODALITIES))
177 if unknown:
178 raise CliError(f"unknown modalities: {', '.join(unknown)}")
179 if not modalities:
180 raise CliError("select at least one modality")
181 if not -(2**31) <= args.seed < 2**31:
182 raise CliError("--seed must be a signed 32-bit integer")
184 columns, rows, event_times = generate_rows(
185 duration=duration,
186 sampling_rate=sampling_rate,
187 heart_rate=heart_rate,
188 respiratory_rate=respiratory_rate,
189 modalities=modalities,
190 seed=args.seed,
191 event_interval=event_interval,
192 )
193 payload = write_csv(
194 args.output,
195 columns,
196 rows,
197 root=args.root,
198 force=args.force,
199 )
200 report: dict[str, Any] = {
201 "artifact": "synthetic_biosignal_fixture",
202 "columns": columns,
203 "duration_s": duration,
204 "event_onsets_s": event_times,
205 "path_redacted": True,
206 "physiological_validation": False,
207 "row_count": len(rows),
208 "sampling_rate_hz": sampling_rate,
209 "seed": args.seed,
210 "sha256": sha256_bytes(payload),
211 "units": {column: UNITS[column] for column in columns},
212 "warning": (
213 "Analytic waveforms are deterministic software fixtures; they do not "
214 "validate a sensor, method, population, diagnosis, or medical device."
215 ),
216 }
217 emit_json(report, output=args.report, root=args.root, force=args.force)
220if __name__ == "__main__":
221 raise SystemExit(run_cli(main))