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

.. only:: html

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

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

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

.. _sphx_glr_auto_examples_multicsc_plot_somato_mu_waves.py:


===========================================================
Extracting :math:`\mu`-wave from the somato-sensory dataset
===========================================================

This example illustrates how to learn rank-1 atoms [1]_ on the multivariate
somato-sensorymotor dataset from :code:`mne`. The displayed results highlight
the presence of :math:`\mu`-waves located in the SI cortex.

.. [1] Dupré La Tour, T., Moreau, T., Jas, M., & Gramfort, A. (2018).
    `Multivariate Convolutional Sparse Coding for Electromagnetic Brain Signals
    <https://arxiv.org/abs/1805.09654v2>`_. Advances in Neural Information
    Processing Systems (NIPS).

.. GENERATED FROM PYTHON SOURCE LINES 15-23

.. code-block:: Python


    # Authors: Thomas Moreau <thomas.moreau@inria.fr>
    #          Mainak Jas <mainak.jas@telecom-paristech.fr>
    #          Tom Dupre La Tour <tom.duprelatour@telecom-paristech.fr>
    #          Alexandre Gramfort <alexandre.gramfort@telecom-paristech.fr>
    #
    # License: BSD (3-clause)








.. GENERATED FROM PYTHON SOURCE LINES 24-25

Let us first define the parameters of our model.

.. GENERATED FROM PYTHON SOURCE LINES 25-32

.. code-block:: Python


    sfreq = 150.

    # Define the shape of the dictionary
    n_atoms = 25
    n_times_atom = int(round(sfreq * 1.0))  # 1000. ms








.. GENERATED FROM PYTHON SOURCE LINES 33-34

Next, we define the parameters for multivariate CSC

.. GENERATED FROM PYTHON SOURCE LINES 34-56

.. code-block:: Python


    from alphacsc import BatchCDL
    cdl = BatchCDL(
        # Shape of the dictionary
        n_atoms=n_atoms,
        n_times_atom=n_times_atom,
        # Request a rank1 dictionary with unit norm temporal and spatial maps
        rank1=True, uv_constraint='separate',
        # Initialize the dictionary with random chunk from the data
        D_init='chunk',
        # rescale the regularization parameter to be 20% of lambda_max
        lmbd_max="scaled", reg=.2,
        # Number of iteration for the alternate minimization and cvg threshold
        n_iter=100, eps=1e-4,
        # solver for the z-step
        solver_z="lgcd", solver_z_kwargs={'tol': 1e-2, 'max_iter': 1000},
        # solver for the d-step
        solver_d='alternate_adaptive', solver_d_kwargs={'max_iter': 300},
        # Technical parameters
        verbose=1, random_state=0, n_jobs=6)









.. GENERATED FROM PYTHON SOURCE LINES 57-60

Here, we load the data from the somato-sensory dataset and preprocess them
in epochs. The epochs are selected around the stim, starting 2 seconds
before and finishing 4 seconds after.

.. GENERATED FROM PYTHON SOURCE LINES 60-66

.. code-block:: Python


    from alphacsc.datasets.mne_data import load_data
    t_lim = (-2, 4)
    X, info = load_data(dataset='somato', epoch=t_lim, sfreq=sfreq)






.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    Using default location ~/mne_data for somato...

      0%|                                               | 0.00/611M [00:00<?, ?B/s]
      1%|▍                                     | 6.32M/611M [00:00<00:09, 63.1MB/s]
      2%|▊                                     | 13.7M/611M [00:00<00:08, 69.4MB/s]
      3%|█▎                                    | 21.1M/611M [00:00<00:08, 71.4MB/s]
      5%|█▊                                    | 28.4M/611M [00:00<00:08, 72.3MB/s]
      6%|██▏                                   | 35.8M/611M [00:00<00:07, 72.7MB/s]
      7%|██▋                                   | 43.1M/611M [00:00<00:07, 72.9MB/s]
      8%|███▏                                  | 50.5M/611M [00:00<00:07, 73.1MB/s]
      9%|███▌                                  | 57.8M/611M [00:00<00:09, 60.3MB/s]
     11%|████                                  | 65.1M/611M [00:00<00:08, 63.9MB/s]
     12%|████▍                                 | 72.2M/611M [00:01<00:08, 65.6MB/s]
     13%|████▉                                 | 78.9M/611M [00:01<00:09, 54.3MB/s]
     14%|█████▎                                | 86.0M/611M [00:01<00:08, 58.5MB/s]
     15%|█████▊                                | 93.4M/611M [00:01<00:08, 62.5MB/s]
     17%|██████▍                                | 101M/611M [00:01<00:07, 65.6MB/s]
     18%|██████▉                                | 108M/611M [00:01<00:07, 67.9MB/s]
     19%|███████▎                               | 115M/611M [00:01<00:12, 41.2MB/s]
     20%|███████▋                               | 121M/611M [00:02<00:12, 39.0MB/s]
     21%|████████▏                              | 127M/611M [00:02<00:10, 44.6MB/s]
     22%|████████▍                              | 133M/611M [00:02<00:16, 29.3MB/s]
     23%|████████▊                              | 138M/611M [00:02<00:14, 32.9MB/s]
     23%|█████████                              | 143M/611M [00:02<00:15, 29.6MB/s]
     24%|█████████▎                             | 147M/611M [00:03<00:20, 23.2MB/s]
     25%|█████████▋                             | 152M/611M [00:03<00:15, 29.0MB/s]
     26%|██████████                             | 158M/611M [00:03<00:13, 33.6MB/s]
     27%|██████████▎                            | 162M/611M [00:03<00:12, 35.6MB/s]
     27%|██████████▋                            | 166M/611M [00:03<00:19, 22.9MB/s]
     28%|██████████▊                            | 170M/611M [00:04<00:24, 18.1MB/s]
     29%|███████████▎                           | 177M/611M [00:04<00:16, 26.3MB/s]
     30%|███████████▌                           | 182M/611M [00:04<00:14, 29.8MB/s]
     30%|███████████▉                           | 186M/611M [00:04<00:18, 22.8MB/s]
     31%|████████████▎                          | 192M/611M [00:04<00:14, 28.9MB/s]
     32%|████████████▌                          | 198M/611M [00:04<00:11, 34.5MB/s]
     33%|████████████▉                          | 202M/611M [00:05<00:13, 30.3MB/s]
     34%|█████████████▏                         | 206M/611M [00:05<00:13, 30.4MB/s]
     35%|█████████████▋                         | 214M/611M [00:05<00:09, 41.8MB/s]
     36%|██████████████▏                        | 223M/611M [00:05<00:07, 51.2MB/s]
     38%|██████████████▊                        | 231M/611M [00:05<00:06, 59.8MB/s]
     39%|███████████████▎                       | 239M/611M [00:05<00:05, 64.4MB/s]
     40%|███████████████▋                       | 246M/611M [00:05<00:07, 49.6MB/s]
     41%|████████████████                       | 252M/611M [00:05<00:07, 47.1MB/s]
     42%|████████████████▌                      | 259M/611M [00:06<00:06, 52.0MB/s]
     44%|████████████████▉                      | 266M/611M [00:06<00:06, 57.2MB/s]
     45%|█████████████████▍                     | 273M/611M [00:06<00:05, 61.4MB/s]
     46%|█████████████████▉                     | 281M/611M [00:06<00:05, 64.8MB/s]
     47%|██████████████████▎                    | 288M/611M [00:06<00:06, 47.7MB/s]
     48%|██████████████████▋                    | 293M/611M [00:06<00:06, 48.2MB/s]
     49%|███████████████████                    | 299M/611M [00:06<00:06, 46.6MB/s]
     50%|███████████████████▍                   | 304M/611M [00:07<00:06, 44.3MB/s]
     51%|███████████████████▉                   | 311M/611M [00:07<00:05, 51.9MB/s]
     52%|████████████████████▍                  | 319M/611M [00:07<00:05, 56.3MB/s]
     53%|████████████████████▊                  | 325M/611M [00:07<00:07, 37.7MB/s]
     54%|█████████████████████                  | 330M/611M [00:07<00:10, 27.0MB/s]
     55%|█████████████████████▍                 | 335M/611M [00:07<00:09, 30.1MB/s]
     55%|█████████████████████▋                 | 339M/611M [00:08<00:13, 20.5MB/s]
     56%|█████████████████████▊                 | 342M/611M [00:08<00:16, 16.4MB/s]
     57%|██████████████████████▏                | 347M/611M [00:08<00:12, 21.2MB/s]
     58%|██████████████████████▌                | 352M/611M [00:08<00:09, 26.5MB/s]
     59%|██████████████████████▉                | 360M/611M [00:09<00:07, 35.4MB/s]
     60%|███████████████████████▍               | 367M/611M [00:09<00:05, 42.8MB/s]
     61%|███████████████████████▉               | 375M/611M [00:09<00:04, 50.5MB/s]
     63%|████████████████████████▍              | 382M/611M [00:09<00:04, 57.0MB/s]
     64%|████████████████████████▉              | 390M/611M [00:09<00:03, 62.4MB/s]
     65%|█████████████████████████▍             | 397M/611M [00:09<00:03, 62.5MB/s]
     66%|█████████████████████████▊             | 404M/611M [00:09<00:06, 33.8MB/s]
     67%|██████████████████████████▎            | 411M/611M [00:10<00:05, 39.7MB/s]
     68%|██████████████████████████▋            | 418M/611M [00:10<00:04, 44.8MB/s]
     69%|███████████████████████████            | 424M/611M [00:10<00:05, 33.0MB/s]
     71%|███████████████████████████▌           | 431M/611M [00:10<00:04, 40.7MB/s]
     72%|███████████████████████████▉           | 438M/611M [00:10<00:03, 45.4MB/s]
     73%|████████████████████████████▎          | 444M/611M [00:10<00:03, 49.1MB/s]
     74%|████████████████████████████▊          | 451M/611M [00:10<00:02, 54.6MB/s]
     75%|█████████████████████████████▎         | 458M/611M [00:10<00:02, 59.1MB/s]
     76%|█████████████████████████████▋         | 465M/611M [00:11<00:02, 62.6MB/s]
     77%|██████████████████████████████▏        | 473M/611M [00:11<00:02, 65.1MB/s]
     79%|██████████████████████████████▋        | 480M/611M [00:11<00:01, 67.0MB/s]
     80%|███████████████████████████████        | 487M/611M [00:11<00:01, 68.7MB/s]
     81%|███████████████████████████████▌       | 494M/611M [00:11<00:01, 68.9MB/s]
     82%|████████████████████████████████       | 501M/611M [00:11<00:01, 67.3MB/s]
     83%|████████████████████████████████▍      | 508M/611M [00:12<00:03, 30.8MB/s]
     84%|████████████████████████████████▊      | 513M/611M [00:12<00:03, 27.6MB/s]
     85%|█████████████████████████████████▏     | 520M/611M [00:12<00:02, 33.3MB/s]
     86%|█████████████████████████████████▋     | 527M/611M [00:12<00:02, 40.5MB/s]
     88%|██████████████████████████████████▏    | 534M/611M [00:12<00:01, 47.3MB/s]
     89%|██████████████████████████████████▌    | 540M/611M [00:12<00:02, 31.7MB/s]
     89%|██████████████████████████████████▊    | 545M/611M [00:13<00:01, 33.4MB/s]
     90%|███████████████████████████████████▏   | 552M/611M [00:13<00:01, 39.3MB/s]
     91%|███████████████████████████████████▌   | 557M/611M [00:13<00:01, 28.7MB/s]
     92%|███████████████████████████████████▉   | 563M/611M [00:13<00:01, 34.6MB/s]
     93%|████████████████████████████████████▍  | 570M/611M [00:13<00:01, 40.7MB/s]
     95%|████████████████████████████████████▊  | 577M/611M [00:13<00:00, 47.5MB/s]
     96%|█████████████████████████████████████▎ | 584M/611M [00:13<00:00, 53.4MB/s]
     97%|█████████████████████████████████████▋ | 590M/611M [00:14<00:00, 32.2MB/s]
     98%|██████████████████████████████████████ | 595M/611M [00:14<00:00, 21.0MB/s]
     98%|██████████████████████████████████████▎| 599M/611M [00:15<00:00, 18.5MB/s]
     99%|██████████████████████████████████████▌| 604M/611M [00:15<00:00, 22.3MB/s]
    100%|██████████████████████████████████████▊| 608M/611M [00:15<00:00, 26.2MB/s]
      0%|                                               | 0.00/611M [00:00<?, ?B/s]
    100%|███████████████████████████████████████| 611M/611M [00:00<00:00, 2.40TB/s]
    Attempting to create new mne-python configuration file:
    /github/home/.mne/mne-python.json
    Download complete in 28s (582.2 MB)
    Opening raw data file /github/home/mne_data/MNE-somato-data/sub-01/meg/sub-01_task-somato_meg.fif...
        Range : 237600 ... 506999 =    791.189 ...  1688.266 secs
    Ready.
    Reading 0 ... 269399  =      0.000 ...   897.077 secs...
    Filtering raw data in 1 contiguous segment
    Setting up band-stop filter

    FIR filter parameters
    ---------------------
    Designing a one-pass, zero-phase, non-causal bandstop filter:
    - Windowed time-domain design (firwin) method
    - Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation
    - Lower transition bandwidth: 0.50 Hz
    - Upper transition bandwidth: 0.50 Hz
    - Filter length: 1983 samples (6.603 s)

    Filtering raw data in 1 contiguous segment
    Setting up high-pass filter at 2 Hz

    FIR filter parameters
    ---------------------
    Designing a one-pass, zero-phase, non-causal highpass filter:
    - Windowed time-domain design (firwin) method
    - Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation
    - Lower passband edge: 2.00
    - Lower transition bandwidth: 2.00 Hz (-6 dB cutoff frequency: 1.00 Hz)
    - Filter length: 497 samples (1.655 s)

    111 events found on stim channel STI 014
    Event IDs: [1]
    Not setting metadata
    111 matching events found
    Setting baseline interval to [-3.9992341833870637, 0.0] s
    Applying baseline correction (mode: mean)
    0 projection items activated
    Using data from preloaded Raw for 111 events and 1202 original time points ...
        Rejecting  epoch based on EOG : ['EOG 061']
        Rejecting  epoch based on EOG : ['EOG 061']
        Rejecting  epoch based on EOG : ['EOG 061']
        Rejecting  epoch based on EOG : ['EOG 061']
        Rejecting  epoch based on EOG : ['EOG 061']
    5 bad epochs dropped
    NOTE: pick_types() is a legacy function. New code should use inst.pick(...).
    /github/workspace/alphacsc/datasets/mne_data.py:98: RuntimeWarning: Something went wrong in the data-driven estimation of the data rank as it exceeds the theoretical rank from the info (204 > 64). Consider setting rank to "auto" or setting it explicitly as an integer.
      cov = mne.compute_covariance(epochs_cov)
    Reducing data rank from 204 -> 204
    Estimating covariance using EMPIRICAL
    Done.
    Number of samples used : 127412
    [done]
    Not setting metadata
    111 matching events found
    Setting baseline interval to [-2.001282051803185, 0.0] s
    Applying baseline correction (mode: mean)
    0 projection items activated
    Using data from preloaded Raw for 111 events and 1803 original time points ...
        Rejecting  epoch based on EOG : ['EOG 061']
        Rejecting  epoch based on EOG : ['EOG 061']
        Rejecting  epoch based on EOG : ['EOG 061']
        Rejecting  epoch based on EOG : ['EOG 061']
        Rejecting  epoch based on EOG : ['EOG 061']
        Rejecting  epoch based on EOG : ['EOG 061']
        Rejecting  epoch based on EOG : ['EOG 061']
        Rejecting  epoch based on EOG : ['EOG 061']
    8 bad epochs dropped
    NOTE: pick_types() is a legacy function. New code should use inst.pick(...).




.. GENERATED FROM PYTHON SOURCE LINES 67-68

Fit the model and learn rank1 atoms

.. GENERATED FROM PYTHON SOURCE LINES 68-70

.. code-block:: Python

    cdl.fit(X)





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    ..............
    [BatchCDL] Converged after 14 iteration, (dz, du) = 8.765e-05, 8.726e-05
    [BatchCDL] Fit in 271.3s

    <alphacsc.convolutional_dictionary_learning.BatchCDL object at 0x7f6121a69040>



.. GENERATED FROM PYTHON SOURCE LINES 71-73

Display the 4-th atom, which displays a :math:`\mu`-waveform in its temporal
pattern.

.. GENERATED FROM PYTHON SOURCE LINES 73-108

.. code-block:: Python


    import mne
    import numpy as np
    import matplotlib.pyplot as plt

    i_atom = 4
    n_plots = 3
    figsize = (n_plots * 5, 5.5)
    fig, axes = plt.subplots(1, n_plots, figsize=figsize, squeeze=False)

    # Plot the spatial map of the learn atom using mne topomap
    ax = axes[0, 0]
    u_hat = cdl.u_hat_[i_atom]
    mne.viz.plot_topomap(u_hat, info, axes=ax, show=False)
    ax.set(title='Learned spatial pattern')

    # Plot the temporal pattern of the learn atom
    ax = axes[0, 1]
    v_hat = cdl.v_hat_[i_atom]
    t = np.arange(v_hat.size) / sfreq
    ax.plot(t, v_hat)
    ax.set(xlabel='Time (sec)', title='Learned temporal waveform')
    ax.grid(True)

    # Plot the psd of the time atom
    ax = axes[0, 2]
    psd = np.abs(np.fft.rfft(v_hat)) ** 2
    frequencies = np.linspace(0, sfreq / 2.0, len(psd))
    ax.semilogy(frequencies, psd)
    ax.set(xlabel='Frequencies (Hz)', title='Power Spectral Density')
    ax.grid(True)
    ax.set_xlim(0, 30)

    plt.tight_layout()
    plt.show()



.. image-sg:: /auto_examples/multicsc/images/sphx_glr_plot_somato_mu_waves_001.png
   :alt: Learned spatial pattern, Learned temporal waveform, Power Spectral Density
   :srcset: /auto_examples/multicsc/images/sphx_glr_plot_somato_mu_waves_001.png
   :class: sphx-glr-single-img






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

   **Total running time of the script:** (5 minutes 9.171 seconds)


.. _sphx_glr_download_auto_examples_multicsc_plot_somato_mu_waves.py:

.. only:: html

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

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

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

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

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

    .. container:: sphx-glr-download sphx-glr-download-zip

      :download:`Download zipped: plot_somato_mu_waves.zip <plot_somato_mu_waves.zip>`


.. only:: html

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

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