function test_suite = test_average_samples
% tests for cosmo_average_samples
%
% # 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_average_samples_
ds=cosmo_synthetic_dataset();
a=cosmo_average_samples(ds);
assertElementsAlmostEqual(sort(a.samples), sort(ds.samples));
assertElementsAlmostEqual(sort(a.samples(:,3)), sort(ds.samples(:,3)));
a=cosmo_average_samples(ds,'ratio',.5);
assertElementsAlmostEqual(sort(a.samples), sort(ds.samples));
assertElementsAlmostEqual(sort(a.samples(:,3)), sort(ds.samples(:,3)));
% check wrong inputs
aet=@(varargin)assertExceptionThrown(@()...
cosmo_average_samples(varargin{:}),'');
aet(ds,'ratio',.1);
aet(ds,'ratio',3);
aet(ds,'ratio',.5,'count',2);
ds.sa.chunks(:)=1;
a=cosmo_average_samples(ds,'ratio',.5);
cosmo_check_dataset(a);
ds=cosmo_slice(ds,3,2);
ns=size(ds.samples,1);
ds.samples=ds.sa.targets*1000+(1:ns)';
a=cosmo_average_samples(ds,'ratio',.5,'nrep',10);
% no mixing of different targets
delta=a.samples/1000-a.sa.targets;
assertTrue(all(.00099<=delta & delta<.05));
assertElementsAlmostEqual(delta*3000,round(delta*3000));
a=cosmo_average_samples(ds,'count',3,'nrep',10);
% no mixing of different targets
delta=a.samples/1000-a.sa.targets;
assertTrue(all(.00099<=delta & delta<.05));
assertElementsAlmostEqual(delta*3000,round(delta*3000));
function test_average_samples_split_by
plural_singular={'targets','targets';...
'chunks','chunks';...
'subjects','subject';...
'modalities','modality';...
};
n_dim=size(plural_singular,1);
combis=cosmo_cartprod(repmat({{true,false}},n_dim,1)');
for k=1:size(combis,1);
combi=cell2mat(combis(k,:));
opt=struct();
opt.seed=0; % truly random data
for j=1:n_dim
count=ceil(rand()*2+1);
opt.(['n' plural_singular{j,1}])=count;
end
ds=cosmo_synthetic_dataset(opt);
values=cell(n_dim,1);
for j=1:n_dim
if combi(j)
values{j}=ds.sa.(plural_singular{j,2});
end
end
values=values(combi);
if any(combi)
[idx,unq_cell]=cosmo_index_unique(values);
else
idx={1:(size(ds.samples,1))};
end
n_avg=numel(idx);
n_features=size(ds.samples,2);
expected_samples=zeros(n_avg,n_features);
for m=1:n_avg
expected_samples(m,:)=mean(ds.samples(idx{m},:),1);
end
result=cosmo_average_samples(ds,...
'split_by',plural_singular(combi,2));
assertEqual(size(result.samples),size(expected_samples));
delta=bsxfun(@minus,result.samples(:,1),expected_samples(:,1)');
mapping=zeros(1,n_avg);
for m=1:n_avg
[mn,mn_idx]=min(abs(delta(m,:)));
assert(mn<1e-5); % deal with rounding
mapping(mn_idx)=m;
end
assertEqual(sort(mapping),1:n_avg);
result_perm=cosmo_slice(result,mapping);
assertElementsAlmostEqual(result_perm.samples,expected_samples);
pos=0;
for j=1:n_dim
if combi(j)
pos=pos+1;
fn=plural_singular{j,2};
assertEqual(unq_cell{pos},result_perm.sa.(fn));
end
end
% check default result
if isequal(plural_singular(combi),{'targets','chunks'});
default_result=cosmo_average_samples(ds);
assertEqual(result,default_result);
end
end
function test_average_samples_split_by_empty()
ds=cosmo_synthetic_dataset('ntargets',ceil(rand()*5+2),...
'nchunks',ceil(rand()*5+2));
result=cosmo_average_samples(ds,'split_by',{});
assertElementsAlmostEqual(result.samples,mean(ds.samples,1));
function test_average_samples_exceptions
aet=@(varargin)assertExceptionThrown(@()...
cosmo_average_samples(varargin{:}),'');
ds=cosmo_synthetic_dataset('nreps',5);
aet([]);
x=struct();
x.samples=randn(4);
aet(x);
% illegal count
aet(ds,'count',6)
aet(ds,'count',[2 2])
aet(ds,'count',3.5);
aet(ds,'count',0);
% illegal ratio
aet(ds,'ratio',1.2)
aet(ds,'ratio',-0.2);
aet(ds,'ratio',[.5 .5]);
% mutually exclusive
aet(ds,'ratio',.5,'count',2);
aet(ds,'repeats',[2 2])
aet(ds,'repeats',-1);
aet(ds,'resamplings',[2 2]);
aet(ds,'resamplings',-1);
aet(ds,'resamplings',1,'repeats',1);
% not existing field
ds_bad=ds;
ds_bad.sa=rmfield(ds_bad.sa,'targets');
aet(ds_bad);
% illegal split-by arguments
aet(ds,'split_by',[]);
aet(ds,'split_by',struct());
aet(ds,'split_by','foo');
aet(ds,'split_by',{1,2});
function test_average_samples_with_repeats
nchunks=ceil(rand()*4+3);
ntargets=ceil(rand()*4+3);
ncombi_max=ceil(rand()*3+4);
max_cyc=5;
ncombi_min=ceil(ncombi_max/2);
ds=cosmo_synthetic_dataset('nchunks',nchunks,...
'ntargets',ntargets,...
'nreps',ncombi_max);
ds.sa=rmfield(ds.sa,'rep');
sp=cosmo_split(ds,{'targets','chunks'});
n_splits=numel(sp);
% select subset of samples, each with at least ncombi_min repeats
combi_count=zeros(nchunks,ntargets);
for k=1:n_splits
if k==1
% ensure at least one with minimum
nkeep=ncombi_min;
else
nkeep=ncombi_min+floor(rand()*(ncombi_max-ncombi_min));
end
ds_k=cosmo_slice(sp{k},1:nkeep);
ds_k.sa.repeats=(1:nkeep)';
combi_count(ds_k.sa.chunks(1),ds_k.sa.targets(1))=nkeep;
sp{k}=ds_k;
end
assert(all(cellfun(@(x)size(x.samples,1),sp)));
ds=cosmo_stack(sp);
[nsamples,nfeatures]=size(ds.samples);
% bit widths for features, chunks, targets, and repeats
bws=[nfeatures,nchunks,ntargets,ceil(log2(max_cyc+1))+ncombi_max];
% encode features, chunks, targets and repeats into single number
dsb=binarize_ds(ds,bws);
% helper function
check_with=@(args,...
count,...
repeats) check_with_helper(dsb,args,count,repeats,...
nchunks,ntargets,...
ncombi_max,combi_count,...
bws);
for repeats=[1,ceil(rand()*ncombi_max)]
for count=[1,ceil(rand()*ncombi_min)];
check_with({'count',count,'repeats',repeats},...
count,repeats);
end
for ratio=[.5,.3+rand()*.7];
count=round(ratio*min(combi_count(:)));
check_with({'ratio',ratio,'repeats',repeats},...
count,repeats);
end
end
for resamplings=[0,1,2+round(rand()*4)]
count=ceil(rand()*ncombi_min);
if resamplings==0
repeats=floor(ncombi_min/count);
args={'count',count};
else
repeats=floor(resamplings*ncombi_min/count);
args={'count',count,'resamplings',resamplings};
end
check_with(args,count,repeats);
end
function check_with_helper(dsb, args, count, repeats,...
nchunks, ntargets, ncombi_max, combi_count, bws)
mu=cosmo_average_samples(dsb,args{:});
[chunks,targets,ids]=unbinarize_ds(mu, bws, count);
nsamples=size(ids,1);
nfeatures=size(dsb.samples,2);
% chunk, target, repeat count
ctr_count=zeros(nchunks,ntargets,ncombi_max);
% keep track of each target and chunk combination
for j=1:nsamples
for k=1:nfeatures
% select same samples for all features
id=ids{j,k};
if k==1
first_id=id;
else
assertEqual(first_id,id);
end
end
% no repeats
id_sorted=sort(id(:));
assert(all(diff(id_sorted)>0));
% count should match
assertEqual(numel(id),count);
ctr_count(chunks(j),targets(j),id)=...
ctr_count(chunks(j),targets(j),id)+1;
end
% ensure each sample selected about equally often
[nchunks,ntargets]=size(combi_count);
for k=1:nchunks
for j=1:ntargets
c=squeeze(ctr_count(k,j,:));
pre=c(1:combi_count(k,j));
assert(max(pre)-min(pre)<=1);
post=c((combi_count(k,j)+1):end);
assert(all(post==0));
end
end
% check each target and chunk combination was used the correct number
% of times to form the average
ct_count=sum(ctr_count,3);
expected_ct_count=count*repeats*ones(nchunks,ntargets);
assert(isequal(ct_count, expected_ct_count));
function [chunks,targets,ids]=unbinarize_ds(ds, bws, counts)
[nsamples, nfeatures]=size(ds.samples);
ids=cell(nsamples,nfeatures);
chunks=zeros(nsamples,1);
targets=zeros(nsamples,1);
for k=1:nsamples
for j=1:nfeatures
% Decode repeats; multiple repeats can be present.
% As there can be multiple repeats, the averaging is undone
% and then each bit represents just one repeat
v_id=quick_dec2bin(mod(ds.samples(k,j)*counts,...
2^bws(end)),...
bws(end));
ids{k,j}=bws(end)-find(v_id)+1;
% decode chunks, targets, ids
v=decode(floor(ds.samples(k,j)/2^bws(end)),bws(1:(end-1)));
assertEqual(log2(v(1))+1,j);
c=log2(v(2))+1;
t=log2(v(3))+1;
if j==1
chunks(k)=c;
targets(k)=t;
else
assertEqual(c,chunks(k));
assertEqual(t,targets(k));
end
end
end
function bds=binarize_ds(ds, bws)
bds=ds;
[nsamples,nfeatures]=size(ds.samples);
for k=1:nsamples
sa=cosmo_slice(ds.sa,k,1,'struct');
for j=1:nfeatures
vs=[j, sa.chunks sa.targets sa.repeats];
bds.samples(k,j)=encode(vs,bws);
end
end
function p=encode(vs, bws)
% encode several decimal numbers in a single one, through
% encode([X1 ... Xn]) = bin2dec([dec2bin(X1) ... dec2bin(Xn)])
% where bws contains the bit width for each number
n=numel(bws);
assert(numel(vs)==n);
bs=cell(1,n);
for k=1:n
bw=bws(k);
bs{k}=zeros(1,bw);
bs{k}(bw-vs(k)+1)=1;
end
p=quick_bin2dec(cat(2,bs{:}));
function vs=decode(p, bws)
% encode single decimal numbers in multiple ones, through
% decode(P) = [bin2dec(PB1) ... bin2dec(PBn)]
% with PBi the binary representation part of P for each binary
% representation part
arr=quick_dec2bin(p,sum(bws));
c=0;
n=numel(bws);
vs=zeros(1,n);
for k=1:n
offset=bws(k);
vs(k)=quick_bin2dec(arr(c+(1:offset)));
c=c+offset;
end
function arr=quick_dec2bin(x,bw)
% converts decimal number x to array with length bw and all
% values in 0 and 1
assert(round(x)==x);
arr=zeros(1,bw);
xbs=dec2bin(x);
arr(bw-numel(xbs)+1:end)=(xbs=='1');
return
function x=quick_bin2dec(arr)
% convert binary array to decimal number
x=sum(2.^((numel(arr)-1):-1:0).*arr);
