Skip to content

Commit

Permalink
NF: cosmo_phase_stat, to compute phase bifurcation index, phase oppos…
Browse files Browse the repository at this point in the history 
…ition sum, or phase opposition product measures
  • Loading branch information
@nno
nno committed Mar 30, 2017
1 parent 656424b commit 4933d5d
Show file tree
Hide file tree
Showing 2 changed files with 377 additions and 0 deletions.
174 changes: 174 additions & 0 deletions mvpa/cosmo_phase_stat.m
View file
@@ -0,0 +1,174 @@
function stat_ds=cosmo_phase_stat(ds,varargin)
% Compute phase perturbation, or opposition sum or product phase statistic
%
% stat_ds=cosmo_phase_stat(ds,...)
%
% Inputs:
% ds dataset struct with fields:
% .samples PxQ complex matrix for P samples (trials,
% observations) and Q features (e.g. combinations
% of time points, frequencies and channels)
% .sa.targets Px1 array with trial conditions.
% There must be exactly two conditions, thus
% .sa.targets must have exactly two unique
% values.
% .sa.chunks Px1 array indicating which samples can be
% considered to be independent. It is required
% that all samples are independent, therefore
% all values in .sa.chunks must be different from
% each other
% .fa } optional feature attributes
% .a } optional sample attributes
% 'output',p Return statistic, one of the following:
% - 'pbi' - phase bifurcation index
% - 'pos' - phase opposition sum
% - 'pop' - phase opposition product
% 'samples_are_unit_length',u (optional)
% If u==true, then all elements in ds.samples
% are assumed to be already of unit length. If
% this is indeed true, this can speed up the
% computation of the output.
% 'sample_balancer', f (optional)
% If targets are unbalanced (there are more
% trials in one condition than in the other one),
% this function is used to balance the targets.
% If the smallest number of targets is S,
% this function must have the signature
% idxs=f(t1,t2,seed)
% where t1 and t2 are vectors with the target
% positions, seed is an optional value that can
% be used for a pseudo-random number generator,
% and idxs is an Sx2 matrix containing the
% balanced indices for the rows in targets for
% the two classes.
% By default a function is used that
% pseudo-randomly selects a subset for each
% class.
% 'check_dataset',c (optional, default=true)
% if c==false, there is no check for consistency
% of the ds input.
%
% Output:
% stat_ds struct with fields
% .samples 1xQ array with 'pbi', 'pos', or 'pop' function
% .a } if present in the input, then the output
% .fa } contains these fields as well
%
%
% # For CoSMoMVPA's copyright information and license terms, #
% # see the COPYING file distributed with CoSMoMVPA. #


default=struct();
default.sample_balancer=@default_sample_balancer;
default.check_dataset=true;
default.seed=1;
opt=cosmo_structjoin(default,varargin{:});

check_inputs(ds,opt);

% balance dataset
balanced_ds=balance_dataset(ds,opt);

% compute inter-trial coherence for the two conditions
itc_ds=cosmo_phase_itc(balanced_ds,opt);

%itc_ds must have three rows in .samples
assert(size(itc_ds.samples,1)==3);
assert(isnan(itc_ds.sa.targets(3)));

% compute PBI, POP or POS
itc1=itc_ds.samples(1,:);
itc2=itc_ds.samples(2,:);
itc_all=itc_ds.samples(3,:);

stat=compute_phase_stat(opt.output,itc1,itc2,itc_all);

% set result
stat_ds=cosmo_slice(itc_ds,1,1,false);
stat_ds.samples=stat;
stat_ds.sa=struct();



function s=compute_phase_stat(name, itc1, itc2, itc_all)
switch 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,'this should not happen');
end




function balanced_ds=balance_dataset(ds,opt)
target_idxs=cosmo_index_unique(ds.sa.targets);

if numel(target_idxs)~=2
error('Need exactly two targets');
end

sample_idxs=opt.sample_balancer(target_idxs{1},target_idxs{2},...
opt.seed);
balanced_ds=cosmo_slice(ds,sample_idxs,1,false);


function idxs=default_sample_balancer(t1,t2,seed)
n1=numel(t1);
n2=numel(t2);

if n1<n2
idxs=[t1,select_subset_from(t2,n1,seed)];
else
idxs=[select_subset_from(t1,n2,seed),t2];
end

function vec_subset=select_subset_from(vec, count, seed)
% pseudo-random selection of subset
idxs=cosmo_randperm(numel(vec),count,'seed',seed);
vec_subset=vec(idxs);


function [t1,t2]=get_two_targets_row(targets)
idxs=cosmo_index_unique(targets);
if numel(idxs)~=2
error('Input must have exactly two unique values in .sa.targets');
end
t1=idxs{1};
t2=idxs{2};

function check_inputs(ds,opt)
if opt.check_dataset
cosmo_check_dataset(ds);
end

if ~(isstruct(ds) ...
&& isfield(ds,'samples') ...
&& isfield(ds,'sa') ...
&& isfield(ds.sa,'targets'))
error(['first input must be struct with fields .samples and '...
'.sa.targets']);
end

if ~isa(opt.sample_balancer,'function_handle')
error('option ''sample_balancer'' must be a function handle');
end

if ~isfield(opt,'output')
error('option ''output'' is required');
end

allowed_values={'pbi','pos','pop'};
if ~cosmo_match({opt.output},allowed_values)
error('option ''output'' must be one of: ''%s''',...
cosmo_strjoin(allowed_values,''', '''));
end

203 changes: 203 additions & 0 deletions tests/test_phase_stat.m
View file
@@ -0,0 +1,203 @@
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_basics
r=randint();
helper_test_phase_stat_with_trial_counts(r,r);
helper_test_phase_stat_with_trial_counts(r,r+randint());
helper_test_phase_stat_with_trial_counts(r+randint(),r);

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;
opt.sample_balancer=@custom_sample_balancer;
opt.seed=randint()*randint()+randint();

result=cosmo_phase_stat(ds,opt);

% compute expected result
samples=ds.samples;

t1=find(ds.sa.targets==1);
t2=find(ds.sa.targets==2);
sample_idxs=opt.sample_balancer(t1,t2,opt.seed);

itc1=compute_itc(samples(sample_idxs(:,1),:));
itc2=compute_itc(samples(sample_idxs(:,2),:));
itc_all=compute_itc(samples(sample_idxs(:),:));

expected_samples=compute_phase_stat(stat_name,itc1,itc2,itc_all);

assertElementsAlmostEqual(result.samples,expected_samples);
assertEqual(ds.a,result.a);
assertEqual(ds.fa,result.fa);
assertEqual(struct(),result.sa);

if numel(t1)==numel(t2)
% balanced case, default balancer should give same reslt
opt=rmfield(opt,'sample_balancer');

result_default_balancer=cosmo_phase_stat(ds,opt);
assertElementsAlmostEqual(result.samples,...
result_default_balancer.samples);
end




function idxs=custom_sample_balancer(t1,t2,seed)

count=min(numel(t1),numel(t2));
[unused,rp]=sort(cosmo_rand(count,1,...
'seed',sum(t1)+sum(t2)+seed));

assert(count>0);

idxs=[t1(rp), t2(rp)];




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(ntrials1,ntrials2)
ds1=cosmo_synthetic_dataset('seed',0,'nchunks',ntrials1);
ds1.sa.targets(:)=1;

ds2=cosmo_synthetic_dataset('seed',0,'nchunks',ntrials2);
ds2.sa.targets(:)=2;

ds=cosmo_stack({ds1,ds2});
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 ok.
bad_ds=ds;
bad_ds.sa.targets(:)=[repmat([1 2],1,5),[1 1]];
cosmo_phase_stat(bad_ds,extra_args{:});

% 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;
aet_arg(bad_ds);
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);





0 comments on commit 4933d5d

Please sign in to comment.