function res=cosmo_distatis(ds, varargin)
% apply DISTATIS measure to each feature
%
% res=cosmo_statis_measure(ds, opt)
%
% Inputs:
% ds dataset struct with dissimilarity values; usually
% the output from @cosmo_dissimilarity_matrix_measure
% applied to each subject followed by cosmo_stack. It
% can also be a cell with datasets (one per subject).
% 'return', d d can be 'distance' (default) or 'crossproduct'.
% 'distance' returns a distance matrix, whereas
% 'crossproduct' returns a crossproduct matrix
% 'split_by', s sample attribute that discriminates chunks
% (participants) (default: 'chunks')
% 'shape', sh shape of output if it were unflattened using
% cosmo_unflatten, either 'square' (default) or
% 'triangle' (which gives the lower diagonal of the
% distance matrix)
%
% Returns:
% res result dataset struct with feature-wise optimal
% compromise distance matrix across subjects
% .samples
%
%
% Example:
% % (This example cannot be documentation tested using Octave,
% % since Octave does not allow for-loops with evalc)
% cosmo_skip_test_if_no_external('matlab');
% %
% ds=cosmo_synthetic_dataset('nsubjects',5,'nchunks',1,'ntargets',4);
% %
% % define neighborhood (here a searchlight with radius of 1 voxel)
% nbrhood=cosmo_spherical_neighborhood(ds,'radius',1,'progress',false);
% %
% % define measure
% measure=@cosmo_dissimilarity_matrix_measure;
% % each subject is a chunk
% ds.sa.chunks=ds.sa.subject;
% % compute DSM for each subject
% sp=cosmo_split(ds,'chunks');
% for k=1:numel(sp)
% sp{k}=cosmo_searchlight(sp{k},nbrhood,measure,'progress',false);
% sp{k}.sa.chunks=ones(6,1)*k;
% end
% % merge results
% dsms=cosmo_stack(sp);
% %
% r=cosmo_distatis(dsms,'return','distance','progress',false);
% cosmo_disp(r);
% %|| .samples
% %|| [ 0 0 0 0 0 0
% %|| 0.818 1.09 0.77 0.653 1.03 0.421
% %|| 0.869 1.3 1.06 1.04 0.932 1.07
% %|| : : : : : :
% %|| 1.16 0.889 0.99 0.631 1.48 0.621
% %|| 0.268 0.952 0.965 0.462 0.943 1.04
% %|| 0 0 0 0 0 0 ]@16x6
% %|| .fa
% %|| .center_ids
% %|| [ 1 2 3 4 5 6 ]
% %|| .i
% %|| [ 1 2 3 1 2 3 ]
% %|| .j
% %|| [ 1 1 1 2 2 2 ]
% %|| .k
% %|| [ 1 1 1 1 1 1 ]
% %|| .nvoxels
% %|| [ 3 4 3 3 4 3 ]
% %|| .radius
% %|| [ 1 1 1 1 1 1 ]
% %|| .quality
% %|| [ 0.685 0.742 0.617 0.648 0.757 0.591 ]
% %|| .nchunks
% %|| [ 5 5 5 5 5 5 ]
% %|| .a
% %|| .fdim
% %|| .labels
% %|| { 'i' 'j' 'k' }
% %|| .values
% %|| { [ 1 2 3 ] [ 1 2 ] [ 1 ] }
% %|| .sdim
% %|| .labels
% %|| { 'targets1' 'targets2' }
% %|| .values
% %|| { [ 1 [ 1
% %|| 2 2
% %|| 3 3
% %|| 4 ] 4 ] }
% %|| .vol
% %|| .mat
% %|| [ 2 0 0 -3
% %|| 0 2 0 -3
% %|| 0 0 2 -3
% %|| 0 0 0 1 ]
% %|| .dim
% %|| [ 3 2 1 ]
% %|| .xform
% %|| 'scanner_anat'
% %|| .sa
% %|| .targets1
% %|| [ 1
% %|| 2
% %|| 3
% %|| :
% %|| 2
% %|| 3
% %|| 4 ]@16x1
% %|| .targets2
% %|| [ 1
% %|| 1
% %|| 1
% %|| :
% %|| 4
% %|| 4
% %|| 4 ]@16x1
%
% Reference:
% - Abdi, H., Valentin, D., O?Toole, A. J., & Edelman, B. (2005).
% DISTATIS: The analysis of multiple distance matrices. In
% Proceedings of the IEEE Computer Society: International conference
% on computer vision and pattern recognition, San Diego, CA, USA
% (pp. 42?47).
%
% Notes:
% - DISTATIS tries to find an optimal compromise distance matrix across
% the different samples (participants)
% - Output can be reshape to matrix or array form using
% cosmo_unflatten(res,1)
%
% # For CoSMoMVPA's copyright information and license terms, #
% # see the COPYING file distributed with CoSMoMVPA. #
cosmo_check_external('distatis');
defaults.return='distance';
defaults.split_by='chunks';
defaults.shape='square';
defaults.mask_output=[];
defaults.progress=100;
defaults.feature_ids=[];
defaults.autoscale=true;
defaults.abs_correlation=false;
defaults.weights='eig';
opt=cosmo_structjoin(defaults,varargin);
subject_cell=get_subject_data(ds,opt);
nsubj=numel(subject_cell);
[dsms,nclasses,dim_labels,dim_values]=get_dsms(subject_cell);
feature_ids=get_feature_ids(size(dsms{1},3),opt);
nfeatures=numel(feature_ids);
quality=zeros(1,nfeatures);
nobservations=zeros(1,nfeatures);
prev_msg='';
clock_start=clock();
show_progress=nfeatures>1 && opt.progress;
for k=1:nfeatures
feature_id=feature_ids(k);
x=zeros(nclasses*nclasses,nsubj);
for j=1:nsubj
dsm=dsms{j}(:,:,feature_id);
x(:,j)=distance2crossproduct(dsm, opt.autoscale);
end
[x,subj_msk]=cosmo_remove_useless_data(x);
nkeep=sum(subj_msk);
% equivalent, but slower:
% [e,v]=eigs(c,1);
[ew,v]=get_weights(x, feature_id, nkeep, opt);
% compute compromise
compromise=x*ew;
result=convert_compromise(compromise, opt);
if feature_id==1
% allocate space
samples=zeros(numel(result),nfeatures);
end
samples(:,k)=result;
quality(:,k)=v/nkeep;
nobservations(:,k)=nkeep;
if show_progress && (k<10 || ...
mod(k, opt.progress)==0 || ...
k==nfeatures)
status=sprintf('quality=%.3f%% (avg)',mean(quality(1:k)));
prev_msg=cosmo_show_progress(clock_start,k/nfeatures,...
status,prev_msg);
end
end
% set output in either triangular or square shape
[res,i,j]=get_samples_in_shape(samples,nclasses,opt.shape);
res=copy_fields(ds,res,{'fa','a'});
% add attributes
res.fa.quality=quality;
res.fa.nchunks=nobservations;
res.a.sdim=struct();
res.a.sdim.labels=dim_labels;
res.a.sdim.values=dim_values;
res.sa.(dim_labels{1})=i(:);
res.sa.(dim_labels{2})=j(:);
cosmo_check_dataset(res);
function [res,i,j]=get_samples_in_shape(samples,nclasses,shape)
res=struct();
switch shape
case 'triangle'
[msk,i,j]=distance_matrix_mask(nclasses);
res.samples=samples(msk(:),:);
case 'square'
res.samples=samples;
[i,j]=find(ones(nclasses));
otherwise
error('unsupported direction %s', shape);
end
function dst=copy_fields(src,dst,keys)
for k=1:numel(keys)
key=keys{k};
if isfield(src,key)
dst.(key)=src.(key);
end
end
function feature_ids=get_feature_ids(nfeatures, opt)
feature_ids=opt.feature_ids;
if isempty(feature_ids);
feature_ids=1:nfeatures;
end
function [ew,v]=get_weights(x, feature_id, nkeep, opt)
switch opt.weights
case 'eig'
[ew,v]=eigen_weights(x, feature_id);
case 'uniform'
% all the same (allowing for comparison with 'eig')
ew=ones(nkeep,1)/nkeep;
v=0;
otherwise
error('illegal weight %s', opt.weights);
end
function subject_cell=get_subject_data(ds,opt)
if isstruct(ds)
subject_cell=cosmo_split(ds,opt.split_by);
else
subject_cell=ds;
end
if numel(subject_cell)==0
error('empty input');
end
function [ew,v]=eigen_weights(x, feature_id)
c=cosmo_corr(x);
negative_c=c<0;
if any(negative_c(:))
[i,j]=find(negative_c);
error(['feature %d has negative correlation between '...
'sample %d and %d, which is not supported by '...
'distatis. DISTATIS assumes that the similarity '...
'data from all samples (typically: participants) '...
'correlate positively. Because that is not the '...
'case, you cannot use DISTATIS analysis on this '...
'data. '],...
feature_id,i(1),j(1));
end
[v,e]=fast_eig1(c);
if all(e<0)
e=-e;
end
assert(all(e>0));
assert(v>0);
% normalize first eigenvector
ew=e/sum(e);
function result=convert_compromise(compromise, opt)
switch opt.return
case 'crossproduct'
result=compromise;
case 'distance'
result=crossproduct2distance(compromise);
otherwise
error('illegal opt.return');
end
function z=crossproduct2distance(x)
n=sqrt(numel(x));
e=ones(n,1);
d=x(1:(n+1):end);
dd=d*e';
ddt=dd';
y=dd(:)+ddt(:)-2*x;
z=ensure_distance_vector(y);
function assert_symmetric(x, tolerance)
if nargin<2, tolerance=1e-8; end
% assert x is a square matrix
sz=size(x);
assert(isequal(sz,sz([2 1])));
xx=x'-x;
msk=xx>tolerance;
if any(msk)
[i,j]=find(msk,1);
error('not symmetric: x(%d,%d)=%d ~= %d=x(%d,%d)',...
i,j,x(i,j),x(j,i),j,i);
end
function z_vec=distance2crossproduct(x, autoscale)
n=size(x,1);
e=ones(n,1);
m=e*(1/n);
ee=eye(n)-e*m';
y=-.5*ee*(x+x')*ee';
if autoscale
z=(1/fast_eig1(y))*y;
else
z=y;
end
assert_symmetric(z);
% equivalent, but slower:
% z=(1/eigs(y,1))*y(:);
z_vec=z(:);
function [lambda,pivot]=fast_eig1(x)
% returns the first eigenvalue in lambda, and the corresponding
% eigenvector in pivot
if cosmo_wtf('is_matlab')
[pivot,lambda]=eigs(x,1);
else
% There seems a bug in Octave for 'eigs',
% so use 'eig' instead.
% http://savannah.gnu.org/bugs/?44004
[e,v]=eig(x);
diag_v=diag(v);
% find largest eigenvalue and eigenvector
[lambda,i]=max(diag_v);
pivot=e(:,i);
end
% The code below is disabled because under certain circumstances
% it would return a near-zero eigenvalue if indeed one eigenvalue (but
% not the largest one) is zero.
% % compute first (largest) eigenvalue and corresponding eigenvector
% % using power iteration method; benchmarking suggests this can be up to
% % five times as fast as using eigs(x,1)
% n=size(x,1);
% pivot=ones(n,1);
% tolerance=1e-8;
% max_iter=1000;
%
% old_lambda=NaN;
% for k=1:max_iter
% z=x*pivot;
% pivot=z / norm(z);
%
% lambda=pivot'*z;
% if abs(lambda-old_lambda)/lambda<tolerance
% z=x*pivot;
% pivot=z / sqrt(sum(z.^2));
%
% lambda=pivot'*z;
% return
% end
% old_lambda=lambda;
% end
%
% % matlab fallback
% [pivot,lambda]=eigs(x,1);
function y=ensure_distance_vector(x)
tolerance=1e-8;
n=sqrt(numel(x));
xsq=reshape(x,n,n);
dx=diag(xsq);
assert(all(dx<tolerance));
xsq=xsq-diag(dx);
delta=xsq-xsq';
assert(all(delta(:)<tolerance));
xsq=.5*(xsq+xsq');
y=xsq(:);
function [dsms,nclasses,dim_labels,dim_values]=get_dsms(data_cell)
nsubj=numel(data_cell);
% allocate
dsms=cell(nsubj,1);
for k=1:numel(data_cell)
data=data_cell{k};
% get data
[dsm,dim_labels,dim_values,is_ds]=get_dsm(data);
% store data
dsms{k}=dsm;
if k==1
nclasses=size(dsm,1);
first_dim_labels=dim_labels;
first_dim_values=dim_values;
data_first=data;
else
if ~isequal(first_dim_labels,dim_labels)
error('dim label mismatch between subject 1 and %d',k);
end
if ~isequal(first_dim_values,dim_values)
error('dim label mismatch between subject 1 and %d',k);
end
% check for compatibility over subjects, raises an error if not
% kosher
if is_ds
cosmo_stack({cosmo_slice(data,1),...
cosmo_slice(data_first,1)},1,'unique');
end
end
end
function [msk,i,j]=distance_matrix_mask(nclasses)
msk=triu(repmat(1:nclasses,nclasses,1),1)'>0;
[i,j]=find(msk);
function [dsm, dim_labels, dim_values, is_ds]=get_dsm(data)
is_ds=isstruct(data);
if is_ds
[dsm,dim_labels,dim_values]=cosmo_unflatten(data,1);
elseif isnumeric(data)
sz=size(data);
if numel(sz)~=2
error('only vectorized distance matrices are supported');
end
[n,nfeatures]=size(data);
side=(1+sqrt(1+8*n))/2; % so that side*(side-1)/2==n
if ~isequal(side, round(side))
error(['size %d of input vector is not correct for '...
'the number of elements below the diagonal of a '...
'square (distance) matrix'], n);
end
[msk,i,j]=distance_matrix_mask(side);
dsm=zeros([side,side,nfeatures]);
assert(numel(i)==n);
for pos=1:n
dsm(i(pos),j(pos),:)=data(pos,:);
end
sq1=cosmo_squareform(data(:,1));
dsm1=dsm(:,:,1);
assert(isequal(sq1,dsm1+dsm1'));
dim_labels={'targets1','targets2'};
dim_values={(1:side)',(1:side)'};
else
error('illegal input: expect dataset struct, or cell with arrays');
end
