function test_suite = test_spherical_neighborhood
% tests for cosmo_spherical_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_simple_neighborhood
ds=cosmo_synthetic_dataset();
nh1=cosmo_spherical_neighborhood(ds,'radius',0,'progress',false);
mp=cosmo_align({ds.fa.i, ds.fa.j, ds.fa.k},...
{nh1.fa.i, nh1.fa.j, nh1.fa.k});
ds_al=cosmo_slice(ds,mp,2);
assertEqual(nh1.a,ds.a);
assertEqual(nh1.fa.i,ds_al.fa.i);
assertEqual(nh1.fa.j,ds_al.fa.j);
assertEqual(nh1.fa.k,ds_al.fa.k);
assertEqual(nh1.fa.nvoxels,ones(1,6));
assertEqual(nh1.fa.radius,zeros(1,6));
assertEqual(nh1.fa.center_ids,mp);
assertEqual(nh1.neighbors,num2cell(mp'));
assertEqual(nh1.origin.fa,ds_al.fa);
assertEqual(nh1.origin.a.fdim,ds_al.a.fdim);
nh2=cosmo_spherical_neighborhood(ds,'radius',1.5,'progress',false);
mp=cosmo_align({ds.fa.i, ds.fa.j, ds.fa.k},...
{nh1.fa.i, nh1.fa.j, nh1.fa.k});
ds_al=cosmo_slice(ds,mp,2);
assertEqual(nh2.a,ds.a);
assertEqual(nh2.fa.i,ds_al.fa.i);
assertEqual(nh2.fa.j,ds_al.fa.j);
assertEqual(nh2.fa.k,ds_al.fa.k);
nvoxels=[4 6 4 4 6 4];
assertEqual(nh2.fa.nvoxels,nvoxels(mp));
assertEqual(nh2.fa.radius,ones(1,6)*1.5);
assertEqual(nh2.fa.center_ids,mp);
assertEqual(numel(nh2.neighbors),6);
nh={ [ 1 4 2 5 ];...
[ 2 1 5 3 4 6 ];...
[ 3 2 6 5 ];...
[ 4 1 5 2 ];...
[ 5 4 2 6 1 3 ];...
[ 6 5 3 2 ] };
for k=1:6
assertEqual(nh2.neighbors{k},nh{mp(k)});
end
nh3=cosmo_spherical_neighborhood(ds,'count',4,'progress',false);
mp=cosmo_align({ds.fa.i, ds.fa.j, ds.fa.k},...
{nh1.fa.i, nh1.fa.j, nh1.fa.k});
ds_al=cosmo_slice(ds,mp,2);
assertEqual(nh3.a,ds.a);
assertEqual(nh3.fa.i,ds_al.fa.i);
assertEqual(nh3.fa.j,ds_al.fa.j);
assertEqual(nh3.fa.k,ds_al.fa.k);
assertEqual(nh3.fa.nvoxels,[4 4 4 4 4 4]);
radii=[sqrt(2) 1 sqrt(2) sqrt(2) 1 sqrt(2)];
assertElementsAlmostEqual(nh3.fa.radius,...
radii(mp),...
'relative',1e-3);
assertEqual(nh3.fa.center_ids,mp);
assertEqual(numel(nh3.neighbors),6);
nh={ [ 1 4 2 5 ];...
[ 2 1 5 3 ];...
[ 3 2 6 5 ];...
[ 4 1 5 2 ];...
[ 5 4 2 6 ];...
[ 6 5 3 2 ] };
for k=1:6
assertEqual(nh3.neighbors{k},nh{mp(k)});
end
function test_exceptions
ds=cosmo_synthetic_dataset();
aet=@(x)assertExceptionThrown(@()...
cosmo_spherical_neighborhood(x{:}),'');
aet({ds});
aet({ds,'foo'});
aet({ds,'foo',1});
aet({ds,'radius',[1 2]});
aet({ds,'count',[1 2]});
aet({ds,'radius',-1});
aet({ds,'count',-1});
aet({ds,'radius',1,'count',1});
aet({ds,'count',7});
aet({'foo','count',7});
function test_sparse_dataset
nfeatures_test=3;
ds=cosmo_synthetic_dataset('size','big');
nf=size(ds.samples,2);
rp=randperm(nf);
nf_full=round(nf*.4);
ids=rp(1:nf_full);
ds_sel=cosmo_slice(ds,ids,2);
ids_full=repmat(ids,1,2);
ds_full=cosmo_slice(ds,ids_full,2);
radius=2+rand();
nh4=cosmo_spherical_neighborhood(ds_full,'radius',radius,...
'progress',false);
assertEqual(numel(nh4.neighbors),numel(ids));
mp=cosmo_align({ds_sel.fa.i, ds_sel.fa.j, ds_sel.fa.k},...
{nh4.fa.i, nh4.fa.j, nh4.fa.k});
ds_al=cosmo_slice(ds_sel,mp,2);
assertEqual(nh4.a,ds.a);
assertEqual(nh4.fa.i,ds_al.fa.i);
assertEqual(nh4.fa.j,ds_al.fa.j);
assertEqual(nh4.fa.k,ds_al.fa.k);
assertEqual(nh4.fa.center_ids,mp);
rp=randperm(size(ds_al.samples,2));
center_ids=rp(1:nfeatures_test);
ijk=[ds_al.fa.i; ds_al.fa.j; ds_al.fa.k];
ijk_full=[ds_full.fa.i; ds_full.fa.j; ds_full.fa.k];
for center_id=center_ids
ijk_center=ijk(:,center_id);
delta=sum(bsxfun(@minus,ijk_center,ijk_full).^2,1).^.5;
nbr_ids=find(delta<=radius);
assertEqual(nbr_ids,sort(nh4.neighbors{center_id}));
end
function test_with_freq_dimension_dataset
ds=cosmo_synthetic_dataset('size','big');
nfeatures=size(ds.samples,2);
freqs=[2 4 6];
nfreqs=numel(freqs);
ds_cell=cell(nfreqs,1);
for k=1:nfreqs
ds_freq=ds;
ds_freq.a.fdim.labels=[{'freq'};ds_freq.a.fdim.labels];
ds_freq.a.fdim.values=[{freqs};ds_freq.a.fdim.values];
ds_freq.fa.freq=ones(1,nfeatures)*k;
ds_cell{k}=ds_freq;
end
ds_full=cosmo_stack(ds_cell,2);
rp=randperm(nfeatures*nfreqs);
id_full=rp(1:((nfreqs-1)*nfeatures));
ds_full=cosmo_slice(ds_full,id_full,2);
ijk_full=[ds_full.fa.i; ds_full.fa.j; ds_full.fa.k];
[ijk_idxs,ijk_unq]=cosmo_index_unique(ijk_full');
ijk_idx=cellfun(@(x)x(1),ijk_idxs);
ds_al=cosmo_slice(ds_full,ijk_idx,2);
radius=5+rand()*3;
nh=cosmo_spherical_neighborhood(ds_full,'radius',radius,...
'progress',false);
assertEqual(nh.fa.i,ds_al.fa.i);
assertEqual(nh.fa.j,ds_al.fa.j);
assertEqual(nh.fa.k,ds_al.fa.k);
assertEqual(nh.a.fdim.labels,ds_full.a.fdim.labels(2:4));
assertEqual(nh.a.fdim.values,ds_full.a.fdim.values(2:4));
assertEqual(numel(nh.neighbors),numel(ijk_idxs));
ijk_full=[ds_full.fa.i; ds_full.fa.j; ds_full.fa.k];
ijk_al=[ds_al.fa.i; ds_al.fa.j; ds_al.fa.k];
% test match for euclidean distance
rp=randperm(numel(ijk_idxs));
rp=rp(1:10);
for r=rp
nbrs=nh.neighbors{r};
ijk_center=ijk_al(:,r);
delta=sum(bsxfun(@minus,ijk_center,ijk_full).^2,1).^.5;
assertEqual(find(delta<=radius), sort(nbrs));
end
% test with transposed dimension values
ds2_full=ds_full;
ds2_full.a.fdim.values=cellfun(@transpose,...
ds2_full.a.fdim.values,...
'UniformOutput',false)';
nh2=cosmo_spherical_neighborhood(ds2_full,'radius',radius,...
'progress',false);
nh2.origin.a.fdim.values=cellfun(@transpose,...
nh2.origin.a.fdim.values,...
'UniformOutput',false)';
assertEqual(nh,nh2);
assertFalse(isfield(nh.fa,'inside'));
function test_meeg_source_dataset
ds=cosmo_synthetic_dataset('type','source','size','normal');
nf=size(ds.samples,2);
[unused,idxs]=sort(cosmo_rand(1,nf*4,'seed',1));
rps=mod(idxs-1,nf)+1;
id_full=rps(round(nf/2)+(1:(3*nf)));
ds_full=cosmo_slice(ds,id_full,2);
pos_full=ds_full.fa.pos;
pos_idxs=cosmo_index_unique(pos_full');
pos_idx=cellfun(@(x)x(1),pos_idxs);
ds_unq=cosmo_slice(ds_full,pos_idx,2);
radius=1.2+.2*rand();
nh=cosmo_spherical_neighborhood(ds_full,'radius',radius,...
'progress',false);
mp=cosmo_align(ds_unq.fa.pos',nh.fa.pos');
ds_al=cosmo_slice(ds_unq,mp,2);
assertEqual(nh.fa.pos,ds_al.fa.pos);
assertEqual(nh.a.fdim.labels{1},ds.a.fdim.labels{1});
assertEqual(nh.a.fdim.values{1},ds.a.fdim.values{1});
assertEqual(numel(nh.neighbors),numel(pos_idx));
count=ceil(4/3*pi*(radius)^3 * .5);
nh2=cosmo_spherical_neighborhood(ds_full,'count',count,...
'progress',false);
assertEqual(nh2.fa.pos,ds_al.fa.pos);
assertEqual(nh2.a.fdim.labels{1},ds.a.fdim.labels{1});
assertEqual(nh2.a.fdim.values{1},ds.a.fdim.values{1});
assertEqual(numel(nh2.neighbors),numel(pos_idx));
pos_al=ds_al.a.fdim.values{1}(:,ds_al.fa.pos);
pos_full=nh.a.fdim.values{1}(:,ds_full.fa.pos);
voxel_size=10;
for r=1:numel(pos_idx)
d=sum(bsxfun(@minus,pos_full,pos_al(:,r)).^2,1).^.5;
idxs=find(d<=(radius*voxel_size));
assertEqual(sort(nh.neighbors{r}),sort(idxs));
end
[p,q]=cosmo_overlap(nh.neighbors,nh2.neighbors);
dp=diag(p);
dq=diag(q);
assertTrue(mean(dp)>.8);
assertTrue(mean(dq)>.3);
function test_small_meeg_source_dataset_without_inside_field
ds=cosmo_synthetic_dataset('type','source','size','small');
voxel_size=10;
idx=cosmo_index_unique(ds.fa.pos');
n_pos=numel(idx);
ds_pos=ds.a.fdim.values{1}(:,ds.fa.pos);
for radius=0:.4:2;
nh=cosmo_spherical_neighborhood(ds,...
'radius',radius,...
'progress',false);
center_pos=nh.a.fdim.values{1}(:,nh.fa.pos);
assertEqual(numel(nh.neighbors),n_pos);
for k=1:n_pos
sel_idx=nh.neighbors{k};
delta=bsxfun(@minus,center_pos(:,k),...
ds_pos);
distance=sqrt(sum((delta/voxel_size).^2,1));
expected_sel_idx=find(distance<=radius);
assertEqual(sort(sel_idx(:)),sort(expected_sel_idx(:)));
end
end
function test_small_meeg_source_dataset_with_inside_field
ds=cosmo_synthetic_dataset('type','source','size','small');
voxel_size=10;
nf=size(ds.samples,2);
idx=cosmo_index_unique(ds.fa.pos');
n_pos=numel(idx);
for keep_ratio=.4:.3:1;
keep=randperm(n_pos);
n_keep=round(n_pos*keep_ratio);
keep=keep(1:n_keep);
ds.fa.inside=false(1,nf);
for k=1:numel(keep)
ds.fa.inside(idx{keep(k)})=true;
end
ds_pos=ds.a.fdim.values{1}(:,ds.fa.pos);
for radius=0:.4:2;
nh=cosmo_spherical_neighborhood(ds,...
'radius',radius,...
'progress',false);
center_pos=nh.a.fdim.values{1}(:,nh.fa.pos);
assertEqual(numel(nh.neighbors),n_keep);
for k=1:n_keep
sel_idx=nh.neighbors{k};
delta=bsxfun(@minus,center_pos(:,k),...
ds_pos);
distance=sqrt(sum((delta/voxel_size).^2,1));
expected_sel_idx=find(distance<=radius & ds.fa.inside);
assertEqual(sort(sel_idx(:)),sort(expected_sel_idx(:)));
end
end
end
function test_meeg_source_illegal_inside()
ds=cosmo_synthetic_dataset('type','source','size','small');
n_features=size(ds.samples,1);
illegal_values={min(max(ds.samples(1,:),0),1),...
repmat({'foo'},1,n_features)};
for k=1:numel(illegal_values)
ds.fa.inside=illegal_values{k};
assertExceptionThrown(@()...
cosmo_spherical_neighborhood(ds,'radius',1),'');
end
function test_meeg_missing_dimension_label()
ds=cosmo_synthetic_dataset();
ds=cosmo_dim_remove(ds,'i');
assertExceptionThrown(@()...
cosmo_spherical_neighborhood(ds,'radius',1),'');
function test_meeg_wrong_dimension_order()
ds=cosmo_synthetic_dataset();
ds.fa.j(:)=1;
ds.fa.k(:)=1;
ds.a.fdim.labels([2,3])=ds.a.fdim.labels([3,2]);
assertExceptionThrown(@()...
cosmo_spherical_neighborhood(ds,'radius',1),'');
function test_meeg_source_sdim_time
ds=cosmo_synthetic_dataset('type','source','size','big');
ds_time=cosmo_dim_transpose(ds,'time',1);
nh=cosmo_spherical_neighborhood(ds_time,'radius',0,'progress',false);
assert(~isfield(nh.a,'sdim'));
assert(~isfield(nh,'sa'));
function test_fmri_fixed_number_of_features()
ds=cosmo_synthetic_dataset('size','normal');
nf=size(ds.samples,2);
[unused,idxs]=sort(cosmo_rand(1,nf*3,'seed',1));
rps=mod(idxs-1,nf)+1;
id_full=rps(round(nf/2)+(1:(2*nf)));
ds_full=cosmo_slice(ds,id_full,2);
id_idxs=cosmo_index_unique(id_full');
id_idx=cellfun(@(x)x(1),id_idxs);
ds_sel=cosmo_slice(ds_full,id_idx,2);
count=20;
nh=cosmo_spherical_neighborhood(ds_full,'count',count,...
'progress',false);
ijk_sel={ds_sel.fa.i; ds_sel.fa.j; ds_sel.fa.k};
ijk_nh={nh.fa.i; nh.fa.j; nh.fa.k};
mp=cosmo_align(ijk_sel,ijk_nh);
ds_al=cosmo_slice(ds_sel,mp,2);
assertEqual(nh.fa.i,ds_al.fa.i);
assertEqual(nh.fa.j,ds_al.fa.j);
assertEqual(nh.fa.k,ds_al.fa.k);
pos_nh=[nh.fa.i; nh.fa.j; nh.fa.k];
pos_full=[ds_full.fa.i;ds_full.fa.j;ds_full.fa.k];
rp=randperm(numel(id_idx));
rp=rp(1:10);
for r=rp
d=sum(bsxfun(@minus,pos_nh(:,r),pos_full).^2,1).^.5;
idxs=nh.neighbors{r};
d_inside=d(idxs);
d_outside=d(setdiff(1:nf,idxs));
assert(max(d_inside)<=min(d_outside));
assertElementsAlmostEqual(max(d_inside),nh.fa.radius(r),...
'absolute',1e-4);
end
