"""
.. _main-tutorial:

============================
About Dipoles in MEG and EEG
============================

For an explanation of what is going on in the demo and background information
about magentoencephalography (MEG) and electroencephalography (EEG) in
general, let's walk through some code. To execute this code, you'll need
to have a working version of python with ``mne`` installed, see the
:ref:`quick-start` documentation for instructions. You'll need the development
version, so you'll need to do
``pip install git+https://github.com/mne-tools/mne-python.git`` You'll also
need to install the requirements such as with
``pip install -r requirements.txt``.
"""
# Author: <NAME> <<EMAIL>>
#
# License: BSD-3-Clause

# %%
# Let's start by importing the dependencies we'll need.

import os.path as op  # comes with python and helps naviagte to files
import numpy as np  # a scientific computing package with arrays
import matplotlib.pyplot as plt  # a plotting library
import mne  # our main analysis software package
from nilearn.plotting import plot_anat  # this package plots brains

# %%
# Background
# ----------
# MEG and EEG researchers record very small electromagentic potentials
# generated by the brain from outside the head. When it comes from the
# recording devices, it looks like this (there are a lot of channels
# so only a subset are shown):

data_path = mne.datasets.sample.data_path()  # get the sample data path
raw = mne.io.read_raw(  # navigate to some raw sample data
    op.join(data_path, 'MEG', 'sample', 'sample_audvis_raw.fif'))
raw_plot = raw.copy()  # make a copy to modify for plotting
raw_plot.pick_channels(raw.ch_names[::10])  # pick only every tenth channel
raw_plot.plot(n_channels=len(raw_plot.ch_names),
              duration=1,  # only a small, 1 second time window
              start=50,  # start part way in
              )

# %%
# The goal of MEG and EEG researchers is to try and understand how activity
# in the brain changes as we respond to stimuli in our environment and
# perform behaviors. To do that, researchers will often use magnetic resonance
# (MR) to create an image of the research subject's brain. These images
# look like this:

# first, get a T1-weighted MR scan file from the MNE example dataset
T1_fname = op.join(data_path, 'subjects', 'sample', 'mri', 'T1.mgz')
plot_anat(T1_fname)  # now we can plot it

# %%
# The T1 MR image can be used to figure out where the surfaces of the
# brain skull and scalp are as well as label the parts of the brain
# in the image using Freesurfer. The command below does this (it takes
# 8 hours so I wouldn't recommend executing it now but it has already
# been done for you in the mne sample data, see `here
# <https://surfer.nmr.mgh.harvard.edu/fswiki/DownloadAndInstall>`_ for
# how to install Freesurfer):
#
# .. code-block:: bash
#
#    recon-all -subjid sample -sd $DATA_PATH/subjects -i $T1_FNAME -all

# %%
# Now let's put it all together and see the problem that MEG and EEG
# researchers face in figuring out what's going on inside the brain from
# electromagnetic potentials on the surface of the scalp. As you can see below,
# there are a lot of MEG and EEG sensors and they cover a large portion of the
# head but its not readily apparent how much of each brain area each sensor
# records from and how to separate the summed activity from all the brain areas
# that is recorded by each sensor into components for each brain area:
#
# .. note::
#
#     The sensor positions don't come aligned to the MR image since they are
#     recorded by a different device so we need to a transformation matrix to
#     transform them from the coordinate frame they are in to MR coordinates.
#     This can be done with :func:`mne.gui.coregistration` to generate the
#     ``trans`` file that is loaded below.

# the subjects_dir is where Freesurfer stored all the surface files
subjects_dir = op.join(data_path, 'subjects')
trans = mne.read_trans(op.join(data_path, 'MEG', 'sample',
                               'sample_audvis_raw-trans.fif'))

# the main plotter for mne, the brain object
brain = mne.viz.Brain(subject_id='sample', hemi='both', surf='pial',
                      subjects_dir=subjects_dir)
brain.add_skull(alpha=0.5)  # alpha sets transparency
brain.add_head(alpha=0.5)
brain.add_sensors(raw.info, trans)
# set a nice view to show
brain.show_view(azimuth=120, elevation=90, distance=500)

# %%
# Making a Source Space and Forward Model
# ---------------------------------------
# First let's setup a space of vertices within the brain that we will consider
# as the sources of signal. In a real brain, there are hundreds of billions
# of cells but we don't have the resolution with only hundreds of sensors to
# determine the activity of each cell, so, instead, we'll choose a regularly
# sampled grid of sources that represent the summed activity of tens of
# thousands of cells. In most analyses in publications, the source space has
# around 8000 vertices, but, for this example, we'll use a smaller source
# space for demonstration.

# %%
# First, we would need to make a boundary element model (BEM) to account for
# differences in conductivity of the brain, skull and scalp. This can be
# done with :func:`mne.make_bem_model` but, in this case, we'll just load
# a pre-computed model. We'll also load the solution to the BEM model for how
# different conductivities of issues effect current dipoles as they pass
# through each of the layers, but this can be computed with
# :func:`mne.make_bem_solution`.

bem_fname = op.join(subjects_dir, 'sample', 'bem',
                    'sample-5120-5120-5120-bem.fif')

# load a pre-computed solution the how the sources within the brain will
# be affected by the different conductivities
bem_sol = op.join(subjects_dir, 'sample', 'bem',
                  'sample-5120-5120-5120-bem-sol.fif')
# plot it, it's saved out in a standard location,
# so we don't have to pass the path
mne.viz.plot_bem(subject='sample', subjects_dir=op.join(data_path, 'subjects'),
                 slices=np.linspace(45, 200, 12).round().astype(int))

# %%
# Making a Dipole
# ---------------
# Now, we're ready to make a dipole and see how its current will be recorded
# at the scalp with MEG and EEG.
#
# .. note::
#
#     You can use ``print(mne.Dipole.__doc__)`` to print the arguments that
#     are required by ``mne.Dipole`` or any other class, method or function.

# make a dipole within the temporal lobe pointing superiorly,
# fake a goodness-of-fit number
dip_pos = [-0.0647572, 0.01315963, 0.07091921]
dip = mne.Dipole(times=[0], pos=[dip_pos], amplitude=[3e-8],
                 ori=[[0, 0, 1]], gof=50)

# plot it!
brain = mne.viz.Brain(subject_id='sample', hemi='both', surf='pial',
                      subjects_dir=subjects_dir, alpha=0.25)
brain.add_dipole(dip, trans, scales=10)
brain.show_view(azimuth=150, elevation=60, distance=500)

# %%
# Simulating Sensor Data
# ----------------------
# We're ready to compute a forward operator using the BEM to make the so-called
# leadfield matrix which multiplies activity at the dipole to give the
# modelled the activity at the sensors. We can then use this to simulate evoked
# data.

fwd, stc = mne.make_forward_dipole(
    dipole=dip, bem=bem_sol, info=raw.info, trans=trans)
# we don't have a few things like the covarience matrix or a number of epochs
# to average so we use these arguments for a reasonable solution
evoked = mne.simulation.simulate_evoked(
    fwd, stc, raw.info, cov=None, nave=np.inf)

# Now we can see what it would look like at the sensors
fig, axes = plt.subplots(1, 3, figsize=(6, 4))  # make a figure with 3 subplots
# use zip to iterate over axes and channel types at the same time
for ax, ch_type in zip(axes, ('grad', 'mag', 'eeg')):
    # we're just looking at the relative pattern so we won't use a colorbar
    evoked.plot_topomap(times=[0], ch_type=ch_type, axes=ax, colorbar=False)
    ax.set_title(ch_type)

# %%
# Wrapping Up
# -----------
# We covered some good intuition but there's lots more to learn! The main thing
# is that MEG and EEG researchers generally don't have the information about
# what's going on inside the brain, that's what they are trying to predict. To
# reverse this process, you need to invert the forward solution (tutorial:
# :ref:`tut-viz-stcs`). There is tons more to explore in the MNE `tutorials
# <https://mne.tools/dev/auto_tutorials/index.html>`_ and `examples
# <https://mne.tools/dev/auto_examples/index.html>`_ pages. Let's leave off by
# setting up a source space of many different dipoles and seeing their
# different activities manifest on the scalp as measured by the sensors.

src = mne.setup_volume_source_space(
    subject='sample', pos=20,  # in mm
    bem=bem_fname, subjects_dir=subjects_dir)

# make the leadfield matrix
fwd = mne.make_forward_solution(
    raw.info, trans=trans, src=src, bem=bem_sol)

# plot our setup
brain = mne.viz.Brain(subject_id='sample', hemi='both', surf='pial',
                      subjects_dir=subjects_dir, alpha=0.25)
brain.add_volume_labels(alpha=0.25, colors='gray')
brain.add_forward(fwd, trans, scale=3)
brain.show_view(azimuth=30, elevation=90, distance=500)

# %%
# Plot the same solution using a source space of dipoles

# take the source space from the forward model because some of the
# vertices are excluded from the vertices in src
n_dipoles = fwd['source_rr'].shape[0]  # rr is the vertex positions
# find the closest dipole to the one we used before (it was in this
# source space) using the euclidean distance (np.linalg.norm)
idx = np.argmin(np.linalg.norm(fwd['source_rr'] - dip_pos, axis=1))
# make an empty matrix of zeros
data = np.zeros((n_dipoles, 3, 1))
data[idx, 2, 0] = 3e-8  # make the same dipole as before
# this is the format that vertiex numbers are stored in
vertices = [fwd['src'][0]['vertno']]
stc = mne.VolVectorSourceEstimate(data, vertices=vertices,
                                  subject='sample', tmin=0, tstep=1)

evoked = mne.simulation.simulate_evoked(
    fwd, stc, raw.info, cov=None, nave=np.inf)

# confirm our replication
fig, axes = plt.subplots(1, 3, figsize=(6, 4))  # make a figure with 3 subplots
for ax, ch_type in zip(axes, ('grad', 'mag', 'eeg')):
    evoked.plot_topomap(times=[0], ch_type=ch_type, axes=ax, colorbar=False)
    ax.set_title(ch_type)

# %%
# Now, go crazy and simulate a bunch of random dipoles

np.random.seed(88)  # always seed random number generation for reproducibility
stc.data = np.random.random(stc.data.shape) * 3e-8 - 1.5e-8
evoked = mne.simulation.simulate_evoked(
    fwd, stc, raw.info, cov=None, nave=np.inf)

# now that's a complicated faked brain pattern, fortunately brain activity
# is much more correlated (neighboring areas have similar activity) which
# makes results a bit easier to interpret
fig, axes = plt.subplots(1, 3, figsize=(6, 4))  # make a figure with 3 subplots
for ax, ch_type in zip(axes, ('grad', 'mag', 'eeg')):
    evoked.plot_topomap(times=[0], ch_type=ch_type, axes=ax, colorbar=False)
    ax.set_title(ch_type)
