function test_suite = test_cluster_neighborhood
% tests for cosmo_cluster_neighborhood
%
% # 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_fmri_cluster_neighborhood
ds = cosmo_synthetic_dataset('type', 'fmri', 'size', 'normal');
nf = size(ds.samples, 2);
imsk = find(rand(1, nf) > .8);
rp = randperm(numel(imsk));
ds = cosmo_slice(ds, [imsk(rp) imsk(rp(end:-1:1))], 2);
nh1 = cosmo_cluster_neighborhood(ds, 'progress', false);
nh2 = cosmo_cluster_neighborhood(ds, 'progress', false, 'fmri', 1);
nh_sph = cosmo_spherical_neighborhood(ds, 'progress', false, ...
'radius', 2.5);
nh3 = cosmo_cluster_neighborhood(ds, 'progress', false, 'fmri', nh_sph);
ds.fa.sizes = ones(1, size(ds.samples, 2));
assertEqual(nh1.a, ds.a);
assertEqual(nh1.fa, ds.fa);
nf = size(ds.samples, 2);
rp = randperm(nf);
nfeatures_test = nf;
ijk = [ds.fa.i; ds.fa.j; ds.fa.k];
feature_ids = rp(1:nfeatures_test);
for j = 1:nfeatures_test
feature_id = feature_ids(j);
delta = sqrt(sum(bsxfun(@minus, ijk(:, feature_id), ijk).^2, 1));
msk1 = delta < sqrt(3) + .001;
assertEqual(sort(nh1.neighbors{feature_id}), find(msk1));
msk2 = delta < sqrt(1) + .001;
assertEqual(sort(nh2.neighbors{feature_id}), find(msk2));
msk3 = delta < 2.5;
assertEqual(sort(nh3.neighbors{feature_id}), find(msk3));
end
function test_tiny_fmri_neighborhood
ds = cosmo_synthetic_dataset('type', 'fmri', 'size', 'normal');
ds = cosmo_slice(ds, [22 22], 2);
% should pass
nh = cosmo_cluster_neighborhood(ds, 'progress', false);
function test_meeg_cluster_neighborhood
if cosmo_skip_test_if_no_external('fieldtrip')
return
end
ds = cosmo_synthetic_dataset('type', 'timefreq', 'size', 'normal');
nf = size(ds.samples, 2);
imsk = find(rand(1, nf) > .4);
rp = randperm(numel(imsk));
ds = cosmo_slice(ds, [imsk(rp) imsk(rp(end:-1:1))], 2);
ds = cosmo_dim_prune(ds);
n = numel(ds.a.fdim.values{1});
ds.a.fdim.values{1} = ds.a.fdim.values{1}(randperm(n));
nf = size(ds.samples, 2);
ds.fa.sizes = ones(1, size(ds.samples, 2));
chan_nbrhood = cosmo_meeg_chan_neighborhood(ds, 'delaunay', true, ...
'chantype', 'all', ...
'label', 'dataset');
assertEqual(ds.a.fdim.values(1), chan_nbrhood.a.fdim.values(1));
ncombi = 7; % all 7 possibilities
test_range = ceil(rand(1, ncombi) * 7);
nfeatures_test = 3;
for i = test_range
use_chan = i <= 4;
use_freq = mod(i, 2) == 1;
use_time = mod(ceil(i / 2), 2) == 1;
use_msk = [use_chan, use_freq, use_time];
labels = {'chan', 'freq', 'time'};
ndim = numel(labels);
args = struct();
for k = 1:ndim
if ~use_msk(k)
label = labels{k};
args.(label) = false;
end
end
cl_nbrhood = cosmo_cluster_neighborhood(ds, args, 'progress', false);
assertEqual(cl_nbrhood.fa, ds.fa);
assertEqual(cl_nbrhood.a, ds.a);
rp = randperm(nf);
for k = 1:nfeatures_test
feature_id = rp(k);
counter = zeros(1, nf);
for j = 1:ndim
label = labels{j};
fa = ds.fa.(label);
if use_msk(j)
if j == 1
% channel
ids = chan_nbrhood.neighbors{fa(feature_id)};
else
% anything else
ids = find(abs(fa - fa(feature_id)) <= 1.5);
end
else
ids = find(fa == fa(feature_id));
end
counter(ids) = counter(ids) + 1;
end
nbrs = find(counter == ndim);
assertEqual(nbrs, cl_nbrhood.neighbors{feature_id});
end
end
function test_tiny_meeg_cluster_neighborhood
if cosmo_skip_test_if_no_external('fieldtrip')
return
end
ds = cosmo_synthetic_dataset('type', 'meeg');
nh = cosmo_cluster_neighborhood(ds, 'progress', false);
ds.fa.sizes = ones(1, 6);
assertEqual(ds.fa, nh.fa);
assertEqual(ds.a, nh.a);
half_neighbors = {[1 4]; [2 3 5 6]; [2 3 5 6]};
assertEqual(nh.neighbors, repmat(half_neighbors, 2, 1));
% transpose of fdim elements should be fine
ds2 = ds;
ds2.a.fdim.values = {ds2.a.fdim.values{1}', ds2.a.fdim.values{2}'};
nh2 = cosmo_cluster_neighborhood(ds2, 'progress', false);
cosmo_check_neighborhood(nh2, ds);
cosmo_check_neighborhood(nh, ds);
assertEqual(nh.a, nh2.a);
assertEqual(nh.fa, nh2.fa);
assertEqual(nh.neighbors, nh2.neighbors);
% different fdim elements i not ok
ds3 = ds2;
ds3.a.fdim.values{1} = ds3.a.fdim.values{1}(end:-1:1);
function test_cluster_neighborhood_transpose
opt = struct();
opt.progress = false;
ds = cosmo_synthetic_dataset('type', 'timefreq', 'size', 'normal');
ds = cosmo_dim_remove(ds, 'chan');
nh = cosmo_cluster_neighborhood(ds, opt);
cp = cosmo_cartprod(repmat({[false, true]}, 4, 1));
n = size(cp, 1);
for k = 1:n
t_label = cp{k, 1};
t_value = cp{k, 2};
t_elem1 = cp{k, 3};
t_elem2 = cp{k, 4};
ds2 = ds;
if t_label
ds2.a.fdim.labels = ds2.a.fdim.labels';
end
if t_value
ds2.a.fdim.values = ds2.a.fdim.values';
end
if t_elem1
ds2.a.fdim.values{1} = ds2.a.fdim.values{1}';
end
if t_elem2
ds2.a.fdim.values{2} = ds2.a.fdim.values{2}';
end
nh2 = cosmo_cluster_neighborhood(ds2, opt);
assertEqual(nh2.a, nh.a);
assertEqual(nh2.fa, nh.fa);
assertEqual(nh2.neighbors, nh.neighbors);
end
function test_cluster_neighborhood_surface
if cosmo_skip_test_if_no_external('surfing')
return
end
ds = cosmo_synthetic_dataset('type', 'surface'); % ,'size','normal');
vertices = [0 0 0 1 1 1
1 2 3 1 2 3
0 0 0 0 0 0]';
faces = [3 2 3 2
2 1 5 4
5 4 6 5]';
opt = struct();
opt.progress = false;
opt.vertices = vertices;
opt.faces = faces;
nh1 = cosmo_cluster_neighborhood(ds, opt);
assertEqual(nh1.neighbors, { [1 2 4]
[1 2 3 4 5]
[2 3 5 6]
[1 2 4 5]
[2 3 4 5 6]
[3 5 6] });
assertEqual(ds.a, nh1.a);
ds.fa.sizes = [1 3 2 2 3 1] / 6;
ds.fa.radius = sqrt([1 2 2 2 2 1]);
ds.fa.area = [6 11 9 9 11 6] / 6;
assertEqual(ds.fa, nh1.fa);
opt.direct = false;
nh2 = cosmo_cluster_neighborhood(ds, opt);
assertEqual(nh2.neighbors, num2cell((1:6)'));
assertEqual(ds.a, nh2.a);
ds.fa.radius(:) = 0;
ds.fa.area = [1 3 2 2 3 1] / 6;
assertEqual(ds.fa, nh2.fa);
% take a subset of features
rp = randperm(6);
nsel = 4;
sel = rp(1:nsel);
ds3 = cosmo_slice(ds, sel, 2);
% re-order node ids
rp2 = randperm(6);
ds.a.fdim.values{1} = ds.a.fdim.values{1}(rp2);
opt.direct = true;
nh3 = cosmo_cluster_neighborhood(ds3, opt);
assertEqual(nh3.a, ds3.a);
assertEqual(nh3.fa.node_indices, ds3.fa.node_indices);
node_ids = ds3.a.fdim.values{1}(ds3.fa.node_indices);
for k = 1:nsel
nbs = nh3.neighbors{k}(:);
node_id = node_ids(k);
nb1 = intersect(nh1.neighbors{node_id}, sel);
assertEqual(sort(nb1'), sort(node_ids(nbs))');
end
function test_cluster_neighborhood_source
ds = cosmo_synthetic_dataset('size', 'normal', 'type', 'source');
nf = size(ds.samples, 2);
[unused, idxs] = sort(cosmo_rand(1, nf * 3));
rps = mod(idxs - 1, nf) + 1;
rp = rps(round(nf / 2) + (1:(2 * nf)));
ds = cosmo_slice(ds, rp, 2);
[unused, rp] = sort(cosmo_rand(1, nf));
rp = rp(1:4);
grid_spacing = 10;
ds_pos = ds.a.fdim.values{1}(:, ds.fa.pos) / grid_spacing;
for connectivity = 0:3
if connectivity == 0
radius = sqrt(3) + .001;
args = {};
else
args = {'source', connectivity};
radius = sqrt(connectivity) + .001;
end
nh = cosmo_cluster_neighborhood(ds, 'progress', false, args);
nh_pos = nh.a.fdim.values{1}(:, nh.fa.pos) / grid_spacing;
for r = rp
idxs = nh.neighbors{r};
d = sum(bsxfun(@minus, nh_pos(:, r), ds_pos).^2, 1).^.5;
d_inside = d(idxs);
outside_mask = true(size(d));
outside_mask(d <= radius) = false;
d_outside = d(outside_mask);
assert(all(d_inside <= radius));
assert(all(d_outside > radius));
end
end
function test_cluster_neighborhood_source_mom
ds = cosmo_synthetic_dataset('size', 'normal', 'type', 'source', ...
'data_field', 'mom');
nf = size(ds.samples, 2);
[unused, idxs] = sort(cosmo_rand(1, nf * 3));
rps = mod(idxs - 1, nf) + 1;
rp = rps(round(nf / 2) + (1:(2 * nf)));
ds = cosmo_slice(ds, rp, 2);
grid_spacing = 10;
for connectivity = 0:3
if connectivity == 0
radius = sqrt(3) + .001;
args = {};
else
args = {'source', connectivity};
radius = sqrt(connectivity) + .001;
end
for has_3d_mom = [false true]
ds2 = ds;
if has_3d_mom
assertExceptionThrown(@() ...
cosmo_cluster_neighborhood(ds, ...
'progress', false, args), '');
continue
end
keep_dim = ceil(rand() * 3);
keep_msk = ds2.fa.mom == keep_dim;
ds2 = cosmo_slice(ds, keep_msk, 2);
ds2.fa.mom(:) = 1;
ds2.a.fdim.values{2} = ds2.a.fdim.values{2}(keep_dim);
ds2_pos = ds2.a.fdim.values{1}(:, ds2.fa.pos) / grid_spacing;
nf = size(ds2.samples, 2);
[unused, rp] = sort(cosmo_rand(1, nf));
rp = rp(1:4);
nh = cosmo_cluster_neighborhood(ds2, 'progress', false, args);
nh_pos = nh.a.fdim.values{1}(:, nh.fa.pos) / grid_spacing;
for r = rp
idxs = nh.neighbors{r};
d = sum(bsxfun(@minus, nh_pos(:, r), ...
ds2_pos).^2, 1).^.5;
d_inside = d(idxs);
outside_mask = true(size(d));
outside_mask(d <= radius) = false;
d_outside = d(outside_mask);
assert(all(d_inside <= radius));
assert(all(d_outside > radius));
end
end
end
function test_cluster_neighborhood_fmri_time
for fa_n_rep = [1, 4]
n_time = ceil(rand() * 3 + 4);
TR = 2;
time_values = TR * (1:n_time);
ds_cell = cell(n_time, fa_n_rep);
for fa_t = 1:n_time
ds = cosmo_synthetic_dataset('size', 'normal', 'seed', 0);
n_features = size(ds.samples, 2);
time_idx = ones(1, n_features) * fa_t;
ds = cosmo_dim_insert(ds, 2, 4, {'time'}, {time_values}, {time_idx});
ds_cell(fa_t, :) = repmat({ds}, 1, fa_n_rep);
end
% select subset
ds_full = cosmo_stack(ds_cell, 2);
n_features = size(ds_full.samples, 2);
rp = randperm(n_features);
n_keep = ceil(n_features / 2);
keep_idxs = rp(1:n_keep);
ds = cosmo_slice(ds_full, keep_idxs, 2);
fa_ijk = [ds.fa.i; ds.fa.j; ds.fa.k];
fa_t = ds.fa.time;
for join_time = -1:1
if join_time == -1
args = {};
else
args = {'time', join_time == 1};
end
nh = cosmo_cluster_neighborhood(ds, args, 'progress', false);
assertEqual(nh.origin.fa, ds.fa);
assertEqual(nh.origin.a, ds.a);
assertEqual(nh.fa.sizes, ones(1, n_keep));
assertEqual(nh.a, ds.a);
time_radius = abs(join_time);
for k = 1:numel(nh.neighbors)
idx = nh.neighbors{k};
delta_xyz = sqrt(sum(bsxfun(@minus, fa_ijk(:, k), fa_ijk).^2, 1));
delta_t = abs(fa_t(k) - fa_t);
msk = delta_xyz <= 1.9 & delta_t <= time_radius;
assertEqual(sort(idx), find(msk));
end
end
end
function test_meeg_cluster_neighborhood_unknown_eeg_channels
if cosmo_skip_test_if_no_external('fieldtrip')
return
end
ds_orig = cosmo_synthetic_dataset('type', 'meeg', ...
'sens', 'eeg1005', ...
'size', 'normal');
for with_unknown_channel = [false, true]
ds = ds_orig;
if with_unknown_channel
nchan = max(ds.fa.chan);
idx = ceil(rand() * nchan);
ds.a.fdim.values{1}{idx} = 'foo';
end
nh = cosmo_cluster_neighborhood(ds, 'progress', false);
if with_unknown_channel
empty_msk = ds.fa.chan == idx;
assert(any(empty_msk));
else
empty_msk = false(size(ds.fa.chan));
end
count = cellfun(@numel, nh.neighbors);
assert(all(count(empty_msk) == 0));
assert(all(count(~empty_msk) > 0));
end
function test_cluster_neighborhood_exceptions
ds = cosmo_synthetic_dataset();
aet = @(varargin)assertExceptionThrown( ...
@()cosmo_cluster_neighborhood(varargin{:}, ...
'progress', false), '');
aet(ds, 'foo');
aet(ds, 'fmri', -1);
aet(ds, 'fmri', true);
aet(ds, 'fmri', [1 1]);
ds.a.fdim.labels{2} = 'foo';
ds.fa.foo = ds.fa.j;
aet(ds);
ds2 = cosmo_synthetic_dataset('type', 'meeg');
ds2.a.fdim.labels{1} = 'freq';
ds2.fa.freq = ones(1, 6);
aet(ds2, 'freq');
aet(ds2, 'freq', 2);
aet(ds2, 'freq', NaN);
aet(ds2, 'freq', [true true]);
