function test_suite=test_phase_stat
% tests for test_phase_stat
%
% # 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 r=randint()
r=ceil(rand()*10+10);
function test_phase_stat_basic_balanced_trials
r=randint();
% test with both balanced and unbalanced number of trials
helper_test_phase_stat_with_trial_counts(r,r);
function test_phase_stat_basic_unbalanced_trials
names={'pos','pop','pbi'};
for k=1:numel(names)
name=names{k};
for delta=[-1,0,1]
r=randint();
ds=generate_phase_dataset(r,r+delta);
opt=struct();
opt.output=name;
func=@()cosmo_phase_stat(ds,opt);
is_balanced=delta==0;
if is_balanced
% should be ok
func();
else
assertExceptionThrown(func,'');
end
end
end
function helper_test_phase_stat_with_trial_counts(ntrials1,ntrials2)
ds=generate_phase_dataset(ntrials1,ntrials2);
names={'pos','pop','pbi'};
for k=1:numel(names)
name=names{k};
helper_test_phase_stat_with_name(ds,name);
end
function helper_test_phase_stat_with_name(ds,stat_name)
opt=struct();
opt.output=stat_name;
% compute result
result=cosmo_phase_stat(ds,opt);
% compute expected result
samples=ds.samples;
t1=find(ds.sa.targets==1);
t2=find(ds.sa.targets==2);
assert(numel(t1)==numel(t2));
itc1=compute_itc(samples(t1,:));
itc2=compute_itc(samples(t2,:));
itc_all=compute_itc(samples);
expected_samples=compute_phase_stat(stat_name,itc1,itc2,itc_all);
% verify output matches expected output
assertElementsAlmostEqual(result.samples,expected_samples);
assertEqual(ds.a,result.a);
assertEqual(ds.fa,result.fa);
assertEqual(struct(),result.sa);
function test_phase_stat_with_signal()
names={'pos','pop','pbi'};
for k=1:numel(names)
name=names{k};
helper_test_phase_stat_with_signal_with_name(name);
end
function helper_test_phase_stat_with_signal_with_name(name)
ds=cosmo_synthetic_dataset('ntargets',2,'nchunks',50);
nsamples=size(ds.samples,1);
ds.sa.chunks(:)=1:nsamples;
ds=cosmo_slice(ds,1,2);
% use small phase angle differences
% generally increase the distance, which should lead to an increase in
% PBI, POS and POP
sd=pi/100;
signals=0:10;
nsignals=numel(signals);
result=zeros(nsignals,1);
for k=1:nsignals
for target=[1,2]
msk=ds.sa.targets==target;
r=rand(sum(msk),size(ds.samples,2));
rng=r;
if target==1
% add increasing difference between two classes
rng=rng+signals(k);
end
angle=2*pi*rng*sd;
x=exp(1i*angle);
ds.samples(msk,:)=x;
end
s=cosmo_phase_stat(ds,'output',name);
result(k)=s.samples;
end
assert(cosmo_corr(result,signals')>.5);
function s=compute_phase_stat(stat_name,itc1,itc2,itc_all)
switch stat_name
case 'pbi'
s=(itc1-itc_all).*(itc2-itc_all);
case 'pop'
s=(itc1.*itc2)-itc_all.^2;
case 'pos'
s=(itc1+itc2)-2*itc_all;
otherwise
assert(false);
end
function itc=compute_itc(samples)
assert(~isreal(samples))
s=samples./abs(samples);
itc=abs(sum(s))/size(s,1);
function idx=select_randomly(targets,value,count)
pos=find(targets(:)==value);
[unused,rp]=sort(rand(numel(pos),1));
idx=pos(rp(1:count));
function ds=generate_phase_dataset(varargin)
ndatasets=numel(varargin);
ds_cell=cell(ndatasets,1);
for k=1:ndatasets
ntrials=varargin{k};
ds_k=cosmo_synthetic_dataset('seed',0,...
'nchunks',ntrials,...
'ntargets',1);
ds_k.sa.targets(:)=k;
ds_cell{k}=ds_k;
end
ds=cosmo_stack(ds_cell);
ds.sa.chunks(:)=1:numel(ds.sa.chunks);
sz=size(ds.samples);
ds.samples=randn(sz)+1i*randn(sz);
function test_phase_stat_exceptions
extra_args={'output','pbi'};
aet=@(varargin)assertExceptionThrown(...
@()cosmo_phase_stat(varargin{:}),'');
aet_arg=@(varargin)aet(varargin,extra_args{:});
ds=cosmo_synthetic_dataset('ntargets',2,'nchunks',6);
nsamples=size(ds.samples,1);
sz=size(ds.samples);
ds.samples=randn(sz)+1i*randn(sz);
ds.sa.chunks(:)=1:nsamples;
cosmo_phase_stat(ds,extra_args{:}); % ok
% with targets being 2 or 3 is also ok
ds.sa.targets=ds.sa.targets+1;
cosmo_phase_stat(ds,extra_args{:}); % ok
% input not imaginary
bad_ds=ds;
bad_ds.samples=randn(sz);
aet_arg(bad_ds);
% chunks not all unique
bad_ds=ds;
bad_ds.sa.chunks(1)=bad_ds.sa.chunks(2);
aet_arg(bad_ds);
% imbalance is not ok.
bad_ds=ds;
bad_ds.sa.targets(:)=[repmat([1 2],1,5),[1 1]];
aet_arg(bad_ds);
% bad values for samples_are_unit_length
bad_samples_are_unit_length_cell={[],'',1,[true false]};
for k=1:numel(bad_samples_are_unit_length_cell)
arg={'samples_are_unit_length',...
bad_samples_are_unit_length_cell{k}};
aet_arg(ds,arg{:});
end
% with samples_are_unit_length=true, raise exception if some values
% are not unit length
aet_arg(ds,'samples_are_unit_length',true);
% number of classes must be exactly 2
for bad_class_count=[1,3,4]
bad_ds=ds;
bad_ds.sa.targets(:)=mod(1:nsamples,bad_class_count)+1;
idx=cosmo_index_unique(bad_ds.sa.targets);
counts=cellfun(@numel,idx);
assert(all(counts==counts(1))); % balanced counts
aet_arg(bad_ds); % yet an error is raised
end
% no samples
bad_ds=cosmo_slice(ds,[],1);
aet_arg(bad_ds);
% single sample
bad_ds=cosmo_slice(ds,1,1);
aet_arg(bad_ds);
% balancer function must be function handle
aet_arg(bad_ds,'balancer_func',struct);
% raise exception when called without the 'output' argument, or wrong
% output
aet(ds);
aet(ds,'output','foo');
aet(ds,'output',1);
