{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n\n# Kernel OPM phantom data\n\nIn this dataset, a Neuromag phantom was placed inside the Kernel OPM helmet and\nstimulated with 7 modules active (121 channels). Here we show some example traces.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Authors: The MNE-Python contributors.\n# License: BSD-3-Clause\n# Copyright the MNE-Python contributors.\n\nimport numpy as np\n\nimport mne\n\ndata_path = mne.datasets.phantom_kernel.data_path()\nfname = data_path / \"phantom_32_100nam_raw.fif\"\nraw = mne.io.read_raw_fif(fname).load_data()\nevents = mne.find_events(raw, stim_channel=\"STI101\")\n\n# Bads identified by inspecting averages\nraw.info[\"bads\"] = [\n    \"RC2.bx.ave\",\n    \"LC3.bx.ave\",\n    \"RC2.by.7\",\n    \"RC2.bz.7\",\n    \"RC2.bx.4\",\n    \"RC2.by.4\",\n    \"LC3.bx.5\",\n]\n# Drop the module-average channels\nraw.drop_channels([ch_name for ch_name in raw.ch_names if \".ave\" in ch_name])\n# Add field correction projectors\nraw.add_proj(mne.preprocessing.compute_proj_hfc(raw.info, order=2))\nraw.pick(\"meg\", exclude=\"bads\")\nraw.filter(0.5, 40)\nepochs = mne.Epochs(\n    raw,\n    events,\n    tmin=-0.1,\n    tmax=0.25,\n    decim=5,\n    preload=True,\n    baseline=(None, 0),\n)\nevoked = epochs[\"17\"].average()  # a high-SNR dipole for these data\nfig = evoked.plot()\nt_peak = 0.016  # based on visual inspection of evoked\nfig.axes[0].axvline(t_peak, color=\"k\", ls=\":\", lw=3, zorder=2)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The data covariance has an interesting structure because of densely packed sensors:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "cov = mne.compute_covariance(epochs, tmax=-0.01)\nmne.viz.plot_cov(cov, raw.info)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "So let's be careful and compute rank ahead of time and regularize:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "rank = mne.compute_rank(epochs, tol=1e-3, tol_kind=\"relative\")\ncov = mne.compute_covariance(epochs, tmax=-0.01, rank=rank, method=\"shrunk\")\nmne.viz.plot_cov(cov, raw.info)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Look at our alignment:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "sphere = mne.make_sphere_model(r0=(0.0, 0.0, 0.0), head_radius=0.08)\ntrans = mne.transforms.Transform(\"head\", \"mri\", np.eye(4))\nalign_kwargs = dict(\n    trans=trans,\n    bem=sphere,\n    surfaces={\"outer_skin\": 0.2},\n    show_axes=True,\n)\nmne.viz.plot_alignment(\n    raw.info,\n    coord_frame=\"meg\",\n    meg=dict(sensors=1.0, helmet=0.05),\n    **align_kwargs,\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's do dipole fits, which are not great because the dev_head_t is approximate and\nthe sensor coverage is sparse:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "data = list()\nfor ii in range(1, 33):\n    evoked = epochs[str(ii)][1:-1].average().crop(t_peak, t_peak)\n    data.append(evoked.data[:, 0])\nevoked = mne.EvokedArray(np.array(data).T, evoked.info, tmin=0.0)\ndel epochs\ndip, residual = mne.fit_dipole(evoked, cov, sphere, n_jobs=None)\nactual_pos, actual_ori = mne.dipole.get_phantom_dipoles()\nactual_amp = np.ones(len(dip))  # fake amp, needed to create Dipole instance\nactual_gof = np.ones(len(dip))  # fake GOF, needed to create Dipole instance\ndip_true = mne.Dipole(dip.times, actual_pos, actual_amp, actual_ori, actual_gof)\n\nfig = mne.viz.plot_alignment(\n    evoked.info, coord_frame=\"head\", meg=\"sensors\", **align_kwargs\n)\nmne.viz.plot_dipole_locations(\n    dipoles=dip_true, mode=\"arrow\", color=(0.0, 0.0, 0.0), fig=fig\n)\nmne.viz.plot_dipole_locations(dipoles=dip, mode=\"arrow\", color=(0.2, 1.0, 0.5), fig=fig)\nmne.viz.set_3d_view(figure=fig, azimuth=30, elevation=70, distance=0.4)"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "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.12.2"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}