function test_suite = test_meeg_io()
% tests for MEEG input/output
%
% # For CoSMoMVPA's copyright information and license terms, #
% # see the COPYING file distributed with CoSMoMVPA. #
try % assignment of 'localfunctions' is necessary in Matlab >= 2016
test_functions = localfunctions();
catch % no problem; early Matlab versions can use initTestSuite fine
end
initTestSuite;
function test_meeg_ft_dataset()
aet = @(varargin)assertExceptionThrown(@() ...
cosmo_meeg_dataset(varargin{:}), '');
dimords = get_ft_dimords();
n = numel(dimords);
for k = 1:n
dimord = dimords{k};
[ft, fdim, data_label] = generate_ft_struct(dimord);
ds = cosmo_meeg_dataset(ft);
assertEqual(ds.a.fdim, fdim);
[nsamples, nfeatures] = size(ds.samples);
% check feature sizes
fdim_sizes = cellfun(@numel, fdim.values);
assertEqual(prod(fdim_sizes), nfeatures);
% check sample size
data = ft.(data_label);
data_size = size(data);
has_rpt = nsamples > 1;
if has_rpt
assertEqual(data_size(2:end)', fdim_sizes);
else
assertEqual(data_size', fdim_sizes);
end
assertElementsAlmostEqual(data(:), ds.samples(:));
ds2 = cosmo_slice(ds, randperm(nfeatures), 2);
ft2 = cosmo_map2meeg(ds2);
if isfield(ft, 'cfg')
ft = rmfield(ft, 'cfg');
end
if isfield(ft, 'avg') && isfield(ft, 'trial')
ft = rmfield(ft, 'avg');
end
assertEqual(ft, ft2);
% wrong size trialinfo should not store trialinfo
assertTrue(isfield(ds2.sa, 'trialinfo'));
ft.trialinfo = [1; 2];
ds3 = cosmo_meeg_dataset(ft);
assertFalse(isfield(ds3.sa, 'trialinfo'));
end
aet(struct());
aet(struct('avg', 1));
aet(struct('avg', 1, 'dimord', 'rpt_foo'));
function test_meeg_ft_dataset_trials()
aet = @(varargin)assertExceptionThrown(@() ...
cosmo_meeg_dataset(varargin{:}), '');
dimords = get_ft_dimords();
n = numel(dimords);
for k = 1:n
dimord = dimords{k};
ft = generate_ft_struct(dimord);
ds = cosmo_meeg_dataset(ft);
% check subset of trials option
ntrials = size(ds.samples, 1);
trial_idx = ceil(rand(1, 2) * ntrials);
ds_single_trial = cosmo_meeg_dataset(ft, ...
'trials', trial_idx);
assertEqual(cosmo_slice(ds, trial_idx), ds_single_trial);
ds_single_trial = cosmo_meeg_dataset(ft, ...
cosmo_structjoin('trials', trial_idx));
assertEqual(cosmo_slice(ds, trial_idx), ds_single_trial);
illegal_args = {ntrials + 1, 0, struct, cell(1, 0), 'foo', true, 1.5};
for j = 1:numel(illegal_args)
arg = illegal_args{j};
aet(ft, 'trials', arg);
end
end
function test_synthetic_meeg_dataset()
combis = cosmo_cartprod({{'timelock', 'timefreq', 'source'}, ...
{'tiny', 'small', 'normal', 'big', 'huge'}});
for k = 1:4:size(combis, 1)
ds = cosmo_synthetic_dataset('type', combis{k, 1}, ...
'size', combis{k, 2});
ft = cosmo_map2meeg(ds);
ds2 = cosmo_meeg_dataset(ft);
assertEqual(ds.samples, ds2.samples);
assertEqual(ds.fa, ds2.fa);
assertEqual(ds.a.meeg.samples_field, ds2.a.meeg.samples_field);
end
ds2 = cosmo_meeg_dataset(ds, 'targets', 1);
assertTrue(all(ds2.sa.targets == 1));
assertExceptionThrown(@()cosmo_meeg_dataset(ds, 'targets', [1 2]), '');
function test_meeg_eeglab_txt_io()
ds = cosmo_synthetic_dataset('type', 'meeg');
tmp_fn = sprintf('_tmp_%06.0f.txt', rand() * 1e5);
file_remover = onCleanup(@()delete(tmp_fn));
fid = fopen(tmp_fn, 'w');
file_closer = onCleanup(@()fclose(fid));
chans = [{' '} ds.a.fdim.values{1}];
fprintf(fid, '%s\t', chans{:});
fprintf(fid, '\n');
times = ds.a.fdim.values{2};
ntime = numel(times);
nsamples = size(ds.samples, 1);
for k = 1:nsamples
for j = 1:ntime
data = ds.samples(k, ds.fa.time == j);
fprintf(fid, '%.3f', times(j));
fprintf(fid, '\t%.4f', data);
fprintf(fid, '\n');
end
end
fclose(fid);
ds2 = cosmo_meeg_dataset(tmp_fn);
assertElementsAlmostEqual(ds.samples, ds2.samples, 'absolute', 1e-4);
assertEqual(ds.a.fdim.values{1}, ds2.a.fdim.values{1});
assertElementsAlmostEqual(ds.a.fdim.values{2}, ...
1000 * ds2.a.fdim.values{2});
assertEqual(ds.a.fdim.labels, ds2.a.fdim.labels);
assertEqual(ds.fa, ds2.fa);
% test trials option
nsamples = size(ds.samples, 1);
trial_idx = ceil(rand(1, 2) * nsamples);
ds_trials = cosmo_meeg_dataset(tmp_fn, 'trials', trial_idx);
ds_expected_trials = cosmo_slice(ds, trial_idx);
assertElementsAlmostEqual(ds_trials.samples, ...
ds_expected_trials.samples, ...
'absolute', 1e-4);
% test illegal options
aet = @(varargin)assertExceptionThrown(@() ...
cosmo_meeg_dataset(varargin{:}), '');
illegal_args = {nsamples + 1, 0, struct, {}, 'foo', true, 1.5};
for j = 1:numel(illegal_args)
arg = illegal_args{j};
aet(tmp_fn, 'trials', arg);
aet(tmp_fn, cosmo_structjoin('trials', arg));
end
% add bogus data, expect exception
fid = fopen(tmp_fn, 'a');
fprintf(fid, '.3');
fclose(fid);
file_closer = [];
aet(tmp_fn);
tmp2_fn = sprintf('_tmp_%06.0f.txt', rand() * 1e5);
file_remover2 = onCleanup(@()delete(tmp2_fn));
cosmo_map2meeg(ds2, tmp2_fn);
ds3 = cosmo_meeg_dataset(tmp2_fn);
assertEqual(ds2, ds3);
function test_meeg_ft_io()
ds = cosmo_synthetic_dataset('type', 'meeg');
tmp_fn = sprintf('_tmp_%06.0f.mat', rand() * 1e5);
file_remover = onCleanup(@()delete(tmp_fn));
cosmo_map2meeg(ds, tmp_fn);
ds2 = cosmo_meeg_dataset(tmp_fn);
assertEqual(ds.a.meeg.samples_field, ds2.a.meeg.samples_field);
ds.a.meeg = [];
ds.sa = struct();
ds2.a.meeg = [];
% deal with rounding errors in Octave
assertElementsAlmostEqual(ds.samples, ds2.samples);
assertElementsAlmostEqual(ds.a.fdim.values{2}, ds2.a.fdim.values{2});
ds3 = cosmo_meeg_dataset(ds2);
assertEqual(ds2, ds3);
ds2.samples = ds.samples;
ds2.a.fdim.values{2} = ds.a.fdim.values{2};
assertEqual(ds, ds2);
function test_meeg_ft_io_exceptions()
aeti = @(varargin)assertExceptionThrown(@() ...
cosmo_meeg_dataset(varargin{:}), '');
aeto = @(varargin)assertExceptionThrown(@() ...
cosmo_map2meeg(varargin{:}), '');
ds = cosmo_synthetic_dataset('type', 'timefreq');
aeti('file_without_extension');
aeti('file.with_unknown_extension');
aeto(ds, 'file_without_extension');
aeto(ds, 'file.with_unknown_extension');
aeto(ds, 'eeglab_timelock.txt'); % not supported
function dimords = get_ft_dimords()
dimords = { 'chan_time', ...
'rpt_chan_time'...
'subj_chan_time'...
'chan_freq', ...
'rpt_chan_freq', ...
'subj_chan_freq', ...
'chan_freq_time', ...
'rpt_chan_freq_time', ...
'subj_chan_freq_time' ...
};
function [ft, fdim, data_label] = generate_ft_struct(dimord)
seed = 1;
fdim = struct();
fdim.values = cell(3, 1);
fdim.labels = cell(3, 1);
ft = struct();
ft.dimord = dimord;
dims = cosmo_strsplit(dimord, '_');
ndim = numel(dims);
sizes = [3 4 5 6];
chan_values = {'MEG0113' 'MEG0112' 'MEG0111' 'MEG0122'...
'MEG0123' 'MEG0121' 'MEG0132'};
freq_values = (2:2:24);
time_values = (-1:.1:2);
data_label = 'avg';
ntrials = 1;
nkeep = 0;
for k = 1:ndim
idxs = 1:sizes(k);
switch dims{k}
case 'rpt'
data_label = 'trial';
ntrials = numel(idxs);
case 'subj'
data_label = 'individual';
ntrials = numel(idxs);
case 'chan'
ft.label = chan_values(idxs);
nkeep = nkeep + 1;
fdim.values{nkeep} = ft.label;
fdim.labels{nkeep} = 'chan';
case 'freq'
ft.freq = freq_values(idxs);
data_label = 'powspctrm';
nkeep = nkeep + 1;
fdim.values{nkeep} = ft.freq;
fdim.labels{nkeep} = 'freq';
case 'time'
ft.time = time_values(idxs);
nkeep = nkeep + 1;
fdim.values{nkeep} = ft.time;
fdim.labels{nkeep} = 'time';
end
end
fdim.values = fdim.values(1:nkeep);
fdim.labels = fdim.labels(1:nkeep);
keep_sizes = sizes(1:k);
ft.(data_label) = cosmo_norminv(cosmo_rand(keep_sizes, 'seed', seed));
ft.cfg = struct();
ft.trialinfo = [(1:ntrials); (ntrials:-1:1)]';
if strcmp(data_label, 'trial')
ft.avg = mean(ft.(data_label), 1);
end
function test_eeglab_io()
datatypes = {'timef', 'erp', 'itc'};
args = cosmo_cartprod({{true, false}, ...
{true, false}, ...
datatypes});
ncombi = size(args, 1);
for k = 1:ncombi
arg = args(k, :);
[s, ds, ext] = build_eeglab_dataset_struct(arg{:});
ds_from_struct = cosmo_meeg_dataset(s);
assertEqual(ds.samples, ds_from_struct.samples);
assertEqual(ds, ds_from_struct);
% store, then read using cosmo_meeg_dataset
fn = sprintf('%s.%s', tempname(), ext);
save(fn, '-mat', '-struct', 's');
cleaner = onCleanup(@()delete(fn));
ds_loaded = cosmo_meeg_dataset(fn);
assertEqual(ds, ds_loaded);
clear cleaner;
s_converted = cosmo_map2meeg(ds, ['-' ext]);
assertEqual(s, s_converted);
% store using cosmo_map2meeg, then read
cosmo_map2meeg(ds, fn);
cleaner = onCleanup(@()delete(fn));
s_loaded = load(fn, '-mat');
assertEqual(s_loaded, s);
assertEqual(s, s_loaded);
clear cleaner;
end
function test_eeglab_io_trials()
% test with loading a subset of trials
datatypes = {'timef', 'erp', 'itc'};
args = cosmo_cartprod({{true, false}, ...
{true, false}, ...
datatypes});
ncombi = size(args, 1);
for k = 1:ncombi
arg = args(k, :);
[s, ds, ext] = build_eeglab_dataset_struct(arg{:});
nsamples = size(ds.samples, 1);
trial_idx = ceil(rand(1, 2) * nsamples);
ds_expected_trials = cosmo_slice(ds, trial_idx);
% with struct input
ds_trials = cosmo_meeg_dataset(s, 'trials', trial_idx);
assertElementsAlmostEqual(ds_trials.samples, ...
ds_expected_trials.samples, ...
'absolute', 1e-4);
% store, then read using cosmo_meeg_dataset
fn = sprintf('%s.%s', tempname(), ext);
save(fn, '-mat', '-struct', 's');
cleaner = onCleanup(@()delete(fn));
ds_loaded = cosmo_meeg_dataset(fn, 'trials', trial_idx);
assertEqual(ds_loaded, ds_expected_trials);
% test illegal options
aet = @(varargin)assertExceptionThrown(@() ...
cosmo_meeg_dataset(varargin{:}), '');
illegal_args = {nsamples + 1, 0, struct, {}, 'foo', true, 1.5};
for j = 1:numel(illegal_args)
arg = illegal_args{j};
aet(s, 'trials', arg);
end
clear cleaner;
end
function test_eeglab_io_ersp()
aet = @(varargin)assertExceptionThrown(@() ...
cosmo_meeg_dataset(varargin{:}), '');
args = cosmo_cartprod({{true, false}, ...
{false}, ...
{'ersp'}, ...
{false, true}});
ncombi = size(args, 1);
for k = 1:ncombi
arg = args(k, 1:3);
[s, ds_cell, ext] = build_eeglab_dataset_struct(arg{:});
% either load baseline data or original data
with_baseline = args{k, 4};
if with_baseline
load_args = {'data_field', 'erspbase'};
ds = ds_cell{2};
else
load_args = {'data_field', 'ersp'};
ds = ds_cell{1};
end
% illegal without arguments or wrong arguments
aet(s);
aet(s, 'data_field', 'foo');
aet(s, 'data_field', false);
% load data with correct arguments
ds_loaded = cosmo_meeg_dataset(s, load_args{:});
assertEqual(ds_loaded, ds);
% store, then read using cosmo_meeg_dataset
fn = sprintf('%s.%s', tempname(), ext);
save(fn, '-mat', '-struct', 's');
cleaner = onCleanup(@()delete(fn));
ds_loaded = cosmo_meeg_dataset(fn, load_args);
% writing the file is not supported
assertExceptionThrown(@()cosmo_map2meeg(ds, fn), '');
clear cleaner;
assertEqual(ds_loaded, ds);
end
function test_eeglab_io_exceptions()
aet_md = @(varargin)assertExceptionThrown(@() ...
cosmo_meeg_dataset(varargin{:}), '');
aet_m2m = @(varargin)assertExceptionThrown(@() ...
cosmo_map2meeg(varargin{:}), '');
s = build_eeglab_dataset_struct(true, true, 'timef');
% bad datatype
s.datatype = 'foo';
aet_md(s);
% output is not a filename
ds = cosmo_synthetic_dataset('type', 'timefreq');
aet_m2m(ds, struct);
% bad fdim
good_labels = {'chan', 'freq', 'time'};
all_bad_labels = {'chan', 'freq', 'time', 'foo'};
for dim = 1:numel(good_labels)
for j = 1:numel(all_bad_labels)
ds_bad_chan_fdim = ds;
bad = all_bad_labels{j};
if ~strcmp(bad, good_labels{dim})
ds_bad_chan_fdim.a.fdim.labels{dim} = bad;
aet_m2m(ds_bad_chan_fdim, '-dattimef');
end
end
end
function [s, ds, ext] = build_eeglab_dataset_struct(has_ica, has_trial, datatype, ...
chan_dim, freq_dim, time_dim)
if nargin < 6
time_dim = randint();
end
if nargin < 5
freq_dim = randint();
end
if nargin < 4
chan_dim = randint();
end
% trial dimension
if has_trial
trial_dim = randint();
else
trial_dim = 1;
end
if strcmp(datatype, 'ersp')
% has baseline corrected data together with baseline data
builder = @build_eeglab_dataset_struct;
args = {chan_dim, freq_dim, time_dim};
[s1, ds1, ext] = builder(has_ica, has_trial, ...
'ersp_baselinecorrected', args{:});
[s2, ds2] = builder(has_ica, has_trial, ...
'erspbase', args{:});
keys = fieldnames(s1);
for k = 1:numel(keys)
key = keys{k};
s2.(key) = s1.(key);
end
s = s2;
s.datatype = upper(datatype);
% make sure parameters are the same
ds1.a.meeg.parameters = s.parameters;
ds2.a.meeg.parameters = s.parameters;
if isfield(s, 'chanlabels')
chan_labels = s.chanlabels;
ds1.a.fdim.values{1} = chan_labels;
ds2.a.fdim.values{1} = chan_labels;
end
ds = {ds1, ds2};
% remove second part from extension
ext = regexprep(ext, '_.*', '');
return
end
% channel / component dimension
if has_ica
chan_prefix = 'comp';
ext_prefix = 'ica';
make_chan_prefix_func = @()chan_prefix;
else
chan_prefix = 'chan';
ext_prefix = 'dat';
make_chan_prefix_func = @randstr;
end
has_freq = 2;
switch datatype
case 'timef'
chan_suffix = '_timef';
case 'erp'
chan_suffix = '';
case 'ersp_baselinecorrected'
chan_suffix = '_ersp';
case 'erspbase'
chan_suffix = '_erspbase';
case 'itc'
chan_suffix = '_itc';
otherwise
assert(false);
end
make_chan_label = @(idx) sprintf('%s%d', make_chan_prefix_func(), idx);
chan_label = {chan_prefix};
chan_value = {arrayfun(make_chan_label, 1:chan_dim, ...
'UniformOutput', false)};
% frequency dimension
switch datatype
case {'timef', 'ersp_baselinecorrected', 'itc', 'erspbase'}
has_freq = true;
case {'erp'}
has_freq = false;
otherwise
assert(false);
end
if has_freq
freq_label = {'freq'};
freq_value = {(1:freq_dim) * 2};
samples_type = 'timefreq';
else
freq_dim = [];
freq_label = {};
freq_value = {};
samples_type = 'timelock';
end
ext_suffix = datatype;
hastime = ~strcmp(datatype, 'erspbase');
if hastime
% include time dimension
time_label = {'time'};
time_value = {(1:time_dim) * .2 - .1};
else
% no time dimension
time_dim = [];
time_label = {};
time_value = {};
end
% data
dim_sizes = [trial_dim, chan_dim, freq_dim, time_dim];
dim_sizes_without_chan = [dim_sizes([1, 3:end]), 1];
data_arr = randn(dim_sizes);
% params
parameters = {randstr(), randstr()};
% make dataset
ds = cosmo_flatten(data_arr, ...
[chan_label, freq_label, time_label], ...
[chan_value, freq_value, time_value]);
ds.sa = struct();
ds.a.meeg.samples_field = 'trial';
ds.a.meeg.samples_type = samples_type;
ds.a.meeg.samples_label = 'rpt';
ds.a.meeg.parameters = parameters;
s = struct();
for k = 1:chan_dim
key = sprintf('%s%d%s', chan_prefix, k, chan_suffix);
value = data_arr(:, k, :);
value_rs = reshape(value, dim_sizes_without_chan);
if has_freq
% it seems that for freq data, single trial data is the last
% dimension, whereas for erp data, single trial data is the first
% dimension.
value_rs = shiftdim(value_rs, 1);
end
s.(key) = value_rs;
end
if ~has_ica
s.chanlabels = chan_value{1};
assert(iscellstr(s.chanlabels));
end
if has_freq
s.freqs = freq_value{1};
end
if hastime
s.times = time_value{1};
end
s.datatype = upper(ext_suffix);
s.parameters = parameters;
ext = [ext_prefix, ext_suffix];
function test_dimord_label()
opt = struct();
opt.samples_label = {'', 'rpt', 'trial'};
opt.nsamples = {1, randint(), 10};
opt.datatype = {'timefreq', 'timelock'};
combis = cosmo_cartprod(opt);
n_combi = numel(combis);
for k = 1:n_combi
c = combis{k};
ds = cosmo_synthetic_dataset('type', c.datatype, ...
'ntargets', 1, ...
'nchunks', c.nsamples, ...
'size', 'big');
assertEqual(size(ds.samples, 1), c.nsamples);
with_samples_label = ~isempty(c.samples_label);
if with_samples_label
ds.a.meeg.samples_label = c.samples_label;
else
ds.a.meeg = rmfield(ds.a.meeg, 'samples_label');
end
data_is_average = c.nsamples == 1 && ~with_samples_label;
switch c.datatype
case 'timefreq'
samples_field = 'powspctrm';
case 'timelock'
if data_is_average
samples_field = 'avg';
else
samples_field = 'trial';
end
otherwise
assert(false);
end
ft = cosmo_map2meeg(ds);
labels = cosmo_strsplit(ft.dimord, '_');
ndim_expected = numel(ds.a.fdim.labels);
if ~data_is_average
ndim_expected = ndim_expected + 1;
end
assertEqual(numel(labels), ndim_expected);
assertEqual(numel(size(ft.(samples_field))), ndim_expected);
end
function test_meeg_source_dataset_pos_dim_inside_fields()
% mapping back and forth should be fine whether or not the
% 'pos', 'dim' and 'inside' fields are there or not
ds_orig = cosmo_synthetic_dataset('type', 'source', 'size', 'huge');
ft_orig = cosmo_map2meeg(ds_orig);
for has_dim = [false, true]
for has_tri = [false, true]
for has_inside = [false, true]
ft = ft_orig;
if has_dim
ft.dim = [2 2 2];
end
n_pos = size(ft.pos, 1);
if has_tri
ft.tri = ceil(rand(5, 3) * max(n_pos));
end
if ~has_inside
ft = rmfield(ft, 'inside');
end
func = @()cosmo_meeg_dataset(ft);
% verify fields are there
ds = func();
assertEqual(has_dim, isfield(ds.a.meeg, 'dim'));
assertEqual(has_tri, isfield(ds.a.meeg, 'tri'));
% map back
ft_back = cosmo_map2meeg(ds);
% cosmo_map2meeg always returns an inside field, so for
% now we remove it if it is present here
if ~has_inside
ft_back = rmfield(ft_back, 'inside');
end
% ensure expected fields are preserved
assertEqual(ft, ft_back);
end
end
end
function x = randint()
x = ceil(rand() * 10 + 5);
function x = randstr()
x = char(rand(1, 10) * 24 + 65);
