Examples/mne

examples/mne/README.md

@@ -0,0 +1,18 @@
+# mne Examples
+
+Each sub-directory contains a self-contained example. The order in
+which the examples are to appear is specified in `order.json` (an
+array of directory names in the expected order).
+
+In each example directory you'll find:
+
+* `config.toml` - must conform to the specification outlined here:
+  https://docs.pyscript.net/latest/user-guide/configuration/ This is
+  parsed and ultimately turned into a JSON representation as part of
+  the package's API object.
+* `setup.py` - Python code for contextual and environmental setup,
+  NOT SEEN BY THE END USER, but is run before the `code.py` code is
+  evaluated. Allows us to create useful (IPython) shims, avoid
+  repeating boilerplate and whatnot.
+* `code.py` - the actual code added to the editor which forms the
+  practical example of using the package.

examples/mne/events_and_epochs/code.py

@@ -0,0 +1,98 @@
+# ---------------------------------------------------------------------
+# Cutting continuous data into epochs around stimulus events and
+# averaging them into an evoked response.
+# ---------------------------------------------------------------------
+
+heading("From continuous data to evoked responses")
+note(
+    "Most M/EEG analysis revolves around three structures: Raw "
+    "(continuous), Epochs (segments around events), and Evoked "
+    "(epoch averages). Here we simulate an experiment with two "
+    "conditions and walk through that pipeline."
+)
+
+# Build continuous data with a brief positive deflection 200 ms after
+# each 'target' event and a smaller deflection after each 'standard'.
+sampling_freq = 500.0
+duration_s = 60.0
+n_samples = int(sampling_freq * duration_s)
+
+# One event every ~2 seconds, alternating standard/target.
+event_samples = np.arange(500, n_samples - 500, 1000)
+event_codes = np.where(np.arange(len(event_samples)) % 2 == 0, 1, 2)
+
+def response_kernel(amplitude):
+    # 300 ms Gaussian bump, peak around sample 100 (= 200 ms post event).
+    t = np.arange(150)
+    return amplitude * np.exp(-((t - 100) ** 2) / (2 * 20 ** 2))
+
+signal = rng.normal(0, 0.3, size=n_samples)
+for sample, code in zip(event_samples, event_codes):
+    amp = 2.0 if code == 2 else 0.6  # targets evoke a bigger response
+    end = sample + 150
+    signal[sample:end] += response_kernel(amp)
+
+# In MNE, an events array has shape (n_events, 3): sample, prev id, id.
+events = np.column_stack([
+    event_samples,
+    np.zeros_like(event_samples),
+    event_codes,
+]).astype(int)
+event_id = {"standard": 1, "target": 2}
+
+info = mne.create_info(ch_names=["Cz"], sfreq=sampling_freq, ch_types="eeg")
+raw = mne.io.RawArray(signal[np.newaxis, :] * 1e-6, info)
+
+note(f"Generated {len(events)} events ({np.sum(event_codes == 2)} targets).")
+
+# `Epochs` slices Raw into windows aligned to events. `tmin`/`tmax` are
+# relative to each event (in seconds). `baseline` subtracts the mean of
+# the pre-stimulus interval from each epoch.
+epochs = mne.Epochs(
+    raw,
+    events,
+    event_id=event_id,
+    tmin=-0.1,
+    tmax=0.5,
+    baseline=(-0.1, 0.0),
+    preload=True,
+)
+
+note("Epochs object summary:")
+display(HTML(f"<pre>{epochs}</pre>"), append=True)
+
+# Indexing by condition name returns the matching subset; `.average()`
+# collapses across trials to produce an Evoked response.
+evoked_standard = epochs["standard"].average()
+evoked_target = epochs["target"].average()
+
+heading("Evoked waveforms")
+note(
+    "Averaging across many noisy trials reveals the stereotyped "
+    "response. Targets show a clear peak ~200 ms after the event; "
+    "standards show a smaller one."
+)
+
+fig, ax = plt.subplots(figsize=(9, 4))
+ax.plot(
+    evoked_standard.times * 1000,
+    evoked_standard.data[0] * 1e6,
+    label="Standard",
+    color="steelblue",
+    linewidth=2,
+)
+ax.plot(
+    evoked_target.times * 1000,
+    evoked_target.data[0] * 1e6,
+    label="Target",
+    color="crimson",
+    linewidth=2,
+)
+ax.axvline(0, color="black", linestyle="--", linewidth=0.8)
+ax.axhline(0, color="gray", linewidth=0.5)
+ax.set_xlabel("Time relative to event (ms)")
+ax.set_ylabel("Amplitude (µV)")
+ax.set_title("Average evoked response at Cz")
+ax.legend()
+fig.tight_layout()
+display(fig, append=True)

examples/mne/events_and_epochs/config.toml

@@ -0,0 +1 @@
+packages = ["mne", "numpy", "matplotlib"]

examples/mne/events_and_epochs/setup.py

@@ -0,0 +1,28 @@
+"""Setup for the epoching example."""
+
+import js
+from pyscript import window, HTML, display as _display
+
+js.alert = window.alert
+
+
+def display(*args, **kwargs):
+    return _display(
+        *args, **kwargs, target=__pyscript_display_target__,
+    )
+
+
+def heading(text, level=2):
+    display(HTML(f"<h{level}>{text}</h{level}>"), append=True)
+
+
+def note(text):
+    display(HTML(f"<p>{text}</p>"), append=True)
+
+
+import numpy as np
+import matplotlib.pyplot as plt
+import mne
+
+mne.set_log_level("WARNING")
+rng = np.random.default_rng(7)

examples/mne/filtering_and_psd/code.py

@@ -0,0 +1,75 @@
+# ---------------------------------------------------------------------
+# Filtering a Raw object and inspecting its frequency content.
+# ---------------------------------------------------------------------
+
+heading("Cleaning a noisy signal with band-pass filtering")
+note(
+    "We'll build a single-channel Raw with a 10 Hz alpha rhythm, "
+    "a 50 Hz line-noise contaminant, and slow drift. Then we apply "
+    "a band-pass filter and compare power spectra before and after."
+)
+
+sampling_freq = 500.0
+duration_s = 30.0
+n_samples = int(sampling_freq * duration_s)
+times = np.arange(n_samples) / sampling_freq
+
+alpha = 1.0 * np.sin(2 * np.pi * 10.0 * times)
+line_noise = 0.8 * np.sin(2 * np.pi * 50.0 * times)
+slow_drift = 1.2 * np.sin(2 * np.pi * 0.3 * times)
+white_noise = rng.normal(0, 0.5, size=n_samples)
+
+# Combined signal in microvolt-range, scaled to volts for MNE.
+signal = (alpha + line_noise + slow_drift + white_noise) * 1e-6
+
+info = mne.create_info(ch_names=["Oz"], sfreq=sampling_freq, ch_types="eeg")
+raw = mne.io.RawArray(signal[np.newaxis, :], info)
+
+# `copy()` keeps the original around so we can compare. `filter` is
+# an in-place band-pass between 1 Hz and 40 Hz, removing both slow
+# drift and the 50 Hz line noise.
+raw_filtered = raw.copy().filter(l_freq=1.0, h_freq=40.0)
+
+heading("Comparing the time series")
+note("First two seconds before and after filtering:")
+
+n_show = int(2 * sampling_freq)
+raw_data, raw_times = raw[0, :n_show]
+filt_data, _ = raw_filtered[0, :n_show]
+
+fig, axes = plt.subplots(2, 1, figsize=(9, 5), sharex=True)
+axes[0].plot(raw_times, raw_data[0] * 1e6, color="gray")
+axes[0].set_title("Raw signal (µV)")
+axes[1].plot(raw_times, filt_data[0] * 1e6, color="darkgreen")
+axes[1].set_title("Band-pass 1-40 Hz (µV)")
+axes[1].set_xlabel("Time (s)")
+fig.tight_layout()
+display(fig, append=True)
+
+heading("Power spectral density")
+note(
+    "MNE's `compute_psd` uses Welch's method by default. The alpha "
+    "peak around 10 Hz survives filtering, while the 50 Hz line and "
+    "low-frequency drift are strongly attenuated."
+)
+
+psd_raw = raw.compute_psd(fmin=0.1, fmax=80.0)
+psd_filt = raw_filtered.compute_psd(fmin=0.1, fmax=80.0)
+
+freqs = psd_raw.freqs
+power_raw = psd_raw.get_data()[0]
+power_filt = psd_filt.get_data()[0]
+
+fig, ax = plt.subplots(figsize=(9, 4))
+ax.semilogy(freqs, power_raw, color="gray", label="Raw")
+ax.semilogy(freqs, power_filt, color="darkgreen", label="Filtered")
+ax.axvline(10, color="steelblue", linestyle="--", linewidth=1,
+           label="Alpha (10 Hz)")
+ax.axvline(50, color="crimson", linestyle="--", linewidth=1,
+           label="Line (50 Hz)")
+ax.set_xlabel("Frequency (Hz)")
+ax.set_ylabel("Power (V²/Hz)")
+ax.set_title("Power spectral density")
+ax.legend()
+fig.tight_layout()
+display(fig, append=True)

examples/mne/filtering_and_psd/config.toml

@@ -0,0 +1 @@
+packages = ["mne", "numpy", "matplotlib"]

examples/mne/filtering_and_psd/setup.py

@@ -0,0 +1,29 @@
+"""Setup for the filtering example. Mirrors the names established in
+the first example without re-registering the IPython shim."""
+
+import js
+from pyscript import window, HTML, display as _display
+
+js.alert = window.alert
+
+
+def display(*args, **kwargs):
+    return _display(
+        *args, **kwargs, target=__pyscript_display_target__,
+    )
+
+
+def heading(text, level=2):
+    display(HTML(f"<h{level}>{text}</h{level}>"), append=True)
+
+
+def note(text):
+    display(HTML(f"<p>{text}</p>"), append=True)
+
+
+import numpy as np
+import matplotlib.pyplot as plt
+import mne
+
+mne.set_log_level("WARNING")
+rng = np.random.default_rng(7)

examples/mne/order.json

@@ -0,0 +1,5 @@
+[
+    "raw_eeg_basics",
+    "filtering_and_psd",
+    "events_and_epochs"
+]

examples/mne/raw_eeg_basics/code.py

@@ -0,0 +1,74 @@
+"""
+A first look at MNE-Python: build a synthetic multi-channel EEG
+recording, wrap it in an MNE Raw object, and inspect it.
+
+MNE-Python (https://mne.tools/) is the de facto toolkit for
+M/EEG analysis. Its central data structures are `Info` (metadata
+about channels and sampling) and `Raw` (continuous time series).
+Here we build both from scratch using NumPy arrays.
+"""
+from IPython.core.display import display, HTML
+import numpy as np
+import matplotlib.pyplot as plt
+import mne
+
+# Quiet MNE's default INFO chatter so the example output stays focused.
+mne.set_log_level("WARNING")
+
+rng = np.random.default_rng(7)
+
+
+heading("1. Building a synthetic EEG recording")
+note(
+    "We'll simulate 60 seconds of data on five EEG channels at "
+    "250 Hz: a noisy background with an alpha-band (10 Hz) "
+    "oscillation that is stronger over the occipital channels."
+)
+
+sampling_freq = 250.0  # Hz
+duration_s = 60.0
+n_samples = int(sampling_freq * duration_s)
+times = np.arange(n_samples) / sampling_freq
+
+channel_names = ["Fz", "Cz", "Pz", "O1", "O2"]
+# Occipital channels (O1, O2) get a stronger 10 Hz alpha rhythm.
+alpha_gain = np.array([0.2, 0.4, 0.8, 1.5, 1.5])
+
+alpha = np.sin(2 * np.pi * 10.0 * times)
+background = rng.normal(0, 1.0, size=(len(channel_names), n_samples))
+# MNE expects EEG signals in volts; scale to microvolt-range.
+data = (alpha_gain[:, None] * alpha + background) * 1e-6
+
+# `create_info` describes the channels; combined with the data array
+# it yields a Raw object that behaves like a real recording.
+info = mne.create_info(
+    ch_names=channel_names,
+    sfreq=sampling_freq,
+    ch_types="eeg",
+)
+raw = mne.io.RawArray(data, info)
+
+note("The Raw object summarizes the recording:")
+display(HTML(f"<pre>{raw}</pre>"), append=True)
+display(HTML(f"<pre>{raw.info}</pre>"), append=True)
+
+heading("2. Plotting the first few seconds")
+note("Pulling out the first 3 seconds and plotting each channel.")
+
+window_data, window_times = raw[:, : int(3 * sampling_freq)]
+
+fig, ax = plt.subplots(figsize=(9, 4))
+offset = 6e-6  # vertical spacing between channels, in volts
+for i, name in enumerate(channel_names):
+    ax.plot(
+        window_times,
+        window_data[i] + i * offset,
+        linewidth=0.8,
+        label=name,
+    )
+ax.set_yticks([i * offset for i in range(len(channel_names))])
+ax.set_yticklabels(channel_names)
+ax.set_xlabel("Time (s)")
+ax.set_title("Synthetic EEG: first 3 seconds")
+fig.tight_layout()
+display(fig, append=True)

examples/mne/raw_eeg_basics/config.toml

@@ -0,0 +1 @@
+packages = ["mne", "numpy", "matplotlib"]

examples/mne/raw_eeg_basics/setup.py

@@ -0,0 +1,42 @@
+"""
+Shim IPython's display API onto PyScript so example code written in a
+Jupyter/IPython idiom runs unmodified in the browser.
+"""
+
+import sys
+import types
+import js
+from pyscript import window, HTML, display as _display
+
+js.alert = window.alert
+
+
+def display(*args, **kwargs):
+    """Wrap pyscript.display so output lands in the example target."""
+    return _display(
+        *args, **kwargs, target=__pyscript_display_target__,
+    )
+
+
+ipython = types.ModuleType("IPython")
+core = types.ModuleType("IPython.core")
+core_display = types.ModuleType("IPython.core.display")
+core_display.display = display
+core_display.HTML = HTML
+ipython.core = core
+core.display = core_display
+ipython.get_ipython = lambda: None
+ipython.display = core_display
+sys.modules["IPython"] = ipython
+sys.modules["IPython.core"] = core
+sys.modules["IPython.core.display"] = core_display
+sys.modules["IPython.display"] = core_display
+
+
+def heading(text, level=2):
+    display(HTML(f"<h{level}>{text}</h{level}>"), append=True)
+
+
+def note(text):
+    display(HTML(f"<p>{text}</p>"), append=True)
+