.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "auto_tutorials/tutorial_04_sliding_win_model_comparison.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        :ref:`Go to the end <sphx_glr_download_auto_tutorials_tutorial_04_sliding_win_model_comparison.py>`
        to download the full example code

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_auto_tutorials_tutorial_04_sliding_win_model_comparison.py:


Sliding time window laminar model comparison using free energy
==============================================================

This tutorial demonstrates how to perform laminar inference of event-related responses in a sliding time window using model comparison based on free energy as a metric of model fit, described in [Bonaiuto et al., 2021, Laminar dynamics of high amplitude beta bursts in human motor cortex](https://doi.org/10.1016/j.neuroimage.2021.118479). A temporal Gaussian function is simulated at a particular cortical location in various layers. Source reconstruction is performed on the whole time window using the Empirical Bayesian Beamformer on the simulated sensor data using a forward model based on the multilayer mesh as a localizer. This is used to select priors on each layer mesh for a sliding time window model comparison using free energy.

.. GENERATED FROM PYTHON SOURCE LINES 9-13

Setting up the simulations
--------------------------

Simulations are based on an existing dataset, which is used to define the sampling rate, number of trials, duration of each trial, and the channel layout.

.. GENERATED FROM PYTHON SOURCE LINES 13-48

.. code-block:: default



    import os
    import shutil
    import numpy as np
    import matplotlib.pyplot as plt
    import tempfile

    from lameg.invert import coregister, invert_ebb
    from lameg.source import load_source_time_series
    from lameg.laminar import sliding_window_model_comparison
    from lameg.simulate import run_current_density_simulation
    from lameg.surf import LayerSurfaceSet
    from lameg.util import get_fiducial_coords
    from lameg.viz import show_surface, color_map
    import spm_standalone

    # Subject information for data to base the simulations on
    subj_id = 'sub-104'
    ses_id = 'ses-01'

    # Fiducial coil coordinates
    fid_coords = get_fiducial_coords(subj_id, '../test_data/participants.tsv')

    # Data file to base simulations on
    data_file = os.path.join(
        '../test_data',
        subj_id,
        'meg',
        ses_id,
        f'spm/pspm-converted_autoreject-{subj_id}-{ses_id}-001-btn_trial-epo.mat'
    )

    spm = spm_standalone.initialize()


.. GENERATED FROM PYTHON SOURCE LINES 49-50

For source reconstructions, we need an MRI and a surface mesh. The simulations will be based on a forward model using the multilayer mesh, and the model comparison will use each layer mesh

.. GENERATED FROM PYTHON SOURCE LINES 50-54

.. code-block:: default


    surf_set_bilam = LayerSurfaceSet(subj_id, 2)
    surf_set = LayerSurfaceSet(subj_id, 11)


.. GENERATED FROM PYTHON SOURCE LINES 55-56

We're going to copy the data file to a temporary directory and direct all output there.

.. GENERATED FROM PYTHON SOURCE LINES 56-77

.. code-block:: default


    # Extract base name and path of data file
    data_path, data_file_name = os.path.split(data_file)
    data_base = os.path.splitext(data_file_name)[0]

    # Where to put simulated data
    tmp_dir = tempfile.mkdtemp()

    # Copy data files to tmp directory
    shutil.copy(
        os.path.join(data_path, f'{data_base}.mat'), 
        os.path.join(tmp_dir, f'{data_base}.mat')
    )
    shutil.copy(
        os.path.join(data_path, f'{data_base}.dat'), 
        os.path.join(tmp_dir, f'{data_base}.dat')
    )

    # Construct base file name for simulations
    base_fname = os.path.join(tmp_dir, f'{data_base}.mat')


.. GENERATED FROM PYTHON SOURCE LINES 78-79

Invert the subject's data using the multilayer mesh. This step only has to be done once - this is just to compute the forward model that will be used in the simulations

.. GENERATED FROM PYTHON SOURCE LINES 79-102

.. code-block:: default


    # Patch size to use for inversion (in this case it matches the simulated patch size)
    patch_size = 5
    # Number of temporal modes to use for EBB inversion
    n_temp_modes = 4

    # Coregister data to multilayer mesh
    coregister(
        fid_coords,
        base_fname,
        surf_set,
        spm_instance=spm
    )

    # Run inversion
    [_,_] = invert_ebb(
        base_fname,
        surf_set,
        patch_size=patch_size, 
        n_temp_modes=n_temp_modes,
        spm_instance=spm
    )


.. GENERATED FROM PYTHON SOURCE LINES 103-106

Simulating a signal on the pial surface
---------------------------------------
We're going to simulate 200ms of a Gaussian with a dipole moment of 5nAm and a width of 25ms

.. GENERATED FROM PYTHON SOURCE LINES 106-122

.. code-block:: default


    # Strength of simulated activity (nAm)
    dipole_moment = 10
    # Temporal width of the simulated Gaussian
    signal_width=.025 # 25ms
    # Sampling rate (must match the data file)
    s_rate = 600

    # Generate 200ms of a Gaussian at a sampling rate of 600Hz (to match the data file)
    time=np.linspace(0,.2,121)
    zero_time=time[int((len(time)-1)/2+1)]
    sim_signal=np.exp(-((time-zero_time)**2)/(2*signal_width**2)).reshape(1,-1)
    plt.plot(time,dipole_moment*sim_signal[0,:])
    plt.xlabel('Time (s)')
    plt.ylabel('Amplitude (nAm)')


.. GENERATED FROM PYTHON SOURCE LINES 123-126

.. image:: ../_static/tutorial_04_sim_signal.png
   :width: 800
   :alt:

.. GENERATED FROM PYTHON SOURCE LINES 128-129

We need to pick a location (mesh vertex) to simulate at

.. GENERATED FROM PYTHON SOURCE LINES 129-143

.. code-block:: default


    # Vertex to simulate activity at
    sim_vertex=10561

    cam_view = [40, -240, 25,
                60, 37, 17,
                0, 0, 1]
    plot = show_surface(
        surf_set,
        marker_vertices=sim_vertex,
        marker_size=5,
        camera_view=cam_view
    )


.. GENERATED FROM PYTHON SOURCE LINES 144-147

.. image:: ../_static/tutorial_04_sim_location.png
   :width: 800
   :alt:

.. GENERATED FROM PYTHON SOURCE LINES 149-150

We'll simulate a 5mm patch of activity with -5 dB SNR at the sensor level. The desired level of SNR is achieved by adding white noise to the projected sensor signals

.. GENERATED FROM PYTHON SOURCE LINES 150-172

.. code-block:: default


    # Simulate at a vertex on the pial surface
    pial_vertex = surf_set.get_multilayer_vertex('pial', sim_vertex)
    prefix = f'sim_{sim_vertex}_pial_'

    # Size of simulated patch of activity (mm)
    sim_patch_size = 5
    # SNR of simulated data (dB)
    SNR = -10

    # Generate simulated data
    pial_sim_fname = run_current_density_simulation(
        base_fname,
        prefix,
        pial_vertex,
        sim_signal,
        dipole_moment,
        sim_patch_size,
        SNR,
        spm_instance=spm
    )


.. GENERATED FROM PYTHON SOURCE LINES 173-176

Localizer inversion
-------------------
Now we'll run a source reconstruction using the multilayer mesh, extract the signal in the pial layer, and select a prior based on the peak.

.. GENERATED FROM PYTHON SOURCE LINES 176-197

.. code-block:: default


    [_,_,MU] = invert_ebb(
        pial_sim_fname,
        surf_set,
        patch_size=patch_size,
        n_temp_modes=n_temp_modes,
        return_mu_matrix=True,
        spm_instance=spm
    )

    verts_per_surf = surf_set.get_vertices_per_layer()
    layer_vertices = np.arange(verts_per_surf)
    layer_ts, time, ch_names = load_source_time_series(
        pial_sim_fname,
        mu_matrix=MU,
        vertices=layer_vertices
    )

    m_layer_max = np.max(np.mean(layer_ts,axis=-1),-1)
    prior = np.argmax(m_layer_max)


.. GENERATED FROM PYTHON SOURCE LINES 198-199

We can see that the prior is the same as the location we simulated at

.. GENERATED FROM PYTHON SOURCE LINES 199-225

.. code-block:: default


    # Plot colors and camera view
    max_abs = np.max(np.abs(m_layer_max))
    c_range = [-max_abs, max_abs]

    # Plot peak
    colors,_ = color_map(
        m_layer_max,
        "RdYlBu_r",
        c_range[0],
        c_range[1]
    )
    thresh_colors=np.ones((colors.shape[0],4))*255
    thresh_colors[:,:3]=colors
    thresh_colors[m_layer_max<np.percentile(m_layer_max,99.9),3]=0

    plot = show_surface(
        surf_set,
        vertex_colors=thresh_colors,
        info=True,
        camera_view=cam_view,
        marker_vertices=prior,
        marker_size=5,
        marker_color=[0,0,255]
    )


.. GENERATED FROM PYTHON SOURCE LINES 226-229

.. image:: ../_static/tutorial_04_localizer.png
   :width: 800
   :alt:

.. GENERATED FROM PYTHON SOURCE LINES 231-234

Sliding time window model comparison (pial - white matter)
----------------------------------------------------------
Now we can run sliding time window model comparison between source models based on the pial and white matter surfaces using free energy. Specifically, we'll look at the difference in free energy between the two models (pial - white matter), in sliding and overlapping windows of 16ms. The free energy difference (pial - white matter) should be positive (more model evidence for the pial surface model) because we simulated activity on the pial surface.

.. GENERATED FROM PYTHON SOURCE LINES 234-257

.. code-block:: default


    # Number of temporal models for sliding time window inversion
    sliding_n_temp_modes = 4
    # Size of sliding window (in ms)
    win_size = 50
    # Whether or not windows overlap
    win_overlap = True

    # Run sliding time window model comparison between the first layer (pial) and the last layer (white matter)
    [Fs,wois] = sliding_window_model_comparison(
        prior, 
        fid_coords,
        pial_sim_fname,
        surf_set_bilam,
        spm_instance=spm,
        invert_kwargs={
            'patch_size': patch_size, 
            'n_temp_modes': sliding_n_temp_modes,
            'win_size': win_size, 
            'win_overlap': win_overlap,
        }
    )


.. GENERATED FROM PYTHON SOURCE LINES 258-259

Plot difference in free energy over time (pial minus white) - this should be positive

.. GENERATED FROM PYTHON SOURCE LINES 259-263

.. code-block:: default

    plt.plot(np.mean(wois,axis=-1), Fs[0,:]-Fs[1,:])
    plt.xlabel('Time (ms)')
    plt.ylabel(r'$\Delta$F')


.. GENERATED FROM PYTHON SOURCE LINES 264-267

.. image:: ../_static/tutorial_04_pial_sim_results.png
   :width: 800
   :alt:

.. GENERATED FROM PYTHON SOURCE LINES 269-272

White matter surface simulation with pial - white matter sliding time window model comparison
---------------------------------------------------------------------------------------------
Let's simulate the same pattern of activity, in the same location, but on the white matter surface. This time, sliding time window model comparison should yield greater model evidence for the white matter surface, and therefore the difference in free energy (pial - white matter) should be negative.

.. GENERATED FROM PYTHON SOURCE LINES 272-304

.. code-block:: default


    # Simulate at the corresponding vertex on the white matter surface
    white_vertex = surf_set.get_multilayer_vertex('white', sim_vertex)
    prefix = f'sim_{sim_vertex}_white_'

    # Generate simulated data
    white_sim_fname = run_current_density_simulation(
        base_fname,
        prefix,
        white_vertex,
        sim_signal,
        dipole_moment,
        sim_patch_size,
        SNR,
        spm_instance=spm
    )

    # Run sliding time window model comparison between the first layer (pial) and the last layer (white matter)
    [Fs,wois] = sliding_window_model_comparison(
        prior,
        fid_coords,
        white_sim_fname,
        surf_set_bilam,
        spm_instance=spm,
        invert_kwargs={
            'patch_size': patch_size, 
            'n_temp_modes': sliding_n_temp_modes,
            'win_size': win_size, 
            'win_overlap': win_overlap,
        }
    )


.. GENERATED FROM PYTHON SOURCE LINES 305-306

Plot difference in free energy over time (pial minus white) - this should be negative

.. GENERATED FROM PYTHON SOURCE LINES 306-310

.. code-block:: default

    plt.plot(np.mean(wois,axis=-1), Fs[0,:]-Fs[1,:])
    plt.xlabel('Time (ms)')
    plt.ylabel(r'$\Delta$F')


.. GENERATED FROM PYTHON SOURCE LINES 311-314

.. image:: ../_static/tutorial_04_white_sim_results.png
   :width: 800
   :alt:

.. GENERATED FROM PYTHON SOURCE LINES 317-320

Simulation in each layer with sliding time window model comparison across layers
--------------------------------------------------------------------------------
That was sliding time window model comparison with two candidate models: one based on the white matter surface, and one on the pial. Let's now simulate on each layer, and for each simulation, run sliding time window model comparison across all layers. We'll turn off SPM visualization here.

.. GENERATED FROM PYTHON SOURCE LINES 320-358

.. code-block:: default


    # Now simulate at the corresponding vertex on each layer, and for each simulation, run sliding window model
    # comparison across all layers
    all_layerF = []
    for l in range(surf_set.n_layers):
        print(f'Simulating in layer {l}')
        l_vertex = surf_set.get_multilayer_vertex(l, sim_vertex)
        prefix = f'sim_{sim_vertex}{l}_'

        l_sim_fname = run_current_density_simulation(
            base_fname,
            prefix,
            l_vertex,
            sim_signal,
            dipole_moment,
            sim_patch_size,
            SNR,
            spm_instance=spm
        )

        [Fs, wois] = sliding_window_model_comparison(
            prior,
            fid_coords,
            l_sim_fname,
            surf_set,
            viz=False,
            spm_instance=spm,
            invert_kwargs={
                'patch_size': patch_size,
                'n_temp_modes': sliding_n_temp_modes,
                'win_size': win_size,
                'win_overlap': win_overlap,
            }
        )

        all_layerF.append(Fs)
    all_layerF = np.squeeze(np.array(all_layerF))


.. GENERATED FROM PYTHON SOURCE LINES 359-360

For each simulation, we can plot the free energy for all models relative to the worst model within a central time window. The layer model with the highest free energy should correspond to the layer that the activity was simulated in.

.. GENERATED FROM PYTHON SOURCE LINES 360-395

.. code-block:: default


    # Average free energy within small time window in center of the epoch
    woi_t = np.mean(wois,axis=-1)
    woi_idx = np.where((woi_t>=-20) & (woi_t<=20))[0]
    m_all_layerF = np.mean(all_layerF[:,:,woi_idx],axis=2)

    col_r = plt.cm.cool(np.linspace(0,1, num=surf_set.n_layers))
    plt.figure(figsize=(10,4))

    # For each simulation, plot the mean free energy of each layer model relative to that of the worst
    # model for that simulation
    plt.subplot(1,2,1)
    for l in range(surf_set.n_layers):
        layerF = m_all_layerF[l,:]
        plt.plot(layerF-np.min(layerF), label=f'{l}', color=col_r[l,:])
    plt.legend()
    plt.xlabel('Eval layer')
    plt.ylabel(r'$\Delta$F')

    # For each simulation, find which layer model had the greatest free energy
    plt.subplot(1,2,2)
    peaks=[]
    for l in range(surf_set.n_layers):
        layerF = m_all_layerF[l,:]
        layerF = layerF-np.min(layerF)
        pk = np.argmax(layerF)
        peaks.append(pk)
    plt.plot(peaks)
    plt.xlim([-0.5,10.5])
    plt.ylim([-0.5,10.5])
    plt.plot([0,10],[0,10],'k--')
    plt.xlabel('Sim layer')
    plt.ylabel(r'Peak $\Delta$F')
    plt.tight_layout()


.. GENERATED FROM PYTHON SOURCE LINES 396-399

.. image:: ../_static/tutorial_04_results.png
   :width: 800
   :alt:

.. GENERATED FROM PYTHON SOURCE LINES 399-420

.. code-block:: default


    # Normalization step
    norm_layerF = np.zeros(m_all_layerF.shape)
    for l in range(surf_set.n_layers):
        norm_layerF[l,:] = m_all_layerF[l,:] - np.min(m_all_layerF[l,:])

    # Transpose for visualization
    im=plt.imshow(norm_layerF.T, cmap='Spectral_r')

    # Find the indices of the max value in each column
    max_indices = np.argmax(norm_layerF, axis=1)

    # Plot an 'X' at the center of the square for each column's maximum
    for idx, max_idx in enumerate(max_indices):
        plt.text(idx, max_idx, 'X', fontsize=12, ha='center', va='center', color='black', weight='bold')

    plt.xlabel('Simulated layer', fontsize=14)
    plt.ylabel('Evaluated layer', fontsize=14)
    cb=plt.colorbar(im)
    cb.set_label(r'$\Delta F$', fontsize=14)


.. GENERATED FROM PYTHON SOURCE LINES 421-424

.. image:: ../_static/tutorial_04_results_matrix.png
   :width: 800
   :alt:

.. GENERATED FROM PYTHON SOURCE LINES 426-431

.. code-block:: default

    spm.terminate()

    # Delete simulation files
    shutil.rmtree(tmp_dir)



.. rst-class:: sphx-glr-timing

   **Total running time of the script:** (0 minutes 0.000 seconds)


.. _sphx_glr_download_auto_tutorials_tutorial_04_sliding_win_model_comparison.py:

.. only:: html

  .. container:: sphx-glr-footer sphx-glr-footer-example




    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: tutorial_04_sliding_win_model_comparison.py <tutorial_04_sliding_win_model_comparison.py>`

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: tutorial_04_sliding_win_model_comparison.ipynb <tutorial_04_sliding_win_model_comparison.ipynb>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_