Demonstrations - full listings

demo_fmri_correlation_searchlight

%% Demo: fMRI searchlights with split-half correlations, classifier, and representational similarity analysis
%
% The data used here is available from http://cosmomvpa.org/datadb.zip
%
% This example uses the following dataset:
% - 'ak6' is based on the following work (please cite if you use it):
%    Connolly et al (2012), Representation of biological classes in the
%    human brain. Journal of Neuroscience,
%    doi 10.1523/JNEUROSCI.5547-11.2012
%
%    Six categories (monkey, lemur, mallard, warbler, ladybug, lunamoth)
%    during ten runs in an fMRI study. Using the General Linear Model
%    response were estimated for each category in each run, resulting
%    in 6*10=60 t-values.
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #



%% Set data paths
% The function cosmo_config() returns a struct containing paths to tutorial
% data. (Alternatively the paths can be set manually without using
% cosmo_config.)
config=cosmo_config();

ak6_study_path=fullfile(config.tutorial_data_path,'ak6');

% show readme information
readme_fn=fullfile(ak6_study_path,'README');
cosmo_type(readme_fn);

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

% set result directory
output_path=config.output_data_path;


%% Example: split-half correlation measure (Haxby 2001-style)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% This example uses the 'ak6' dataset

% define data filenames & load data from even and odd runs

data_path=fullfile(ak6_study_path,'s01'); % data from subject s01
mask_fn=fullfile(data_path, 'brain_mask.nii'); % whole brain mask

data_odd_fn=fullfile(data_path,'glm_T_stats_odd.nii');
ds_odd=cosmo_fmri_dataset(data_odd_fn,'mask',mask_fn,...
                            'targets',1:6,'chunks',1);


data_even_fn=fullfile(data_path,'glm_T_stats_even.nii');
ds_even=cosmo_fmri_dataset(data_even_fn,'mask',mask_fn,...
                            'targets',1:6,'chunks',2);

% Combine even and odd runs
ds_odd_even=cosmo_stack({ds_odd, ds_even});

% print dataset
fprintf('Dataset input:\n');
cosmo_disp(ds_odd_even);

% Use cosmo_correlation_measure.
% This measure returns, by default, a split-half correlation measure
% based on the difference of mean correlations for matching and
% non-matching conditions (a la Haxby 2001).
measure=@cosmo_correlation_measure;

% define spherical neighborhood with radius of 3 voxels
radius=3; % voxels
nbrhood=cosmo_spherical_neighborhood(ds_odd_even,'radius',radius);

% Run the searchlight with a 3 voxel radius
corr_results=cosmo_searchlight(ds_odd_even,nbrhood,measure);

% print output
fprintf('Dataset output:\n');
cosmo_disp(corr_results);


% Plot the output
cosmo_plot_slices(corr_results);

% Define output location
output_fn=fullfile(output_path,'corr_searchlight.nii');

% Store results to disc
cosmo_map2fmri(corr_results, output_fn);

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

demo_fmri_distatis

%% Demo: DISTATIS
%
% The data used here is available from http://cosmomvpa.org/datadb.zip
%
% It is based on the following work:
% * Connolly et al (2012), Representation of biological classes in the human
%   brain. Journal of Neuroscience, doi 10.1523/JNEUROSCI.5547-11.2012
%
% Six categories (monkey, lemur, mallard, warbler, ladybug, lunamoth)
% during ten runs in an fMRI study.
%
% This example shows the application of DISTATIS, which tries to find an
% optimal 'compromise' dissimilarity matrix across a set of observations
% (participants)
%
% Reference:
%   - Abdi, H., Valentin, D., O?Toole, A. J., & Edelman, B. (2005).
%     DISTATIS: The analysis of multiple distance matrices. In
%     Proceedings of the IEEE Computer Society: International conference
%     on computer vision and pattern recognition, San Diego, CA, USA
%     (pp. 42?47).
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #

%% Set data paths
% The function cosmo_config() returns a struct containing paths to tutorial
% data. (Alternatively the paths can be set manually without using
% cosmo_config.)
config=cosmo_config();
study_path=fullfile(config.tutorial_data_path,'ak6');
output_path=config.output_data_path;

readme_fn=fullfile(study_path,'README');
cosmo_type(readme_fn);

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

%% Preprocessing for DISTATIS: RSM analysis
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subject_ids={'s01','s02','s03','s04','s05','s06','s07','s08'};
nsubjects=numel(subject_ids);

mask_label='vt_mask';

ds_rsms=cell(nsubjects,1); % allocate space for output
for subject_num=1:nsubjects
    subject_id=subject_ids{subject_num};

    % Code from here is pretty much identical to that above >>>

    % set path for this subject
    data_path=fullfile(study_path,subject_id);

    % Define data locations and load data from even and odd runs
    mask_fn=fullfile(data_path, [mask_label '.nii']); % vt mask

    % Use odd runs only
    data_fn=fullfile(data_path,'glm_T_stats_odd.nii');
    ds=cosmo_fmri_dataset(data_fn,'mask',mask_fn,...
                            'targets',1:6,'chunks',1);

    ds_rsm=cosmo_dissimilarity_matrix_measure(ds);

    % set chunks (one chunk per subject)
    ds_rsm.sa.chunks=subject_num*ones(size(ds_rsm.samples,1),1);
    ds_rsms{subject_num}=ds_rsm;
end

% combine data from all subjects
all_ds=cosmo_stack(ds_rsms);

%% Run DISTATIS
distatis=cosmo_distatis(all_ds);

%% show comprimise distance matrix
[compromise_matrix,dim_labels,values]=cosmo_unflatten(distatis,1);

labels={'monkey', 'lemur', 'mallard', 'warbler', 'ladybug', 'lunamoth'};
n_labels=numel(labels);
figure();
imagesc(compromise_matrix)
title('DSM');
set(gca,'YTick',1:n_labels,'YTickLabel',labels);
set(gca,'XTick',1:n_labels,'XTickLabel',labels);
ylabel(dim_labels{1});
xlabel(dim_labels{2});
colorbar

% skip if stats toolbox is not present
if cosmo_check_external('@stats',false)
    figure();
    hclus = linkage(compromise_matrix);
    dendrogram(hclus,'labels',labels,'orientation','left');
    title('dendrogram');

    figure();
    F = cmdscale(squareform(compromise_matrix));
    text(F(:,1), F(:,2), labels);
    title('2D MDS plot');
    mx = max(abs(F(:)));
    xlim([-mx mx]); ylim([-mx mx]);
end


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

demo_fmri_rois

%% Demo: fMRI Region-Of-Interest (ROI) analyses
%
% The data used here is available from http://cosmomvpa.org/datadb.zip
%
% It is based on the following work:
% * Connolly et al (2012), Representation of biological classes in the human
%   brain. Journal of Neuroscience, doi 10.1523/JNEUROSCI.5547-11.2012
%
% Six categories (monkey, lemur, mallard, warbler, ladybug, lunamoth)
% during ten runs in an fMRI study. Using the General Linear Model response
% were estimated for each category in each run, resulting in 6*10=60
% t-values.
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #

%% Set data paths
% The function cosmo_config() returns a struct containing paths to tutorial
% data. (Alternatively the paths can be set manually without using
% cosmo_config.)
config=cosmo_config();
study_path=fullfile(config.tutorial_data_path,'ak6');
output_path=config.output_data_path;

readme_fn=fullfile(study_path,'README');
cosmo_type(readme_fn);

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

%% Example 1: split-half correlation measure (Haxby 2001-style)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subject_id='s01';
mask_label='ev_mask';
data_path=fullfile(study_path,subject_id); % data from subject s01

% Define data locations and load data from even and odd runs
mask_fn=fullfile(data_path, [mask_label '.nii']); % whole brain

data_odd_fn=fullfile(data_path,'glm_T_stats_odd.nii');
ds_odd=cosmo_fmri_dataset(data_odd_fn,'mask',mask_fn,...
                            'targets',1:6,'chunks',1);


data_even_fn=fullfile(data_path,'glm_T_stats_even.nii');
ds_even=cosmo_fmri_dataset(data_even_fn,'mask',mask_fn,...
                            'targets',1:6,'chunks',2);

% Combine even and odd runs
ds_odd_even=cosmo_stack({ds_odd, ds_even});

% remove constant features
ds_odd_even=cosmo_remove_useless_data(ds_odd_even);

% print dataset
fprintf('Dataset input:\n');
cosmo_disp(ds_odd_even);

% compute correlations
ds_corr=cosmo_correlation_measure(ds_odd_even);

% show result
fprintf(['Average correlation difference between matching and '...
            'non-matching categories in %s for %s is %.3f\n'],...
            mask_label, subject_id, ds_corr.samples);


%% Example 2: split-half correlation measure with group analysis
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

subject_ids={'s01','s02','s03','s04','s05','s06','s07','s08'};
nsubjects=numel(subject_ids);

mask_label='vt_mask';

ds_corrs=cell(nsubjects,1); % allocate space for output
for subject_num=1:nsubjects
    subject_id=subject_ids{subject_num};

    % Code from here is pretty much identical to that above >>>

    % set path for this subject
    data_path=fullfile(study_path,subject_id);

    % Define data locations and load data from even and odd runs
    mask_fn=fullfile(data_path, [mask_label '.nii']); % whole brain

    data_odd_fn=fullfile(data_path,'glm_T_stats_odd.nii');
    ds_odd=cosmo_fmri_dataset(data_odd_fn,'mask',mask_fn,...
                                'targets',1:6,'chunks',1);


    data_even_fn=fullfile(data_path,'glm_T_stats_even.nii');
    ds_even=cosmo_fmri_dataset(data_even_fn,'mask',mask_fn,...
                                'targets',1:6,'chunks',2);

    % Combine even and odd runs
    ds_odd_even=cosmo_stack({ds_odd, ds_even});

    % remove constant features
    ds_odd_even=cosmo_remove_useless_data(ds_odd_even);

    ds_corr=cosmo_correlation_measure(ds_odd_even);

    % <<< identical up to here

    % set targets and chunks for the output, so that cosmo_stat can be used
    % below
    ds_corr.sa.targets=1;
    ds_corr.sa.chunks=subject_num;

    ds_corrs{subject_num}=ds_corr;
end

% combine the data from all subjects
ds_all=cosmo_stack(ds_corrs);

%%
samples=ds_all.samples; % get the correlations for all subjects

% run one-sample t-test again zero

% Using cosmo_stats
ds_t=cosmo_stat(ds_all,'t');     % t-test against zero
ds_p=cosmo_stat(ds_all,'t','p'); % convert to p-value

fprintf(['correlation difference in %s at group level: '...
           '%.3f +/- %.3f, %s=%.3f, p=%.5f (using cosmo_stat)\n'],...
            mask_label,mean(samples),std(samples),...
            ds_t.sa.stats{1},ds_t.samples,ds_p.samples);

% Using matlab's stat toolbox (if present)
if cosmo_check_external('@stats',false)
    [h,p,ci,stats]=ttest(samples);
    fprintf(['Correlation difference in %s at group level: '...
            '%.3f +/- %.3f, t_%d=%.3f, p=%.5f (using matlab stats '...
            'toolbox)\n'],...
            mask_label,mean(samples),std(samples),stats.df,stats.tstat,p);
else
    fprintf('Matlab stats toolbox not found\n');
end


%% Example 3: comparison of four classifiers in two regions of interest
config=cosmo_config();
data_path=fullfile(config.tutorial_data_path,'ak6','s01');
data_fn=fullfile(data_path,'glm_T_stats_perrun.nii');

% Define classifiers and mask labels
classifiers={@cosmo_classify_nn,...
             @cosmo_classify_naive_bayes,...
             @cosmo_classify_lda};

% Use svm classifiers, if present
svm_name2func=struct();
svm_name2func.matlabsvm=@cosmo_classify_matlabsvm;
svm_name2func.libsvm=@cosmo_classify_libsvm;
svm_name2func.svm=@cosmo_classify_svm;
svm_names=fieldnames(svm_name2func);
for k=1:numel(svm_names)
    svm_name=svm_names{k};
    if cosmo_check_external(svm_name,false)
        classifiers{end+1}=svm_name2func.(svm_name);
    else
        warning('Classifier %s skipped because not available', svm_name);
    end
end

mask_labels={'vt_mask','ev_mask'};

%
nclassifiers=numel(classifiers);
nmasks=numel(mask_labels);

labels={'monkey', 'lemur', 'mallard', 'warbler', 'ladybug', 'lunamoth'};
nlabels=numel(labels);

% little helper function to replace underscores by spaces
underscore2space=@(x) strrep(x,'_',' ');

for j=1:nmasks
    mask_label=mask_labels{j};
    mask_fn=fullfile(data_path,[mask_label '.nii']);
    ds=cosmo_fmri_dataset(data_fn,'mask',mask_fn,...
                        'targets',repmat(1:6,1,10),...
                        'chunks',floor(((1:60)-1)/6)+1);

    % remove constant features
    ds=cosmo_remove_useless_data(ds);

    % print dataset
    fprintf('Dataset input:\n');
    cosmo_disp(ds);

    % Define partitions
    partitions=cosmo_nfold_partitioner(ds);

    % print dataset
    fprintf('Partitions:\n');
    cosmo_disp(partitions);

    % show result for each classifier
    for k=1:nclassifiers
        classifier=classifiers{k};

        % get predictions for each fold
        [pred,accuracy]=cosmo_crossvalidate(ds, classifier, partitions);

        % get confusion matrix for each fold
        confusion_matrix_folds=cosmo_confusion_matrix(ds.sa.targets,pred);

        % sum confusion for each ground-truth target and prediction,
        % resulting in an nclasses x nclasses matrix
        confusion_matrix=sum(confusion_matrix_folds,3);
        figure
        imagesc(confusion_matrix,[0 10])
        cfy_label=underscore2space(func2str(classifier));
        title_=sprintf('%s using %s: accuracy=%.3f', ...
                        underscore2space(mask_label), cfy_label, accuracy);
        title(title_)
        set(gca,'XTick',1:nlabels,'XTickLabel',labels);
        set(gca,'YTick',1:nlabels,'YTickLabel',labels);
        ylabel('target');
        xlabel('predicted');
    end
end

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

demo_fmri_searchlight_lda

%% Demo: fMRI searchlights LDA classifier
%
% The data used here is available from http://cosmomvpa.org/datadb.zip
%
% This example uses the following dataset:
% + 'digit'
%    A participant made finger pressed with the index and middle finger of
%    the right hand during 4 runs in an fMRI study. Each run was divided in
%    4 blocks with presses of each finger and analyzed with the GLM,
%    resulting in 2*4*4=32 t-values
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #


%% Set data paths
% The function cosmo_config() returns a struct containing paths to tutorial
% data. (Alternatively the paths can be set manually without using
% cosmo_config.)
config=cosmo_config();

digit_study_path=fullfile(config.tutorial_data_path,'digit');
readme_fn=fullfile(digit_study_path,'README');
cosmo_type(readme_fn);

output_path=config.output_data_path;

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

%% LDA classifier searchlight analysis
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This analysis identified brain regions where the categories can be
% distinguished using an odd-even partitioning scheme and a Linear
% Discriminant Analysis (LDA) classifier.

data_path=digit_study_path;
data_fn=fullfile(data_path,'glm_T_stats_perblock+orig.HEAD');
mask_fn=fullfile(data_path,'brain_mask+orig.HEAD');

targets=repmat(1:2,1,16)';    % class labels: 1 2 1 2 1 2 1 2 1 2 ... 1 2
chunks=floor(((1:32)-1)/8)+1; % run labels:   1 1 1 1 1 1 1 1 2 2 ... 4 4

ds_per_run = cosmo_fmri_dataset(data_fn, 'mask', mask_fn,...
                                'targets',targets,'chunks',chunks);

% print dataset
fprintf('Dataset input:\n');
cosmo_disp(ds_per_run);


% Use the cosmo_cross_validation_measure and set its parameters
% (classifier and partitions) in a measure_args struct.
measure = @cosmo_crossvalidation_measure;
measure_args = struct();

% Define which classifier to use, using a function handle.
% Alternatives are @cosmo_classify_{svm,matlabsvm,libsvm,nn,naive_bayes}
measure_args.classifier = @cosmo_classify_lda;

% Set partition scheme. odd_even is fast; for publication-quality analysis
% nfold_partitioner is recommended.
% Alternatives are:
% - cosmo_nfold_partitioner    (take-one-chunk-out crossvalidation)
% - cosmo_nchoosek_partitioner (take-K-chunks-out  "             ").
measure_args.partitions = cosmo_oddeven_partitioner(ds_per_run);

% print measure and arguments
fprintf('Searchlight measure:\n');
cosmo_disp(measure);
fprintf('Searchlight measure arguments:\n');
cosmo_disp(measure_args);

% Define a neighborhood with approximately 100 voxels in each searchlight.
nvoxels_per_searchlight=100;
nbrhood=cosmo_spherical_neighborhood(ds_per_run,...
                        'count',nvoxels_per_searchlight);


% Run the searchlight
lda_results = cosmo_searchlight(ds_per_run,nbrhood,measure,measure_args);

% print output dataset
fprintf('Dataset output:\n');
cosmo_disp(lda_results);

% Plot the output
cosmo_plot_slices(lda_results);

% Define output location
output_fn=fullfile(output_path,'lda_searchlight+orig');

% Store results to disc
cosmo_map2fmri(lda_results, output_fn);

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

demo_fmri_searchlight_naive_bayes

%% Demo: fMRI searchlight with naive bayes classifier
%
% The data used here is available from http://cosmomvpa.org/datadb.zip
%
% This example uses the following dataset:
% + 'digit'
%    A participant made finger pressed with the index and middle finger of
%    the right hand during 4 runs in an fMRI study. Each run was divided in
%    4 blocks with presses of each finger and analyzed with the GLM,
%    resulting in 2*4*4=32 t-values
%
% This example uses the cosmo_naive_bayes_classifier_searchlight, which is
% a fast alternative to using the regular searchlight with a
% crossvalidation measure and a classifier
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #



%% Set data paths
% The function cosmo_config() returns a struct containing paths to tutorial
% data. (Alternatively the paths can be set manually without using
% cosmo_config.)
config=cosmo_config();

digit_study_path=fullfile(config.tutorial_data_path,'digit');
readme_fn=fullfile(digit_study_path,'README');
cosmo_type(readme_fn);

output_path=config.output_data_path;

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

%% LDA classifier searchlight analysis
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This analysis identified brain regions where the categories can be
% distinguished using an odd-even partitioning scheme and a Naive Bayes
% classifier.

data_path=digit_study_path;
data_fn=fullfile(data_path,'glm_T_stats_perblock+orig.HEAD');
mask_fn=fullfile(data_path,'brain_mask+orig.HEAD');

targets=repmat(1:2,1,16)';    % class labels: 1 2 1 2 1 2 1 2 1 2 ... 1 2
chunks=floor(((1:32)-1)/8)+1; % run labels:   1 1 1 1 1 1 1 1 2 2 ... 4 4

ds_per_run = cosmo_fmri_dataset(data_fn, 'mask', mask_fn,...
                                'targets',targets,'chunks',chunks);

% print dataset
fprintf('Dataset input:\n');
cosmo_disp(ds_per_run);


% set parameters for naive bayes searchlight (partitions) in a
% measure_args struct.
measure_args = struct();

% Set partition scheme. odd_even is fast; for publication-quality analysis
% nfold_partitioner is recommended.
% Alternatives are:
% - cosmo_nfold_partitioner    (take-one-chunk-out crossvalidation)
% - cosmo_nchoosek_partitioner (take-K-chunks-out  "             ").
measure_args.partitions = cosmo_oddeven_partitioner(ds_per_run);

% print measure and arguments
fprintf('Searchlight measure arguments:\n');
cosmo_disp(measure_args);

% Define a neighborhood with approximately 100 voxels in each searchlight.
nvoxels_per_searchlight=100;
nbrhood=cosmo_spherical_neighborhood(ds_per_run,...
                        'count',nvoxels_per_searchlight);


%% Run the searchlight
%
nb_results = cosmo_naive_bayes_classifier_searchlight(ds_per_run,...
                                                nbrhood,measure_args);
%% Show results
% print output dataset
fprintf('Dataset output:\n');
cosmo_disp(nb_results);

% Plot the output
cosmo_plot_slices(nb_results);

% Define output location
output_fn=fullfile(output_path,'naive_bayes_searchlight+orig');

% Store results to disc
cosmo_map2fmri(nb_results, output_fn);

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

demo_fmri_searchlight_rsm

%% Demo: fMRI searchlights with representational similarity analysis
%
% The data used here is available from http://cosmomvpa.org/datadb.zip
%
% This example uses the following dataset:
% - 'ak6' is based on the following work (please cite if you use it):
%    Connolly et al (2012), Representation of biological classes in the
%    human brain. Journal of Neuroscience,
%    doi 10.1523/JNEUROSCI.5547-11.2012
%
%    Six categories (monkey, lemur, mallard, warbler, ladybug, lunamoth)
%    during ten runs in an fMRI study. Using the General Linear Model
%    response were estimated for each category in each run, resulting
%    in 6*10=60 t-values.
%
% The example shows a searchlight analysis matching local neural similarity
% patterns to three different target similarity matrices
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #


%% Set data paths
% The function cosmo_config() returns a struct containing paths to tutorial
% data. (Alternatively the paths can be set manually without using
% cosmo_config.)
config=cosmo_config();

ak6_study_path=fullfile(config.tutorial_data_path,'ak6');

% show readme information
readme_fn=fullfile(ak6_study_path,'README');
cosmo_type(readme_fn);

output_path=config.output_data_path;

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


%% Load data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% This example uses the 'ak6' dataset
% In this example only one sample (response estimate) per condition (class)
% per feature (voxel) is used. Here this is done using t-stats from odd
% runs. One could also use output from a GLM based on an entire
% scanning session experiment.
%

% define data filenames & load data from even and odd runs

data_path=fullfile(ak6_study_path,'s01'); % data from subject s01
mask_fn=fullfile(data_path, 'brain_mask.nii'); % whole brain mask

data_fn=fullfile(data_path,'glm_T_stats_odd.nii');
ds=cosmo_fmri_dataset(data_fn,'mask',mask_fn,...
                            'targets',1:6,'chunks',1);

%% Set animal species & class
ds.sa.labels={'monkey','lemur','mallard','warbler','ladybug','lunamoth'}';
ds.sa.animal_class=[1 1 2 2 3 3]';

% simple sanity check to ensure all attributes are set properly
cosmo_check_dataset(ds);

% print dataset
fprintf('Dataset input:\n');
cosmo_disp(ds);

%% Define feature neighorhoods
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% For the searchlight, define neighborhood for each feature (voxel).
nvoxels_per_searchlight=100;

% The neighborhood defined here is used three times (one for each target
% similarity matrix), so it is not recomputed for every searchlight call.
fprintf('Defining neighborhood for each feature\n');
nbrhood=cosmo_spherical_neighborhood(ds,'count',nvoxels_per_searchlight);

% print neighborhood
fprintf('Searchlight neighborhood definition:\n');
cosmo_disp(nbrhood);


%% Simple RSM searchlight
nsamples=size(ds.samples,1);
target_dsm=zeros(nsamples);

% define 'simple' target structure where primates (monkey, lemur),
% birds (mallard, warbler) and insects (ladybug, lunamoth) are the same
% (distance=0), all other pairs are equally dissimilar (distance=1).
% Given the ds.sa.targets assignment, pairs (1,2), (3,4) and (5,6) have
% distance 0, all others distance 1.
animal_class=ds.sa.animal_class;
for row=1:nsamples
    for col=1:nsamples
        same_animal_class=animal_class(row)==animal_class(col);

        if same_animal_class
            target_dsm(row,col)=0;
        else
            target_dsm(row,col)=1;
        end
    end
end

fprintf('Using the following target dsm\n');
disp(target_dsm);
imagesc(target_dsm)
set(gca,'XTick',1:nsamples,'XTickLabel',ds.sa.labels,...
        'YTick',1:nsamples,'YTickLabel',ds.sa.labels)

% set measure
measure=@cosmo_target_dsm_corr_measure;
measure_args=struct();
measure_args.target_dsm=target_dsm;

% print measure and arguments
fprintf('Searchlight measure:\n');
cosmo_disp(measure);
fprintf('Searchlight measure arguments:\n');
cosmo_disp(measure_args);

% run searchlight
ds_rsm_binary=cosmo_searchlight(ds,nbrhood,measure,measure_args);

% Note: when these results are used for group analysis across multiple
% participants, it may be good to Fisher-transform the correlation values,
% so that they are more normally distributed. This can be done by:
%
% ds_rsm_binary.samples=atanh(ds_rsm_binary.samples);


% show results
cosmo_plot_slices(ds_rsm_binary);

% store results
output_fn=fullfile(output_path,'rsm_binary.nii');
cosmo_map2fmri(ds_rsm_binary,output_fn);


%% Using another RSM

% This example is very similar to the previous example.
% - This example uses a different target representational similarity
%   matrix. The code below allows for identifying regions that show a
%   linear dissimilarity across animal class, with primates of distance 1
%   from birdsand distance 2 from insects, and insects distance 1 from
%   birds. Graphically:
%
%            +------------------+------------------+
%        primates             birds             insects
%     {monkey,lemur}     {mallard,warbler}   {ladybug,lunamoth}
%
% - It uses a Spearman rather than Pearson correlation measure to match
%   the neural similarity to the target similarity measure

animal_class=ds.sa.animal_class;

% compute absolute difference between each pair of samples
target_dsm=abs(bsxfun(@minus,animal_class,animal_class'));

fprintf('Using the following target dsm\n');
disp(target_dsm);
imagesc(target_dsm)
set(gca,'XTick',1:nsamples,'XTickLabel',ds.sa.labels,...
        'YTick',1:nsamples,'YTickLabel',ds.sa.labels)

% set measure
measure=@cosmo_target_dsm_corr_measure;
measure_args=struct();
measure_args.target_dsm=target_dsm;

% 'Spearman' requires  matlab with stats toolbox; if not present use
% 'Pearson'
if cosmo_check_external('@stats',false)
    measure_args.type='Spearman';
else
    measure_args.type='Pearson';
    fprintf('Matlab stats toolbox not present, using %s correlation\n',...
                measure_args.type)
end

% run searchlight
ds_rsm_linear=cosmo_searchlight(ds,nbrhood,measure,measure_args);

% Note: when these results are used for group analysis across multiple
% participants, it may be good to Fisher-transform the correlation values,
% so that they are more normally distributed. This can be done by:
%
% ds_rsm_linear.samples=atanh(ds_rsm_linear.samples);

% show results
cosmo_plot_slices(ds_rsm_linear);

% store results
output_fn=fullfile(output_path,'rsm_linear.nii');
cosmo_map2fmri(ds_rsm_linear,output_fn);


%% Using a behavioural RSM
% This example is very similar to the one above, but now the target
% similarity structure is based on behavioural similarity ratings

% load behavioural similarity matrix from disc
behav_model_fn=fullfile(ak6_study_path,'models','behav_sim.mat');
behav_model=importdata(behav_model_fn);

target_dsm=behav_model;

fprintf('Using the following target dsm\n');
disp(target_dsm);
imagesc(target_dsm)
set(gca,'XTick',1:nsamples,'XTickLabel',ds.sa.labels,...
        'YTick',1:nsamples,'YTickLabel',ds.sa.labels)

% set measure
measure=@cosmo_target_dsm_corr_measure;
measure_args=struct();
measure_args.target_dsm=target_dsm;

% run searchlight
ds_rsm_behav=cosmo_searchlight(ds,nbrhood,measure,measure_args);

% Note: when these results are used for group analysis across multiple
% participants, it may be good to Fisher-transform the correlation values,
% so that they are more normally distributed. This can be done by:
%
% ds_rsm_behav.samples=atanh(ds_rsm_behav.samples);

% show results
cosmo_plot_slices(ds_rsm_behav);

% store results
output_fn=fullfile(output_path,'rsm_behav.nii');
cosmo_map2fmri(ds_rsm_behav,output_fn);

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

demo_meeg_dataset_operations

%% MEEG dataset operations
%
% This example shows MVPA analyses performed on MEEG data.
%
% The input dataset involved a paradigm with electrical median nerve
% stimulation for durations of 2s at 20Hz.
%
%
% The code presented here can be adapted for other MEEG analyses, but
% there are a few potential caveats:
% * assignment of targets (labels of conditions) is based here on
%   stimulation periods versus pre-stimulation periods. In typical
%   analyses the targets should be based on different trial conditions, for
%   example as set a FieldTrip .trialinfo field.
% * assignment of chunks (parts of the data that are assumed to be
%   independent) is done randomly. In typical analyses they can be based on
%   different runs, with one chunk value per run.
% * the time window used for analyses is rather small. This means that in
%   particular for time-freq analysis a lot of data is missing, especially
%   for early and late timepoints in the lower frequency bands. For typical
%   analyses it may be preferred to use a wider time window.
% * 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_20hz');

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

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

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

% convert to cosmomvpa struct
ds_tl=cosmo_meeg_dataset(data_tl);

% set the target (trial condition)
ds_tl.sa.targets=ds_tl.sa.trialinfo(:,1); % 1=pre, 2=post

% in addition give a label to each trial
index2label={'pre','post'}; % 1=pre, 2=peri/post
ds_tl.sa.labels=cellfun(@(x)index2label(x),num2cell(ds_tl.sa.targets));

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

%% Print some information about dataset

% print the size. Note that the number of features is the number of time
% points times the number of channels.
[ns,nf]=size(ds_tl.samples);
fprintf('There are %d samples and %d features\n', ns, nf);

fdim_labels=ds_tl.a.fdim.labels;
fdim_values=ds_tl.a.fdim.values;
ndim=numel(fdim_labels);
for dim=1:ndim
    dim_label=fdim_labels{dim};
    nvalues=numel(fdim_values{dim});
    fprintf('Dimension %d (%s) has %d values\n', dim, dim_label, nvalues);
end

%% Simple data manipulations
% select five trials from post-stim period
post_mask=cosmo_match(ds_tl.sa.labels,'post');
ds_tl_post=cosmo_slice(ds_tl,post_mask);
select_trials=2:2:18;

% the last argument indicates to slice along the first dimension.
% The value 1 is the default and can be omitted
ds_trial=cosmo_slice(ds_tl_post, select_trials, 1);

% select a single channel (multiple channels are also possible but
% makes plotting the samples directly a bit more tricky)
chan={'MEG1042'};
chan_msk=cosmo_dim_match(ds_trial,'chan', chan);

% slice again, this time along the feature dimension (dimension 2),
% to get a new dataset with just one channel
ds_trial_chan=cosmo_slice(ds_trial, chan_msk, 2);

% plot five individual trials
figure();
subplot(2,1,1);
plot(data_tl.time, ds_trial_chan.samples');
title(sprintf('%d trials of sensor %s', numel(select_trials), chan{1}));

% use original full data to select data in just one channel
chan_msk=cosmo_dim_match(ds_tl,'chan', chan);
ds_chan=cosmo_slice(ds_tl, chan_msk, 2);

ds_mean=cosmo_fx(ds_chan,@(x)mean(x,1),'targets');
% plot signal time course averaged over all samples
% note: the 'pre' data has been shifted in time to match the 'post' data
subplot(2,1,2);
[time_dim, time_index]=cosmo_dim_find(ds_mean,'time');
assert(time_dim==2); % it is a feature dimension
plot(ds_mean.a.fdim.values{time_index},ds_mean.samples);
title(sprintf('average of sensor %s', chan{1}));
legend(ds_mean.sa.labels);

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

demo_meeg_timefreq_searchlight

%% MEEG time-frequency searchlight
%
% This example shows MVPA analyses performed on MEEG data, using a
% searchlight across the time, frequency and channel dimensions
%
% The input dataset involved a paradigm with electrical median nerve
% stimulation for durations of 2s at 20Hz.
%
% The code presented here can be adapted for other MEEG analyses, but
% there are a few potential caveats:
% * assignment of targets (labels of conditions) is based here on
%   stimulation periods versus pre-stimulation periods. In typical
%   analyses the targets should be based on different trial conditions, for
%   example as set a FieldTrip .trialinfo field.
% * assignment of chunks (parts of the data that are assumed to be
%   independent) is based on a trial-by-trial basis. For cross-validation,
%   the number of chunks is reduced to two to speed up the analysis.
% * the time window used for analyses is rather small. This means that in
%   particular for time-freq analysis a lot of data is missing, especially
%   for early and late timepoints in the lower frequency bands. For typical
%   analyses it may be preferred to use a wider time window.
% * 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_20hz');

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

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

% load data
data_fn=fullfile(data_path,'subj102_B01_20Hz_timefreq.mat');
data_tf=load(data_fn);

% convert to cosmomvpa struct
ds_tf=cosmo_meeg_dataset(data_tf);

% set the target (trial condition)
ds_tf.sa.targets=ds_tf.sa.trialinfo(:,1); % 1=pre, 2=post

% set the chunks (independent measurements)
% in this dataset, the first half of the samples (in order)
% are the post-trials;
% the second half the pre-trials
ds_tf.sa.chunks=[(1:145) (1:145)]';


% in addition give a label to each trial
index2label={'pre','post'}; % 1=pre, 2=peri/post
ds_tf.sa.labels=cellfun(@(x)index2label(x),num2cell(ds_tf.sa.targets));

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

% get rid of features with at least one NaN value across samples
fa_nan_mask=sum(isnan(ds_tf.samples),1)>0;
fprintf('%d / %d features have NaN\n', ...
            sum(fa_nan_mask), numel(fa_nan_mask));
ds_tf=cosmo_slice(ds_tf, ~fa_nan_mask, 2);


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


% set chunks
% again for speed just two chunks
% (targets were already set above)
nchunks=2;
ds_tf.sa.chunks=cosmo_chunkize(ds_tf, nchunks);

% define neighborhood parameters for each dimension

% channel neighborhood uses meg_combined_from_planar, which means that the
% input are planar channels but the output has combined-planar channels.
% to use the magnetometers, use 'meg_axial'
chan_type='meg_combined_from_planar';
chan_count=10;        % use 10 channel locations (relative to the combined
                      % planar channels)
                      % as we use meg_combined_from_planar there are
                      % 20 channels in each searchlight because
                      % gradiometers are paired
time_radius=2; % 2*2+1=5 time bines
freq_radius=4; % 4*2+1=9 freq bins


% define the neighborhood for each dimensions
chan_nbrhood=cosmo_meeg_chan_neighborhood(ds_tf, 'count', chan_count, ...
                                                'chantype', chan_type);
freq_nbrhood=cosmo_interval_neighborhood(ds_tf,'freq',...
                                            'radius',freq_radius);
time_nbrhood=cosmo_interval_neighborhood(ds_tf,'time',...
                                            'radius',time_radius);

% cross neighborhoods for chan-time-freq searchlight
nbrhood=cosmo_cross_neighborhood(ds_tf,{chan_nbrhood,...
                                        freq_nbrhood,...
                                        time_nbrhood});

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

% only keep features with at least 10 neighbors
% (some have zero neighbors - in particular, those with low frequencies
% early or late in time)
center_ids=find(nbrhood_nfeatures>10);

% for illustration purposes use the split-half measure because it is
% relatively fast - but clasifiers can also be used
measure=@cosmo_correlation_measure;

% split-half, as there are just two chunks
% (when using a classifier, do not use 'half' but the number of chunks to
% leave out for testing, e.g. 1).
measure_args=struct();
measure_args.partitions=cosmo_nchoosek_partitioner(ds_tf,'half');


%% run searchlight
sl_tf_ds=cosmo_searchlight(ds_tf,nbrhood,measure,measure_args,...
                                      'center_ids',center_ids);
%% visualize results

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

% map to FT struct for visualization
sl_tf_ft=cosmo_map2meeg(sl_tf_ds);

% show figure
figure()
cfg = [];
if cosmo_wtf('is_octave')
    % GNU Octave does not show data when labels are shown
    cfg.interactive='no';
    cfg.showlabels='no';
else
    % Matlab supports interactive viewing and labels
    cfg.interactive = 'yes';
    cfg.showlabels = 'yes';
end
cfg.zlim='maxabs';
cfg.layout       = layout;
ft_multiplotTFR(cfg, sl_tf_ft);

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

demo_meeg_timelock_searchlight

%% MEEG time-lock searchlight
%
% This example shows MVPA analyses performed on MEEG data.
%
% The input dataset involved a paradigm with electrical median nerve
% stimulation for durations of 2s at 20Hz.
%
% Using a time-channel neighborhood, a searchlight map is computed
% indicating in time and space (channel) the pre-stimulation period can be
% distinguished from the peri/post stimulation period.
%
% The code presented here can be adapted for other MEEG analyses, but
% there are a few potential caveats:
% * assignment of targets (labels of conditions) is based here on
%   stimulation periods versus pre-stimulation periods. In typical
%   analyses the targets should be based on different trial conditions, for
%   example as set a FieldTrip .trialinfo field.
% * assignment of chunks (parts of the data that are assumed to be
%   independent) is based on a trial-by-trial basis. For cross-validation,
%   the number of chunks is reduced to two to speed up the analysis.
% * 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_20hz');

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

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

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

% convert to cosmomvpa struct
ds_tl=cosmo_meeg_dataset(data_tl);

% set the target (trial condition)
ds_tl.sa.targets=ds_tl.sa.trialinfo(:,1); % 1=pre, 2=post

% set the chunks (independent measurements)
% in this dataset, the first half of the samples (in order)
% are the post-trials;
% the second half the pre-trials
ds_tl.sa.chunks=[(1:145) (1:145)]';


% in addition give a label to each trial
index2label={'pre','post'}; % 1=pre, 2=peri/post
ds_tl.sa.labels=cellfun(@(x)index2label(x),num2cell(ds_tl.sa.targets));

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

% get rid of features with at least one NaN value across samples
fa_nan_mask=sum(isnan(ds_tl.samples),1)>0;
fprintf('%d / %d features have NaN\n', ...
            sum(fa_nan_mask), numel(fa_nan_mask));
ds_tl=cosmo_slice(ds_tl, ~fa_nan_mask, 2);


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


% set chunks
% again for speed just two chunks
% (targets were already set above)
nchunks=2;
ds_tl.sa.chunks=cosmo_chunkize(ds_tl, nchunks);

% define neighborhood parameters for each dimension

% channel neighborhood uses meg_combined_from_planar, which means that the
% input are planar channels but the output has combined-planar channels.
% to use the magnetometers, use 'meg_axial'
chan_type='meg_combined_from_planar';
chan_count=10;        % use 10 channel locations (relative to the combined
                      % planar channels)
                      % as we use meg_combined_from_planar there are
                      % 20 channels in each searchlight because
                      % gradiometers are paired
time_radius=2; % 2*2+1=5 time bines

% define the neighborhood for each dimensions
chan_nbrhood=cosmo_meeg_chan_neighborhood(ds_tl, 'count', chan_count, ...
                                                'chantype', chan_type);
time_nbrhood=cosmo_interval_neighborhood(ds_tl,'time',...
                                            'radius',time_radius);

% cross neighborhoods for chan-time searchlight
nbrhood=cosmo_cross_neighborhood(ds_tl,{chan_nbrhood,...
                                        time_nbrhood});

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

% only keep features with at least 10 neighbors
center_ids=find(nbrhood_nfeatures>10);

% for illustration purposes use the split-half measure because it is
% relatively fast - but clasifiers can also be used
measure=@cosmo_correlation_measure;

% split-half, as there are just two chunks
% (when using a classifier, do not use 'half' but the number of chunks to
% leave out for testing, e.g. 1).
measure_args=struct();
measure_args.partitions=cosmo_nchoosek_partitioner(ds_tl,'half');


%% run searchlight
sl_tl_ds=cosmo_searchlight(ds_tl,nbrhood,measure,measure_args,...
                                      'center_ids',center_ids);
%% visualize timeseries results

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

% map to FT struct for visualization
sl_tl_ft=cosmo_map2meeg(sl_tl_ds);

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

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

%% visualize topology results
% show figure with topology for 0 to 600ms after stimulus onset in bins of
% 100 ms
figure();
cfg.xlim=-0.1:0.1:0.5;
ft_topoplotER(cfg, sl_tl_ft);

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

demo_meeg_timeseries_classification

%% MEEG timeseries classification
%
% This example shows MVPA analyses performed on MEEG data.
%
% The input dataset involved a paradigm with electrical median nerve
% stimulation for durations of 2s at 20Hz.
%
% The code presented here can be adapted for other MEEG analyses, but
% there are a few potential caveats:
%
% * assignment of targets (labels of conditions) is based here on
%   stimulation periods versus pre-stimulation periods. In typical
%   analyses the targets should be based on different trial conditions, for
%   example as set a FieldTrip .trialinfo field.
% * assignment of chunks (parts of the data that are assumed to be
%   independent) is based on a trial-by-trial basis. For cross-validation,
%   the number of chunks is reduced to four to speed up the analysis.
% * the time window used for analyses is rather small. This means that in
%   particular for time-freq analysis a lot of data is missing, especially
%   for early and late timepoints in the lower frequency bands. For typical
%   analyses it may be preferred to use a wider time window.
% * 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_20hz');

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

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

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

% convert to cosmomvpa struct
ds_tl=cosmo_meeg_dataset(data_tl);

% set the target (trial condition)
ds_tl.sa.targets=ds_tl.sa.trialinfo(:,1); % 1=pre, 2=post

% set the chunks (independent measurements)
% in this dataset, the first half of the samples (in order)
% are the post-trials;
% the second half the pre-trials
ds_tl.sa.chunks=[(1:145) (1:145)]';


% in addition give a label to each trial
index2label={'pre','post'}; % 1=pre, 2=peri/post
ds_tl.sa.labels=cellfun(@(x)index2label(x),num2cell(ds_tl.sa.targets));

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

%% Prepare MVPA
% reset chunks: use four chunks
nchunks=4;
ds_tl.sa.chunks=cosmo_chunkize(ds_tl,nchunks);

% do a take-one-fold out cross validation.
% except when using a splithalf correlation measure it is important that
% the partitions are *balanced*, i.e. each target (or class) is presented
% equally often in each chunk
partitions=cosmo_nchoosek_partitioner(ds_tl,1);
partitions=cosmo_balance_partitions(partitions, ds_tl);

npartitions=numel(partitions);
fprintf('There are %d partitions\n', numel(partitions.train_indices));
fprintf('# train samples:%s\n', sprintf(' %d', cellfun(@numel, ...
                                        partitions.train_indices)));
fprintf('# test samples:%s\n', sprintf(' %d', cellfun(@numel, ...
                                        partitions.test_indices)));




%% Run time-series searchlight on magneto- and gradio-meters seperately

% try two different classification approaches:
% 1) without averaging the samples in the train set
% 2) by averaging 5 samples at the time in the train set, and re-using
%    every sample 3 times.
% (Note: As of July 2015, there is no clear indication in the literature
%  which approach is 'better'. These two approaches are used here to
%  illustrate how they can be used with a searchlight).
average_train_args_cell={{},...                                  % {  1
                         {'average_train_count',5,...            %  { 2
                                'average_train_resamplings',3}}; %  { 2
n_average_train_args=numel(average_train_args_cell);

% compute and plot accuracies for magnetometers and gradiometers separately
chantypes={'meg_axial','meg_planar'};

% in the time searchlight analysis, select the time-point itself, the two
% timepoints after it, and the two timepoints before it
time_radius=2;
nchantypes=numel(chantypes);

ds_chantypes=cosmo_meeg_chantype(ds_tl);
plot_counter=0;
for j=1:n_average_train_args
    % define the measure and its argument.
    % here a simple naive baysian classifier is used.
    % Alternative are @cosmo_classify_{svm,nn,lda}.
    measure=@cosmo_crossvalidation_measure;
    measure_args=struct();
    measure_args.classifier=@cosmo_classify_naive_bayes;
    measure_args.partitions=partitions;

    % add the options to average samples to the measure arguments.
    % (if no averaging is desired, this step can be left out.)
    average_train_args=average_train_args_cell{j};
    measure_args=cosmo_structjoin(measure_args, average_train_args);


    for k=1:nchantypes
        parent_type=chantypes{k};

        % find feature indices of channels matching the parent_type
        chantype_idxs=find(cosmo_match(ds_chantypes,parent_type));

        % define mask with channels matching those feature indices
        chan_msk=cosmo_match(ds_tl.fa.chan,chantype_idxs);

        % slice the dataset to select only the channels matching the channel
        % types
        ds_tl_sel=cosmo_dim_slice(ds_tl, chan_msk, 2);
        ds_tl_sel=cosmo_dim_prune(ds_tl_sel); % remove non-indexed channels

        % define neighborhood over time; for each time point the time
        % point itself is included, as well as the two time points before
        % and the two time points after it
        nbrhood=cosmo_interval_neighborhood(ds_tl_sel,'time',...
                                                'radius',time_radius);

        % run the searchlight using the measure, measure arguments, and
        % neighborhood defined above.
        % Note that while the input has both 'chan' and 'time' as feature
        % dimensions, the output only has 'time' as the feature dimension
        sl_map=cosmo_searchlight(ds_tl_sel,nbrhood,measure,measure_args);
        fprintf('The output has feature dimensions: %s\n', ...
                        cosmo_strjoin(sl_map.a.fdim.labels,', '));

        plot_counter=plot_counter+1;
        subplot(n_average_train_args,nchantypes,plot_counter);

        time_values=sl_map.a.fdim.values{1}; % first dim (channels got nuked)
        plot(time_values,sl_map.samples);

        ylim([.4 .8])
        xlim([min(time_values),max(time_values)]);
        ylabel('classification accuracy (chance=.5)');
        xlabel('time');

        if isempty(average_train_args)
            postfix=' no averaging';
        else
            postfix=' with averaging';
        end

        descr=sprintf('%s - %s', strrep(parent_type,'_',' '), postfix);
        title(descr);
    end
end
% Show citation information
cosmo_check_external('-cite');

demo_meeg_timeseries_generalization

%% MEEG time-by-time transfer classification
%
% This example shows MVPA generalization across time. The input is a
% time-locked dataset (chan x time), the output is either:
%
% * (train_time x test_time) indicating for each combination of time points
%   in the training and test set how similar the patterns are
% * (train_time x test_time x chan), which uses a searchlight and
%   indicates for each combination of time points in the training and
%   test set how similar the patterns are in the neighborhood of each
%   channel.
%
% Results can be visualized in FieldTrip.
%
% The input dataset involved a paradigm with electrical median nerve
% stimulation for durations of 2s at 20Hz.
%
% The code presented here can be adapted for other MEEG analyses, but
% there are a few potential caveats:
%
% * assignment of targets (labels of conditions) is based here on
%   stimulation periods versus pre-stimulation periods. In typical
%   analyses the targets should be based on different trial conditions, for
%   example as set a FieldTrip .trialinfo field.
% * the time window used for analyses is rather small. This means that in
%   particular for time-freq analysis a lot of data is missing, especially
%   for early and late timepoints in the lower frequency bands. For typical
%   analyses it may be preferred to use a wider time window.
% * 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.           #

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

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

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

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

% convert to cosmomvpa struct
ds_tl=cosmo_meeg_dataset(data_tl);

% set the target (trial condition)
ds_tl.sa.targets=ds_tl.sa.trialinfo(:,1); % 1=pre, 2=peri/post

% set the chunks (independent measurements)
% in this dataset, the first half of the samples (in order)
% are the post-trials;
% the second half the pre-trials
ds_tl.sa.chunks=[(1:145) (1:145)]';

% in addition give a label to each trial
index2label={'pre','post'}; % 1=pre, 2=peri/post
ds_tl.sa.labels=cellfun(@(x)index2label(x),num2cell(ds_tl.sa.targets));

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

%% reduce the number of time points
% to allow this example to run relatively fast, only use the 'even'
% timepoints and only use timepoints between 0 and 300 ms relative to
% stimulus onset.
% for a publication-quality analysis this step could be omitted.
msk=cosmo_dim_match(ds_tl,'time',@(x) 0<x &x<.3) & ...
                        mod(ds_tl.fa.time,2)==0;
ds=cosmo_slice(ds_tl,msk,2);
ds=cosmo_dim_prune(ds);


%% Set up chunks
% use half of the data for training and half for testing. Here both chunks
% come from the same session, but it is also possible to use different
% sessions for training and testing.
%
% For example, the train data can contain responses to auditory stimuli
% with the words 'face' and 'house', while the test session measures
% responses to visual stimuli. In that case the chunks should be assigned
% 'manually' (without cosmo_chunkize)
nchunks=2; % two chunks are required for this analysis
ds.sa.chunks=cosmo_chunkize(ds,nchunks);


%% time-by-time generalization on magneto- and gradio-meters seperately
% compute and plot accuracies for magnetometers and gradiometers separately
chan_types={'meg_axial','meg_planar'};
nchan_types=numel(chan_types);

% set arguments for the cosmo_dim_generalization_measure
measure_args=struct();

% the cosmo_dim_generalization_measure requires that another
% measure (here: the crossvalidation measure) is specified. The
% specified measure is applied for each combination of time points
measure_args.measure=@cosmo_crossvalidation_measure;

% When used ordinary, the cosmo_crossvalidation_measure itself
% requires two arguments:
% - classifier (here: LDA)
% - partitions
% However, because the cosmo_dim_generalization_measure defines
% the partitions itself, they are not set here.
measure_args.classifier=@cosmo_classify_lda;

% define the dimension over which generalization takes place
measure_args.dimension='time';

% define the radius for the time dimension. Here not just a single
% time-point is used, but also the time-point before it and the time-point
% after it.
measure_args.radius=1;

% try both channel types (axial and planar)
for k=1:nchan_types
    use_chan_type=chan_types{k};

    % get layout for mag or planar labels, and select only those channels
    ds_chan_types=cosmo_meeg_chantype(ds);
    ds_chan_idxs=find(cosmo_match(ds_chan_types,use_chan_type));
    feature_msk=cosmo_match(ds.fa.chan,ds_chan_idxs);
    ds_sel=cosmo_slice(ds, feature_msk, 2);

    % make 'time' a sample dimension
    % (this necessary for cosmo_dim_generalization_measure)
    ds_time=cosmo_dim_transpose(ds_sel,'time',1);

    fprintf('The input for channel type %s is:\n', use_chan_type)
    cosmo_disp(ds_time);

    % run transfer across time with the searchlight neighborhood
    %cdt_ds=cosmo_cartesian_dim_transfer(ds_time,'time',measure,...
    %                        'args',measure_args);
    cdt_ds=cosmo_dim_generalization_measure(ds_time,measure_args);

    fprintf('The output is:\n')
    cosmo_disp(cdt_ds);

    % unflatten the data to get train_time x test_time matrix
    [data, labels, values]=cosmo_unflatten(cdt_ds,1);

    % show the results
    figure()
    imagesc(data, [.3 .7]);
    title(sprintf('classification accuracy for %s', use_chan_type));
    colorbar();
    nticks=5;

    ytick=round(linspace(1, numel(values{1}), nticks));
    ylabel(strrep(labels{1},'_',' '));
    set(gca,'Ytick',ytick,'YTickLabel',values{1}(ytick));

    xtick=round(linspace(1, numel(values{2}), nticks));
    xlabel(strrep(labels{2},'_',' '));
    set(gca,'Xtick',xtick,'XTickLabel',values{2}(xtick));

    colorbar();
end


%% time-by-time generalization using a searchlight with correlation measure
% how many channel locations in each searchlight
nchannel_locations_searchlight=10;


% make 'time' a sample dimension (necessary for cartesian_dim_transfer)
ds_time=cosmo_dim_transpose(ds_sel,'time',1);


% use planar channels as input; output is
% like combine-planar
use_chan_type='meg_combined_from_planar';

% define searchlight neighborhood
nbrhood=cosmo_meeg_chan_neighborhood(ds_time,...
                        'count',nchannel_locations_searchlight,...
                        'chantype',use_chan_type);

fprintf('The input is:\n')
cosmo_disp(ds_time);

fprintf('The neighborhood is:\n');
cosmo_disp(nbrhood);

% set the measure to be the dim generalization measure
measure=@cosmo_dim_generalization_measure;

% define the arguments for the measure. One of them is called 'measure',
% because that measure is applied to combinations of samples in the
% train and test set
measure_args=struct();
measure_args.measure=@cosmo_correlation_measure;
measure_args.radius=1;
measure_args.dimension='time';


% run transfer across time with the searchlight neighborhood
cdt_ds=cosmo_searchlight(ds_time,nbrhood,measure,measure_args);

fprintf('The output is:\n')
cosmo_disp(cdt_ds)

% move {train,test}_time from being sample dimensions to feature
% dimensions, so they can be mapped to a fieldtrip struct
cdt_tf_ds=cosmo_dim_transpose(cdt_ds,{'train_time','test_time'},2);

fprintf('The output after transposing is:\n')
cosmo_disp(cdt_ds)

% trick fieldtrip into thinking this is a time-freq-chan dataset, by
% renaming train_time and test_time to freq and time, respectively
cdt_tf_ds=cosmo_dim_rename(cdt_tf_ds,'train_time','freq');
cdt_tf_ds=cosmo_dim_rename(cdt_tf_ds,'test_time','time');

%% Show correlation time-by-time generalization searchlight results

% determine layout
layout=cosmo_meeg_find_layout(cdt_tf_ds);


% convert to fieldtrip format
ft=cosmo_map2meeg(cdt_tf_ds);

% show train_time x test_time x chan figure
%
% * train_time is on the vertical ('freq') axis
% * test_time on the horizontal ('time') axis)
figure();
cfg = [];
cfg.interactive = 'yes';
cfg.showlabels = 'yes';
cfg.zlim=[-1 1];
cfg.layout       = layout;
ft_multiplotTFR(cfg, ft);

% show train_time x test_time figure for selected channels
figure();
cfg = [];
cfg.interactive = 'yes';
cfg.showlabels = 'yes';
cfg.zlim=[-1 1];
cfg.layout       = layout;
cfg.channel={'MEG0322+0323', 'MEG0332+0333', 'MEG0412+0413', ...
            'MEG0422+0423', 'MEG0632+0633', 'MEG0642+0643'};

ft_singleplotTFR(cfg,ft);

% show channel topology for timewindow over train_time x test_time
figure();
cfg = [];
cfg.interactive = 'yes';
cfg.xlim=[.2 .25];
cfg.ylim=[.2 .25];
cfg.zlim=[-1 1];
cfg.layout       = layout;
ft_topoplotTFR(cfg, ft);


%% show citation information

cosmo_check_external('-cite');

demo_surface_searchlight_lda

%% Demo: fMRI surface-based searchlights with LDA classifier
%
% The data used here is available from http://cosmomvpa.org/datadb.zip
%
% This example uses the following dataset:
% + 'digit'
%    A participant made finger pressed with the index and middle finger of
%    the right hand during 4 runs in an fMRI study. Each run was divided in
%    4 blocks with presses of each finger and analyzed with the GLM,
%    resulting in 2*4*4=32 t-values
%
% The example shows four possible searchlight analyses, covering typical
% use cases:
%   1) single or twin surfaces
%      + Caret and BrainVoyager use a single surface; a parameter 'offsets'
%        is used to define which voxels are considered to be in the "grey
%        matter" (but this may be not so precise)
%      + FreeSurfer uses twin surfaces (pial and white), and voxels in
%        between or on them are considered to be in the grey matter
%   2) lower resolution output map
%      + in the canonical surface-based searchlight, each node on the input
%        surface(s) is assigned a measure value (accuracy, in this example)
%      + it is also possible to have output in a lower resolution version
%        than the input surfaces; this reduces both the execution time
%        (a Good Thing) and spatial precision (a Bad Thing). Two approaches
%        are illustrated to use a lower resolution surface for output:
%        1) from MapIcosahedron, with a lower value for the number of
%          divisions of the triangles
%        2) using a surface subsampling approach, implemented by
%           surfing_subsample_surface
%
% In all cases a searchlight is run with a 100 voxel searchlight, using a
% disc for which the metric radius varies from node to node. For a fixed
% metric radius of the disc, use a positive value for 'radius' below.
% Distances are measured across the cortical surface using a geodesic
% distance metric.
%
% This example requires the surfing toolbox, github.com/nno/surfing
%
% This example may take quite some time to run. For faster execution, set
% ld=16 (instead of ld=64) below
%
% If you use this code for a publication, please cite:
% Oosterhof, N.N., Wiestler, T, Downing, P.E., & Diedrichsen, J. (2011)
% A comparison of volume-based and surface-based information mapping.
% Neuroimage. DOI:10.1016/j.neuroimage.2010.04.270
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #

%% Check externals
cosmo_check_external('surfing');

%% Set data paths
% The function cosmo_config() returns a struct containing paths to tutorial
% data. (Alternatively the paths can be set manually without using
% cosmo_config.)
config=cosmo_config();

digit_study_path=fullfile(config.tutorial_data_path,'digit');
readme_fn=fullfile(digit_study_path,'README');
cosmo_type(readme_fn);

output_path=config.output_data_path;

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

% resolution parameter for input surfaces
% 64 is for high-quality results; use 16 for fast execution
ld=64;

% Twin surface case (FS)
pial_fn=fullfile(digit_study_path,...
                            sprintf('ico%d_mh.pial_al.asc', ld));
white_fn=fullfile(digit_study_path,...
                            sprintf('ico%d_mh.smoothwm_al.asc', ld));

% Single surface case (Caret/BV)
intermediate_fn=fullfile(digit_study_path,...
                            sprintf('ico%d_mh.intermediate_al.asc', ld));

% Used for visualization
inflated_fn=fullfile(digit_study_path,...
                         sprintf('ico%d_mh.inflated_alCoMmedial.asc', ld));

%%
% Set parameters

% Searchlight radius: select 100 features in each searchlight
% (to use a fixed radius of 8mm, use:
%    cosmo_surficial_neighborhood(...,'radius',8)
% in the code below)

feature_count=100;

% Single surface case: select voxels that are 3 mm or closer to the surface
% on the white-matter side, up to voxels that are 2 mm from the surface on
% the pial matter side
single_surf_offsets=[-2 3];

% Single surface case: number of iterations to downsample surface
lowres_output_onesurf_niter=10;

% Twin surface case: number of linear divisions from MapIcosahedron
lowres_output_twosurf_icold=16;
lowres_intermediate_fn=fullfile(digit_study_path,...
                                sprintf('ico%d_mh.intermediate_al.asc',...
                                        lowres_output_twosurf_icold));


% Use the cosmo_cross_validation_measure and set its parameters
% (classifier and partitions) in a measure_args struct.
measure = @cosmo_crossvalidation_measure;
measure_args = struct();

% Define which classifier to use, using a function handle.
% Alternatives are @cosmo_classify_{svm,nn,naive_bayes}
measure_args.classifier = @cosmo_classify_lda;



%% Load functional data
data_path=digit_study_path;
data_fn=fullfile(data_path,'glm_T_stats_perblock+orig');

targets=repmat(1:2,1,16)';    % class labels: 1 2 1 2 1 2 1 2 1 2 ... 1 2
chunks=floor(((1:32)-1)/8)+1; % run labels:   1 1 1 1 1 1 1 1 2 2 ... 4 4

ds = cosmo_fmri_dataset(data_fn,'targets',targets,'chunks',chunks);

% remove zero elements
zero_msk=all(ds.samples==0,1);
ds = cosmo_slice(ds, ~zero_msk, 2);

fprintf('Dataset has %d samples and %d features\n', size(ds.samples));

% print dataset
fprintf('Dataset input:\n');
cosmo_disp(ds);

%% Set partition scheme. odd_even is fast; for publication-quality analysis
% nfold_partitioner is recommended.
% Alternatives are:
% - cosmo_nfold_partitioner    (take-one-chunk-out crossvalidation)
% - cosmo_nchoosek_partitioner (take-K-chunks-out  "             ").
measure_args.partitions = cosmo_oddeven_partitioner(ds);

% print measure and arguments
fprintf('Searchlight measure:\n');
cosmo_disp(measure);
fprintf('Searchlight measure arguments:\n');
cosmo_disp(measure_args);

%% Read inflated surface
[v_inf,f_inf]=surfing_read(inflated_fn);
fprintf('The inflated surface has %d vertices, %d faces\n',...
            size(v_inf,1), size(f_inf,1))

%% Run four types of searchlights
for one_surf=[true,false]
    if one_surf
        desc='1surf';
    else
        desc='2surfs';
    end

    for lowres_output=[false,true]
        if lowres_output
            desc=sprintf('%s_lowres', desc);
        end
        fprintf('\n\n *** Starting analysis with %s *** \n\n\n', desc)

        % define searchlight surface paramters for each type of analysis
        if one_surf && lowres_output

            % single surface (Caret/BV) with lower-res output
            surf_def={intermediate_fn,single_surf_offsets,...
                            lowres_output_onesurf_niter};

        elseif one_surf && ~lowres_output

            % single surface (Caret/BV) with original-res output
            surf_def={intermediate_fn,single_surf_offsets};

        elseif ~one_surf && lowres_output

            % single surface (FS) with lower-res output
            surf_def={white_fn,pial_fn,lowres_intermediate_fn};

        elseif ~one_surf && ~lowres_output

            % single surface (FS) with original-res output
            surf_def={white_fn,pial_fn};

        else
            assert(false); % should never get here
        end

        % Define the feature neighborhood for each node on the surface
        % - nbrhood has the neighborhood information
        % - vo and fo are vertices and faces of the output surface
        % - out2in is the mapping from output to input surface
        fprintf('Defining neighborhood with %s\n', desc);
        [nbrhood,vo,fo,out2in]=cosmo_surficial_neighborhood(ds,surf_def,...
                                                    'count',feature_count);

        % print neighborhood
        fprintf('Searchlight neighborhood definition:\n');
        cosmo_disp(nbrhood);


        fprintf('The output surface has %d vertices, %d nodes\n',...
                        size(vo,1), size(fo,1));



        % Run the searchlight
        lda_results = cosmo_searchlight(ds,nbrhood,measure,measure_args);


        % print searchlight output
        fprintf('Dataset output:\n');
        cosmo_disp(lda_results);

        % Apply the node mapping from the surifical neighborhood
        % to the high-res inflated surface.
        % (This example shows how such a mapping can be applied to new
        % surfaces)
        if lowres_output
            v_inf_out=v_inf(out2in,:);
            f_inf_out=fo;
        else
            v_inf_out=v_inf;
            f_inf_out=f_inf;
        end

        % visualize the surfaces, if the afni matlab toolbox is present
        if cosmo_check_external('afni',false)
            nvertices=size(v_inf_out,1);

            opt=struct();

            for show_edge=[false, true]
                opt.ShowEdge=show_edge;

                if show_edge
                    t='with edges';
                else
                    t='without edges';
                end

                header=strrep([desc ' ' t],'_',' ');


                DispIVSurf(vo,fo,1:nvertices,lda_results.samples',0,opt);
                title(sprintf('Original %s', header));

                DispIVSurf(v_inf_out,f_inf_out,1:nvertices,...
                                        lda_results.samples',0,opt);
                title(sprintf('Inflated %s', header));
            end
        else
            fprintf('skip surface display; no afni matlab toolbox\n');
        end

        if lowres_output && one_surf
            % in this example only this case a new surface was generated.
            % To aid visualization using external tools, store it to disc.

            % The surface is stored in ASCII, GIFTI and BV SRF
            % formats, if the required externals are present
            surf_output_fn=fullfile(output_path,['inflated_' desc]);

            % AFNI/SUMA ASC
            surfing_write([surf_output_fn '.asc'],v_inf_out,f_inf_out);

            % BV SRF
            if cosmo_check_external('neuroelf',false)
                surfing_write([surf_output_fn '.srf'],v_inf_out,f_inf_out);
            end

            % GIFTI
            if cosmo_check_external('gifti',false)
                surfing_write([surf_output_fn '.gii'],v_inf_out,f_inf_out);
            end
        end

        % store searchlight results
        data_output_fn=fullfile(output_path,['lda_' desc]);

        if cosmo_check_external('afni',false)
            cosmo_map2surface(lda_results, [data_output_fn '.niml.dset']);
        end

        if cosmo_check_external('neuroelf',false)
            cosmo_map2surface(lda_results, [data_output_fn '.smp']);
        end

        % store voxel counts (how often each voxel is in a neighborhood)
        % take a random sample (the first one) from the input dataset
        % and count how often each voxel was selected.
        % If everything works, then voxels in the grey matter have high
        % voxel counts but voxels outside it low or zero counts.
        % Thus, this can be used as a sanity check that can be visualized
        % easily.

        vox_count_ds=cosmo_slice(ds,1);
        vox_count_ds.samples(:)=0;

        ncenters=numel(nbrhood.neighbors);
        for k=1:ncenters
            idxs=nbrhood.neighbors{k}; % feature indices in neigborhood
            vox_count_ds.samples(idxs)=vox_count_ds.samples(idxs)+1;
        end

        vox_count_output_fn=fullfile(output_path,['vox_count_' desc]);

        % store voxel count results
        cosmo_map2fmri(vox_count_ds, [vox_count_output_fn '.nii']);

        if cosmo_check_external('afni',false)
            cosmo_map2fmri(vox_count_ds, [vox_count_output_fn '+orig']);
        end
    end
end

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

demo_surface_tfce

%% Demo: Threshold-Free Cluster Enhancement (TFCE) on surface dataset
%
% The data used here is available from http://cosmomvpa.org/datadb.zip
%
% This example uses the following dataset:
% + 'digit'
%    A participant made finger pressed with the index and middle finger of
%    the right hand during 4 runs in an fMRI study. Each run was divided in
%    4 blocks with presses of each finger and analyzed with the GLM,
%    resulting in 2*4*4=32 t-values
%
% This example illustrates the use of Threshold-Free Cluster Enhancement
% with a permutation test to correct for multiple comparisons.
%
% TFCE reference: Stephen M. Smith, Thomas E. Nichols, Threshold-free
% cluster enhancement: Addressing problems of smoothing, threshold
% dependence and localisation in cluster inference, NeuroImage, Volume 44,
% Issue 1, 1 January 2009, Pages 83-98.
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #

%% Check externals
cosmo_check_external({'surfing','afni'});

%% Set data paths
% The function cosmo_config() returns a struct containing paths to tutorial
% data. (Alternatively the paths can be set manually without using
% cosmo_config.)
config=cosmo_config();

digit_study_path=fullfile(config.tutorial_data_path,'digit');
readme_fn=fullfile(digit_study_path,'README');
cosmo_type(readme_fn);

output_path=config.output_data_path;

% resolution parameter for input surfaces
% 64 is for high-quality results; use 16 for fast execution
ld=16;

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

% load single surface
intermediate_fn=fullfile(digit_study_path,...
                            sprintf('ico%d_mh.intermediate_al.asc', ld));
[vertices,faces]=surfing_read(intermediate_fn);


%% Load functional data
data_path=digit_study_path;
data_fn=fullfile(data_path,'glm_T_stats_perblock+orig');

targets=repmat(1:2,1,16)';    % class labels:  1 2 1 2 1 2 1 2 1 2 ... 1 2
chunks=floor(((1:32)-1)/4)+1; % half-run:      1 1 1 1 2 2 2 2 3 3 ... 8 8

vol_ds = cosmo_fmri_dataset(data_fn,'targets',targets,'chunks',chunks);

%% Map univariate response data to surface

% this measure averages the data near each node to get a surface dataset
radius=0;
surf_band_range=[-2 2]; % get voxel data within 2mm from surface
surf_def={vertices,faces,[-2 2]};
nbrhood=cosmo_surficial_neighborhood(vol_ds,surf_def,'radius',radius);

measure=@(x,opt) cosmo_structjoin('samples',mean(x.samples,2),'sa',x.sa);

surf_ds=cosmo_searchlight(vol_ds,nbrhood,measure);

fprintf('Univariate surface data:\n');
cosmo_disp(surf_ds);

%% Average data in each chunk
% for this example only consider the samples in the first condition
% (targets==1), and average the samples in each chunk
%
% for group analysis: set chunks to (1:nsubj)', assuming each sample is
% data from a single participant
surf_ds=cosmo_slice(surf_ds,surf_ds.sa.targets==1);
surf_ds=cosmo_average_samples(surf_ds);

fn_surf_ds=fullfile(output_path, 'digit_target1.niml.dset');

% save to disc
cosmo_map2surface(surf_ds, fn_surf_ds);
fprintf('Input data saved to %s\n', fn_surf_ds);

%% Run Threshold-Free Cluster Enhancement (TFCE)

% All data is prepared; surf_ds has 8 samples and 5124 nodes. We want to
% see if there are clusters that show a significant difference from zero in
% their response. Thus, .sa.targets is set to all ones (the same
% condition), whereas .sa.chunks is set to (1:8)', indicating that all
% samples are assumed to be independent.
%
% (While this is a within-subject analysis, exactly the same logic can be
% applied to a group-level analysis)

% define neighborhood for each feature
% (cosmo_cluster_neighborhood can be used also for meeg or volumetric
% fmri datasets)
cluster_nbrhood=cosmo_cluster_neighborhood(surf_ds,...
                                        'vertices',vertices,'faces',faces);

fprintf('Cluster neighborhood:\n');
cosmo_disp(cluster_nbrhood);

opt=struct();

% number of null iterations. for publication-quality, use >=1000;
% 10000 is even better
opt.niter=250;

% in this case we run a one-sample test against a mean of 0, and it is
% necessary to specify the mean under the null hypothesis
% (when testing classification accuracies, h0_mean should be set to chance
% level, assuming a balanced crossvalidation scheme was used)
opt.h0_mean=0;

% this example uses the data itself (with resampling) to obtain cluster
% statistcs under the null hypothesis. This is (in this case) somewhat
% conservative due to how the resampling is performed.
% Alternatively, and for better estimates (at the cost of computational
% cost), one can generate a set of (say, 50) datasets using permuted data
% e.g. using cosmo_randomize_targets), put them in a cell and provide
% them as the null argument.
opt.null=[];

fprintf('Running multiple-comparison correction with these options:\n');
cosmo_disp(opt);

% Run TFCE-based cluster correction for multiple comparisons.
% The output has z-scores for each node indicating the probablity to find
% the same, or higher, TFCE value under the null hypothesis
tfce_ds=cosmo_montecarlo_cluster_stat(surf_ds,cluster_nbrhood,opt);

%% Show results

fprintf('TFCE z-score dataset\n');
cosmo_disp(tfce_ds);

nfeatures=size(tfce_ds.samples,2);
percentiles=(1:nfeatures)/nfeatures*100;
plot(percentiles,sort(tfce_ds.samples))
title('sorted TFCE z-scores');
xlabel('feature percentile');
ylabel('z-score');


nvertices=size(vertices,1);
disp_opt=struct();
disp_opt.DataRange=[-2 2];

DispIVSurf(vertices,faces,1:nvertices,tfce_ds.samples',0,disp_opt);

% store results
fn_tfce_ds=fullfile(output_path, 'digit_target1_tfce.niml.dset');
cosmo_map2surface(tfce_ds, fn_tfce_ds);

surf_fn=fullfile(output_path, 'digit_intermediate.asc');
surfing_write(surf_fn,vertices,faces);


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