tutorial_circuits_analysis.m

% TUTORIAL_CIRCUITS_ANALYSIS
% Estimates transfer entropy (TE) between voltage measurements obtained
% from two coupled Mackey-Glass electronic circuits.
%
% results = tutorial_circuits_analysis(pset)
%
% where
%
% PSET is a point-set cell array generated by function generate_pset().
%
% RESULTS is a struct containing the results of the analysis.
%
% See also: generate_pset, tutorial_circuits_figures

% Description: Time-varying transfer entropy versus lag 
% Documentation: tutorial_circuits.txt
% Author: German Gomez-Herrero

function results = tutorial_circuits_analysis(pset)

OUT_FILE = 'tutorial_circuits_analysis';
DIM_SINK = 1;        % embedding dimension for the sink of a flow
DIM_SOURCE = 1;      % embedding dimension for the source of a flow
DIM_FUTURE = 1;      % embed. dim. for the future of the sink of a low
STEP = 1;            % embedding delay
K = 20;              % number of nearest neighbors
WINDOW_RADIUS = 5;   % time-window radius for the temporal TE estimates
LAGS = 5:30;         % lag values
FILTER_ORDER = 20;   % order of the post-processing moving average filter
ALPHA = 0.05;         % significance level

% Function handle to the transfer entropy estimator
estimator = @(x) transfer_entropy_t(x(1,:), x(2,:), x(3,:), ...
    WINDOW_RADIUS, 'yLag', LAGS, 'k', K);
estimator_smooth = @(x) filter((1/FILTER_ORDER)*ones(1,FILTER_ORDER),1, ...
    estimator(x),[],2);

% Assess information flow in the direction  2 <- 1
Ssink = delay_embed(pset(2,:), DIM_SINK, STEP);
Ssource = delay_embed(pset(1,:), DIM_SOURCE, STEP);
W = delay_embed_future(Ssink, STEP);
te21 = estimator_smooth([Ssink;Ssource;W]);

% Significance threshold for flow 2 <- 1
te21sig = permutation_test([Ssink; Ssource;W], [1 2 1], ...
    estimator_smooth, ALPHA);

% Assess information flow in the direction  1 <- 2
Ssink = delay_embed(pset(1,:), DIM_SINK, STEP);
Ssource = delay_embed(pset(2,:), DIM_SOURCE, STEP);
W = delay_embed_future(Ssink, STEP);
te12 = estimator_smooth([Ssink;Ssource;W]);

% Significance threshold for flow 1 <- 2
te12sig = permutation_test([Ssink; Ssource;W], [1 2 1], ...
    estimator_smooth, ALPHA);

results = struct('te21', te21, 'te21sig', te21sig, ...
    'te12', te12, 'te12sig', te12sig, ...
    'windowRadius', WINDOW_RADIUS, 'filterOrder', FILTER_ORDER, ...
    'k', K, 'alpha', ALPHA, 'step', STEP, 'dimSink', DIM_SINK, ...
    'dimSource', DIM_SOURCE, 'dimFuture', DIM_FUTURE, 'lags', LAGS);

save([OUT_FILE '.mat'], 'results');