run meeg timelock measures

run_up:

run_meeg_timelock_measures

%% MEEG time-lock searchlight
%
% This example shows MVPA analyses performed on MEEG data.
%
% The input dataset involved a paradigm where a participant saw
% images of six object categories.
%
%
% The code presented here can be adapted for other MEEG analyses, but
% there please note:
% * the current examples do not perform baseline corrections or signal
%   normalizations, which may reduce discriminatory power.
%
% Note: running this code requires FieldTrip.
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #

%% get timelock data in CoSMoMVPA format

% set configuration
config = cosmo_config();
data_path = fullfile(config.tutorial_data_path, 'meg_obj6');

% show dataset information
readme_fn = fullfile(data_path, 'README');
cosmo_type(readme_fn);

% reset citation list
cosmo_check_external('-tic');

% load preprocessed data
data_fn = fullfile(data_path, 'meg_obj6_s00.mat');
data_tl = load(data_fn);

% Show data_tl
cosmo_disp(data_tl);

%%

% convert to cosmomvpa struct, using cosmo_meeg_dataset
ds = cosmo_meeg_dataset(data_tl);

% show the dataset
cosmo_disp(ds);

%%

% set the targets in ds.sa.targets (trial condition)
% Hint: use the first column from ds.sa.trialinfo
ds.sa.targets = ds.sa.trialinfo(:, 1); % 6 categories

% set the chunks in ds.sa.chunks (independent measurements)
% all trials are here considered to be independent, so the chunks
% must all have a different value
nsamples = size(ds.samples, 1);
ds.sa.chunks = (1:nsamples)';

% in addition give a label to each trial
index2label = {'body', 'car', 'face', 'flower', 'insect', 'scene'};
ds.sa.labels = cellfun(@(x)index2label(x), num2cell(ds.sa.targets));

% just to check everything is ok
cosmo_check_dataset(ds);

%% Count number of channels, time points and trials
fprintf('There are %d channels, %d time points and %d trials\n', ...
        numel(unique(ds.fa.chan)), numel(unique(ds.fa.time)), ...
        size(ds.samples, 1));
%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Part I: compute difference between faces and scenes
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% For each time point and sensor; then visualize the results
% for the magnetometer (meg_axial) sensors.

% slice 'ds' twice to get 'ds_face' and 'ds_scene', each with only trials
% from the face and scene categories

ds_face = cosmo_slice(ds, cosmo_match(ds.sa.labels, 'face'));
ds_scene = cosmo_slice(ds, cosmo_match(ds.sa.labels, 'scene'));

% prepare dataset for output
ds_faceVSscene = cosmo_slice(ds_face, 1);
ds_faceVSscene.sa = struct(); % destroy sample attributes

% Compute difference between average of faces versus average of scenes;
% store the result in the samples field of ds_faceVSscene
ds_faceVSscene.samples = mean(ds_face.samples) - mean(ds_scene.samples);

% Convert ds_faceVSscene to a fieldtrip structure and convert
ft_faceVSscene = cosmo_map2meeg(ds_faceVSscene);
%%
% Use FieldTrip to visualize the face versus house contrast
chantype = 'meg_axial';
layout = cosmo_meeg_find_layout(ds_faceVSscene, 'chantype', chantype);

figure();
cfg = struct();
cfg.interactive = 'yes';
cfg.zlim = [-1 1];
cfg.layout = layout;

% show figure with plots for each sensor
ft_multiplotER(cfg, ft_faceVSscene);
%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Part 2: run searchlight over time
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% set MVPA parameters
fprintf('The input has feature dimensions %s\n', ...
        cosmo_strjoin(ds.a.fdim.labels, ', '));

% only select relevant time period and sensors
sensor_posterior_axial = {'MEG1631', 'MEG1641', 'MEG1731', 'MEG1841', ...
                          'MEG1911', 'MEG1921', 'MEG1941', 'MEG2231', ...
                          'MEG2311', 'MEG2321', 'MEG2341', 'MEG2431', ...
                          'MEG2441', 'MEG2511', 'MEG2531'};

% define the mask
msk = cosmo_dim_match(ds, 'time', @(t) t >= -.1 & t <= .4, ...
                      'chan', sensor_posterior_axial);

% first slice the dataset, then use cosmo_dim_prune to avoid using
% non-selected data. Assign the result to ds_sel
%%%% >>> Your code here <<< %%%%
ds_sel_orig = cosmo_slice(ds, msk, 2);
cosmo_disp(ds_sel_orig);

ft_orig = cosmo_map2meeg(ds_sel_orig);
%%
ds_sel = cosmo_dim_prune(ds_sel_orig);
cosmo_disp(ds_sel);
ft_sel = cosmo_map2meeg(ds_sel);

%%

% define the neighborhood for time with a time radius of 2 time points,
% and assign to time_nbrhood.
% Hint: use cosmo_interval_neighborhood,

time_nbrhood = cosmo_interval_neighborhood(ds_sel, 'time', 'radius', 2);
cosmo_disp(time_nbrhood);

%%

%%%% >>> Your code here <<< %%%%

% Define the measure to be cosmo_crossvalidation_measure,
% and assign to measure
measure = @cosmo_crossvalidation_measure;

%%%% >>> Your code here <<< %%%%

% nsamples=numel(ds_sel.samples);
% rp=cosmo_randperm(nsamples);
% ds_sel.sa.targets=ds_sel.sa.targets(rp);

partitions = cosmo_independent_samples_partitioner(ds_sel, ...
                                                   'fold_count', 5, ...
                                                   'test_ratio', 0.2);

cosmo_disp(partitions);
%%
% Define the partitioning scheme using
% cosmo_independent_samples_partitioner, and assign to partitions.
% Use 'fold_count',5 to use 5 folds,
% and use 'test_ratio',.2 to use 20% of the data for testing (and 80% for
% training) in each fold.

%%%% >>> Your code here <<< %%%%

% Use the LDA classifier and the partitions just defined,
measure_args = struct();
measure_args.partitions = partitions;
measure_args.classifier = @cosmo_classify_lda;

ds_sl = cosmo_searchlight(ds_sel, time_nbrhood, measure, measure_args);

plot(ds_sl.a.fdim.values{1}, ds_sl.samples);
xlabel('time');
ylabel('classification accuracy');

%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Part III: channel-time searchlight
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% select limited time period
msk = cosmo_dim_match(ds, 'time', @(t) t >= -.1 & t <= .4);

% first slice the dataset using the mask and assign to 'ds_sel'
ds_sel = cosmo_slice(ds, msk, 2);

% Now use cosmo_dim_prune to avoid attempts to using
% non-selected data in the searchlight
ds_sel = cosmo_dim_prune(ds_sel);

% Set the seachlight parameters
chan_count = 10; % 10 center channels in each searchlight
time_radius = 2; % 2*2+1=5 time bins
chan_type = 'meg_combined_from_planar';

% define the neighborhood for each dimensions
% First, set 'chan_nbrhood' using cosmo_meeg_chan_neighborhood,
% and use the 'chantype' and 'count' parameters to set the channel type
% and the number of sensors in each searchlight.
chan_nbrhood = cosmo_meeg_chan_neighborhood(ds_sel, 'chantype', chan_type, ...
                                            'count', chan_count);

cosmo_disp(chan_nbrhood);
%%

%%%% >>> Your code here <<< %%%%

% Second, set 'time_nbrhood' using cosmo_interval_neighborhood,
% using the 'time' dimension

time_nbrhood = cosmo_interval_neighborhood(ds_sel, 'time', 'radius', 2);

%%%% >>> Your code here <<< %%%%

% cross neighborhoods for chan-time searchlight
% Hint: use cosmo_cross_neighborhood, and use chan_nbrhood and time_nbrhood
% (in that order) in a cell as the second argument
nbrhood = cosmo_cross_neighborhood(ds_sel, {chan_nbrhood, time_nbrhood});

%%%% >>> Your code here <<< %%%%

% print how many neighbors features have on average
nbrhood_nfeatures = cellfun(@numel, nbrhood.neighbors);
fprintf('Features have on average %.1f +/- %.1f neighbors\n', ...
        mean(nbrhood_nfeatures), std(nbrhood_nfeatures));

%%

% set the 'measure' variable to a function handle to the
% split-half correlation measure
measure = @cosmo_correlation_measure;

% Define the partitioning scheme using
% cosmo_independent_samples_partitioner.
% Use 'fold_count',1 to use 1 folds,
% and use 'test_ratio',.5 to use 50% of the data for testing (and 50% for
% training) in the single fold.
partitions = cosmo_independent_samples_partitioner(ds_sel, 'fold_count', 1, ...
                                                   'test_ratio', 0.5);

%%%% >>> Your code here <<< %%%%
measure_args = struct();
measure_args.partitions = partitions;

%% run searchlight
% run the searchlight using the parameters above, and assign the result
% to a variable 'ds_sl'

ds_sl = cosmo_searchlight(ds_sel, nbrhood, measure, measure_args);

%%%% >>> Your code here <<< %%%%

%% visualize timeseries results

% deduce layout from output
layout = cosmo_meeg_find_layout(ds_sl);
fprintf('The output uses layout %s\n', layout.name);

% map ds_sl to a FieldTrip structure. Assign the result to 'sl_ft'
sl_ft = cosmo_map2meeg(ds_sl);

%%%% >>> Your code here <<< %%%%

figure();
cfg = [];
cfg.interactive = 'yes';
cfg.zlim = [-1 1];
cfg.layout       = layout;

% show figure with fisher-transformed correlations for each sensor
ft_multiplotER(cfg, sl_ft);

%% visualize topology results
% show figure with topology for 100 before to 400ms after stimulus onset
% in bins of 50 ms
figure();
cfg.xlim = -0.1:0.05:0.4;
ft_topoplotER(cfg, sl_ft);

%% Show citation information
cosmo_check_external('-cite');