.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "auto_examples/multicsc/plot_sample_evoked_response.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_sample_evoked_response.py>`
        to download the full example code.

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

.. _sphx_glr_auto_examples_multicsc_plot_sample_evoked_response.py:


=====================================================================
Extracting artifact and evoked response atoms from the sample dataset
=====================================================================

This example illustrates how to learn rank-1 [1]_ atoms on the multivariate
sample dataset from :code:`mne`. We display a selection of atoms, featuring
heartbeat and eyeblink artifacts, two atoms of evoked responses, and a
non-sinusoidal oscillation.

.. [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 17-26

.. 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 27-28

Let us first define the parameters of our model.

.. GENERATED FROM PYTHON SOURCE LINES 28-49

.. code-block:: Python


    # Sampling frequency the signal will be resampled to
    sfreq = 150.

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

    # Regularization parameter which control sparsity
    reg = 0.1

    # Number of processors for parallel computing
    n_jobs = 5

    # To accelerate the run time of this example, we split the signal in n_slits.
    # The number of splits should actually be the smallest possible to avoid
    # introducing border artifacts in the learned atoms and it should be not much
    # larger than n_jobs.
    n_splits = 10









.. GENERATED FROM PYTHON SOURCE LINES 50-51

Next, we define the parameters for multivariate CSC

.. GENERATED FROM PYTHON SOURCE LINES 51-87

.. code-block:: Python


    from alphacsc import GreedyCDL
    cdl = GreedyCDL(
        # 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',
        # apply a temporal window reparametrization
        window=True,
        # at the end, refit the activations with fixed support and no reg to unbias
        unbiased_z_hat=True,
        # Initialize the dictionary with random chunk from the data
        D_init='chunk',
        # rescale the regularization parameter to be a percentage of lambda_max
        lmbd_max="scaled",
        reg=reg,
        # 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': 10000},
        # solver for the d-step
        solver_d='alternate_adaptive',
        solver_d_kwargs={'max_iter': 300},
        # sort atoms by explained variances
        sort_atoms=True,
        # Technical parameters
        verbose=1,
        random_state=0,
        n_jobs=n_jobs)









.. GENERATED FROM PYTHON SOURCE LINES 88-92

Load the sample data from MNE-python and select the gradiometer channels.
The MNE sample data contains MEG recordings of a subject with visual and
auditory stimuli. We load the data using utilities from MNE-python as a Raw
object and select the gradiometers from the signal.

.. GENERATED FROM PYTHON SOURCE LINES 92-107

.. code-block:: Python


    import os
    import mne
    import numpy as np

    print("Loading the data...", end='', flush=True)
    data_path = mne.datasets.sample.data_path()
    subjects_dir = os.path.join(data_path, "subjects")
    data_dir = os.path.join(data_path, 'MEG', 'sample')
    file_name = os.path.join(data_dir, 'sample_audvis_raw.fif')
    raw = mne.io.read_raw_fif(file_name, preload=True, verbose=False)
    raw.pick_types(meg='grad', eeg=False, eog=False, stim=True)
    print('done')






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

 .. code-block:: none

    Loading the data...Using default location ~/mne_data for sample...

      0%|                                              | 0.00/1.65G [00:00<?, ?B/s]
      0%|                                     | 1.48M/1.65G [00:00<02:04, 13.3MB/s]
      0%|                                     | 2.81M/1.65G [00:00<02:34, 10.7MB/s]
      0%|▏                                    | 8.01M/1.65G [00:00<01:01, 26.8MB/s]
      1%|▎                                    | 13.8M/1.65G [00:00<00:43, 38.0MB/s]
      1%|▍                                    | 20.2M/1.65G [00:00<00:34, 46.7MB/s]
      2%|▌                                    | 27.8M/1.65G [00:00<00:28, 56.4MB/s]
      2%|▊                                    | 34.0M/1.65G [00:00<00:27, 58.0MB/s]
      2%|▉                                    | 39.9M/1.65G [00:00<00:28, 56.3MB/s]
      3%|█                                    | 46.7M/1.65G [00:00<00:26, 59.8MB/s]
      3%|█▏                                   | 53.7M/1.65G [00:01<00:25, 62.8MB/s]
      4%|█▎                                   | 60.8M/1.65G [00:01<00:24, 65.3MB/s]
      4%|█▌                                   | 68.0M/1.65G [00:01<00:23, 67.3MB/s]
      5%|█▋                                   | 75.1M/1.65G [00:01<00:23, 68.5MB/s]
      5%|█▊                                   | 82.3M/1.65G [00:01<00:22, 69.5MB/s]
      5%|██                                   | 89.8M/1.65G [00:01<00:21, 71.1MB/s]
      6%|██▏                                  | 96.9M/1.65G [00:01<00:31, 49.9MB/s]
      6%|██▎                                   | 103M/1.65G [00:01<00:30, 50.7MB/s]
      7%|██▌                                   | 110M/1.65G [00:02<00:28, 54.8MB/s]
      7%|██▋                                   | 118M/1.65G [00:02<00:25, 61.2MB/s]
      8%|██▊                                   | 124M/1.65G [00:02<00:35, 43.1MB/s]
      8%|██▉                                   | 129M/1.65G [00:02<00:34, 43.7MB/s]
      8%|███▏                                  | 138M/1.65G [00:02<00:29, 52.2MB/s]
      9%|███▎                                  | 146M/1.65G [00:02<00:25, 59.5MB/s]
      9%|███▌                                  | 153M/1.65G [00:02<00:24, 62.4MB/s]
     10%|███▋                                  | 160M/1.65G [00:02<00:23, 63.9MB/s]
     10%|███▊                                  | 167M/1.65G [00:03<00:24, 61.8MB/s]
     10%|███▉                                  | 173M/1.65G [00:03<00:56, 26.3MB/s]
     11%|████                                  | 178M/1.65G [00:03<00:56, 26.1MB/s]
     11%|████▏                                 | 184M/1.65G [00:03<00:46, 31.7MB/s]
     11%|████▎                                 | 189M/1.65G [00:04<00:55, 26.6MB/s]
     12%|████▌                                 | 197M/1.65G [00:04<00:40, 35.8MB/s]
     12%|████▋                                 | 206M/1.65G [00:04<00:32, 44.9MB/s]
     13%|████▉                                 | 214M/1.65G [00:04<00:26, 53.4MB/s]
     13%|█████                                 | 222M/1.65G [00:04<00:23, 60.5MB/s]
     14%|█████▎                                | 230M/1.65G [00:04<00:21, 64.8MB/s]
     14%|█████▍                                | 238M/1.65G [00:04<00:23, 59.2MB/s]
     15%|█████▌                                | 244M/1.65G [00:05<00:41, 33.7MB/s]
     15%|█████▋                                | 250M/1.65G [00:05<00:38, 36.6MB/s]
     15%|█████▊                                | 255M/1.65G [00:05<00:43, 31.9MB/s]
     16%|█████▉                                | 260M/1.65G [00:05<00:38, 36.0MB/s]
     16%|██████                                | 266M/1.65G [00:05<00:35, 39.6MB/s]
     16%|██████▏                               | 270M/1.65G [00:05<00:33, 41.3MB/s]
     17%|██████▍                               | 278M/1.65G [00:05<00:28, 48.7MB/s]
     17%|██████▌                               | 285M/1.65G [00:06<00:25, 52.8MB/s]
     18%|██████▋                               | 291M/1.65G [00:06<00:24, 54.9MB/s]
     18%|██████▊                               | 297M/1.65G [00:06<00:37, 35.9MB/s]
     18%|██████▉                               | 301M/1.65G [00:06<00:41, 32.8MB/s]
     19%|███████                               | 309M/1.65G [00:06<00:32, 41.8MB/s]
     19%|███████▏                              | 315M/1.65G [00:06<00:29, 45.9MB/s]
     19%|███████▍                              | 321M/1.65G [00:06<00:26, 49.3MB/s]
     20%|███████▌                              | 327M/1.65G [00:07<00:25, 51.4MB/s]
     20%|███████▋                              | 333M/1.65G [00:07<00:24, 53.8MB/s]
     21%|███████▊                              | 339M/1.65G [00:07<00:24, 53.0MB/s]
     21%|███████▉                              | 345M/1.65G [00:07<00:36, 35.4MB/s]
     21%|████████                              | 349M/1.65G [00:07<00:38, 34.3MB/s]
     21%|████████                              | 353M/1.65G [00:07<00:38, 33.5MB/s]
     22%|████████▏                             | 357M/1.65G [00:07<00:37, 34.5MB/s]
     22%|████████▎                             | 363M/1.65G [00:08<00:31, 41.2MB/s]
     22%|████████▌                             | 371M/1.65G [00:08<00:25, 49.5MB/s]
     23%|████████▋                             | 379M/1.65G [00:08<00:21, 58.2MB/s]
     23%|████████▉                             | 387M/1.65G [00:08<00:19, 64.8MB/s]
     24%|█████████                             | 394M/1.65G [00:08<00:18, 66.6MB/s]
     24%|█████████▎                            | 402M/1.65G [00:08<00:17, 71.4MB/s]
     25%|█████████▍                            | 411M/1.65G [00:08<00:16, 75.4MB/s]
     25%|█████████▋                            | 419M/1.65G [00:08<00:15, 77.6MB/s]
     26%|█████████▊                            | 428M/1.65G [00:08<00:15, 79.9MB/s]
     26%|██████████                            | 436M/1.65G [00:09<00:21, 56.1MB/s]
     27%|██████████▏                           | 443M/1.65G [00:09<00:21, 57.6MB/s]
     27%|██████████▎                           | 449M/1.65G [00:09<00:19, 60.5MB/s]
     28%|██████████▌                           | 457M/1.65G [00:09<00:18, 64.8MB/s]
     28%|██████████▋                           | 464M/1.65G [00:09<00:18, 63.7MB/s]
     29%|██████████▊                           | 472M/1.65G [00:09<00:17, 66.9MB/s]
     29%|███████████                           | 479M/1.65G [00:09<00:17, 67.9MB/s]
     29%|███████████▏                          | 486M/1.65G [00:09<00:18, 63.6MB/s]
     30%|███████████▎                          | 493M/1.65G [00:09<00:17, 65.5MB/s]
     30%|███████████▍                          | 500M/1.65G [00:10<00:17, 66.9MB/s]
     31%|███████████▋                          | 508M/1.65G [00:10<00:15, 71.6MB/s]
     31%|███████████▊                          | 515M/1.65G [00:10<00:25, 44.0MB/s]
     32%|███████████▉                          | 521M/1.65G [00:10<00:33, 34.3MB/s]
     32%|████████████                          | 526M/1.65G [00:11<00:46, 24.1MB/s]
     32%|████████████▏                         | 530M/1.65G [00:11<00:42, 26.7MB/s]
     32%|████████████▎                         | 537M/1.65G [00:11<00:33, 33.7MB/s]
     33%|████████████▍                         | 543M/1.65G [00:11<00:28, 38.9MB/s]
     33%|████████████▋                         | 551M/1.65G [00:11<00:22, 48.5MB/s]
     34%|████████████▊                         | 560M/1.65G [00:11<00:19, 57.3MB/s]
     34%|█████████████                         | 567M/1.65G [00:12<00:45, 23.9MB/s]
     35%|█████████████▏                        | 572M/1.65G [00:12<00:46, 23.3MB/s]
     35%|█████████████▏                        | 576M/1.65G [00:12<00:51, 20.8MB/s]
     35%|█████████████▎                        | 580M/1.65G [00:12<00:46, 23.2MB/s]
     36%|█████████████▌                        | 588M/1.65G [00:13<00:32, 32.6MB/s]
     36%|█████████████▋                        | 595M/1.65G [00:13<00:27, 39.0MB/s]
     36%|█████████████▊                        | 600M/1.65G [00:13<00:35, 30.0MB/s]
     37%|█████████████▉                        | 606M/1.65G [00:13<00:29, 35.1MB/s]
     37%|██████████████                        | 613M/1.65G [00:13<00:25, 40.6MB/s]
     37%|██████████████▏                       | 619M/1.65G [00:13<00:22, 46.5MB/s]
     38%|██████████████▎                       | 625M/1.65G [00:13<00:22, 46.3MB/s]
     38%|██████████████▍                       | 630M/1.65G [00:14<00:30, 33.5MB/s]
     38%|██████████████▌                       | 635M/1.65G [00:14<00:27, 36.7MB/s]
     39%|██████████████▋                       | 641M/1.65G [00:14<00:25, 40.4MB/s]
     39%|██████████████▊                       | 647M/1.65G [00:14<00:22, 45.3MB/s]
     39%|██████████████▉                       | 652M/1.65G [00:14<00:21, 46.4MB/s]
     40%|███████████████                       | 657M/1.65G [00:14<00:31, 31.9MB/s]
     40%|███████████████▏                      | 661M/1.65G [00:14<00:29, 34.1MB/s]
     40%|███████████████▎                      | 665M/1.65G [00:15<00:28, 34.8MB/s]
     41%|███████████████▍                      | 672M/1.65G [00:15<00:22, 42.7MB/s]
     41%|███████████████▌                      | 679M/1.65G [00:15<00:20, 48.5MB/s]
     41%|███████████████▊                      | 686M/1.65G [00:15<00:17, 54.2MB/s]
     42%|███████████████▉                      | 692M/1.65G [00:15<00:18, 52.8MB/s]
     42%|████████████████                      | 698M/1.65G [00:15<00:17, 55.8MB/s]
     43%|████████████████▏                     | 704M/1.65G [00:15<00:16, 57.7MB/s]
     43%|████████████████▎                     | 711M/1.65G [00:15<00:15, 60.0MB/s]
     43%|████████████████▌                     | 718M/1.65G [00:15<00:14, 63.2MB/s]
     44%|████████████████▋                     | 725M/1.65G [00:15<00:14, 65.9MB/s]
     44%|████████████████▊                     | 733M/1.65G [00:16<00:13, 68.5MB/s]
     45%|█████████████████                     | 740M/1.65G [00:16<00:13, 70.1MB/s]
     45%|█████████████████▏                    | 747M/1.65G [00:16<00:35, 25.7MB/s]
     46%|█████████████████▎                    | 753M/1.65G [00:16<00:29, 30.5MB/s]
     46%|█████████████████▍                    | 760M/1.65G [00:17<00:24, 35.8MB/s]
     46%|█████████████████▋                    | 767M/1.65G [00:17<00:20, 42.7MB/s]
     47%|█████████████████▊                    | 775M/1.65G [00:17<00:17, 49.3MB/s]
     47%|█████████████████▉                    | 781M/1.65G [00:17<00:19, 45.3MB/s]
     48%|██████████████████                    | 787M/1.65G [00:17<00:20, 43.1MB/s]
     48%|██████████████████▏                   | 793M/1.65G [00:17<00:18, 46.9MB/s]
     48%|██████████████████▍                   | 800M/1.65G [00:17<00:16, 53.2MB/s]
     49%|██████████████████▌                   | 808M/1.65G [00:17<00:14, 58.6MB/s]
     49%|██████████████████▋                   | 815M/1.65G [00:18<00:13, 62.4MB/s]
     50%|██████████████████▉                   | 822M/1.65G [00:18<00:12, 65.4MB/s]
     50%|███████████████████                   | 829M/1.65G [00:18<00:31, 26.3MB/s]
     50%|███████████████████▏                  | 834M/1.65G [00:18<00:29, 27.3MB/s]
     51%|███████████████████▎                  | 841M/1.65G [00:19<00:24, 33.4MB/s]
     51%|███████████████████▍                  | 848M/1.65G [00:19<00:20, 39.5MB/s]
     52%|███████████████████▋                  | 854M/1.65G [00:19<00:18, 44.3MB/s]
     52%|███████████████████▊                  | 860M/1.65G [00:19<00:23, 33.2MB/s]
     52%|███████████████████▉                  | 865M/1.65G [00:19<00:22, 35.7MB/s]
     53%|████████████████████                  | 870M/1.65G [00:19<00:19, 40.0MB/s]
     53%|████████████████████▏                 | 877M/1.65G [00:19<00:16, 45.7MB/s]
     53%|████████████████████▎                 | 882M/1.65G [00:20<00:22, 34.4MB/s]
     54%|████████████████████▍                 | 889M/1.65G [00:20<00:18, 40.5MB/s]
     54%|████████████████████▌                 | 894M/1.65G [00:20<00:17, 42.4MB/s]
     54%|████████████████████▋                 | 899M/1.65G [00:20<00:16, 45.5MB/s]
     55%|████████████████████▊                 | 904M/1.65G [00:20<00:24, 30.7MB/s]
     55%|████████████████████▉                 | 909M/1.65G [00:20<00:21, 33.9MB/s]
     55%|█████████████████████                 | 913M/1.65G [00:21<00:30, 24.2MB/s]
     55%|█████████████████████                 | 917M/1.65G [00:21<00:29, 24.6MB/s]
     56%|█████████████████████▏                | 920M/1.65G [00:21<00:28, 25.9MB/s]
     56%|█████████████████████▎                | 925M/1.65G [00:21<00:24, 30.0MB/s]
     56%|█████████████████████▎                | 928M/1.65G [00:21<00:30, 24.0MB/s]
     57%|█████████████████████▍                | 934M/1.65G [00:21<00:22, 31.6MB/s]
     57%|█████████████████████▋                | 942M/1.65G [00:21<00:17, 41.1MB/s]
     57%|█████████████████████▊                | 950M/1.65G [00:21<00:13, 50.5MB/s]
     58%|██████████████████████                | 958M/1.65G [00:22<00:11, 59.0MB/s]
     58%|██████████████████████▏               | 966M/1.65G [00:22<00:10, 66.0MB/s]
     59%|██████████████████████▍               | 975M/1.65G [00:22<00:09, 71.4MB/s]
     59%|██████████████████████▌               | 983M/1.65G [00:22<00:18, 37.1MB/s]
     60%|██████████████████████▋               | 988M/1.65G [00:23<00:26, 24.8MB/s]
     60%|██████████████████████▊               | 995M/1.65G [00:23<00:22, 29.9MB/s]
     61%|██████████████████████▍              | 1.00G/1.65G [00:23<00:18, 36.1MB/s]
     61%|██████████████████████▌              | 1.01G/1.65G [00:23<00:14, 45.2MB/s]
     62%|██████████████████████▊              | 1.02G/1.65G [00:23<00:12, 49.6MB/s]
     62%|██████████████████████▉              | 1.02G/1.65G [00:23<00:11, 54.7MB/s]
     62%|███████████████████████              | 1.03G/1.65G [00:23<00:10, 61.4MB/s]
     63%|███████████████████████▎             | 1.04G/1.65G [00:23<00:09, 67.5MB/s]
     63%|███████████████████████▍             | 1.05G/1.65G [00:24<00:20, 29.1MB/s]
     64%|███████████████████████▌             | 1.05G/1.65G [00:25<00:35, 16.8MB/s]
     64%|███████████████████████▋             | 1.06G/1.65G [00:25<00:33, 17.6MB/s]
     64%|███████████████████████▊             | 1.06G/1.65G [00:25<00:27, 21.1MB/s]
     65%|███████████████████████▉             | 1.07G/1.65G [00:25<00:21, 26.8MB/s]
     65%|████████████████████████             | 1.07G/1.65G [00:25<00:21, 26.3MB/s]
     65%|████████████████████████▏            | 1.08G/1.65G [00:26<00:22, 25.9MB/s]
     66%|████████████████████████▎            | 1.09G/1.65G [00:26<00:16, 33.6MB/s]
     66%|████████████████████████▍            | 1.09G/1.65G [00:26<00:16, 35.1MB/s]
     66%|████████████████████████▍            | 1.09G/1.65G [00:26<00:15, 37.0MB/s]
     66%|████████████████████████▌            | 1.10G/1.65G [00:26<00:14, 38.9MB/s]
     67%|████████████████████████▋            | 1.10G/1.65G [00:26<00:24, 22.1MB/s]
     67%|████████████████████████▊            | 1.11G/1.65G [00:26<00:18, 28.9MB/s]
     68%|█████████████████████████            | 1.12G/1.65G [00:27<00:14, 37.6MB/s]
     68%|█████████████████████████▏           | 1.12G/1.65G [00:27<00:13, 38.9MB/s]
     68%|█████████████████████████▏           | 1.13G/1.65G [00:27<00:17, 30.8MB/s]
     68%|█████████████████████████▎           | 1.13G/1.65G [00:27<00:18, 28.2MB/s]
     69%|█████████████████████████▍           | 1.13G/1.65G [00:27<00:19, 26.4MB/s]
     69%|█████████████████████████▍           | 1.14G/1.65G [00:27<00:19, 27.0MB/s]
     69%|█████████████████████████▌           | 1.14G/1.65G [00:27<00:16, 31.7MB/s]
     69%|█████████████████████████▋           | 1.15G/1.65G [00:28<00:15, 31.7MB/s]
     70%|█████████████████████████▊           | 1.15G/1.65G [00:28<00:19, 25.9MB/s]
     70%|█████████████████████████▉           | 1.16G/1.65G [00:28<00:14, 33.2MB/s]
     70%|█████████████████████████▉           | 1.16G/1.65G [00:28<00:14, 34.6MB/s]
     70%|██████████████████████████           | 1.16G/1.65G [00:28<00:20, 24.3MB/s]
     71%|██████████████████████████▏          | 1.17G/1.65G [00:28<00:19, 25.1MB/s]
     71%|██████████████████████████▎          | 1.17G/1.65G [00:29<00:14, 33.5MB/s]
     71%|██████████████████████████▍          | 1.18G/1.65G [00:29<00:10, 44.6MB/s]
     72%|██████████████████████████▌          | 1.19G/1.65G [00:29<00:13, 34.5MB/s]
     72%|██████████████████████████▋          | 1.19G/1.65G [00:29<00:12, 36.9MB/s]
     72%|██████████████████████████▊          | 1.20G/1.65G [00:29<00:11, 41.1MB/s]
     73%|██████████████████████████▉          | 1.20G/1.65G [00:29<00:13, 32.6MB/s]
     73%|███████████████████████████          | 1.21G/1.65G [00:29<00:10, 40.8MB/s]
     74%|███████████████████████████▏         | 1.22G/1.65G [00:29<00:08, 48.5MB/s]
     74%|███████████████████████████▍         | 1.22G/1.65G [00:30<00:07, 57.3MB/s]
     74%|███████████████████████████▌         | 1.23G/1.65G [00:30<00:07, 56.4MB/s]
     75%|███████████████████████████▋         | 1.24G/1.65G [00:30<00:10, 38.8MB/s]
     75%|███████████████████████████▊         | 1.24G/1.65G [00:30<00:14, 27.5MB/s]
     75%|███████████████████████████▉         | 1.25G/1.65G [00:31<00:16, 25.2MB/s]
     76%|████████████████████████████         | 1.25G/1.65G [00:31<00:13, 29.2MB/s]
     76%|████████████████████████████         | 1.25G/1.65G [00:31<00:16, 23.6MB/s]
     76%|████████████████████████████▏        | 1.26G/1.65G [00:31<00:16, 24.0MB/s]
     77%|████████████████████████████▎        | 1.26G/1.65G [00:31<00:11, 32.4MB/s]
     77%|████████████████████████████▍        | 1.27G/1.65G [00:31<00:15, 25.5MB/s]
     77%|████████████████████████████▍        | 1.27G/1.65G [00:31<00:14, 26.7MB/s]
     77%|████████████████████████████▌        | 1.28G/1.65G [00:32<00:16, 23.5MB/s]
     77%|████████████████████████████▋        | 1.28G/1.65G [00:32<00:12, 29.9MB/s]
     78%|████████████████████████████▊        | 1.28G/1.65G [00:32<00:12, 28.6MB/s]
     78%|████████████████████████████▉        | 1.29G/1.65G [00:32<00:09, 37.0MB/s]
     78%|█████████████████████████████        | 1.30G/1.65G [00:32<00:08, 41.8MB/s]
     79%|█████████████████████████████▏       | 1.30G/1.65G [00:32<00:07, 48.9MB/s]
     79%|█████████████████████████████▎       | 1.31G/1.65G [00:32<00:06, 53.4MB/s]
     80%|█████████████████████████████▍       | 1.32G/1.65G [00:32<00:06, 51.7MB/s]
     80%|█████████████████████████████▋       | 1.32G/1.65G [00:33<00:05, 57.5MB/s]
     80%|█████████████████████████████▊       | 1.33G/1.65G [00:33<00:05, 56.9MB/s]
     81%|█████████████████████████████▉       | 1.34G/1.65G [00:33<00:05, 60.7MB/s]
     81%|██████████████████████████████       | 1.34G/1.65G [00:33<00:04, 65.5MB/s]
     82%|██████████████████████████████▏      | 1.35G/1.65G [00:33<00:06, 48.9MB/s]
     82%|██████████████████████████████▎      | 1.36G/1.65G [00:33<00:07, 37.4MB/s]
     82%|██████████████████████████████▍      | 1.36G/1.65G [00:33<00:07, 41.4MB/s]
     83%|██████████████████████████████▌      | 1.37G/1.65G [00:34<00:06, 43.1MB/s]
     83%|██████████████████████████████▋      | 1.37G/1.65G [00:34<00:06, 45.1MB/s]
     83%|██████████████████████████████▊      | 1.38G/1.65G [00:34<00:08, 31.8MB/s]
     84%|██████████████████████████████▉      | 1.38G/1.65G [00:34<00:08, 33.3MB/s]
     84%|███████████████████████████████      | 1.39G/1.65G [00:34<00:13, 19.8MB/s]
     84%|███████████████████████████████      | 1.39G/1.65G [00:35<00:15, 16.5MB/s]
     84%|███████████████████████████████▏     | 1.39G/1.65G [00:35<00:16, 16.1MB/s]
     85%|███████████████████████████████▎     | 1.40G/1.65G [00:35<00:10, 24.9MB/s]
     85%|███████████████████████████████▍     | 1.41G/1.65G [00:35<00:07, 33.2MB/s]
     85%|███████████████████████████████▌     | 1.41G/1.65G [00:36<00:11, 21.0MB/s]
     86%|███████████████████████████████▋     | 1.42G/1.65G [00:36<00:08, 27.4MB/s]
     86%|███████████████████████████████▊     | 1.42G/1.65G [00:36<00:10, 22.4MB/s]
     86%|███████████████████████████████▉     | 1.43G/1.65G [00:36<00:09, 24.3MB/s]
     87%|████████████████████████████████     | 1.43G/1.65G [00:36<00:07, 30.0MB/s]
     87%|████████████████████████████████▏    | 1.44G/1.65G [00:36<00:06, 35.7MB/s]
     87%|████████████████████████████████▎    | 1.44G/1.65G [00:36<00:04, 43.3MB/s]
     88%|████████████████████████████████▍    | 1.45G/1.65G [00:37<00:04, 48.3MB/s]
     88%|████████████████████████████████▋    | 1.46G/1.65G [00:37<00:03, 55.6MB/s]
     89%|████████████████████████████████▊    | 1.47G/1.65G [00:37<00:03, 60.4MB/s]
     89%|████████████████████████████████▉    | 1.47G/1.65G [00:37<00:02, 65.4MB/s]
     90%|█████████████████████████████████▏   | 1.48G/1.65G [00:37<00:03, 44.1MB/s]
     90%|█████████████████████████████████▎   | 1.49G/1.65G [00:37<00:03, 45.5MB/s]
     90%|█████████████████████████████████▍   | 1.49G/1.65G [00:37<00:03, 50.8MB/s]
     91%|█████████████████████████████████▌   | 1.50G/1.65G [00:37<00:02, 53.4MB/s]
     91%|█████████████████████████████████▋   | 1.51G/1.65G [00:38<00:02, 59.2MB/s]
     92%|█████████████████████████████████▉   | 1.51G/1.65G [00:38<00:02, 63.7MB/s]
     92%|██████████████████████████████████   | 1.52G/1.65G [00:38<00:04, 32.7MB/s]
     92%|██████████████████████████████████▏  | 1.53G/1.65G [00:39<00:05, 21.7MB/s]
     93%|██████████████████████████████████▎  | 1.53G/1.65G [00:39<00:04, 26.3MB/s]
     93%|██████████████████████████████████▍  | 1.54G/1.65G [00:39<00:03, 32.7MB/s]
     94%|██████████████████████████████████▌  | 1.55G/1.65G [00:39<00:02, 40.0MB/s]
     94%|██████████████████████████████████▋  | 1.55G/1.65G [00:39<00:03, 28.6MB/s]
     94%|██████████████████████████████████▊  | 1.56G/1.65G [00:39<00:02, 32.3MB/s]
     95%|██████████████████████████████████▉  | 1.56G/1.65G [00:39<00:02, 37.3MB/s]
     95%|███████████████████████████████████  | 1.57G/1.65G [00:40<00:02, 40.4MB/s]
     95%|███████████████████████████████████▎ | 1.57G/1.65G [00:40<00:01, 47.4MB/s]
     96%|███████████████████████████████████▍ | 1.58G/1.65G [00:40<00:01, 49.4MB/s]
     96%|███████████████████████████████████▌ | 1.59G/1.65G [00:40<00:02, 29.3MB/s]
     96%|███████████████████████████████████▋ | 1.59G/1.65G [00:40<00:01, 36.0MB/s]
     97%|███████████████████████████████████▊ | 1.60G/1.65G [00:40<00:01, 27.7MB/s]
     97%|███████████████████████████████████▊ | 1.60G/1.65G [00:41<00:02, 22.1MB/s]
     97%|███████████████████████████████████▉ | 1.61G/1.65G [00:41<00:02, 23.2MB/s]
     97%|████████████████████████████████████ | 1.61G/1.65G [00:41<00:01, 26.5MB/s]
     98%|████████████████████████████████████▏| 1.61G/1.65G [00:41<00:01, 31.2MB/s]
     98%|████████████████████████████████████▏| 1.62G/1.65G [00:41<00:01, 22.7MB/s]
     98%|████████████████████████████████████▎| 1.62G/1.65G [00:42<00:00, 28.8MB/s]
     99%|████████████████████████████████████▍| 1.63G/1.65G [00:42<00:01, 19.1MB/s]
     99%|████████████████████████████████████▌| 1.63G/1.65G [00:42<00:01, 15.3MB/s]
     99%|████████████████████████████████████▋| 1.64G/1.65G [00:42<00:00, 20.3MB/s]
     99%|████████████████████████████████████▋| 1.64G/1.65G [00:42<00:00, 23.5MB/s]
     99%|████████████████████████████████████▊| 1.64G/1.65G [00:43<00:00, 17.8MB/s]
    100%|████████████████████████████████████▉| 1.65G/1.65G [00:43<00:00, 24.0MB/s]
      0%|                                              | 0.00/1.65G [00:00<?, ?B/s]
    100%|█████████████████████████████████████| 1.65G/1.65G [00:00<00:00, 9.00TB/s]
    Download complete in 01m15s (1576.2 MB)
    NOTE: pick_types() is a legacy function. New code should use inst.pick(...).
    done




.. GENERATED FROM PYTHON SOURCE LINES 108-111

Then, we remove the powerline artifacts and high-pass filter to remove the
drift which can impact the CSC technique. The signal is also resampled to
150 Hz to reduce the computationnal burden.

.. GENERATED FROM PYTHON SOURCE LINES 111-119

.. code-block:: Python


    print("Preprocessing the data...", end='', flush=True)
    raw.notch_filter(np.arange(60, 181, 60), n_jobs=n_jobs, verbose=False)
    raw.filter(2, None, n_jobs=n_jobs, verbose=False)
    raw = raw.resample(sfreq, npad='auto', n_jobs=n_jobs, verbose=False)
    print('done')






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

 .. code-block:: none

    Preprocessing the data.../github/workspace/examples/multicsc/plot_sample_evoked_response.py:115: RuntimeWarning: Resampling of the stim channels caused event information to become unreliable. Consider finding events on the original data and passing the event matrix as a parameter.
      raw = raw.resample(sfreq, npad='auto', n_jobs=n_jobs, verbose=False)
    done




.. GENERATED FROM PYTHON SOURCE LINES 120-124

Load the data as an array and split it in chunks to allow parallel processing
during the model fit. Each split is considered as independent.
To reduce the impact of border artifacts, we use `apply_window=True`
which scales down the border of each split with a tukey window.

.. GENERATED FROM PYTHON SOURCE LINES 124-132

.. code-block:: Python


    from alphacsc.utils.signal import split_signal

    X = raw.get_data(picks=['meg'])
    info = raw.copy().pick_types(meg=True).info  # info of the loaded channels
    X_split = split_signal(X, n_splits=n_splits, apply_window=True)






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

 .. code-block:: none

    NOTE: pick_types() is a legacy function. New code should use inst.pick(...).




.. GENERATED FROM PYTHON SOURCE LINES 133-134

Fit the model and learn rank1 atoms

.. GENERATED FROM PYTHON SOURCE LINES 134-138

.. code-block:: Python


    cdl.fit(X_split)






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

 .. code-block:: none

    .................................................+
    ......
    [GreedyCDL] Converged after 56 iteration, (dz, du) = 6.867e-05, 9.388e-05
    [GreedyCDL] Fit in 1366.8s

    <alphacsc.convolutional_dictionary_learning.GreedyCDL object at 0x7f6123e55550>



.. GENERATED FROM PYTHON SOURCE LINES 139-143

Then we call the `transform` method, which returns the sparse codes
associated with X, without changing the dictionary learned during the `fit`.
Note that we transform on the *unsplit* data so that the sparse codes
reflect the original data and not the windowed data.

.. GENERATED FROM PYTHON SOURCE LINES 143-146

.. code-block:: Python

    z_hat = cdl.transform(X[None, :])






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

 .. code-block:: none

    Refitting the activation to avoid amplitude bias...done




.. GENERATED FROM PYTHON SOURCE LINES 147-152

Display a selection of atoms
----------------------------

We recognize a heartbeat artifact, an eyeblink artifact, two atoms of evoked
responses, and a non-sinusoidal oscillation.

.. GENERATED FROM PYTHON SOURCE LINES 152-201

.. code-block:: Python


    import matplotlib.pyplot as plt

    # preselected atoms of interest
    plotted_atoms = [0, 1, 2, 6, 4]

    n_plots = 3  # number of plots by atom
    n_columns = min(6, len(plotted_atoms))
    split = int(np.ceil(len(plotted_atoms) / n_columns))
    figsize = (4 * n_columns, 3 * n_plots * split)
    fig, axes = plt.subplots(n_plots * split, n_columns, figsize=figsize)
    for ii, kk in enumerate(plotted_atoms):

        # Select the axes to display the current atom
        print("\rDisplaying {}-th atom".format(kk), end='', flush=True)
        i_row, i_col = ii // n_columns, ii % n_columns
        it_axes = iter(axes[i_row * n_plots:(i_row + 1) * n_plots, i_col])

        # Select the current atom
        u_k = cdl.u_hat_[kk]
        v_k = cdl.v_hat_[kk]

        # Plot the spatial map of the atom using mne topomap
        ax = next(it_axes)
        mne.viz.plot_topomap(u_k, info, axes=ax, show=False)
        ax.set(title="Spatial pattern %d" % (kk, ))

        # Plot the temporal pattern of the atom
        ax = next(it_axes)
        t = np.arange(n_times_atom) / sfreq
        ax.plot(t, v_k)
        ax.set_xlim(0, n_times_atom / sfreq)
        ax.set(xlabel='Time (sec)', title="Temporal pattern %d" % kk)

        # Plot the power spectral density (PSD)
        ax = next(it_axes)
        psd = np.abs(np.fft.rfft(v_k, n=256)) ** 2
        frequencies = np.linspace(0, sfreq / 2.0, len(psd))
        ax.semilogy(frequencies, psd, label='PSD', color='k')
        ax.set(xlabel='Frequencies (Hz)', title="Power spectral density %d" % kk)
        ax.grid(True)
        ax.set_xlim(0, 30)
        ax.set_ylim(1e-4, 1e2)
        ax.legend()
    print("\rDisplayed {} atoms".format(len(plotted_atoms)).rjust(40))

    fig.tight_layout()





.. image-sg:: /auto_examples/multicsc/images/sphx_glr_plot_sample_evoked_response_001.png
   :alt: Spatial pattern 0, Spatial pattern 1, Spatial pattern 2, Spatial pattern 6, Spatial pattern 4, Temporal pattern 0, Temporal pattern 1, Temporal pattern 2, Temporal pattern 6, Temporal pattern 4, Power spectral density 0, Power spectral density 1, Power spectral density 2, Power spectral density 6, Power spectral density 4
   :srcset: /auto_examples/multicsc/images/sphx_glr_plot_sample_evoked_response_001.png
   :class: sphx-glr-single-img


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

 .. code-block:: none


    Displaying 0-th atom
    Displaying 1-th atom
    Displaying 2-th atom
    Displaying 6-th atom
    Displaying 4-th atom                      
    Displayed 5 atoms




.. GENERATED FROM PYTHON SOURCE LINES 202-210

Display the evoked reconstructed envelope
-----------------------------------------

The MNE sample data contains data for auditory (event_id=1 and 2) and
visual stimuli (event_id=3 and 4). We extract the events now so that we can
later identify the atoms related to different events. Note that the
convolutional sparse coding method does not need to know the events for
learning atoms.

.. GENERATED FROM PYTHON SOURCE LINES 210-217

.. code-block:: Python


    event_id = [1, 2, 3, 4]
    events = mne.find_events(raw, stim_channel='STI 014')
    events = mne.pick_events(events, include=event_id)
    events[:, 0] -= raw.first_samp






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

 .. code-block:: none

    319 events found on stim channel STI 014
    Event IDs: [ 1  2  3  4  5 32]




.. GENERATED FROM PYTHON SOURCE LINES 218-223

For each atom (columns), and for each event (rows), we compute the envelope
of the reconstructed signal, align it with respect to the event onsets, and
take the average. For some atoms, the activations are correlated with the
events, leading to a large evoked envelope. The gray area corresponds to
not statistically significant values, computing with sampling.

.. GENERATED FROM PYTHON SOURCE LINES 223-268

.. code-block:: Python


    from alphacsc.utils.signal import fast_hilbert
    from alphacsc.viz.epoch import plot_evoked_surrogates
    from alphacsc.utils.convolution import construct_X_multi

    # time window around the events. Note that for the sample datasets, the time
    # inter-event is around 0.5s
    t_lim = (-0.1, 0.5)

    n_plots = len(event_id)
    n_columns = min(6, len(plotted_atoms))
    split = int(np.ceil(len(plotted_atoms) / n_columns))
    figsize = (4 * n_columns, 3 * n_plots * split)
    fig, axes = plt.subplots(n_plots * split, n_columns, figsize=figsize)

    for ii, kk in enumerate(plotted_atoms):

        # Select the axes to display the current atom
        print("\rDisplaying {}-th atom envelope".format(kk), end='', flush=True)
        i_row, i_col = ii // n_columns, ii % n_columns
        it_axes = iter(axes[i_row * n_plots:(i_row + 1) * n_plots, i_col])

        # Select the current atom
        v_k = cdl.v_hat_[kk]
        v_k_1 = np.r_[[1], v_k][None]
        z_k = z_hat[:, kk:kk + 1]
        X_k = construct_X_multi(z_k, v_k_1, n_channels=1)[0, 0]

        # compute the 'envelope' of the reconstructed signal X_k
        correlation = np.abs(fast_hilbert(X_k))

        # loop over all events IDs
        for this_event_id in event_id:
            this_events = events[events[:, 2] == this_event_id]
            # plotting function
            ax = next(it_axes)
            this_info = info.copy()
            event_info = dict(event_id = this_event_id, events=events)
            this_info['temp'] = event_info
            plot_evoked_surrogates(correlation, info=this_info, t_lim=t_lim, ax=ax,
                                   n_jobs=n_jobs, label='event %d' % this_event_id)
            ax.set(xlabel='Time (sec)', title="Evoked envelope %d" % kk)
    print("\rDisplayed {} atoms".format(len(plotted_atoms)).rjust(40))
    fig.tight_layout()




.. image-sg:: /auto_examples/multicsc/images/sphx_glr_plot_sample_evoked_response_002.png
   :alt: Evoked envelope 0, Evoked envelope 1, Evoked envelope 2, Evoked envelope 6, Evoked envelope 4, Evoked envelope 0, Evoked envelope 1, Evoked envelope 2, Evoked envelope 6, Evoked envelope 4, Evoked envelope 0, Evoked envelope 1, Evoked envelope 2, Evoked envelope 6, Evoked envelope 4, Evoked envelope 0, Evoked envelope 1, Evoked envelope 2, Evoked envelope 6, Evoked envelope 4
   :srcset: /auto_examples/multicsc/images/sphx_glr_plot_sample_evoked_response_002.png
   :class: sphx-glr-single-img


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

 .. code-block:: none


    Displaying 0-th atom envelopeUsing data from preloaded Raw for 72 events and 91 original time points ...
    0 bad epochs dropped
    Using data from preloaded Raw for 73 events and 91 original time points ...
    0 bad epochs dropped
    Using data from preloaded Raw for 73 events and 91 original time points ...
    0 bad epochs dropped
    Using data from preloaded Raw for 70 events and 91 original time points ...
    0 bad epochs dropped

    Displaying 1-th atom envelopeUsing data from preloaded Raw for 72 events and 91 original time points ...
    0 bad epochs dropped
    Using data from preloaded Raw for 73 events and 91 original time points ...
    0 bad epochs dropped
    Using data from preloaded Raw for 73 events and 91 original time points ...
    0 bad epochs dropped
    Using data from preloaded Raw for 70 events and 91 original time points ...
    0 bad epochs dropped

    Displaying 2-th atom envelopeUsing data from preloaded Raw for 72 events and 91 original time points ...
    0 bad epochs dropped
    Using data from preloaded Raw for 73 events and 91 original time points ...
    0 bad epochs dropped
    Using data from preloaded Raw for 73 events and 91 original time points ...
    0 bad epochs dropped
    Using data from preloaded Raw for 70 events and 91 original time points ...
    0 bad epochs dropped

    Displaying 6-th atom envelopeUsing data from preloaded Raw for 72 events and 91 original time points ...
    0 bad epochs dropped
    Using data from preloaded Raw for 73 events and 91 original time points ...
    0 bad epochs dropped
    Using data from preloaded Raw for 73 events and 91 original time points ...
    0 bad epochs dropped
    Using data from preloaded Raw for 70 events and 91 original time points ...
    0 bad epochs dropped

    Displaying 4-th atom envelopeUsing data from preloaded Raw for 72 events and 91 original time points ...
    0 bad epochs dropped
    Using data from preloaded Raw for 73 events and 91 original time points ...
    0 bad epochs dropped
    Using data from preloaded Raw for 73 events and 91 original time points ...
    0 bad epochs dropped
    Using data from preloaded Raw for 70 events and 91 original time points ...
    0 bad epochs dropped
                      
    Displayed 5 atoms




.. GENERATED FROM PYTHON SOURCE LINES 269-286

Display the equivalent dipole for a learned topomap
---------------------------------------------------

Finally, let us fit a dipole to one of the atoms. To fit a dipole,
we need the following:

* BEM solution: Obtained by running the cortical reconstruction pipeline
  of Freesurfer and describes the conductivity of different tissues in
  the head.
* Trans: An affine transformation matrix needed to bring the data
  from sensor space to head space. This is usually done by coregistering
  the fiducials with the MRI.
* Noise covariance matrix: To whiten the data so that the assumption
  of Gaussian noise model with identity covariance matrix is satisfied.

We recommend users to consult the MNE documentation for further information.


.. GENERATED FROM PYTHON SOURCE LINES 286-294

.. code-block:: Python

    subjects_dir = os.path.join(data_path, 'subjects')
    fname_bem = os.path.join(subjects_dir, 'sample', 'bem',
                             'sample-5120-bem-sol.fif')
    fname_trans = os.path.join(data_path, 'MEG', 'sample',
                               'sample_audvis_raw-trans.fif')
    fname_cov = os.path.join(data_path, 'MEG', 'sample', 'sample_audvis-cov.fif')









.. GENERATED FROM PYTHON SOURCE LINES 295-298

Let us construct an evoked object for MNE with the spatial pattern of the
atoms.


.. GENERATED FROM PYTHON SOURCE LINES 298-301

.. code-block:: Python

    evoked = mne.EvokedArray(cdl.u_hat_.T, info)









.. GENERATED FROM PYTHON SOURCE LINES 302-304

Fit a dipole to each of the atoms.


.. GENERATED FROM PYTHON SOURCE LINES 304-308

.. code-block:: Python

    dip = mne.fit_dipole(evoked, fname_cov, fname_bem, fname_trans,
                         n_jobs=n_jobs, verbose=False)[0]









.. GENERATED FROM PYTHON SOURCE LINES 309-312

Plot the dipole fit from the 3rd atom, linked to mu-wave and display the
goodness of fit.


.. GENERATED FROM PYTHON SOURCE LINES 312-340

.. code-block:: Python


    atom_dipole_idx = 4

    from mpl_toolkits.mplot3d import Axes3D

    fig = plt.figure(figsize=(10, 4))

    # Display the dipole fit
    ax = fig.add_subplot(1, 3, 1, projection='3d')
    dip.plot_locations(fname_trans, 'sample', subjects_dir, idx=atom_dipole_idx,
                       ax=ax)
    ax.set_title('Atom #{} (GOF {:.2f}%)'.format(atom_dipole_idx,
                                                 dip.gof[atom_dipole_idx]))

    # Plot the spatial map
    ax = fig.add_subplot(1, 3, 2)
    mne.viz.plot_topomap(cdl.u_hat_[atom_dipole_idx], info, axes=ax)

    # Plot the temporal atom
    ax = fig.add_subplot(1, 3, 3)
    t = np.arange(n_times_atom) / sfreq
    ax.plot(t, cdl.v_hat_[atom_dipole_idx])
    ax.set_xlim(0, n_times_atom / sfreq)
    ax.set(xlabel='Time (sec)', title="Temporal pattern {}"
           .format(atom_dipole_idx))

    fig.suptitle('')
    fig.tight_layout()



.. image-sg:: /auto_examples/multicsc/images/sphx_glr_plot_sample_evoked_response_003.png
   :alt: , Atom #4 (GOF 47.27%), Temporal pattern 4
   :srcset: /auto_examples/multicsc/images/sphx_glr_plot_sample_evoked_response_003.png
   :class: sphx-glr-single-img






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

   **Total running time of the script:** (28 minutes 6.999 seconds)


.. _sphx_glr_download_auto_examples_multicsc_plot_sample_evoked_response.py:

.. only:: html

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

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

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

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

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

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

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


.. only:: html

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

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