%% 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
%%%% >>> Your code here <<< %%%%
%%
% convert to cosmomvpa struct, using cosmo_meeg_dataset
%%%% >>> Your code here <<< %%%%
% show the dataset
%%%% >>> Your code here <<< %%%%
%%
% set the targets in ds.sa.targets (trial condition)
% Hint: use the first column from ds.sa.trialinfo
%%%% >>> Your code here <<< %%%%
% 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
%%%% >>> Your code here <<< %%%%
% 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
%%%% >>> Your code here <<< %%%%
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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
%%%% >>> Your code here <<< %%%%
% 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
%%%% >>> Your code here <<< %%%%
% 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 varibale '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');