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

.. only:: html

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

        Click :ref:`here <sphx_glr_download_auto_paper_cluster_01_generate_seeds_tc.py>`
        to download the full example code

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

.. _sphx_glr_auto_paper_cluster_01_generate_seeds_tc.py:


Part 2: Generate seeds time coursed and locations
=================================================

.. GENERATED FROM PYTHON SOURCE LINES 8-20

.. code-block:: default


    import mne
    import numpy as np
    import os.path as op
    import os
    from scipy import signal
    import functions_code as myf
    import sys
    target = '.'
    noct = '6'



.. GENERATED FROM PYTHON SOURCE LINES 21-22

Loading data

.. GENERATED FROM PYTHON SOURCE LINES 22-46

.. code-block:: default


    # data path
    file_fwd = op.join(target, 'data', 'oct'+noct+'_fwd.fif')

    # load fwd
    fwd = mne.read_forward_solution(file_fwd, verbose=False)
    fwd = mne.convert_forward_solution(
        fwd, surf_ori=True, force_fixed=True, use_cps=True, verbose=False)
    fwd = mne.pick_types_forward(fwd, meg='mag', eeg=False, ref_meg=False)

    G = fwd['sol']['data']  # leadfield matrix
    G = 10**5*G  # rescale to avoid small numbers


    # dipole position
    dip_pos = fwd['source_rr']

    # dipole orientations
    dip_or = fwd['source_nn']

    # load cortico-cortical distance matrix
    cortico_dist_file = op.join(target, 'data', 'cortico_dist_oct'+noct+'.npy')
    cortico_dist = np.load(cortico_dist_file)


.. GENERATED FROM PYTHON SOURCE LINES 47-48

define features

.. GENERATED FROM PYTHON SOURCE LINES 48-71

.. code-block:: default


    N_mod = int(sys.argv[1])  # Number of simulated AR models (with connections)
    N_loc = int(sys.argv[2])  # Number of different connected pairs of locations
    # Standard deviation of the entries of the AR model
    alpha = np.random.rand(N_mod)*0.9+0.1
    alpha.sort()

    T = int(10000)   # Number of time points
    fs = int(128)    # Samples frequency
    delta_t = 1/fs
    fmin = 8
    fmax = 12

    M = G.shape[0]  # Number of sensor
    N_dense = G.shape[1]  # Number of sources in the dense source space
    N_act = int(2)  # Number of active patches
    P = int(5)   # MVAR order
    ratio_max = 1.5  # Ratio max between the intensities of the seed sources

    # store relevant features
    features = {'N_mod': N_mod, 'N_loc': N_loc, 'T': T,
                'fs': fs, 'fmin': fmin, 'fmax': fmax, 'N_act': N_act}


.. GENERATED FROM PYTHON SOURCE LINES 72-73

generate N_mod MVAR models of dimension 2 (pairs of time courses)

.. GENERATED FROM PYTHON SOURCE LINES 73-110

.. code-block:: default


    # Initialize the time courses of the seeds of the active patches
    seed_tcs = np.zeros((N_act, T, N_mod))

    nperseg = 256  # length of the window for the fourier transform
    nfft = nperseg  # number of frequencies
    i_mod = 0
    while i_mod < N_mod:
        ratio = np.inf
        while ratio > ratio_max:
            AR_mod = myf.gen_ar_model(N_act, P, alpha[i_mod])
            X_act, AR_mod = myf.gen_ar_series(AR_mod, T)
            norm_X = np.linalg.norm(X_act, axis=1)
            norm_X.sort()
            # retain pairs of time courses whose intensities are close
            # (int_max/int_min < ratio_max)
            ratio = norm_X[-1]/norm_X[0]
        norm_const = np.mean(np.std(X_act, axis=1))
        # normalize the time series so that different pairs have similar intensity
        X_act = X_act/norm_const
        # compute power spectrum in the freq range of interest
        f, Pwe = signal.welch(X_act, fs=fs, window='hann', nperseg=nperseg,
                              noverlap=nperseg // 2, nfft=nfft, detrend='constant',
                              return_onesided=True, scaling='density', axis=-1)
        P_tot = Pwe[0, :]+Pwe[1, :]
        f_in = np.intersect1d(np.where(f >= fmin)[0], np.where(f <= fmax)[0])
        # retain only time courses with sufficiently high power in the frequency
        # range of interest
        if np.sum(P_tot[f_in])/len(f_in) > 1.2*np.sum(P_tot)/len(f):
            b, a = signal.butter(3, np.array(
                [8, 12]), btype='bandpass', analog=False, output='ba', fs=fs)
            X_act = signal.filtfilt(b, a, X_act, axis=- 1, padtype='odd',
                                    padlen=None, method='pad', irlen=None)
            seed_tcs[:, :, i_mod] = X_act
            i_mod = i_mod + 1



.. GENERATED FROM PYTHON SOURCE LINES 111-112

generate patches sources

.. GENERATED FROM PYTHON SOURCE LINES 112-117

.. code-block:: default


    # select pairs of active sources (seeds) in the source space
    seed_locs = myf.select_sources(G, dip_pos, N_loc)



.. GENERATED FROM PYTHON SOURCE LINES 118-119

store the generated data in a dictionary and save it

.. GENERATED FROM PYTHON SOURCE LINES 119-128

.. code-block:: default


    seed_tc_loc = {'features': features,
                   'seed_tcs': seed_tcs, 'seed_locs': seed_locs}

    run_data_path = op.join('.', 'run_data')
    if not op.exists(run_data_path):
        os.mkdir(run_data_path)

    np.save(op.join(run_data_path, 'seed_tc_loc.npy'), seed_tc_loc)


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

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


.. _sphx_glr_download_auto_paper_cluster_01_generate_seeds_tc.py:

.. only:: html

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


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

      :download:`Download Python source code: 01_generate_seeds_tc.py <01_generate_seeds_tc.py>`

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

      :download:`Download Jupyter notebook: 01_generate_seeds_tc.ipynb <01_generate_seeds_tc.ipynb>`


.. only:: html

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

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