diff --git a/workflows/eeg-viewer/dev/230327_MNE_Qt_EEG_browser.ipynb b/workflows/eeg-viewer/dev/230327_MNE_Qt_EEG_browser.ipynb index e0c78c1..4eebeca 100644 --- a/workflows/eeg-viewer/dev/230327_MNE_Qt_EEG_browser.ipynb +++ b/workflows/eeg-viewer/dev/230327_MNE_Qt_EEG_browser.ipynb @@ -19,7 +19,7 @@ }, { "cell_type": "code", - "execution_count": 1, + "execution_count": 2, "id": "e848b88f-271f-4608-a1d3-9a76046e57b7", "metadata": { "tags": [] @@ -40,7 +40,17 @@ }, { "cell_type": "code", - "execution_count": 2, + "execution_count": 3, + "id": "f64dc29b-9963-4acb-8927-e112de1a0af7", + "metadata": {}, + "outputs": [], + "source": [ + "mne.set_config('MNE_DATASETS_SAMPLE_PATH', '/Users/droumis/data/mne_data/')" + ] + }, + { + "cell_type": "code", + "execution_count": 4, "id": "ac5dac4a-c73f-4307-b083-11cea9b7f147", "metadata": { "tags": [] @@ -49,10 +59,10 @@ { "data": { "text/plain": [ - "PosixPath('/Users/droumis/data/MNE-sample-data/MNE-sample-data')" + "PosixPath('/Users/droumis/data/mne_data/MNE-sample-data')" ] }, - "execution_count": 2, + "execution_count": 4, "metadata": {}, "output_type": "execute_result" } @@ -72,7 +82,7 @@ }, { "cell_type": "code", - "execution_count": 3, + "execution_count": 12, "id": "d2c0b381-c224-483d-8c02-48b9f722b1f4", "metadata": { "tags": [] @@ -81,10 +91,10 @@ { "data": { "text/plain": [ - "PosixPath('/Users/droumis/data/MNE-sample-data/MNE-sample-data/MEG/sample/sample_audvis_raw.fif')" + "PosixPath('/Users/droumis/data/mne_data/MNE-sample-data/MEG/sample/sample_audvis_raw.fif')" ] }, - "execution_count": 3, + "execution_count": 12, "metadata": {}, "output_type": "execute_result" } @@ -98,7 +108,7 @@ }, { "cell_type": "code", - "execution_count": 4, + "execution_count": 14, "id": "b00c6d20-9239-4102-9935-ac7a8726beb9", "metadata": { "tags": [] @@ -108,7 +118,7 @@ "name": "stdout", "output_type": "stream", "text": [ - "Opening raw data file /Users/droumis/data/MNE-sample-data/MNE-sample-data/MEG/sample/sample_audvis_raw.fif...\n", + "Opening raw data file /Users/droumis/data/mne_data/MNE-sample-data/MEG/sample/sample_audvis_raw.fif...\n", " Read a total of 3 projection items:\n", " PCA-v1 (1 x 102) idle\n", " PCA-v2 (1 x 102) idle\n", @@ -124,7 +134,50 @@ }, { "cell_type": "code", - "execution_count": 5, + "execution_count": 6, + "id": "16fadf8f-ad0f-4f6a-bb7b-9624e0b288fe", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "'/Users/droumis/data/mne_data/MNE-testing-data/EDF/test_edf_stim_resamp.edf'" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "sample_data_raw_file = ('/Users/droumis/data/mne_data/MNE-testing-data/EDF/test_edf_stim_resamp.edf')\n", + "sample_data_raw_file" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "id": "7ba5e861-c715-4c33-ac78-751354032b76", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Extracting EDF parameters from /Users/droumis/data/mne_data/S001R04.edf...\n", + "EDF file detected\n", + "Setting channel info structure...\n", + "Creating raw.info structure...\n" + ] + } + ], + "source": [ + "raw = mne.io.read_raw_edf('/Users/droumis/data/mne_data/S001R04.edf')" + ] + }, + { + "cell_type": "code", + "execution_count": 15, "id": "3c91f539-941f-4b96-824e-0b2ac76058dc", "metadata": { "tags": [] @@ -225,7 +278,7 @@ ">" ] }, - "execution_count": 5, + "execution_count": 15, "metadata": {}, "output_type": "execute_result" } @@ -244,7 +297,7 @@ }, { "cell_type": "code", - "execution_count": 6, + "execution_count": 16, "id": "d3e794f9-47cf-4872-a267-a47296519da4", "metadata": { "tags": [] @@ -254,17 +307,16 @@ "name": "stdout", "output_type": "stream", "text": [ - "Using qt as 2D backend.\n", "Using pyopengl with version 3.1.6\n" ] }, { "data": { "text/plain": [ - "" + "" ] }, - "execution_count": 6, + "execution_count": 16, "metadata": {}, "output_type": "execute_result" }, @@ -278,7 +330,7 @@ } ], "source": [ - "raw.plot(duration=10, n_channels=40) # launches GUI" + "raw.plot(duration=10, n_channels=40, clipping=None) # launches GUI" ] }, { diff --git a/workflows/eeg-viewer/dev/230802_subcoordinate_y.ipynb b/workflows/eeg-viewer/dev/230802_subcoordinate_y.ipynb new file mode 100644 index 0000000..9709b69 --- /dev/null +++ b/workflows/eeg-viewer/dev/230802_subcoordinate_y.ipynb @@ -0,0 +1,816 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Try out [the new subplot approach](https://github.com/holoviz/holoviews/pull/5840)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%load_ext autoreload\n", + "%autoreload 2\n", + "import numpy as np\n", + "\n", + "import holoviews as hv; hv.extension('bokeh')\n", + "from holoviews.plotting.links import RangeToolLink\n", + "from holoviews.operation.datashader import rasterize\n", + "from holoviews import Dataset\n", + "\n", + "from neurodatagen.eeg import generate_eeg_powerlaw\n", + "from hvneuro import download_file\n", + "\n", + "from bokeh.models import HoverTool\n", + "from scipy.stats import zscore\n", + "import panel as pn; pn.extension(template='material')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "n_channels = 10\n", + "n_seconds = 5\n", + "fs = 512\n", + "\n", + "data, time, channels = generate_eeg_powerlaw(n_channels, n_seconds, fs)\n", + "\n", + "print(f'shape: {data.shape} (n_channels, samples) ')\n", + "data" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Visualize synthetic data" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "channel_curves = []\n", + "ch_subcoord = 1./n_channels\n", + "for i, channel_data in enumerate(data):\n", + " \n", + " channel_curves.append(hv.Curve(channel_data).opts(\n", + " subcoordinate_y=(i*ch_subcoord, (i+1)*ch_subcoord), color=\"black\", line_width=1, \n", + " tools=['hover', 'xwheel_zoom'], shared_axes=False))\n", + "\n", + "eeg_viewer = hv.Overlay(channel_curves, kdims=\"Channel\").opts(width=800, height=500)\n", + "eeg_viewer" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## my shared version" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import pandas as pd\n", + "import holoviews as hv; hv.extension('bokeh')\n", + "from holoviews import opts\n", + "from holoviews import Dataset\n", + "from holoviews.plotting.links import RangeToolLink\n", + "from bokeh.models import HoverTool\n", + "import colorcet as cc\n", + "import panel as pn; pn.extension()\n", + "from scipy.stats import zscore\n", + "\n", + "n_channels = 10\n", + "n_seconds = 5\n", + "fs = 512\n", + "max_ch_disp = 5 # max channels to initially display\n", + "max_t_disp = 3 # max time in seconds to initially display\n", + "\n", + "total_samples = n_seconds * fs\n", + "time = np.linspace(0, n_seconds, total_samples)\n", + "data = np.random.randn(n_channels, total_samples).cumsum(axis=1)\n", + "channels = ['EEG {}'.format(i) for i in range(n_channels)]\n", + "\n", + "channel_curves = []\n", + "ch_subcoord = 1./n_channels\n", + "\n", + "hover = HoverTool(tooltips=[\n", + " (\"Channel\", \"@channel\"),\n", + " (\"Time\", \"$x s\"),\n", + " (\"Amplitude\", \"$y µV\")])\n", + "\n", + "for i, channel_data in enumerate(data):\n", + " ds = Dataset((time, channel_data, channels[i]), [\"Time\", \"Amplitude\", \"channel\"])\n", + " curve = hv.Curve(ds, \"Time\", [\"Amplitude\", \"channel\"], label=channels[i]).opts(\n", + " subcoordinate_y=(i*ch_subcoord, (i+1)*ch_subcoord), color=\"black\", line_width=1, \n", + " tools=[hover, 'xwheel_zoom'], shared_axes=False)\n", + " channel_curves.append(curve)\n", + "\n", + "eeg_viewer = hv.Overlay(channel_curves, kdims=\"Channel\").opts(\n", + " padding=0, xlabel=\"Time (s)\", ylabel=\"Channel\",\n", + " show_legend=False, aspect=1.5, responsive=True, # yticks=yticks,\n", + " shared_axes=False, backend_opts={\n", + " \"x_range.bounds\": (time.min(), time.max()),\n", + " \"y_range.bounds\": (data.min(), data.max())})\n", + "\n", + "# # Create events\n", + "# events = [\n", + "# {'start': 1, 'end': 1.5, 'label': 'event1'},\n", + "# {'start': 2, 'end': 2.5, 'label': 'event1'},\n", + "# {'start': 3, 'end': 3.5, 'label': 'event1'},\n", + "# {'start': 1.5, 'end': 2, 'label': 'event2'},\n", + "# {'start': 2.5, 'end': 3, 'label': 'event2'},\n", + "# ]\n", + "\n", + "# labels = list(set([event['label'] for event in events]))\n", + "# colors = cc.glasbey_bw_minc_20_minl_30[:len(labels)]\n", + "# label_colors = dict(zip(labels, colors))\n", + "\n", + "# for event in events:\n", + "# vspan = hv.VSpan(event['start'], event['end']).opts(fill_alpha=0.5, fill_color=label_colors[event['label']])\n", + "# eeg_viewer *= vspan\n", + "\n", + "yticks = [( (i*ch_subcoord + (i+1)*ch_subcoord) / 2, ich) for i, ich in enumerate(channels)]\n", + "y_positions, _ = zip(*yticks)\n", + "\n", + "z_data = zscore(data, axis=1)\n", + "\n", + "minimap = hv.Image((time, y_positions, z_data), [\"Time (s)\", \"Channel\"], \"Amplitude (uV)\")\n", + "minimap = minimap.opts(\n", + " cmap=\"RdBu_r\", colorbar=False, xlabel='', alpha=.5, yticks=[yticks[0], yticks[-1]],\n", + " height=100, responsive=True, default_tools=[''], shared_axes=False, clim=(-z_data.std()*2.5, z_data.std()*2.5))\n", + "\n", + "\n", + "if len(channels) < max_ch_disp:\n", + " max_ch_disp = len(channels)\n", + "max_y_disp = (max_ch_disp+2)*ch_subcoord\n", + "\n", + "time_s = len(time)/fs\n", + "if time_s < max_t_disp:\n", + " max_t_disp = time_s\n", + " \n", + "RangeToolLink(minimap, list(eeg_viewer.values())[0], axes=[\"x\", \"y\"],\n", + " boundsx=(None, max_t_disp),\n", + " boundsy=(None, max_y_disp))\n", + "\n", + "# eeg_app = pn.Column(pn.Row(eeg_viewer, min_height=500), minimap).servable(target='main')\n", + "eeg_app = pn.Column((eeg_viewer + minimap).cols(1), min_height=650).servable(target='main')\n", + "eeg_app" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## philipps version after pr work" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import pandas as pd\n", + "import holoviews as hv; hv.extension('bokeh')\n", + "from holoviews import opts\n", + "from holoviews import Dataset\n", + "from holoviews.plotting.links import RangeToolLink\n", + "from bokeh.models import HoverTool\n", + "import panel as pn; pn.extension()\n", + "from scipy.stats import zscore\n", + "\n", + "n_channels = 10\n", + "n_seconds = 5\n", + "fs = 512\n", + "max_ch_disp = 5 # max channels to initially display\n", + "max_t_disp = 3 # max time in seconds to initially display\n", + "\n", + "total_samples = n_seconds * fs\n", + "time = np.linspace(0, n_seconds, total_samples)\n", + "data = np.random.randn(n_channels, total_samples).cumsum(axis=1)\n", + "channels = ['EEG {}'.format(i) for i in range(n_channels)]\n", + "\n", + "channel_curves = []\n", + "ch_subcoord = 1./n_channels\n", + "\n", + "hover = HoverTool(tooltips=[\n", + " (\"Channel\", \"@channel\"),\n", + " (\"Time\", \"$x s\"),\n", + " (\"Amplitude\", \"$y µV\")])\n", + "\n", + "for i, channel_data in enumerate(data):\n", + " ds = Dataset((time, channel_data, channels[i]), [\"Time\", \"Amplitude\", \"channel\"])\n", + " channel_curves.append(hv.Curve(ds, \"Time\", [\"Amplitude\", \"channel\"], label=channels[i]).opts(\n", + " subcoordinate_y=(i*ch_subcoord, (i+1)*ch_subcoord), color=\"black\", line_width=1, \n", + " tools=[hover, 'xwheel_zoom'], shared_axes=False))\n", + "\n", + "annotation = hv.VSpan(0.3, 0.5)\n", + "eeg_viewer = (hv.Overlay(channel_curves, kdims=\"Channel\") * annotation).opts(\n", + " padding=0, xlabel=\"Time (s)\", ylabel=\"Channel\",\n", + " show_legend=False, aspect=1.5, min_height=500, responsive=True,\n", + " shared_axes=False, backend_opts={\n", + " \"x_range.bounds\": (time.min(), time.max()),\n", + " \"y_range.bounds\": (0, 1)}) \n", + "\n", + "yticks = [( (i*ch_subcoord + (i+1)*ch_subcoord) / 2, ich) for i, ich in enumerate(channels)]\n", + "y_positions, _ = zip(*yticks)\n", + "\n", + "z_data = zscore(data, axis=1)\n", + "\n", + "minimap = hv.Image((time, y_positions, z_data), [\"Time (s)\", \"Channel\"], \"Amplitude (uV)\")\n", + "minimap = minimap.opts(\n", + " cmap=\"RdBu_r\", colorbar=False, xlabel='', alpha=.5, yticks=[yticks[0], yticks[-1]],\n", + " height=100, responsive=True, default_tools=[''], shared_axes=False, clim=(-z_data.std()*2.5, z_data.std()*2.5))\n", + "\n", + "if len(channels) < max_ch_disp:\n", + " max_ch_disp = len(channels)\n", + "max_y_disp = (max_ch_disp+2)*ch_subcoord\n", + "\n", + "time_s = len(time)/fs\n", + "if time_s < max_t_disp:\n", + " max_t_disp = time_s\n", + " \n", + "RangeToolLink(minimap, eeg_viewer, axes=[\"x\", \"y\"],\n", + " boundsx=(None, max_t_disp),\n", + " boundsy=(None, max_y_disp))\n", + "\n", + "eeg_app = pn.Column((eeg_viewer + minimap * annotation).cols(1), min_height=650).servable(target='main')\n", + "eeg_app" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Bokeh" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from bokeh.models import ColumnDataSource, Range1d, RangeTool\n", + "from bokeh.plotting import figure\n", + "from bokeh.layouts import column\n", + "from bokeh.palettes import Viridis10\n", + "from bokeh.transform import linear_cmap\n", + "from bokeh.io import output_notebook, show; output_notebook()\n", + "import numpy as np\n", + "from scipy.stats import zscore\n", + "\n", + "n_channels = 10\n", + "n_seconds = 5\n", + "fs = 512\n", + "total_samples = n_seconds * fs\n", + "time = np.linspace(0, n_seconds, total_samples)\n", + "data = np.random.randn(n_channels, total_samples).cumsum(axis=1)\n", + "channels = ['EEG {}'.format(i) for i in range(n_channels)]\n", + "\n", + "p = figure(height=500, width=900, tools=\"xwheel_zoom,xpan,xbox_zoom,reset\", active_scroll=\"xwheel_zoom\")\n", + "\n", + "for i, channel_data in enumerate(data):\n", + " y_range = Range1d(start=i, end=i+1)\n", + " subplot = p.subplot(x_source=p.x_range, x_target=p.x_range, y_source=p.y_range, y_target=y_range)\n", + " source = ColumnDataSource(data=dict(x=time, y=channel_data, channel=[channels[i]]*len(time)))\n", + " subplot.line('x', 'y', source=source, line_color=\"black\", line_width=1, alpha=0.7)\n", + "\n", + "# p.yaxis.ticker = list(range(n_channels))\n", + "# p.yaxis.major_label_overrides = {i: channel for i, channel in enumerate(channels)}\n", + "p.xaxis.axis_label = 'Time (s)'\n", + "p.yaxis.axis_label = 'Channel'\n", + "p.grid.visible = False\n", + "\n", + "show(p)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Original" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "max_ch_disp = 10 # max channels to initially display\n", + "max_t_disp = 4 # max time in seconds to initially display\n", + "\n", + "spacing = 2.5 # Spacing between channels\n", + "offset = np.std(data) * spacing\n", + "\n", + "# Create a hv.Curve element per chan.\n", + "# Note: alternative is to call hv.Path once on offset-adjusted data, but \n", + "# then we couldn't independently apply formating to the channels (which \n", + "# we aren't doing yet, but we likely will in the future)\n", + "channel_curves = []\n", + "max_data = data.max()\n", + "hover = HoverTool(tooltips=[\n", + " (\"Channel\", \"@channel\"),\n", + " (\"Time\", \"$x s\"),\n", + " (\"Amplitude\", \"@original_amplitude µV\")])\n", + "for i, channel_data in enumerate(data):\n", + " offset_data = channel_data + (i * offset)\n", + " max_data = max(offset_data.max(), max_data) # update max\n", + " ds = Dataset((time, offset_data, channel_data, channels[i]), [\"Time\", \"Amplitude\", \"original_amplitude\", \"channel\"])\n", + " channel_curves.append(\n", + " hv.Curve(ds, \"Time\", [\"Amplitude\", \"original_amplitude\", \"channel\"]).opts(\n", + " color=\"black\", line_width=1,\n", + " tools=[hover, 'xwheel_zoom'], shared_axes=False))\n", + "\n", + "# Create mapping from yaxis location to ytick for each channel\n", + "# so we can have categorical-style labeling on a continuous axis.\n", + "# Note: this would/should change when we implement independent \n", + "# coordinates.\n", + "yticks = [(i * offset, ich) for i, ich in enumerate(channels)]\n", + "\n", + "\n", + "# Create an overlay of curves\n", + "# TODO.. setting x/y_range bounds does not yet restrict the RangeTool from going beyond these limits\n", + "# TODO.. the zoom out will stop when it hits any single bound, and not continue zooming out in other directions/dims\n", + "eeg_viewer = hv.Overlay(channel_curves, kdims=\"Channel\").opts(\n", + " padding=0, xlabel=\"Time (s)\", ylabel=\"Channel\", #default_tools=['hover', 'pan', 'box_zoom', 'save', 'reset'],\n", + " yticks=yticks, show_legend=False, aspect=1.5, responsive=True,\n", + " shared_axes=False, backend_opts={\n", + " \"x_range.bounds\": (time.min(), time.max()),\n", + " \"y_range.bounds\": (data.min(), max_data)})\n", + "\n", + "\n", + "# # Create events\n", + "# events = [\n", + "# {'start': 1, 'end': 1.5, 'label': 'event1'},\n", + "# {'start': 2, 'end': 2.5, 'label': 'event1'},\n", + "# {'start': 3, 'end': 3.5, 'label': 'event1'},\n", + "# {'start': 1.5, 'end': 2, 'label': 'event2'},\n", + "# {'start': 2.5, 'end': 3, 'label': 'event2'},\n", + "# ]\n", + "\n", + "# labels = list(set([event['label'] for event in events]))\n", + "# colors = cc.glasbey_category10[:len(labels)]\n", + "# label_colors = dict(zip(labels, colors))\n", + "\n", + "# for event in events:\n", + "# vspan = hv.VSpan(event['start'], event['end']).opts(responsive=True, fill_alpha=0.5, fill_color=tuple(label_colors[event['label']]))\n", + "# eeg_viewer *= vspan\n", + "\n", + "# Get the y positions of the yticks to use as yaxis of minimap image\n", + "y_positions, _ = zip(*yticks)\n", + "\n", + "# Compute z-scores across time for each channel\n", + "z_data = zscore(data, axis=1)\n", + "\n", + "# Generate the zscored image for the minimap using the y tiack positions from the eeg_viewer\n", + "minimap = hv.Image((time, y_positions, z_data), [\"Time (s)\", \"Channel\"], \"Amplitude (uV)\")\n", + "\n", + "# RangeTool doesn't work with rasterized object? TODO: file issue\n", + "# minimap = rasterize(minimap)\n", + "\n", + "# maybe I should datashade/2d-bin the data before creating the hv.Image\n", + "# I could use lttb (1d so per channel) or ResampleOperation2D (but I think that applies to the entire nb)\n", + "# or some operation from datashader to return the 2d hist\n", + "\n", + "\n", + "# Style the minimap \n", + "minimap = minimap.opts(\n", + " cmap=\"RdBu_r\", colorbar=False, xlabel='', alpha=.5, yticks=[yticks[0], yticks[-1]],\n", + " height=100, responsive=True, default_tools=[''], shared_axes=False, clim=(-z_data.std()*2.5, z_data.std()*2.5))\n", + " \n", + "# Create RangeToolLink between the minimap and the main EEG viewer \n", + "# (quirk: apply to just one eeg trace and it will apply to all. see HoloViews #4472)\n", + "max_y_disp = np.max(data[max_ch_disp-1,:] + ((max_ch_disp-1) * offset))\n", + "RangeToolLink(minimap, list(eeg_viewer.values())[0], axes=[\"x\", \"y\"],\n", + " boundsx=(None, max_t_disp),\n", + " boundsy=(None, max_y_disp))\n", + "\n", + "# Display vertically\n", + "# layout = (eeg_viewer + minimap).cols(1).opts(shared_axes=False, merge_tools=False)\n", + "# eeg_app = pn.Row(layout).servable() # too much spacing between plots in served app\n", + "# eeg_app = pn.Column(pn.Row(eeg_viewer, min_height=500, sizing_mode='stretch_both'), minimap, sizing_mode='stretch_both')#.servable()#target='main') # BUG Panel #5315: rangetool is variably active in the bokeh toolbar on eeg viewer plot.. not respecting shared_axes=False\n", + "# eeg_app\n", + "eeg_app = pn.Column((eeg_viewer + minimap).cols(1), min_height=650).servable(target='main')\n", + "eeg_app" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Real data pipeline" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import mne" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Intake data" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# This dataset is 2.6 MB on disk\n", + "url = \"https://physionet.org/files/eegmmidb/1.0.0/S001/S001R04.edf?download\"\n", + "local_data_path = \"~/data/mne_data\"\n", + "\n", + "\n", + "# Will not download if already present at local_data_path\n", + "local_file_path = download_file(url, local_data_path)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "raw = mne.io.read_raw_edf(local_file_path, preload=True)\n", + "raw.info" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "# preview the channel names, types, signal ranges, and uncompressed size\n", + "raw.describe()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Clean channel names, set sensor positions, and reference data" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# clean up the channel names\n", + "raw.rename_channels(lambda s: s.strip(\".\"));" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "# preview available montages that are shipped with MNE\n", + "# mne.channels.get_builtin_montages(descriptions=True)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Let's use the standard 10-20\n", + "montage = mne.channels.make_standard_montage(\"standard_1020\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# plot the assigned positions of our data channels\n", + "raw.set_montage(montage, match_case=False)\n", + "sphere=(0, 0.015, 0, 0.099) #manually adjust the y origin coord and radius a bit\n", + "raw.plot_sensors(show_names=True, sphere=sphere);" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# re-reference EEG data to the average over all recording channels\n", + "raw.set_eeg_reference(\"average\");" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Gather the data for plotting into simple arrays" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "time = raw.times\n", + "channels = raw.ch_names\n", + "\n", + "# get the EEG data (for this data set, all channels are EEG anyways)\n", + "eeg_indices = mne.pick_types(raw.info, eeg=True)\n", + "data = raw.get_data(picks=eeg_indices, units={\"eeg\":\"uV\"})" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# data = data[:5,600:1600]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Visualize real data" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Philipp's new version with real data" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import pandas as pd\n", + "import holoviews as hv; hv.extension('bokeh')\n", + "from holoviews import opts\n", + "from holoviews import Dataset\n", + "from holoviews.plotting.links import RangeToolLink\n", + "from bokeh.models import HoverTool\n", + "import panel as pn; pn.extension()\n", + "from scipy.stats import zscore\n", + "\n", + "n_channels = len(channels)\n", + "# n_seconds = 5\n", + "fs = raw.info['sfreq']\n", + "max_ch_disp = 5 # max channels to initially display\n", + "max_t_disp = 3 # max time in seconds to initially display\n", + "\n", + "# total_samples = n_seconds * fs\n", + "# time = np.linspace(0, n_seconds, total_samples)\n", + "# data = np.random.randn(n_channels, total_samples).cumsum(axis=1)\n", + "# channels = ['EEG {}'.format(i) for i in range(n_channels)]\n", + "\n", + "channel_curves = []\n", + "ch_subcoord = 1./n_channels\n", + "\n", + "hover = HoverTool(tooltips=[\n", + " (\"Channel\", \"@channel\"),\n", + " (\"Time\", \"$x s\"),\n", + " (\"Amplitude\", \"$y µV\")])\n", + "\n", + "for i, channel_data in enumerate(data):\n", + " ds = Dataset((time, channel_data, channels[i]), [\"Time\", \"Amplitude\", \"channel\"])\n", + " channel_curves.append(hv.Curve(ds, \"Time\", [\"Amplitude\", \"channel\"], label=channels[i]).opts(\n", + " subcoordinate_y=((i)*(1./n_channels), (i+10)*(1./n_channels)), color=\"black\", line_width=1, \n", + " tools=[hover, 'xwheel_zoom'], shared_axes=False))\n", + "\n", + "annotation = hv.VSpan(0.3, 0.5)\n", + "eeg_viewer = (hv.Overlay(channel_curves, kdims=\"Channel\") * annotation).opts(\n", + " padding=0, xlabel=\"Time (s)\", ylabel=\"Channel\",\n", + " show_legend=False, aspect=1.5, min_height=500, responsive=True,\n", + " shared_axes=False, backend_opts={\n", + " \"x_range.bounds\": (time.min(), time.max()),\n", + " \"y_range.bounds\": (0, 1)}) \n", + "\n", + "yticks = [( (i*ch_subcoord + (i+1)*ch_subcoord) / 2, ich) for i, ich in enumerate(channels)]\n", + "y_positions, _ = zip(*yticks)\n", + "\n", + "z_data = zscore(data, axis=1)\n", + "\n", + "minimap = hv.Image((time, y_positions, z_data), [\"Time (s)\", \"Channel\"], \"Amplitude (uV)\")\n", + "minimap = minimap.opts(\n", + " cmap=\"RdBu_r\", colorbar=False, xlabel='', alpha=.5, yticks=[yticks[0], yticks[-1]],\n", + " height=100, responsive=True, default_tools=[''], shared_axes=False, clim=(-z_data.std()*2.5, z_data.std()*2.5))\n", + "\n", + "minimap = rasterize(minimap)\n", + "\n", + "if len(channels) < max_ch_disp:\n", + " max_ch_disp = len(channels)\n", + "max_y_disp = (max_ch_disp+2)*ch_subcoord\n", + "\n", + "time_s = len(time)/fs\n", + "if time_s < max_t_disp:\n", + " max_t_disp = time_s\n", + " \n", + "RangeToolLink(minimap, eeg_viewer, axes=[\"x\", \"y\"],\n", + " boundsx=(None, max_t_disp),\n", + " boundsy=(None, max_y_disp))\n", + "\n", + "eeg_app = pn.Column((eeg_viewer + minimap * annotation).cols(1), min_height=650).servable(target='main')\n", + "eeg_app" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "max_ch_disp = 10 # max channels to initially display\n", + "max_t_disp = 5 # max time in seconds to initially display\n", + "\n", + "n_channels = data.shape[0]\n", + "channel_curves = []\n", + "ch_subcoord = 1./n_channels\n", + "\n", + "hover = HoverTool(tooltips=[\n", + " (\"Channel\", \"@channel\"),\n", + " (\"Time\", \"$x s\"),\n", + " (\"Amplitude\", \"$y µV\")])\n", + "\n", + "for i, channel_data in enumerate(data):\n", + " ds = Dataset((time, channel_data, channels[i]), [\"Time\", \"Amplitude\", \"channel\"])\n", + " channel_curves.append(hv.Curve(ds, \"Time\", [\"Amplitude\", \"channel\"], label=channels[i]).opts(\n", + " subcoordinate_y=((i-1)*ch_subcoord, (i+2)*ch_subcoord), color=\"black\", line_width=1, \n", + " tools=[hover, 'xwheel_zoom'], shared_axes=False))\n", + "\n", + "eeg_viewer = hv.Overlay(channel_curves, kdims=\"Channel\").opts(responsive=True, ylabel=\"Channel\")\n", + "\n", + "yticks = [(((i+2)*ch_subcoord)+((i-1)*ch_subcoord/2), ich) for i, ich in enumerate(channels)]\n", + "y_positions, _ = zip(*yticks)\n", + "\n", + "z_data = zscore(data, axis=1)\n", + "\n", + "minimap = hv.Image((time, y_positions, z_data), [\"Time (s)\", \"Channel\"], \"Amplitude (uV)\")\n", + "minimap = minimap.opts(\n", + " cmap=\"RdBu_r\", colorbar=False, xlabel='', alpha=.5, yticks=[yticks[0], yticks[-1]],\n", + " height=100, responsive=True, default_tools=[''], shared_axes=False, clim=(-z_data.std()*2.5, z_data.std()*2.5))\n", + "\n", + "\n", + "if len(channels) < max_ch_disp:\n", + " max_ch_disp = len(channels)\n", + "max_y_disp = (max_ch_disp+2)*ch_subcoord\n", + "\n", + "time_s = len(time)/raw.info['sfreq']\n", + "if time_s < max_t_disp:\n", + " max_t_disp = time_s\n", + " \n", + "RangeToolLink(minimap, list(eeg_viewer.values())[0], axes=[\"x\", \"y\"],\n", + " boundsx=(None, max_t_disp),\n", + " boundsy=(None, max_y_disp))\n", + "\n", + "eeg_app = pn.Column(pn.Row(eeg_viewer, min_height=500), minimap).servable(target='main')\n", + "eeg_app" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "max_ch_disp = 10 # max channels to initially display\n", + "max_t_disp = 5 # max time in seconds to initially display\n", + "\n", + "spacing = 2.5 # Spacing between channels\n", + "offset = np.std(data) * spacing\n", + "\n", + "# Create a hv.Curve element per channel\n", + "channel_curves = []\n", + "max_data = data.max()\n", + "hover = HoverTool(tooltips=[\n", + " (\"Channel\", \"@channel\"),\n", + " (\"Time\", \"$x s\"),\n", + " (\"Amplitude\", \"@original_amplitude µV\")])\n", + "for i, channel_data in enumerate(data):\n", + " offset_data = channel_data + (i * offset)\n", + " max_data = max(offset_data.max(), max_data) # update max\n", + " ds = Dataset((time, offset_data, channel_data, channels[i]), [\"Time\", \"Amplitude\", \"original_amplitude\", \"channel\"])\n", + " channel_curves.append(\n", + " hv.Curve(ds, \"Time\", [\"Amplitude\", \"original_amplitude\", \"channel\"]).opts(\n", + " color=\"black\", line_width=1,\n", + " tools=[hover, 'xwheel_zoom'], shared_axes=False))\n", + "\n", + "# Create mapping from yaxis location to ytick for each channel\n", + "yticks = [(i * offset, ich) for i, ich in enumerate(channels)]\n", + "\n", + "# Create an overlay of curves\n", + "eeg_viewer = hv.Overlay(channel_curves, kdims=\"Channel\").opts(\n", + " padding=0, xlabel=\"Time (s)\", ylabel=\"Channel\", yticks=yticks, show_legend=False, aspect=1.5, responsive=True,\n", + " shared_axes=False, backend_opts={\n", + " \"x_range.bounds\": (time.min(), time.max()),\n", + " \"y_range.bounds\": (data.min(), max_data)})\n", + "\n", + "# Get the y positions of the yticks to use as yaxis of minimap image\n", + "y_positions, _ = zip(*yticks)\n", + "\n", + "# Compute z-scores across time for each channel\n", + "z_data = zscore(data, axis=1)\n", + "\n", + "# Generate the zscored image for the minimap using the y tick positions from the eeg_viewer\n", + "minimap = hv.Image((time, y_positions, z_data), [\"Time (s)\", \"Channel\"], \"Amplitude (uV)\")\n", + "minimap = minimap.opts(\n", + " cmap=\"RdBu_r\", colorbar=False, xlabel='', alpha=.5, yticks=[yticks[0], yticks[-1]],\n", + " height=100, responsive=True, default_tools=[''], shared_axes=False, clim=(-z_data.std()*2.5, z_data.std()*2.5))\n", + "\n", + "# Create RangeToolLink between the minimap and the main EEG viewer \n", + "if len(channels) < max_ch_disp:\n", + " max_ch_disp = len(channels)\n", + "max_y_disp = np.max(data[max_ch_disp-1,:] + ((max_ch_disp-1) * offset))\n", + "\n", + "time_s = len(time)/raw.info['sfreq']\n", + "if time_s < max_t_disp:\n", + " max_t_disp = time_s\n", + " \n", + "RangeToolLink(minimap, list(eeg_viewer.values())[0], axes=[\"x\", \"y\"],\n", + " boundsx=(None, max_t_disp),\n", + " boundsy=(None, max_y_disp))\n", + "\n", + "eeg_app = pn.Column(pn.Row(eeg_viewer, min_height=500), minimap).servable(target='main')\n", + "eeg_app" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.9.16" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/workflows/eeg-viewer/workflow_eeg-viewer.ipynb b/workflows/eeg-viewer/workflow_eeg-viewer.ipynb index c42ba26..fab2639 100644 --- a/workflows/eeg-viewer/workflow_eeg-viewer.ipynb +++ b/workflows/eeg-viewer/workflow_eeg-viewer.ipynb @@ -124,16 +124,17 @@ "spacing = 5.5 # Spacing between channels\n", "offset = np.std(data) * spacing\n", "\n", - "# Create a hv.Curve element per chan.\n", - "# Note: alternative is to call hv.Path once on offset-adjusted data, but \n", - "# then we couldn't independently apply formating to the channels (which \n", - "# we aren't doing yet, but we likely will in the future)\n", + "annotation = hv.VSpan(1, 2) # example annotation (start, end) time\n", + "\n", + "# Create a hv.Curve element per chan\n", "channel_curves = []\n", "max_data = data.max()\n", + " \n", "hover = HoverTool(tooltips=[\n", " (\"Channel\", \"@channel\"),\n", " (\"Time\", \"$x s\"),\n", " (\"Amplitude\", \"@original_amplitude µV\")])\n", + "\n", "for i, channel_data in enumerate(data):\n", " offset_data = channel_data + (i * offset)\n", " max_data = max(offset_data.max(), max_data) # update max\n", @@ -149,11 +150,10 @@ "# coordinates.\n", "yticks = [(i * offset, ich) for i, ich in enumerate(channels)]\n", "\n", - "\n", "# Create an overlay of curves\n", "# TODO.. setting x/y_range bounds does not yet restrict the RangeTool from going beyond these limits\n", "# TODO.. the zoom out will stop when it hits any single bound, and not continue zooming out in other directions/dims\n", - "eeg_viewer = hv.Overlay(channel_curves, kdims=\"Channel\").opts(\n", + "eeg_viewer = (annotation * hv.Overlay(channel_curves, kdims=\"Channel\")).opts(\n", " padding=0, xlabel=\"Time (s)\", ylabel=\"Channel\", #default_tools=['hover', 'pan', 'box_zoom', 'save', 'reset'],\n", " yticks=yticks, show_legend=False, aspect=1.5, responsive=True,\n", " shared_axes=False, backend_opts={\n", @@ -167,20 +167,13 @@ "z_data = zscore(data, axis=1)\n", "\n", "# Generate the zscored image for the minimap using the y tiack positions from the eeg_viewer\n", - "minimap = hv.Image((time, y_positions, z_data), [\"Time (s)\", \"Channel\"], \"Amplitude (uV)\")\n", - "\n", - "# RangeTool doesn't work with rasterized object? TODO: file issue\n", - "# minimap = rasterize(minimap)\n", - "\n", - "# maybe I should datashade/2d-bin the data before creating the hv.Image\n", - "# I could use lttb (1d so per channel) or ResampleOperation2D (but I think that applies to the entire nb)\n", - "# or some operation from datashader to return the 2d hist\n", - "\n", + "minimap = rasterize(hv.Image((time, y_positions, z_data), [\"Time (s)\", \"Channel\"], \"Amplitude (uV)\"))\n", "\n", "# Style the minimap \n", + "clim_mul = 1.2\n", "minimap = minimap.opts(\n", " cmap=\"RdBu_r\", colorbar=False, xlabel='', alpha=.5, yticks=[yticks[0], yticks[-1]],\n", - " height=100, responsive=True, default_tools=[''], shared_axes=False, clim=(-z_data.std()*2.5, z_data.std()*2.5))\n", + " height=100, responsive=True, default_tools=[''], shared_axes=False, clim=(-z_data.std()*clim_mul, z_data.std()*clim_mul))\n", " \n", "# Create RangeToolLink between the minimap and the main EEG viewer \n", "# (quirk: apply to just one eeg trace and it will apply to all. see HoloViews #4472)\n", @@ -189,10 +182,13 @@ " boundsx=(None, max_t_disp),\n", " boundsy=(None, max_y_disp))\n", "\n", - "# Display vertically\n", + "\n", "# layout = (eeg_viewer + minimap).cols(1).opts(shared_axes=False, merge_tools=False)\n", "# eeg_app = pn.Row(layout).servable() # too much spacing between plots in served app\n", - "eeg_app = pn.Column(pn.Row(eeg_viewer, min_height=500, sizing_mode='stretch_both'), minimap, sizing_mode='stretch_both')#.servable()#target='main') # BUG Panel #5315: rangetool is variably active in the bokeh toolbar on eeg viewer plot.. not respecting shared_axes=False\n", + "# eeg_app = pn.Column(pn.Row(eeg_viewer, min_height=500, sizing_mode='stretch_both'), minimap, sizing_mode='stretch_both')#.servable()#target='main') # BUG Panel #5315: rangetool is variably active in the bokeh toolbar on eeg viewer plot.. not respecting shared_axes=False\n", + "\n", + "# reverting approach because of the rangetool bug.. will deal with the spacing in the served app later\n", + "eeg_app = pn.Column((eeg_viewer + minimap * annotation).cols(1), min_height=650).servable(target='main', title='EEG Viewer with HoloViz and Bokeh')\n", "eeg_app" ] }, @@ -216,7 +212,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "### Download, read real dataset" + "### Intake data" ] }, { @@ -280,7 +276,7 @@ }, "outputs": [], "source": [ - "# preview available montages that are shipped with MNE\n", + "# # preview available montages that are shipped with MNE\n", "# mne.channels.get_builtin_montages(descriptions=True)" ] }, @@ -290,8 +286,8 @@ "metadata": {}, "outputs": [], "source": [ - "# Let's use the standard 10-20\n", - "montage = mne.channels.make_standard_montage(\"standard_1020\")" + "# # Let's use the standard 10-20\n", + "# montage = mne.channels.make_standard_montage(\"standard_1020\")" ] }, { @@ -300,10 +296,10 @@ "metadata": {}, "outputs": [], "source": [ - "# plot the assigned positions of our data channels\n", - "raw.set_montage(montage, match_case=False)\n", - "sphere=(0, 0.015, 0, 0.099) #manually adjust the y origin coord and radius a bit\n", - "raw.plot_sensors(show_names=True, sphere=sphere);" + "# # plot the assigned positions of our data channels\n", + "# raw.set_montage(montage, match_case=False)\n", + "# sphere=(0, 0.015, 0, 0.099) #manually adjust the y origin coord and radius a bit\n", + "# raw.plot_sensors(show_names=True, sphere=sphere);" ] }, { @@ -389,7 +385,7 @@ "z_data = zscore(data, axis=1)\n", "\n", "# Generate the zscored image for the minimap using the y tick positions from the eeg_viewer\n", - "minimap = hv.Image((time, y_positions, z_data), [\"Time (s)\", \"Channel\"], \"Amplitude (uV)\")\n", + "minimap = rasterize(hv.Image((time, y_positions, z_data), [\"Time (s)\", \"Channel\"], \"Amplitude (uV)\"))\n", "minimap = minimap.opts(\n", " cmap=\"RdBu_r\", colorbar=False, xlabel='', alpha=.5, yticks=[yticks[0], yticks[-1]],\n", " height=100, responsive=True, default_tools=[''], shared_axes=False, clim=(-z_data.std()*2.5, z_data.std()*2.5))\n", @@ -407,7 +403,7 @@ " boundsx=(None, max_t_disp),\n", " boundsy=(None, max_y_disp))\n", "\n", - "eeg_app = pn.Column(pn.Row(eeg_viewer, min_height=500), minimap).servable(target='main')\n", + "eeg_app = pn.Column((eeg_viewer + minimap).cols(1), min_height=650) #.servable(target='main')\n", "eeg_app" ] },