Add implicit frustration demo

README.md

@@ -7,6 +7,7 @@ Currently, you can find:
 - a [hyperscanning](../../tree/main/hyperscanning/) demos using Jinja and environment variables for customization ;
 - a fully-configurable [oddball protocol](../../tree/main/oddball) with visual, auditory and haptic support ;
 - a simple neurofeedback application based on [frequency bands](../../tree/main/neurofeedback/bands/).
+- a [implicit frustration](implicit/frustration) application with classification scripts.
 
 More demos will be added soon.
 Have fun!

implicit/frustration/.gitignore

@@ -0,0 +1 @@
+data

implicit/frustration/LICENSE

@@ -0,0 +1,21 @@
+MIT License
+
+Copyright (c) 2021
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.

implicit/frustration/README.md

@@ -0,0 +1,11 @@
+Implicit Frustration application based on [Neuroadaptive technology enables implicit cursor
+control based on medial prefrontal cortex activity](https://www.pnas.org/doi/abs/10.1073/pnas.1605155114).
+
+The application is already implemented for two eeg headsets:
+- the muse 2016 BLED by Muse in LSL
+    - stream the LSL using muse_stream.sh
+    - launch the application using timeflux : timeflux muse.yaml
+- the unicorn by G.TEC through Brainflow
+    - just launch the application using timeflux : timeflux unicorn.yaml
+
+See [GUI](gui) for application details.

implicit/frustration/analysis/visualization.py

@@ -0,0 +1,241 @@
+import pandas as pd
+import numpy as np
+from numpy import array, arccos
+from numpy.linalg import norm
+from math import degrees
+import mne
+import json
+import matplotlib.pyplot as plt
+
+fnames = [
+    "../exp_data/20220525-120400-S1-Acq2.hdf5",
+    "../exp_data/20220525-125204-S2-Acq2.hdf5",
+]
+snames = ["S1", "S2"]
+
+tmin = -0.2
+baseline = (None, 0)
+l_freq, h_freq = 0.1, 25
+
+
+def plot_evoked_5_dir(epochs, suffix=""):
+    partition = {"0-10": [], "10-45": [], "45-90": [], "90-135": [], "135-180": []}
+    for ei in epochs.event_id:
+        if ei.isalpha():
+            continue
+        elif 0 <= int(ei) <= 10 and ei not in partition["0-10"]:
+            partition["0-10"].append(ei)
+        elif 10 < int(ei) <= 45 and ei not in partition["10-45"]:
+            partition["10-45"].append(ei)
+        elif 45 < int(ei) <= 90 and ei not in partition["45-90"]:
+            partition["45-90"].append(ei)
+        elif 90 < int(ei) <= 135 and ei not in partition["90-135"]:
+            partition["90-135"].append(ei)
+        elif 135 < int(ei) <= 180 and ei not in partition["135-180"]:
+            partition["135-180"].append(ei)
+    old_event_ids = [
+        partition["0-10"],
+        partition["10-45"],
+        partition["45-90"],
+        partition["90-135"],
+        partition["135-180"],
+    ]
+    new_event_id = [
+        {"<0": 100},
+        {"<45": 101},
+        {"<90": 102},
+        {"<135": 103},
+        {"<180": 104},
+    ]
+    for old, new in zip(old_event_ids, new_event_id):
+        epochs = mne.epochs.combine_event_ids(epochs, old, new)
+
+    evokeds = {
+        "dir_0": epochs["<0"].average(),
+        "dir_45": epochs["<45"].average(),
+        "dir_90": epochs["<90"].average(),
+        "dir_135": epochs["<135"].average(),
+        "dir_180": epochs["<180"].average(),
+    }
+    for elec in ["AF7", "AF8", "TP9", "TP10"]:
+        nave = ", ".join([str(e.nave) for e in evokeds.values()])
+        mne.viz.plot_compare_evokeds(
+            evokeds,
+            picks=[elec],
+            cmap="viridis",
+            title=f"{elec} - nave: {nave}",
+            show=False,
+        )
+        plt.savefig(f"{subject}-{elec}{suffix}.png")
+        plt.close()
+    return epochs
+
+
+def plot_evoked_3_dir(epochs, suffix="alt"):
+    partition = {"0-60": [], "60-120": [], "120-180": []}
+    for ei in epochs.event_id:
+        if ei.isalpha():
+            continue
+        elif 0 <= int(ei) <= 60 and ei not in partition["0-60"]:
+            partition["0-60"].append(ei)
+        elif 60 < int(ei) <= 120 and ei not in partition["60-120"]:
+            partition["60-120"].append(ei)
+        elif 120 < int(ei) <= 180 and ei not in partition["120-180"]:
+            partition["120-180"].append(ei)
+    old_event_ids = [partition["0-60"], partition["60-120"], partition["120-180"]]
+    new_event_id = [{"<60": 100}, {"<120": 101}, {"<180": 102}]
+    for old, new in zip(old_event_ids, new_event_id):
+        epochs = mne.epochs.combine_event_ids(epochs, old, new)
+
+    evokeds = {
+        "dir_60": epochs["<60"].average(),
+        "dir_120": epochs["<120"].average(),
+        "dir_180": epochs["<180"].average(),
+    }
+    for elec in ["AF7", "AF8", "TP9", "TP10"]:
+        nave = ", ".join([str(e.nave) for e in evokeds.values()])
+        mne.viz.plot_compare_evokeds(
+            evokeds,
+            picks=[elec],
+            cmap="viridis",
+            title=f"{elec} - nave: {nave}",
+            show=False,
+        )
+        plt.savefig(f"{subject}-{elec}-{suffix}.png")
+        plt.close()
+    return epochs
+
+
+def plot_ep_img_3_dir(epochs):
+    elecs = ["AF7", "AF8", "TP9", "TP10"]
+    for evo in ["<60", "<120", "<180"]:
+        fig = mne.viz.plot_epochs_image(epochs[evo], vmin=0, vmax=20, show=False)
+        fig[0].savefig(f"{subject}-{evo}-erp-gfp.png")
+        fig = mne.viz.plot_epochs_image(
+            epochs[evo], vmin=-10, vmax=10, show=False, picks=elecs, combine=None
+        )
+        for e, f in zip(elecs, fig):
+            f.savefig(f"{subject}-{evo}-{e}-erp-.png")
+    plt.close("all")
+
+
+for fname, subject in zip(fnames, snames):
+    store = pd.HDFStore(fname, "r")
+
+    # for key in store.keys():
+    #     print(key)
+
+    events = store.select("events")
+    steps = events.loc[(events["label"] == "step") | (events["label"] == "key")]
+    onsets = steps.index.values.astype(np.int64) * 1e-9
+    angles = []
+    for l, d in zip(events.label, events.data):
+        if l == "grid":
+            d = json.loads(d)
+            rows, target, cursor_curr = d["rows"], d["target"], d["cursor"]
+            continue
+        elif l == "step":
+            cursor_prev, cursor_curr = cursor_curr, json.loads(d)["cursor"]
+        elif l == "key":
+            d = json.loads(d)
+            if d["key"]:
+                angles.append("dirOk")  # , subject typed stim go toward target
+            else:
+                angles.append("dirNOk")  # , subject typed stim go away
+            continue
+        else:
+            continue
+
+        get_x = lambda i: (i - 1) % rows
+        get_y = lambda i: i // rows
+        get_point = lambda i: array([get_x(i), get_y(i)])
+        get_vec = lambda i, j: get_point(j) - get_point(i)
+        get_ang = lambda i, j, k: degrees(
+            arccos(
+                (get_vec(i, j) @ get_vec(i, k))
+                / (norm(get_vec(i, j)) * norm(get_vec(i, k)))
+            )
+        )
+
+        theta = get_ang(cursor_prev, target, cursor_curr)
+        angles.append(f"{theta:.0f}")
+    annotations = mne.Annotations(onsets, 0, angles, orig_time=0)
+
+    signal = store.select("eeg")
+    channels = list(signal.columns.values)
+    times = signal.index.values.astype(np.int64) * 1e-9
+    data = signal.values.T * 1e-6
+    rate = store.get_node("eeg")._v_attrs["meta"]["rate"]
+    info = mne.create_info(
+        ch_names=channels, sfreq=rate, ch_types=4 * ["eeg"] + ["misc"]
+    )
+    info.set_montage("standard_1020")
+    raw = mne.io.RawArray(data, info)
+    raw._filenames = [
+        fname
+    ]  # Hotfix: see https://github.com/mne-tools/mne-python/issues/9385
+    raw.set_meas_date(times[0])
+    raw.set_annotations(annotations)
+
+    raw.notch_filter(50, method="iir")
+    raw.filter(l_freq, h_freq, method="iir")
+    events, event_id = mne.events_from_annotations(raw)
+    epochs = mne.Epochs(
+        raw,
+        events=events,
+        event_id=event_id,
+        tmin=tmin,
+        tmax=0.5,
+        baseline=baseline,
+        reject=None,  # {'eeg': 90e-6},
+        preload=True,
+        verbose=False,
+    )
+    print(f"Dropped: ", (1 - len(epochs.events) / len(events)) * 100, "%")
+    epochs.plot_drop_log()
+
+    # remove EOG with an offline method: WARNING, won't work online!
+    # WARNING 2: frontal and temporal electrodes are interverted, removing EOG with TP10
+    _, betas = mne.preprocessing.regress_artifact(
+        epochs.copy().subtract_evoked(),
+        picks=["AF7", "TP9", "AF8"],
+        picks_artifact=["TP10"],
+    )
+    epochs, _ = mne.preprocessing.regress_artifact(
+        epochs, betas=betas, picks=["AF7", "TP9", "AF8"], picks_artifact=["TP10"]
+    )
+
+    epochs_orig = epochs.copy()
+    _ = plot_evoked_5_dir(epochs)
+
+    epochs = epochs_orig.copy()
+    epochs = plot_evoked_3_dir(epochs)
+    plot_ep_img_3_dir(epochs)
+
+    epochs = epochs_orig.copy().pick_types(eeg=True)
+    # epochs = epochs[['dirNOk', 'dirOk']]
+    xd = mne.preprocessing.Xdawn(n_components=2, reg="ledoit_wolf")
+    y = [
+        1 if e == epochs.event_id["dirOk"] else 0
+        for e in epochs["dirNOk", "dirOk"].events[:, 2]
+    ]
+    xd.fit(epochs["dirNOk", "dirOk"], y)
+    ep_xd = xd.apply(epochs)
+    _ = plot_evoked_3_dir(ep_xd["dirOk"], suffix="xdok")
+    _ = plot_evoked_3_dir(ep_xd["dirNOk"], suffix="xdnok")
+
+    # epochs.pick_types(eeg=True)
+    # xd = mne.preprocessing.Xdawn(n_components=2, reg="ledoit_wolf")
+    # y = [1 if e == 20 else 0 for e in epochs["<60", "<120"].events[:, 2]]
+    # xd.fit(epochs["<60", "<120"], y)
+    # ep_xd = xd.apply(epochs)
+    # elecs = ["AF7", "AF8", "TP9", "TP10"]
+    # for evo in ["<60", "<120"]:
+    #     fig = mne.viz.plot_epochs_image(ep_xd[evo], vmin=0, vmax=20, show=False)
+    #     fig[0].savefig(f"{subject}-{evo}-xdawn-erp-gfp.png")
+    #     fig = mne.viz.plot_epochs_image(
+    #         epochs[evo], vmin=-10, vmax=10, show=False, picks=elecs, combine=None
+    #     )
+    #     for e, f in zip(elecs, fig):
+    #         f.savefig(f"{subject}-{evo}-xdawn-erp-{e}.png")
+    # plt.close("all")

implicit/frustration/gui/README.md

@@ -0,0 +1,15 @@
+## Events sent to stream `events`
+
+| Label        | Data                                                                                    |
+| ------------ | --------------------------------------------------------------------------------------- |
+| `calibrated` | `null`                                                                                  |
+| `start`      | `null`                                                                                  |
+| `end`        | `null`                                                                                  |
+| `suspend`    | `null`                                                                                  |
+| `running`    | `null`                                                                                  |
+| `grid`       | `{"rows": <number>, "columns": <number>, "target": <number id>, "cursor": <number id>}` |
+| `step`       | `{"angle": <number>, "closer": <boolean>, "cursor": <number id>, "win": <bool>}`        |
+
+## `Cells IDs` on the grid
+
+![alt text](./img/cells_ids.png)

implicit/frustration/gui/assets/css/app.css

@@ -0,0 +1,33 @@
+#symbols {
+    display: flex;
+    align-items: center;
+    justify-content: center;
+}
+
+#footer {
+    align-items: center;
+    display: flex;
+    justify-content: center;
+    flex-direction: column;
+}
+
+#instructions {
+    font-size: 20px;
+}
+
+#button {
+    height: 10%;
+    width: 10%;
+    background-color: black;
+    border: none;
+    color: white;
+    padding: 16px 32px;
+    text-align: center;
+    font-size: 19px;
+    margin: 4px 2px;
+    opacity: 0.8;
+    transition: 0.3s;
+    display: inline-block;
+    text-decoration: none;
+    cursor: pointer;
+}

implicit/frustration/gui/assets/js/app.js

@@ -0,0 +1,48 @@
+"use strict";
+
+load_settings().then((settings) => {
+  let engine = new Engine(settings);
+  engine.io.subscribe("model");
+  engine.io.on("model", on_model);
+
+  let instructions = document.getElementById("instructions");
+  let button = document.getElementById("button");
+  button.addEventListener("click", on_click);
+  instructions.innerHTML =
+    "Press the button when the eeg headset is calibrated";
+  button.innerHTML = "Calibrated";
+
+  async function on_click() {
+    switch (engine.status) {
+      case "ready":
+        instructions.innerHTML = "Press start when you are ready";
+        button.innerHTML = "Start";
+        await engine.calibrated();
+        break;
+      case "calibrated":
+        button.innerHTML = "Suspend";
+        engine.run(instructions);
+        break;
+      case "running":
+        instructions.innerHTML = "Suspend...";
+        button.innerHTML = "Continue";
+        engine.suspend();
+        break;
+      case "suspend":
+        instructions.innerHTML = "Press Green or Red";
+        button.innerHTML = "Suspend";
+        engine.continue();
+        break;
+    }
+  }
+
+  async function on_model(data, meta) {
+    for (let row of Object.values(data)) {
+      switch (
+        row.label
+        // pass
+      ) {
+      }
+    }
+  }
+});