Examples

(tested under Python 3.6 and 3.8)

Loading data

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

# 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

# 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

# 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')