.. _tutorial: Examples ================== *(tested under Python 3.6 and 3.8)* Loading data ######## .. code-block:: python import connectivipy as cp # remember that data should be in a shape (k, N, R), # where k - number of channels, N - data points, R - number of trials # for numpy.array simply put that array as a first argument # when initializing Data class # fs means sampling frequency # chan_names is a list with channel names (length of list must be # the same as first dimension of data) # data_info - additional infromation about the data dt = cp.Data(numpy_array_data, fs=32., chan_names=['Fp1','O1'], data_info='sml') # Matlab data we can read giving a path to a matlab file # and in data_info we put Matlab variable name as a string dd = cp.Data('adata.mat', data_info='bdata') # similarly for SignalML data, but in data_info you need to point out # that you want to read 'sml' data from *.raw and *.xml files with the # same name dt = cp.Data('cdata.raw', data_info='sml') Data class example ######## .. code-block:: python # Example 1 import numpy as np import connectivipy as cp from connectivipy import mvar_gen ### MVAR model coefficients """ MVAR parameters taken from: Sameshima K. & Baccala L. A., Partial directed coherence : a new concept in neural structure determination. Biol. Cybern. (2001) You can compare results with Fig. 3. from that article. """ # let's build mvar model matrix A = np.zeros((2, 5, 5)) # 2 - first dimension is model order # 5 - second and third dimensions mean number of channels A[0, 0, 0] = 0.95 * 2**0.5 A[1, 0, 0] = -0.9025 A[0, 1, 0] = -0.5 A[1, 2, 1] = 0.4 A[0, 3, 2] = -0.5 A[0, 3, 3] = 0.25 * 2**0.5 A[0, 3, 4] = 0.25 * 2**0.5 A[0, 4, 3] = -0.25 * 2**0.5 A[0, 4, 4] = 0.25 * 2**0.5 # multitrial signal generation from a matrix above # let's generate 5-channel signal with 1000 data points # and 5 trials using function mvar_gen ysig = np.zeros((5, 10**3, 5)) ysig[:, :, 0] = mvar_gen(A, 10**3) ysig[:, :, 1] = mvar_gen(A, 10**3) ysig[:, :, 2] = mvar_gen(A, 10**3) ysig[:, :, 3] = mvar_gen(A, 10**3) ysig[:, :, 4] = mvar_gen(A, 10**3) #### connectivity analysis data = cp.Data(ysig, 128, ["Fp1", "Fp2", "Cz", "O1", "O2"]) # you may want to plot data (in multitrial case only one trial is shown) data.plot_data() # fit mvar using Yule-Walker algorithm and order 2 data.fit_mvar(2, 'yw') # you can capture fitted parameters and residual matrix ar, vr = data.mvar_coefficients # now we investigate connectivity using gDTF gdtf_values = data.conn('gdtf') gdtf_significance = data.significance(Nrep=200, alpha=0.05) data.plot_conn('gDTF') # short time version with default parameters pdc_shorttime = data.short_time_conn('pdc', nfft=100, no=10) data.plot_short_time_conn("PDC") How to use specific classes ######## .. code-block:: python # Example 2 import numpy as np import matplotlib.pyplot as plt import connectivipy as cp from connectivipy import mvar_gen """ In this example we don't use Data class """ fs = 256. acf = np.zeros((3, 3, 3)) # matrix shape meaning # (p,k,k) k - number of channels, # p - order of mvar parameters acf[0, 0, 0] = 0.3 acf[0, 1, 0] = 0.6 acf[1, 0, 0] = 0.1 acf[1, 1, 1] = 0.2 acf[1, 2, 0] = 0.6 acf[2, 2, 2] = 0.2 acf[2, 1, 0] = 0.4 # generate 3-channel signal from matrix above y = mvar_gen(acf, int(10e4)) # assign static class cp.Mvar to variable mv mv = cp.Mvar # find best model order using Vieira-Morf algorithm best, crit = mv.order_akaike(y, 15, 'vm') plt.plot(1+np.arange(len(crit)), crit, 'g') plt.show() print(best) # here we know that this is 3 but in real-life cases # we are always uncertain about it # now let's fit parameters to the signal av, vf = mv.fit(y, best, 'vm') # and check whether values are correct +/- 0.01 print(np.allclose(acf, av, 0.01, 0.01)) # now we can calculate Directed Transfer Function from the data dtf = cp.conn.DTF() dtfval = dtf.calculate(av, vf, 128) # all possible methods are visible in that dictionary: print(cp.conn.conn_estim_dc.keys()) cp.plot_conn(dtfval, 'DTF values', fs) Instantaneous ######## .. code-block:: python # Example 3 import numpy as np import matplotlib.pyplot as plt import connectivipy as cp """ This example reproduce simulation from article: Erla S et all (2009) "Multivariate autoregressive model with instantaneous effects to improve brain connectivity estimation" """ # let's make a matrix from original article bcf = np.zeros((4, 5, 5)) # matrix shape meaning (k, k, p) k - number of channels, # p - order of mvar parameters bcf[1, 0, 0] = 1.58 bcf[2, 0, 0] = -0.81 bcf[0, 1, 0] = 0.9 bcf[2, 1, 1] = -0.01 bcf[3, 1, 4] = -0.6 bcf[1, 2, 1] = 0.3 bcf[1, 2, 2] = 0.8 bcf[2, 2, 1] = 0.3 bcf[2, 2, 2] = -0.25 bcf[3, 2, 1] = 0.3 bcf[0, 3, 1] = 0.9 bcf[1, 3, 1] = -0.6 bcf[3, 3, 1] = 0.3 bcf[1, 4, 3] = -0.3 bcf[2, 4, 0] = 0.9 bcf[2, 4, 3] = -0.3 bcf[3, 4, 2] = 0.6 # now we build a corresponding MVAR process without instantenous effect L = np.linalg.inv(np.eye(5)-bcf[0]) acf = np.zeros((3, 5, 5)) for i in range(3): acf[i] = np.dot(L, bcf[i+1]) # generate 5-channel signals from matrix above signal_inst = cp.mvar_gen_inst(bcf, int(10e4)) signal = cp.mvar_gen(acf, int(10e4)) # fit MVAR parameters bv, vfb = cp.Mvar.fit(signal_inst, 3, 'yw') av, vfa = cp.Mvar.fit(signal, 3, 'yw') # use connectivity estimators ipdc = cp.conn.iPDC() ipdcval = ipdc.calculate(bv, vfb, 1.) pdc = cp.conn.PDC() pdcval = pdc.calculate(av, vfa, 1.) def plot_double_conn(values_a, values_b, name='', fs=1, ylim=None, xlim=None, show=True): "function to plot two sets of connectivity values" fq, k, k = values_a.shape fig, axes = plt.subplots(k, k) freqs = np.linspace(0, fs*0.5, fq) if not xlim: xlim = [0, np.max(freqs)] if not ylim: ylim = [0, 1] for i in range(k): for j in range(k): axes[i, j].fill_between(freqs, values_b[:, i, j], 0, facecolor='red', alpha=0.5) axes[i, j].fill_between(freqs, values_a[:, i, j], 0, facecolor='black', alpha=0.5) axes[i, j].set_xlim(xlim) axes[i, j].set_ylim(ylim) plt.suptitle(name,y=0.98) plt.tight_layout() plt.subplots_adjust(top=0.92) if show: plt.show() plot_double_conn(pdcval**2, ipdcval**2, 'PDC / iPDC')