sails documentation¶

Contents:

Tutorials¶

TODO: Write this

Tutorial 1 - A pink noise system¶

In this tutorial we will demonstrate how to set up a simple univariate AR model which models a pink noise process. We will use the model to demonstrate how to extract the transfer function using both a Fourier and Modal estimator.

We start by importing the routines we will need from the numpy, sails and matplotlib. We set the ggplot style for all plots in the tutorials.

import numpy as np
from sails import generate_pink_roots, root_plot, model_from_roots
import matplotlib.pyplot as plt

plt.style.use('ggplot')

We now define the meta-parameters of our system. We will arbitrarily pretend that we are using a sampling rate of 100Hz. From this, we can compute our nyquist frequency, and a vector of frequencies at which we will evaluate measures such as our Fourier estimator. The variable freq_vect will come up often during our tutorials. In this case, we will estimate 64 frequencies which are linearly spaced between 0Hz and Nyquist. You can alter freq_vect and examine how this affects the results of the later estimations.

sample_rate = 100
nyq = sample_rate / 2.

freq_vect = np.linspace(0, nyq, 64)

As mentioned above, for our first system we will emulate a system which behaves as a pink noise process. Whereas we would normally learn our model from some data, for this example we will construct a model from the polynomial form of the AR model. A function (generate_pink_roots()) has been provided which calculates a set of roots which will generate the appropriate roots:

roots = generate_pink_roots(1)

Given the roots of the model, we can plot these in a Z-plane representation to examine them. A function (root_plot()) is provided to make this straightforward.

ax = root_plot(roots)
ax.figure.show()

The plot shows frequency increasing counterclockwise from the x-axis. Nyquist frequency is on the negative x-axis. For real systems, the structure of the pole plot will be mirrored across the x-axis.

As mentioned above, we would normally set up a model by learning the parameters from data. During this exploratory tutorial, we can analytically create the model from the roots. The model_from_roots() function will perform this task for you.

m = model_from_roots(roots)

No matter how we create our model, it will be a subclass of AbstractLinearModel. All subclasses of AbstractLinearModel provide a number of guarantees so that they can all be used with the various estimator routines we discuss below. We will discuss fitting a model from data in a later tutorial.

Now that we have our model, we can extract the transfer function (H) using two different methods. The first method extracts the transfer function using a Fourier transform of the parameters.

from sails import FourierMvarMetrics

F = FourierMvarMetrics.initialise(m, sample_rate, freq_vect)
F_H = F.H

The second method using a modal decomposition of the parameter matrix to break the parameters apart into, hopefully, interpretable components. We first of all create a ModalMvarMetrics and then use the MvarModalDecomposition object inside of this to extract the transfer function for each mode. Each mode will consist of either a pair of poles in the complex plane or a single real mode. The per_mode_transfer_function() will take this into account for you.

from sails import ModalMvarMetrics

M = ModalMvarMetrics.initialise(m, sample_rate, freq_vect)
M_H = M.modes.per_mode_transfer_function(sample_rate, freq_vect)

Finally, we can plot our modes in both forms to examine the relationship between them.

# Plot our fourier and modal spectra

f2 = plt.figure()

plt.plot(freq_vect, np.abs(F_H).squeeze(), 'o');
plt.plot(freq_vect, np.abs(M_H).squeeze());

plt.xlabel('Frequency (Hz)')
plt.ylabel('Frequency Response')

f2.show()

Tutorial 2 - A more interesting oscillation¶

In this tutorial, we will extend our example to include oscillations on top of the background pink noise. We begin in the same way as with the previous tutorial.

import numpy as np
from sails import generate_pink_roots, root_plot, model_from_roots
import matplotlib.pyplot as plt

plt.style.use('ggplot')

sample_rate = 100
nyq = sample_rate / 2.

freq_vect = np.linspace(0, nyq, 64)

roots = generate_pink_roots(1)

We now add in energy at two frequencies and plot the poles.

roots[3:5] = roots[3:5]*1.1

roots[7:9] = roots[7:9]*1.1

ax = root_plot(roots)

ax.figure.show()

In this plot we can see that in addition to our pink noise mode on the x-axis (compare back to the equivalent plot in the first tutorial), we have two poles which now sit closer to the unit circle - these correspond to the additional energy which we have placed in the signal.

We can now continue to set up our model from our roots, generate the Fourier and Modal transfer functions and plot them as before:

from sails import FourierMvarMetrics, ModalMvarMetrics

m = model_from_roots(roots)

F = FourierMvarMetrics.initialise(m, sample_rate, freq_vect)
F_H = F.H

M = ModalMvarMetrics.initialise(m, sample_rate, freq_vect)
M_H = M.modes.transfer_function(sample_rate, freq_vect, sum_modes=True)

# Plot our fourier and modal spectra
f2 = plt.figure()

plt.plot(freq_vect, np.abs(F_H).squeeze(), 'o');
plt.plot(freq_vect, np.abs(M_H).squeeze());

plt.xlabel('Frequency (Hz)')
plt.ylabel('Frequency Response')

f2.show()

In comparison to our plot in the first tutorial, there are (as expected) two peaks in the Fourier transfer function at around 8Hz and 17Hz. We can also see that each of these two peaks has been extracted into a separate mode in the modal decomposition. This can be seen slightly more clearly by using a log scale for the y axis:

f3 = plt.figure()

plt.semilogy(freq_vect, np.abs(F_H).squeeze(), 'o');

plt.semilogy(freq_vect, np.abs(M_H).squeeze());

plt.xlabel('Frequency (Hz)')

plt.ylabel('Frequency Response')

f3.show()

We can also go on to extract the magnitude of the eigenvalues for each of the poles:

ev = np.abs(M.modes.evals)

Note that we have an eigenvalue for each of the poles, not each of the modes. As the eigenvalue will be identical for both of the poles in a given mode (where the mode consists of a pole pair), we can reduce down to examining one pole for each mode. To get the relevant indices, we can extract the information from the mode_indices property on the MvarModalDecomposition object. We can then extract other interesting pieces of information such as the period of the oscillation in the same way:

idx = [i[0] for i in M.modes.mode_indices]

print(ev[idx])

[[ 0.96393242]
 [ 0.83308007]
 [ 0.88077379]
 [ 0.78215297]
 [ 0.84650308]
 [ 0.7603617 ]
 [ 0.75345934]
 [ 0.74824023]
 [ 0.74435422]
 [ 0.73893252]
 [ 0.73980532]
 [ 0.74158676]]

print(M.modes.peak_frequency[idx])

[[  0.        ]
 [  4.37353106]
 [  8.83055516]
 [ 13.2106547 ]
 [ 17.56247546]
 [ 21.90065104]
 [ 26.23121568]
 [ 30.55715939]
 [ 34.88016321]
 [ 47.84046803]
 [ 43.52119577]
 [ 39.20127164]]

From this, we can see that the modes which fit our two oscillations are found at 8.8Hz and 17.6Hz.

Tutorial 3 - Fitting real univariate data¶

In the previous two tutorials we set up our system using the polynomial representation. In most cases, we will want to learn the parameters of our model from real data. In this section, we cover how this is done.

We will be using some data from a single MEG Virtual Electrode reconstruction. The data has been downsampled from the original sampling rate for the purposes of demonstration and ease of analysis.

We start with some imports and finding the location of our example data:

from os.path import join
import numpy as np
import matplotlib.pyplot as plt

from sails import find_example_path, root_plot, model_from_roots

plt.style.use('ggplot')

ex_dir = find_example_path()

The original MEG data was sampled at 2034.51Hz and has been downsampled by 24 times giving us a sampling rate of just under 85Hz. The data has also been Z-scored. The data is stored in an HDF5 file in the example data directory as meg_single.hdf5; the data is stored as the dataset X.

import h5py

sample_rate = 2034.51 / 24

nyq = sample_rate / 2.

freq_vect = np.linspace(0, nyq, 64)

X = h5py.File(join(ex_dir, 'meg_single.hdf5'), 'r')['X'][...]

print(X.shape)

(1, 30157, 1)

The form of the data for the fitting routines is (nsignals, nsamples, ntrials). In the current case, we are performing a univariate analysis so we only have one signal. We have just over 30,000 data points and one one trial.

We are now in a position to set up our model. We will use the Vieira-Morf algorithm to fit our model, and fit a model of order 19 using the delay_vect argument. We will discuss model order selection in a later tutorial.

from sails import VieiraMorfLinearModel

delay_vect = np.arange(20)

m = VieiraMorfLinearModel.fit_model(X, delay_vect)

print(m.order)

We can now compute some diagnostics in order to check that our model looks sensible. We can start by looking at the R-squared value:

diag = m.compute_diagnostics(X)

print(diag.R_square)

0.41148872414

We can see that we are explaining just over 40% of the variance with our model which, given that we are modelling human MEG data collected over roughly 6 minutes, is reasonable.

The diagnostics class also gives us access to various other measures:

print(diag.AIC)

-25563.4038876

print(diag.BIC)

-25396.883104

print(diag.LL)

-25603.4038876

print(diag.DW)

[ 2.00103186]

print(diag.SI)

0.976853801393

print(diag.SR)

0.0663411856574

In turn, these are the Akaike Information Criterion (AIC), Bayesian Information Criterion (BIC), Log Likelihood (LL), Durbin-Watson coefficient (DW), Stability Index (SI) and Stability Ratio (SR). It is also possible to access the Percent Consistency (PC), although this is not computed by default due to it being memory intensive - you can compute this using the percent_consistency() routine.

As in our previous examples, we can extract our metrics and plot the transfer functions using both the Fourier and Modal methods as we have previously done:

from sails import FourierMvarMetrics, ModalMvarMetrics

F = FourierMvarMetrics.initialise(m, sample_rate, freq_vect)
F_H = F.H

M = ModalMvarMetrics.initialise(m, sample_rate, freq_vect)
M_H = M.modes.per_mode_transfer_function(sample_rate, freq_vect)

# Plot our fourier and modal spectra

f2 = plt.figure()

plt.plot(freq_vect, np.abs(F_H).squeeze(), 'o');
plt.plot(freq_vect, np.abs(M_H).squeeze());

plt.xlabel('Frequency (Hz)')
plt.ylabel('Frequency Response')

f2.show()

f3 = plt.figure()

plt.semilogy(freq_vect, np.abs(F_H).squeeze(), 'o');
plt.semilogy(freq_vect, np.abs(M_H).squeeze());

plt.xlabel('Frequency (Hz)')
plt.ylabel('Frequency Response')

f3.show()

In our previous examples, the model was defined by the structure of the polynomial, and we could analytically write down the form of the poles. In this example, we have learned the parameters from data. We may want to look at the plot of the system roots:

ax = M.modes.pole_plot()
ax.figure.show()

As previously, we can also go on to extract the magnitude of the eigenvalues and period of each of the modes:

ev = np.abs(M.modes.evals)

idx = [i[0] for i in M.modes.mode_indices]

print(ev[idx])

[[ 0.9768538 ]
 [ 0.88790649]
 [ 0.83392372]
 [ 0.82257801]
 [ 0.7869429 ]
 [ 0.79447331]
 [ 0.80151847]
 [ 0.81061091]
 [ 0.81709057]
 [ 0.80730369]]

print(M.modes.peak_frequency[idx])

[[  0.        ]
 [  8.88121272]
 [  4.51372792]
 [ 13.89286835]
 [ 18.15521077]
 [ 22.01471814]
 [ 26.47871083]
 [ 30.9324283 ]
 [ 35.59229698]
 [ 40.1740829 ]]

From this, we can see that the mode which primarily fits this data is an alpha mode at 8.9Hz.

Tutorial 4 - Exploring model order¶

In this tutorial, we will examine how to look at different model orders in order to determine the most appropriate model.

The first section of this tutorial should be familiar to you from previous tutorials. We start by loading in some data:

from os.path import join

import h5py

import numpy as np
import matplotlib.pyplot as plt

from sails import find_example_path, VieiraMorfLinearModel

plt.style.use('ggplot')

ex_dir = find_example_path()

sample_rate = 2034.51 / 24
nyq = sample_rate / 2.

freq_vect = np.linspace(0, nyq, 64)

X = h5py.File(join(ex_dir, 'meg_single.hdf5'), 'r')['X'][...]

We will use the AIC values to examine how appropriate our model is. We therefore set up an empty list to store our AIC values in and fit a model for a range of different orders. We can then plot up the values of AIC.

AICs = []

orders = []

for delays in range(2, 35):
    delay_vect = np.arange(delays)
    m = VieiraMorfLinearModel.fit_model(X, delay_vect)
    diag = m.compute_diagnostics(X)
    orders.append(m.order)
    AICs.append(diag.AIC)

f = plt.figure()

plt.plot(orders, AICs, 'o');

plt.xlabel('Model Order')
plt.ylabel('AIC')

f.show()

The same principle can be used for any of the other measures such as BIC and R-squared. You should always be aware that the change the number of parameters in the model for different orders means that certain measures (such as R-squared) will always increase for increasing model order.

R_squares = []

orders = []

for delays in range(2, 35):
    delay_vect = np.arange(delays)
    m = VieiraMorfLinearModel.fit_model(X, delay_vect)
    diag = m.compute_diagnostics(X)
    orders.append(m.order)
    R_squares.append(diag.R_square)

f2 = plt.figure()

plt.plot(orders, R_squares, 'o');

plt.xlabel('Model Order')
plt.ylabel('$R^2$')

f2.show()

Tutorial 5 - MVAR Connectivity Estimation¶

In this tutorial, we will explore a range of connectivity estimators in a simulated network.

We start by importing sails and defining some meta-parameters as we did in previous tutorials.

import numpy as np

import sails

sample_rate = 128

We will use an example network from Baccala & Sameshima 2001. This is defined within the simulation module in Sails. We can import a signal generator based on figure 2 from Baccala & Sameshima 2001.

siggen = sails.simulate.Baccala2001_fig2()

The siggen object contains the autoregressive parameters defining the network and can create a model containg these paratmeters which we can use in further analyses. Here we create a model containing the ‘true’ autoregressive parameters and use it to compute the connectivity metrics.

m = siggen.generate_model()

freq_vect = np.linspace(0,sample_rate/2,36)

F = sails.FourierMvarMetrics.initialise(m, sample_rate, freq_vect)

F is an object containing methods to compute a range of frequency domain connectivity metrics. Each metric is evaluated across the range of frequencies_defined in freq_vect and can be plotted using the plot_vector function in sails. Next we plot the Magnitude Squared Coherence estimated from our simulation paramter.

import matplotlib.pyplot as plt

sails.plotting.plot_vector( np.abs(F.S), freq_vect, diag=True )

plt.show()

This will generate a matrix of plots, each plot represents the coherence as a function of frequency between the node specified in the column and row labels. In this case, we find that all our nodes are strongly coherence with each other. This is as the coherence does not distinguish between direct and indirect connections. For example, nodes 1 and 5 are only connected through node 4 yet the coherence still shows a connection. The coherence is also symmetric, that is the connection from node 1->2 is the same as node 2->1.

Next we plot the Directed Transfer Function which is a directed measure that is able to show when connections are not symmetrical

sails.plotting.plot_vector( F.directed_transfer_function, freq_vect, diag=True )

plt.show()

The Directed Transfer Function shows far fewer connections than the Magnitude Squared Coherence. We can now see that the connections between node 1 and the rest of the nodes are asymmetrical. This means that node 1 is driving the others. The interaction between nodes 4 and 5 is now also isolated. A remaining issue is that the Directed Transfer Function is still sensitive to indirect connections, as we can see by the power in the subplot between node 1 and 5. The Partial Directed Coherence aims to address this problem.

sails.plotting.plot_vector( F.partial_directed_coherence, freq_vect, diag=True )

plt.show()

The Partial Directed Coherence now shows only the direct connections within our network. We retain our frequency resolution and the sensitivity to asymmetrical connections. There are many other MVAR derived connectivity metrics available within sails with different properties and sensitivities, these include:

Coherency (sails.mvar_metrics.coherency(), sails.mvar_metrics.AbstractMVARMetrics.coherency)
Imaginary Coherence (sails.mvar_metrics.AbstractMVARMetrics.imaginary_coherence)
Phase Coherence (sails.mvar_metrics.AbstractMVARMetrics.phase_coherence)
Magnitude Squared Coherence (sails.mvar_metrics.AbstractMVARMetrics.magnitude_squared_coherence)
Partial Coherence (sails.mvar_metrics.partial_coherence(), sails.mvar_metrics.AbstractMVARMetrics.partial_coherence)
Directed Transfer Function (sails.mvar_metrics.directed_transfer_function(), sails.mvar_metrics.AbstractMVARMetrics.directed_transfer_function)
Full Frequency Directed Transfer Function (sails.mvar_metrics.AbstractMVARMetrics.ff_directed_transfer_function)
Directed Directed Transfer Function (sails.mvar_metrics.AbstractMVARMetrics.d_directed_transfer_function)
Partial Directed Coherence (sails.mvar_metrics.partial_directed_coherence(), sails.mvar_metrics.AbstractMVARMetrics.partial_directed_coherence)
Isolated Effective Coherence (sails.mvar_metrics.isolated_effective_coherence(), sails.mvar_metrics.AbstractMVARMetrics.isolated_effective_coherence)

In the second part of this tutorial we will look at fitting and MVAR model and the Partial Directed Coherence to simulated data, rather than from the ‘true’ model.

We can generate data from our simulated model using the sails.simulate:AbstractSigGen.generate_signal() method and specifying the sample_rate and number of samples to generate

X = siggen.generate_signal(sample_rate=128,num_samples=640)

X is a (nchannels x nsamples) array containing our simulated data. We can plot X using matplotlib

plt.figure()

for ii in range(5):
    plt.plot(X[ii,:]+(ii*10))

plt.show()

We now have a figure containing 5 time-series from our simulation. We can see there is an oscillation by eye and that some of the time-series vary together more than others.

We can fit a model to the simulated data and compute connectivity metrics as we did in previous tutorials.

delay_vect = np.arange(4)

m = sails.VieiraMorfLinearModel.fit_model(X, delay_vect)

F = sails.FourierMvarMetrics.initialise(m, sample_rate, freq_vect)

diag = m.compute_diagnostics( X )

We check that our model is fitting well by interrogating the diagnostics. Here we see that we are explaining around 56% of the total variance in the signal and that our model is stable (diag.SI = .91).

Let’s compare the Partial Directed Coherence from our fitted model to the Partial Directed Coherence from the ‘true’ model.

m0 = siggen.generate_model() # This is our true model

F0 = sails.FourierMvarMetrics.initialise( m0, sample_rate, freq_vect )

pdc = np.concatenate( (F0.partial_directed_coherence,F.partial_directed_coherence), axis=3 )

sails.plotting.plot_vector( pdc, freq_vect, diag=True,line_labels=['True','Fitted'])

plt.show()

The resulting figure shows the nodes by nodes matrix of subplots containing the PDC estimates. We can see that our model is doing a pretty good job approximating the true pattern of connectivity. There may be some false-positive connections which show power for the fitted model but not for the true model.

Try re-running the simulation with a higher or lower number of samples in the time series. You should see that the estimation starts to really break down (lots of false positives and a distorted spectrum shape) when we have too few samples (e.g. num_samples = 128) and becomes nearly perfect when we have a very long time-series (e.g. num_sample = 2048)

Tutorial 6 - Sliding window univariate model estimation¶

In this tutorial, we will examine fitting multiple models using a sliding window approach.

For this tutorial, we will use the same MEG example data which we have used in previous tutorials.

We start by importing our modules and finding and loading the example data.

from os.path import join

import h5py

import numpy as np
import matplotlib.pyplot as plt

from sails import find_example_path

plt.style.use('ggplot')

ex_dir = find_example_path()

sample_rate = 2034.51 / 24
nyq = sample_rate / 2.

freq_vect = np.linspace(0, nyq, 64)

X = h5py.File(join(ex_dir, 'meg_single.hdf5'), 'r')['X'][:, 0:11000, :]

As a reminder, our data is (nsignals, nsamples, ntrials):

print(X.shape)

(1, 30517, 1)

The idea behind a sliding window analysis is to fit a separate MVAR model to short windows of the data. For the purposes of this example, we will fit consecutive windows which are 200 samples long. We can do this by simply extracting the relevant data and fitting the model in the same way that we have done so before. We are using 10 delays in each of our models. After computing the model, we can then calculate some of our diagnostics and examine, for example, the \(R^2\) value. We could also choose to examine our Fourier MVAR metrics.

from sails import VieiraMorfLinearModel, FourierMvarMetrics

delay_vect = np.arange(10)

m1 = VieiraMorfLinearModel.fit_model(X[:, 0:200, :], delay_vect)

d1 = m1.compute_diagnostics(X[:, 0:200, :])

print(d1.R_square)

0.170369632033

We can then repeat this procedure for the next, overlapping, window:

delay_vect = np.arange(10)

m2 = VieiraMorfLinearModel.fit_model(X[:, 1:201, :], delay_vect)

d2 = m1.compute_diagnostics(X[:, 1:201, :])

print(d2.R_square)

0.172925597904

This is obviously a time-consuming task, so we provide a helper function which will take a set of time-series data and compute a model and diagnostic information for each window. The sails.modelfit.sliding_window_fit() function takes at least five parameters:

the class to use for fitting (here, we use sails.modelfit.VieiraMorfLinearModel)
the data to which the models are being fitted
the delay vector to use for each model
the length of each window in samples
the step between windows. For example, if this is set to 1 sample, consecutive overlapping windows will be used. If this is set to 10, each consecutive window will start 10 samples apart. If this is equal to the length of the window or greater, the windows will not overlap.

The sliding_window_fit routine returns specially constructed model in which the parameters matrix is set up as (nsources, nsources, ndelays, nwindows). In other words, there is a model (indexed on the final dimension) for each of the windows. A ModelDiagnostics class is also returned which contains the diagnostic information for all of the windows. Use of the sliding_window_fit() routine is straightforward:

from sails import sliding_window_fit

M, D = sliding_window_fit(VieiraMorfLinearModel, X, delay_vect, 200, 1)

Once we have constructed our models, we can go ahead and construct our FourierMvarMetrics class which will allow us to, for example, extract the transfer function for each window:

F = FourierMvarMetrics.initialise(M, sample_rate, freq_vect)

We will now produce a composite plot which illustrates how the behaviour of the system evolves over time. We are going to limit ourselves to the first 10000 data points and windows.

f1, axes = plt.subplots(nrows=5, ncols=1, figsize=(12, 8))

On the top row, we plot our data.

plt.subplot(5, 1, 1)

plt.plot(X[0,:,0])

plt.xlim(0,10000); plt.grid(True)

plt.ylabel('Amplitude')

On the next two rows, we plot our transfer function for the first 10000 windows:

plt.subplot(5, 1, (2, 3))

plt.contourf(M.time_vect, F.freq_vect[3:], F.H[0, 0, 3:, :])

plt.xlim(0,10000); plt.grid(True)

plt.ylabel('Frequency (Hz)')

Underneath the transfer function, we examine the stability index (in red) and the \(R^2\) value (in blue):

plt.subplot(5, 1, 4)

plt.plot(D.SI)

plt.plot(D.R_square)

plt.xlim(0, 10000); plt.grid(True)

plt.ylabel('SI / $R^2$')

On the bottom row, we plot the AIC value for each of the models:

plt.subplot(5, 1, 5)

plt.plot(D.AIC)

plt.xlim(0,10000); plt.grid(True)

plt.xlabel('Time (samples)')
plt.ylabel('AIC')

Finally, we can look at our overall figure:

f1.tight_layout()

plt.show()

Tutorial 7 - Plotting Helpers¶

In this tutorial, we will explore the plotting helper functions which are available for use in sails; primarily for plotting netmats and circular connectivity pots.

For this tutorial, we will use some example plot definition files which are provided with sails.

We start by importing our modules and finding and finding our example path

from os.path import join

from sails import find_example_path

ex_dir = find_example_path()

group_csv = join(ex_dir, 'aal_cortical_plotting_groups.csv')
region_csv = join(ex_dir, 'aal_cortical_plotting_regions.csv')

The two files which are imported above describe the layout of the plots which we will create. The example files given are for a specific cortical subset of the AAL atlas which we have used in some of our own work; this parcellation contains 78 regions; 39 per hemisphere.

The groups.csv file contains a list of the groups of regions which we will plot. In this case, we have divided the 78 regions into 10 groups:

Frontal left / right

Medial left / right

Temporal left / right

Parietal left / right

Occipital left / right

The regions.csv file describes each ROI and places it within these groups. It also provides the indexing information so that we can match between our netmats and the plots.

More details of both of these files can be found in the docstring for the from_csv_files function of the sails.CircosHandler class.

We will now build a connectivity diagram using this structure.

Setting up the handler¶

We start by setting up a handler which parses the CSV files:

from sails import CircosHandler

c = CircosHandler.from_csv_files(group_csv, region_csv)

There is one additional argument available to the from_csv_files function; this argument is analysis_column. By default, the order in which our analysed data is stored comes from the AnalysisID column of the regions.csv spreadsheet. You can alter which column is used for this by setting the analysis_column argument to the function. This is useful in cases where you have data in different orders from different analyses but where you want to produce plots with the same ordering with little effort.

Extracting Region Ordering and Indices¶

We can use use this handler to extract various orderings of our ROIs and the indices needed to re-order data into them.

# Show the order in circular form (clockwise) of the groups
print(c.circular_groups)

['frontal_L', 'medial_L', 'temporal_L', 'parietal_L', 'occipital_L', 'occipital_R', 'parietal_R', 'temporal_R', 'medial_R', 'frontal_R']

# Show the order in netmat form (top to bottom) of the groups
print(c.netmat_groups)

['frontal_L', 'medial_L', 'temporal_L', 'parietal_L', 'occipital_L', 'frontal_R', 'medial_R', 'temporal_R', 'parietal_R', 'occipital_R']

# Show the last 10 regions in circular order and the corresponding
# indices into our data
print(c.circular_regions[-10:])
print(c.circular_indices[-10:])

['Frontal_Sup_Medial_R', 'Frontal_Inf_Tri_R', 'Frontal_Inf_Oper_R', 'Frontal_Mid_R', 'Frontal_Sup_R', 'Frontal_Inf_Orb_R', 'Frontal_Mid_Orb_R', 'Frontal_Med_Orb_R', 'Frontal_Sup_Orb_R', 'Rectus_R']
[28, 19, 15, 25, 31, 17, 24, 21, 30, 61]

# Show the last 10 regions in netmat order and the corresponding
# indices into our data
print(c.netmat_regions[-10:])
print(c.netmat_indices[-10:])

['SupraMarginal_R', 'Rolandic_Oper_R', 'Precuneus_R', 'Occipital_Sup_R', 'Occipital_Mid_R', 'Occipital_Inf_R', 'Calcarine_R', 'Cuneus_R', 'Lingual_R', 'Fusiform_R']
[67, 63, 59, 45, 43, 41, 5, 13, 39, 33]

Circular Plots¶

To generate circular plots, we use the program circos. You will need to have the program installed on your computer in order to generate such plots. circos is available from http://circos.ca/software/. On Debian and similar, you can simply apt install circos circos-data libsvg-perl.

Note that if you are creating plots for publication using Circos, you should cite the relevant publication: Krzywinski, M. et al. Circos: an Information Aesthetic for Comparative Genomics. Genome Res (2009) 19:1639-1645.

Generating a karyotype¶

Circos uses the term karyotype to describe the ordering of chromosomes and bands within them. In our case, we are using this to describe groups of regions and individual ROIs respectively. We only need to generate a karyotype file once for each layout; we do not need a different karyotype file for each individual plot.

To generate a karytotype file, we open a file and write the contents of the karyotype() function into it.

f = open('aal_karyotype.txt', 'w')
f.write(c.karyotype())
f.close()

In our example case, the start of the file will look like this:

chr - frontal_L 1 0 12000 red
chr - medial_L 2 0 6000 yellow
chr - temporal_L 3 0 6000 purple
chr - parietal_L 4 0 8000 green
chr - occipital_L 5 0 7000 blue
chr - occipital_R 6 0 7000 vvdblue
chr - parietal_R 7 0 8000 vvdgreen
chr - temporal_R 8 0 6000 vvdpurple
chr - medial_R 9 0 6000 vvdyellow
chr - frontal_R 10 0 12000 vvdred
band frontal_L Rectus_L Rectus_L 0 1000 grey
band frontal_L Frontal_Sup_Orb_L Frontal_Sup_Orb_L 1000 2000 grey

We will use this file in the next section to make our plots.

Generating a connectivity plot based on the karyotype¶

We start by generating some controlled data with connections between only a few regions

import numpy as np

data = np.zeros((78, 78)) * np.nan

# Add a strong positive connection between Amygdala_L and Insula_R
data[0, 37] = 10.0
# Add a weaker negative connection between Cuneus_R and Heschl_R
data[13, 35] = -4.0

Note that when the data values are used in the circular plotting routines, they will be used as line widths in pixels.

We now generate a set of circos configuration files from this data. We also need to pass it our karyotype file name and the output base name.

c.gen_circos_files(data, 'aal_karyotype.txt', 'aal_test_out')

The code above assumes that you have the circos config files installed in /etc/circos. If you have them in another location, pass the circos_path variable to the gen_circos_files routine, e.g.:

c.gen_circos_files(data, 'aal_karyotype.txt', 'aal_test_out',
                   circos_path='/usr/local/etc/circos')

This will generate two files: aal_test_out.conf and aal_test_out.txt. The former is the circos configuration file and the latter is the file which contains the actual information about the connections. In our case, we can see that the latter only contains two lines; one for each of our two connections.

Manually Generating the plot¶

To manually generate the plot from the configuration files, we use a normal shell and run the circos command.

circos -conf aal_test_out.conf

You may get an error which contains the following:

*** CIRCOS ERROR ***

cwd: /tmp

command: /usr/bin/circos -conf aal_test_out.conf

CONFIGURATION FILE ERROR

This is because of a problem with circos finding some of its configuration files. You can fix this by running the following commands from your shell (whilst in the directory containing the files). This assumes that the circos config files are in /etc/circos:

mkdir etc
    ln -s /etc/circos/tracks etc/

You should then re-run the circos command.

Two files will be created: test_out.png and test_out.svg. The image should look like this:

Automatically Generating the plot¶

Instead of going to the effort of the above, there is a helper routine gen_circos_plot which will create a temporary directory, generate all of the configuration files, run circos and then copy the plots to your final destination. It can be used directly in place of gen_circos_files. If the circos call fails, an exception will be raised which will include the contents of stdout and stderr so that you can investigate the problem.

Note that you do not pass a karyotype file name to the gen_circos_plot routine as it will generate a karyotype file in the temporary directory for you.

c.gen_circos_plot(data, 'aal_test_out',
                  circos_path='/usr/local/etc/circos')

Modifying link colours¶

By default, links with a negative value will be shown in blue, links with a positive value in red and links with a strictly zero value in black. You can modify this in two ways:

By setting a tuple for (zero_colour, negative_colour, positive_colour)

By passing a numpy matrix with dtype object containing a string for each connection. In this case, you can also pass a defaultcolor which will be used if a matrix entry is None

An example of using the tuple syntax:

c.gen_circos_files(data, 'aal_karyotype.txt', 'aal_test_out_2',
                   linkcolors=('black', 'green', 'purple'))

and an example of using the matrix syntax:

# As an example here, we set the colour of the Amygdala/Insula link
# explicitly and set the other link using the default color syntax
colors = np.empty((78, 78), dtype=object)

colors[0, 37] = 'yellow'

c.gen_circos_files(data, 'aal_karyotype.txt', 'aal_test_out_3',
                   linkcolors=colors,
                   defaultlinkcolor='orange')

Netmat Plots¶

In this section, we will discuss how to use the netmat plotting routines which come as part of the CircosHandler above. Further down there is documentation on how to use the raw plot_netmat routine.

We assume that we are in the same session as above and that c still represents a CircosHandler object.

# Set up some netmat data
netdata = np.zeros((78, 78))

# Link from Calcarine_L (4) to Heschl_R (35) with a positive score
netdata[4, 35] = 1.0

# Link from Calcarine_L (4) to Heschl_L (34) with a negative score
netdata[4, 34] = -1.0

# Link from Frontal_Mid_L (22) to Fusiform_L (32) with a double positive socre
netdata[22, 32] = 2.0

# Create our plot - a figure and axes will be created as we have
# not supplied any

# Note that, as demonstrated here, you can use named arguments to pass
# extra options into the pcolormesh call.  See the docstring of
# the CircosHandler.plot_netmat function for more help
ax = c.plot_netmat(netdata, label_fontsize=4, cmap='RdBu_r', vmin=-2.0, vmax=2.0)

Raw Netmat Plots¶

If you do not wish to set up a full set of CSV files etc to plot your netmats, you can call the sails.plotting.plot_netmat routine having prepared your arguments yourself.

You will need to set up colour mappings and lists as well as (optionally) labels for your regions. An example is shown below which uses the same data as above but changes the colour layout slightly:

from sails import plot_netmat

    # These original colours were taken as RGB tuples from circos
orig_colors = {'red':        (217, 120, 99),
               'orange':     (217, 144, 89),
               'yellow':     (220, 213, 103),
               'purple':     (155, 152, 183),
               'green':      (127, 180, 128),
               'blue':       (121, 166, 193),
               'vvdred':     (152, 49, 58),
               'vvdorange':  (143, 79, 52),
               'vvdyellow':  (178, 170, 49),
               'vvdpurple':  (99, 62, 139),
               'vvdgreen':   (49, 109, 82),
               'vvdblue':    (54, 95, 148)}

# We need to convert them to to matplotlib compatible RGB 0-1 tuples
colors = {}

    for cname, cval in orig_colors.items():
        colors[cname] = [(x / 255.0) for x in cval]

    # This is the layout of our AAL parcellation division
    cnames = ['red']       * 12 + \
         ['yellow']    * 6 + \
         ['purple']    * 6 + \
         ['green']     * 8 + \
         ['blue']      * 7 + \
         ['vvdred']    * 12 + \
         ['vvdyellow'] * 6 + \
         ['vvdpurple'] * 6 + \
         ['vvdgreen']  * 8 + \
         ['vvdblue']   * 7

tick_pos = [12, 14, 18, 24, 32, 39, 51, 53, 57, 63, 71]

    # Some temporary labels for testing
    labels = ['R{}'.format(x) for x in range(1, 79)]

    # Note that in reality, we would have to re-order our data to be
# in the correct order for the netmat here (in the above example,
# this is handled by the CircosHandler class).
# Here we leave it unordered, which is why the plot ends up
# different (and wrong given what we claimed we were doing with
# the data above)

    plot_netmat(netdata, colors=colors, cnames=cnames, tick_pos=tick_pos,
        labels=labels, label_fontsize=4, cmap='RdBu_r', vmin=-2.0, vmax=2.0)

Gallery¶

Simulation Gallery¶

This page contains a set of figures re-created around existing papers which used simulated or analytic systems. These graphs are provided for comparison with the original papers.

The code for generating these figures can be found in sim_gallery.py.

Baccala and Sameshima 2001¶

We re-create figures 1 to 5 of this paper. Each figure compares PDC and DTF connectivity estimators.

Pascual-Marqui et al 2014¶

We re-create figure 3 of this paper using IEC and PDC. PDC and DTF connectivity estimators.

Module reference¶

Model fitting and Diagnostics¶

Model fitting classes¶

class sails.modelfit.VieiraMorfLinearModel[source]¶

TODO: Description

classmethod fit_model(data, delay_vect)[source]¶

Estimates the multichannel autoregressive spectral estimators using a Vieira-Morf algorithm. Equations are referenced to Marple, appendix 15.B.

This is the multitrial versions of the algorithm, using the AMVAR method outlined in Ding 2000.

Parameters:

data (numpy.ndarray) – array of shape [nsignals, nsamples, ntrials]
delay_vect (numpy.ndarray) – Vector containing evenly spaced delays to be assessed in the model. Delays are represented in samples. Must start with zero.
assess_order (bool) – Optional True/False value indicating whether additional information about each iteration of the fit should be returned. Defaults to False

Returns:

A populated object containing the fitted forward and backwards coefficients and several other useful variables and methods.

get_residuals(data, forward_parameters=True)[source]¶

Returns the prediction error from a fitted model. This is a wrapper function for get_residuals()

Parameters:	forward_parameters (bool) – If True, use forward parameters, otherwise use backward parameters

class sails.modelfit.OLSLinearModel[source]¶

TODO: Description

classmethod fit_model(data, delay_vect)[source]¶: TODO: Description

Model fitting helper functions¶

sails.modelfit.sliding_window_fit(model_class, data, delay_vect, win_len_samp, win_step_samp, compute_diagnostics=True)[source]¶

Model diagnostics¶

class sails.diags.DelayDiagnostics[source]¶

Class which computes the mutual information as a function of lag from zero lag to the first zero crossing of the autocorrelation function.

classmethod delay_search(data, maxdelay, step, sample_rate, constant_window=True)[source]¶

Compute MI as a function of lag from zero lag to the first zero crossing of the autocorrelation function

If NafOptions().verbose is True, progress will be reported to stdout.

Parameters:	data (numpy.ndarray) – array of [signals, samples, trials] maxdelay (int) – maximum delay to consider step (int) – step to increment the delay by sample_rate (float) – sample rate of data constant_window (bool) – Flag indicating that the same number of datapoints should be included at each delay. Default is True

class sails.diags.ModelDiagnostics[source]¶: TODO: Description

sails.modelfit.get_residuals(data, parameters, delay_vect, backwards=False)[source]¶

TODO: Description

Parameters:	data – TODO parameters – TODO delay_vect – TODO backwards – TODO
Returns:	TODO

Model Decomposition¶

class sails.modal.MvarModalDecomposition[source]¶

Object which calculates the modal decomposition of a fitted model derived from the AbstractLinearModel class.

compute_modal_parameters()[source]¶: Compute the order 1 or 2 autoregressive coefficients associated with each pole or conjugate pair.

compute_residues()[source]¶: Compute the residue matrix from the eigenvectors associated with each pole.

excitation(resid_cov)[source]¶: http://climate-dynamics.org/wp-content/uploads/2015/08/arfit.pdf

classmethod initialise(model, sample_rate=None, normalise=False)[source]¶

Compute the modal pole-residue decomposition of the A matrix from a fitted MVAR model.

Currently only implemented for a single realisation (A.ndim == 3)

Parameters:	model – Fitted model object. Must be derived from AbstractLinearModel. sample_rate – sample rate of system in Hz normalise – Whether to adjust the phase of the eigenvectors
Returns:	Modal decomposition object

modal_freqz()[source]¶: Compute filter characteristics for each pole

modal_timecourse(data)[source]¶

Compute the modal time-course from a dataset and a fitted model. This is the tranformation of a set of [nsignals x time] data observations into modal coordinates of shape [nmodes x time].

Parameters:	data – TODO

modal_transfer_function(sample_rate, modes=None)[source]¶

Compute a transfer function matrix based on the excitation period for each mode

Parameters:	sample_rate – sample rate of system in Hz modes – List of mode indices to evaluate over (optional)

mode_indices¶

Returns a list of tuples of mode indices where each tuple contains the indices of either a pole-pair or an unpaired pole

Returns:	list of tuples containing indices into the modes

per_mode_transfer_function(sample_rate, freq_vect)[source]¶

Extracts the transfer function for each mode by summing the transfer function across pole-pairs where necessary

Parameters:	sample_rate – sample rate of system in Hz freq_vect – Vector of frequencies at which to evaluate function

period_hz¶: This is an old and deprecated way of accessing peak_frequency

pole_plot(ax=None, normscale=50, plottype='unitcircle')[source]¶

Plot the poles of the current system.

Parameters:	ax – Optional axis on which to plot normscale – Scaling factor for points on plot plottyle – String defining which plot to make, can be ‘unitcircle’ or ‘eigenspectrum’
Returns:	Reference to axes on which plot was drawn

residue_norm¶: TODO

transfer_function(sample_rate, freq_vect, modes=None, sum_modes=True)[source]¶

Compute the transfer function in pole-residue form, splitting the system into modes with separate transfer functions. The full system transfer function is then a linear sum of each modal transfer function.

Parameters:	sample_rate – sample rate of system in Hz freq_vect – Vector of frequencies at which to evaluate function modes – List of mode indices to evaluate over (optional) sum_modes – Boolean indicating whether to sum modes or return all

Model metric estimatation¶

class sails.mvar_metrics.FourierMvarMetrics[source]¶

TODO: Document

classmethod initialise(model, sample_rate, freq_vect, nmodes=None)[source]¶: Compute some basic values from A across our frequency vector

class sails.mvar_metrics.ModalMvarMetrics[source]¶

classmethod initialise(model, sample_rate, freq_vect, nmodes=None, sum_modes=True)[source]¶

Compute some basic values from the pole-residue decomposition of the A matrix

Currently only implemented for a single realisation (A.ndim == 3)

Parameters:	model – TODO sample_rate – TODO freq_vect – TODO nmodes – TODO

sails.mvar_metrics.modal_transfer_function(evals, evecl, evecr, nchannels, sample_rate=None, freq_vect=None)[source]¶: Comute the transfer function in pole-residue form, splitting the system into modes with separate transfer functions. The full system transfer function is then a linear sum of each modal transfer function.

Plotting¶

sails.plotting.root_plot(rts, ax=None, figargs={}, plotargs={})[source]¶

Plot a set of roots (complex numbers).

Parameters:	rts – Roots to plot ax – Optional Axes on which to place plot.
Returns:	Axes object on which plot was drawn

sails.plotting.plot_vector(metric, x_vect, y_vect=None, x_label=None, y_label=None, title=None, labels=None, line_labels=None, F=None, cmap=<matplotlib.colors.LinearSegmentedColormap object>, triangle=None, diag=False, thresh=None, two_tail_thresh=False, font_size=10, use_latex=False)[source]¶

Function for plotting frequency domain connectivity at a single time-point.

Params metric:	matrix containing connectivity values [nsignals x signals x frequencies x participants] in which the first dimension refers to source nodes and the second dimension refers to target nodes
Params x_vect:	vector of frequencies to label the x axis
Params y_vect:	y_vect: vector containing the values for the y-axis
Parameters:	x_label (string [optional]) – label for the x axis y_label (string [optional]) – label for the y axis title (string [optional]) – title for the figure labels (list) – list of node labels for columns and vectors line_labels (list) – list of labels for each separate line (participant dimension in metric) F (figurehandle [optional]) – handle of existing figure to plot within triangle (string [optional]) – string to indicate whether only the ‘upper’ or ‘lower’ triangle of the matrix should be plotted diag (bool [optional]) – flag to indicate whether the diagonal elements should be plotted thresh (ndarray [optional]) – matrix containing thesholds to be plotted alongside connectivity values [nsignals x nsignals x frequencies] two_tailed_thresh (bool [optional]) – flag to indicate whether both signs (+/-) of the threshold should be plotted font_size (int [optional]) – override the default font size

sails.plotting.plot_matrix(metric, x_vect, y_vect, x_label=None, y_label=None, z_vect=None, title=None, labels=None, F=None, cmap=<matplotlib.colors.LinearSegmentedColormap object>, font_size=8, use_latex=False, diag=True)[source]¶

Function for plotting frequency domain connectivity over many time points

Parameters:

metric (ndarray) – matrix containing connectivity values [nsignals x signals x frequencies x participants] in which the first dimension refers to source nodes and the second dimension refers to target nodes
x_vect (1d array) – vector of frequencies to label the x axis
y_vect (1d array [optional]) – vector containing the values for the y-axis
z_vect (1d array [optional]) – vector containing values for the colour scale
x_label (string [optional]) – label for the x axis
y_label (string [optional]) – label for the y axis
title (string [optional]) – title for the figure
labels (list) – list of node labels for columns and vectors
F (figurehandle [optional]) – handle of existing figure to plot within
cmap (matplotlib colormap [optional]) – matplotlib.cm.<colormapname> to use for colourscale (redundant for plot vector??)
font_size (int [optional]) – override the default font size

Simulation¶

class sails.simulate.AbstractSigGen(extraargs=None)[source]¶

Abstract base class for signal generation

generate_signal(meg_reader)[source]¶

Generate a custom signal for embedding in a dataset

Parameters:	meg_reader (naf.meg.abstractdata.AbstractMEGReader) – MEG data for which to generate signal
Return type:	numpy.ndarray
Returns:	numpy array of shape (samples). Note that samples must be equal to the value of meg_reader.num_samples

get_num_sources()[source]¶

Return the number of sources which this class provides.

If not override, returns 1

class sails.simulate.AbstractLagSigGen(extraargs=None)[source]¶

Generic class for making lagged signals

generate_basesignal(nsamples)[source]¶: Function defining the starting signal which will be lagged throughout the network

generate_noise(nsamples)[source]¶: Function defining the additive noise added at each lag

generate_signal(sample_rate=None, num_samples=None, noise_dB=None)[source]¶

Method for generating a set of time courses based on either the structure of a meg_reader object or passed parameters and the self.generate_parameters method defined in this class

Returns a timeseries of dimension [samples x signals]

get_lags()[source]¶: Method returning the lagging of the signal through the system. Should return a [nsignals x nsignals x nlags] matrix in which the first dimension refers to the source node and the second dimension the target node. Each non-zero value indicates that the signal in the first dimension will be added into the signal in the second dimension at the specified weight.

get_num_sources()[source]¶

Return the number of sources which this class provides.

If not override, returns 1

class sails.simulate.AbstractMVARSigGen(extraargs=None)[source]¶

Generic class for making MVAR signals

check_stability()[source]¶: Checks the stablility of the parameter matrix by the magnitude of the largest eigenvalue of the first lag in the parameter matrix. Lutkepohl 2005 ch 1 for more detail

generate_model()[source]¶: Returns a LinearModel containing containing the true model parameters.

generate_parameters()[source]¶: Method to generate the parameters of the MVAR system. Should return a [nsignals x nsignals x nlags] parameter matrix in which the first dimension refers to the source node and the second dimension the target node

generate_signal(sample_rate=None, num_samples=None, noise_dB=None, num_realisations=1)[source]¶

Method for generating a set of time courses based on either the structure of a meg_reader object or passed parameters and the self.generate_parameters method defined in this class

Returns a timeseries of dimension [samples x signals]

get_num_sources()[source]¶

Return the number of sources which this class provides.

If not override, returns 1

class sails.simulate.Baccala2001_fig1(extraargs=None)[source]¶

https://doi.org/10.1007/PL00007990

generate_parameters()[source]¶: Method to generate the parameters of the MVAR system. Should return a [nsignals x nsignals x nlags] parameter matrix in which the first dimension refers to the source node and the second dimension the target node

get_num_sources()[source]¶

Return the number of sources which this class provides.

If not override, returns 1

class sails.simulate.Baccala2001_fig2(extraargs=None)[source]¶

https://doi.org/10.1007/PL00007990

generate_parameters()[source]¶: Method to generate the parameters of the MVAR system. Should return a [nsignals x nsignals x nlags] parameter matrix in which the first dimension refers to the source node and the second dimension the target node

get_num_sources()[source]¶

Return the number of sources which this class provides.

If not override, returns 1

class sails.simulate.Baccala2001_fig3(extraargs=None)[source]¶

https://doi.org/10.1007/PL00007990

generate_parameters()[source]¶: Method to generate the parameters of the MVAR system. Should return a [nsignals x nsignals x nlags] parameter matrix in which the first dimension refers to the source node and the second dimension the target node

get_num_sources()[source]¶

Return the number of sources which this class provides.

If not override, returns 1

class sails.simulate.Baccala2001_fig4(extraargs=None)[source]¶

https://doi.org/10.1007/PL00007990

generate_parameters()[source]¶: Method to generate the parameters of the MVAR system. Should return a [nsignals x nsignals x nlags] parameter matrix in which the first dimension refers to the source node and the second dimension the target node

get_num_sources()[source]¶

Return the number of sources which this class provides.

If not override, returns 1

class sails.simulate.Baccala2001_fig5(extraargs=None)[source]¶

https://doi.org/10.1007/PL00007990

generate_parameters()[source]¶: Method to generate the parameters of the MVAR system. Should return a [nsignals x nsignals x nlags] parameter matrix in which the first dimension refers to the source node and the second dimension the target node

get_num_sources()[source]¶

Return the number of sources which this class provides.

If not override, returns 1

class sails.simulate.Fasoula2013_eqn26(extraargs=None)[source]¶

https://doi.org/10.1016/j.jneumeth.2013.02.021

generate_parameters()[source]¶: Method to generate the parameters of the MVAR system. Should return a [nsignals x nsignals x nlags] parameter matrix in which the first dimension refers to the source node and the second dimension the target node

get_num_sources()[source]¶

Return the number of sources which this class provides.

If not override, returns 1

class sails.simulate.Schelter2006_fig1(extraargs=None)[source]¶

https://doi.org/10.1016/j.jneumeth.2005.09.001

generate_parameters()[source]¶: Method to generate the parameters of the MVAR system. Should return a [nsignals x nsignals x nlags] parameter matrix in which the first dimension refers to the source node and the second dimension the target node

get_num_sources()[source]¶

Return the number of sources which this class provides.

If not override, returns 1

class sails.simulate.PascualMarqui2014_fig3(extraargs=None)[source]¶

https://dx.doi.org/10.3389/fnhum.2014.00448

generate_parameters()[source]¶: Method to generate the parameters of the MVAR system. Should return a [nsignals x nsignals x nlags] parameter matrix in which the first dimension refers to the source node and the second dimension the target node

get_num_sources()[source]¶

Return the number of sources which this class provides.

If not override, returns 1

class sails.simulate.Korzeniewska2003Lag(extraargs=None)[source]¶

Generate realisations of the time-lagged network defined in Korzeniewska 2003 figure 1.

https://doi.org/10.1016/S0165-0270(03)00052-9

generate_basesignal(nsamples)[source]¶

Generate oscillatory signal to be propogated through the netowrk. This creates a single resonance at one-half of the Nyqvist frequency via direct pole placement.

Note: the original base-signal in the paper is a realisation from an AR model. I couldn’t find the parameters for this mode, so use pole placement here. The spectral profile will be different to the paper, but the connectivity patterns are the same.

get_lags()[source]¶: Method returning the lagging of the signal through the system. Should return a [nsignals x nsignals x nlags] matrix in which the first dimension refers to the source node and the second dimension the target node. Each non-zero value indicates that the signal in the first dimension will be added into the signal in the second dimension at the specified weight.

get_num_sources()[source]¶

Return the number of sources which this class provides.

If not override, returns 1

class sails.simulate.SailsTutorialExample(extraargs=None)[source]¶

generate_basesignal(f, r, sample_rate, num_samples)[source]¶: Generate a simple signal by pole placement

generate_signal(f1, f2, sample_rate, num_samples, noise_dB=None, num_realisations=1, magnitude_noise=0.02)[source]¶

Method for generating a set of time courses based on either the structure of a meg_reader object or passed parameters and the self.generate_parameters method defined in this class

Returns a timeseries of dimension [samples x signals x realisations]

get_num_sources()[source]¶

Return the number of sources which this class provides.

If not override, returns 1

Support routines¶

Abstract Classes¶

class sails.modelfit.AbstractLinearModel[source]¶

TODO: Description

companion¶: Parameter matrix in companion form

compute_diagnostics(data, compute_pc=False)[source]¶

Compute several diagnostic metrics from a fitted MVAR model and a dataset.

Return type:	diags.ModelDiagnostics
Returns:	ModelDiagnostics object containing LL, AIC, BIC, stability ratio, durbin-watson, R squared and percent consistency measures

get_residuals(data)[source]¶

Returns the prediction error from a fitted model. This is a wrapper function for get_residuals()

Parameters:	data – TODO

nsignals¶: Number of signals in fitted model

order¶: Order of fitted model

class sails.mvar_metrics.AbstractMVARMetrics[source]¶

TODO: Description

PSD¶: Power Spectral Density matrix

S¶: Spectral Matrix of shape [nsignals, nsignals, order, epochs]

coherency¶: TODO: Document

d_directed_transfer_function¶: TODO: Document

directed_transfer_function¶: TODO: Document

ff_directed_transfer_function¶: TODO: Document

geweke_granger_causality¶: TODO: Document

imaginary_coherence¶: TODO: Document

inv_S¶: TODO: Document

isolated_effective_coherence¶: TODO: Document

magnitude_squared_coherence¶: TODO: Document

partial_coherence¶: TODO: Document

partial_directed_coherence¶: TODO: Document

phase_coherence¶: TODO: Document

Statistical Functions¶

sails.mvar_metrics.sdf_spectrum(A, sigma, sample_rate, freq_vect)[source]¶

Estimate of the Spectral Density as found on wikipedia and Quirk & Liu 1983

This assumes that the spectral representation of A is invertable

sails.mvar_metrics.psd_spectrum(A, sigma, sample_rate, freq_vect)[source]¶

Estimate the PSD representation of a set of MVAR coefficients as stated in Penny 2000 Signal processing course array.ps 7.41

This assumes that the spectral representation of A is invertable

TODO: This doesn’t behave as expected for univariate?, sdf_spectrum recommended

Parameters:	A – TODO sigma – TODO sample_rate – TODO freq_vect – TODO
Returns:	TODO

sails.mvar_metrics.ar_spectrum(A, sigma, sample_rate, freq_vect)[source]¶

Estimate the spectral representation of a set of MVAR coefficients as suggested by Baccala & Sameshima 2001, the equation without a number just below equation 13.

Parameters:	A – TODO sigma – TODO sample_rate – TODO freq_vect – TODO
Returns:	TODO

sails.mvar_metrics.transfer_function(Af)[source]¶

Function for computing the transfer function of a system

Parameters:	Af – Frequency domain version of parameters (can be calculated using ar_spectrum function) [nchannels x nchannels x nfreqs x nepochs]
Returns:	numpy array of dimensions [nsignals x nsignals x frequency x epochs]

sails.mvar_metrics.spectral_matrix(H, noise_cov)[source]¶

Function for computing the spectral matrix, the matrix of spectra and cross-spectra.

Parameters:	H – The transfer matrix [nchannels x nchannels x nfreqs x nepochs] noise_cov – The noise covariance matrix of the system [nchannels x nchannels x nepochs]
Returns:	numpy array of dimensions [nsignals x nsignals x frequency x epochs]

sails.mvar_metrics.coherency(S)[source]¶

Method for computing the Coherency. This is the complex form of coherence, from which metrics such as magnitude squared coherence can be derived

Params S:	The spectral matrix: 3D or 4D [nchannels x nchannels x nfreqs x nepochs]
Returns:	numpy array of dimensions [nsignals x nsignals x frequency x epochs]

sails.mvar_metrics.partial_coherence(inv_S)[source]¶

Method for computing the Partial Coherence.

Params inv_S:	The inverse spectral matrix [nchannels x nchannels x nfreqs x nepochs]
Returns:	numpy array of dimensions [nsignals x nsignals x frequency x epochs]

sails.mvar_metrics.partial_directed_coherence(Af)[source]¶

Function to estimate the partial directed coherence from a set of multivariate parameters.

Parameters:	Af – Frequency domain version of parameters (can be calculated using ar_spectrum function)
Returns:	PDC : ndarray containing the PDC estimates from the parameters of dimensions [nchannels x nchannels x order]

sails.mvar_metrics.isolated_effective_coherence(Af, noise_cov)[source]¶

Function for estimating the Isolated Effective Coherence as defined in Pascual-Marqui et al 2014

Parameters:	Af – Frequency domain version of parameters (can be calculated using ar_spectrum function) [nchannels x nchannels x nfreqs x nepochs] noise_cov – The noise covariance matrix of the system [nchannels x nchannels x nepochs]
Returns:	iec: numpy array of dimensions [nsignals x nsignals x frequency x epochs]

sails.mvar_metrics.directed_transfer_function(H)[source]¶

Method for computing the Directed Transfer Function as defined in Kaminski & Blinowska 1991

Parameters:	H – The transfer matrix [nchannels x nchannels x nfreqs x nepochs]
Returns:	numpy array of dimensions [nsignals x nsignals x frequency x epochs]

sails.modal.adjust_phase(x)[source]¶

Method for normalising eigenvector matrices to meet the conditions in eqn 21 of Neumaier & Schneider 2001.

Parameters:	x – Vector of complex numbers to adjust the phase of
Returns:	Vector of numbers with adjusted phase

sails.stats.find_cov_multitrial(a, b, ddof=1)[source]¶

Estimate the average covariance between two datasets across many trials

Parameters:	a (ndarray) – a multivariate time series of shape nsignals x nsamples x ntrials b (ndarray) – a multivariate time series of shape nsignals x nsamples x ntrials ddof (None orint) – if None, no normalisation will be performed. If an integer, covariance will be normalised by (nsamples - ddof)
Return type:	ndarray
Returns:	A covariance matrix between a and b of shape nsignals x nsignals

sails.stats.find_cov(a, b, ddof=1)[source]¶

Estimate the covariance between two datasets.

Parameters:	a (ndarray) – a multivariate time series of shape nsignals x nsamples b (ndarray) – a multivariate time series of shape nsignals x nsamples ddof (int) – if None, no normalisation will be performed. If an integer, covariance will be normalised by (nsamples - ddof)
Return type:	ndarray
Returns:	A covariance matrix between a and b of shape nsignals x nsignals

sails.stats.mutual_information(A, B, r=None, nbins=21, base=<ufunc 'log2'>)[source]¶

Calculate the mutual information between two signals.

Parameters:	A (ndarray) – First signal, of shape [channels x samples] B (ndarray) – Second signal, of shape [channels x samples] r (tuple (float, float)) – Range of the random variables used for bin calculations. If not set, takes the minimum and maximum values from the combination of A and B nbins (int) – Number of bins to compute the frequencies within base (function) – Base to use for MI calculation - defaults to np.log2
Return type:	ndarray
Returns:	the mutual information within each signal passed

sails.stats.durbin_watson(residuals, step=1)[source]¶

Test for autocorrelation in the residuals of a model.

The result is a value between 0 and 4. A value of 0 indicates a positive autocorrelaiton, 2 indicates no correlation and a value of 4 indicates a negative autocorrelation.

If many trials are found the statistic is estimated for each trial separately.

Parameters:	residuals (ndarray) – Residuals of shape [channels, samples, ntrials] If only two dimensions, will be treated as [channels, samples] If only one dimension, will be treated as [samples] step (int) – step-size to use through residuals when calculating DW. Defaults to 1.
Return type:	ndarray
Returns:	durbin-watson index

sails.stats.percent_consistency(data, residuals)[source]¶

Estimate how well a model retains the correlational structure of a dataset.

This uses the percent consistency measure outlined in Ding et al 2000. All auto and cross correlations within a multivariate data set are computed and compared to the auto and cross correlations of the data as predicted by a model (in this case the data - model residuals, as these are easily obtained).

Parameters:	data (numpy.ndarray) – An array containing the observed data, of shape [nsignals x nsamples x ntrials] residuals – An array containing the remaining variance of data after model fitting, of shape [nsignals x nsamples x ntrials]
Rtype float:
Returns:	A value between 0 and 100 indicating the percentage of the autocorrelation in the data that is retained in the model.

sails.stats.rsquared_from_residuals(data, residuals, per_signal=False)[source]¶

Estimate the variance explained by a model from the original data and the model residuals

If data and residuals are 2D, they will be treated as [channels, nsamples]

If data and residuals are 1D, they will be treated as [nsamples]

Parameters:	data (ndarray) – Array containing the original data, of shape [channels, nsamples, ntrials] residuals (ndarray) – Array containing the model residuals, of shape [channels, nsamples, ntrials] per_signal (bool) – Boolean indicating whether the variance explained should be computed separately per channel
Return type:	float
Returns:	Value between 0 and 1 indicating the proportion of variance explained by the model.

sails.stats.mahalanobis(X, Y=None, sigma=None)[source]¶

Estimate the mahalanobis distance of each point in an array across the samples dimension.

Parameters:	X (ndarray) – Array containing data of points: shape [channels x samples] Y (ndarray or None) – Optional Array containing data of second points: shape [channels x samples]. If not provided, the origin will be used as the point for comparison. sigma (ndarray or None) – A precomputed covariance matrix to use in the calculation, of shape [channels x channels]. If not provided, sigma will be calculated from X
Return type:	ndarray
Returns:	The mahalanobis distance for each sample, of shape [samples]

sails.stats.profile_likelihood(X)[source]¶

Find the ‘elbow’ point in a curve.

This function uses the automatic dimensionality selection method from Zhu et al 2006 to estimate the elbow or inflection point in a given curve.

Method can be found in Mu Zhu, Ali Ghodsi. Computational Statistics & Data Analysis 51 (2006) 918 - 930

Parameters:	X (ndarray) – 1d array containing the curve to be analysed, of shape [samples]
Return type:	ndarray
Returns:	The log likelihood that each point on the curve separates the distributions to the left and right, of shape [samples]

sails.stats.stability_index(A)[source]¶

Compute the stability index as defined in Lutkepohl 2006 pg 11.

This is a proxy for the stationarity assumption of MVAR models, if the magnitude of the largest principal component of the model parameters is less than 1, the model is stable and thus, stationary.

Parameters:	a (numpy.ndarray) – The MVAR parameter matrix
Return type:	float
Returns:	The stability index

sails.stats.stability_ratio(A)[source]¶

Compute the stability ratio as a measure of stationarity

A stronger test for stationarity than the stability index. Computes the ratio between the largest eigenvalue of the parameters and all the others. If the largest parameter is larger than all other parmeters (indicated by a value < 1) we have a super-stable system.

Parameters:	a (numpy.ndarray) – The MVAR parameter matrix
Rtype float:
Returns:	The stability ratio

sails.modelvalidation.approx_model_selection_criteria(A, EF, method='ss')[source]¶

Estimate Akaike’s Information Criterion and Bayes Information Criterion from a model fit using an approximation of the loglikelihood.

Two approximations are available:: ‘ss’: The sum square of the residuals. This approximation is outlined in Burnham & Anderson 2004 and Chatfield 2003. ‘det’: The determinant of the residual covariance matrix.

The sum squares method is used as default.

Parameters:	A – Array containing the parameters of a fitted model, of shape [channels x channels x maxorder] EF – Array containing the residuals of the fitted model, of shape [channels x samples x trials] method – string indicating which method should be used to approximate the loglikelihood (‘ss’ or ‘det’)
Returns:	tuple of (aLL, AIC, BIC) aLL: Approximation of the loglikelihood (determined by the method defined above) AIC: Estimate of Akaike’s Information Criterion BIC: Estimate of the Bayesian Information Criterion

Tutorial helper functions¶

sails.tutorial_utils.generate_pink_roots(power, order=24)[source]¶

Generates the poles of a system with a pink spectrum.

Parameters:	power – Power of pole order – Order of system to generate
Returns:	Roots of the system (length of order)

sails.tutorial_utils.model_from_roots(roots)[source]¶: Sets up a univariate model based on the roots of a polynomial

sails.tutorial_utils.find_example_path()[source]¶: Returns the path to the example data directory

sails documentation¶

Tutorials¶

Tutorial 1 - A pink noise system¶

Tutorial 2 - A more interesting oscillation¶

Tutorial 3 - Fitting real univariate data¶

Tutorial 4 - Exploring model order¶

Tutorial 5 - MVAR Connectivity Estimation¶

Tutorial 6 - Sliding window univariate model estimation¶

Tutorial 7 - Plotting Helpers¶

Setting up the handler¶

Extracting Region Ordering and Indices¶

Circular Plots¶

Generating a karyotype¶

Generating a connectivity plot based on the karyotype¶

Manually Generating the plot¶

Automatically Generating the plot¶

Modifying link colours¶

Netmat Plots¶

Raw Netmat Plots¶

Gallery¶

Simulation Gallery¶

Baccala and Sameshima 2001¶

Pascual-Marqui et al 2014¶

Module reference¶

Model fitting and Diagnostics¶

Model fitting classes¶

Model fitting helper functions¶

Model diagnostics¶

Model Decomposition¶

Model metric estimatation¶

Plotting¶

Simulation¶

Support routines¶

Abstract Classes¶

Statistical Functions¶

Tutorial helper functions¶

Indices and tables¶