function stat_repr=cosmo_statcode(ds, output_format)
% Convert statcode for different analysis packages
%
% stat_repr=cosmo_statcode(ds[, output_format])
%
% Inputs:
% ds dataset struct with N samples, or Nx1 cell with strings
% with statistic labels, or AFNI, BV or NIFTI header
% struct for N samples.
% output_format Optional; one of 'afni', 'bv', or 'nifti'; or empty.
%
% Returns:
% stat_repr - If output_format is empty or omitted: an Nx1 cell
% with a string representation of the statistic in each
% sample, e.g. 'Ftest(123,2)' or 'Zscore()' or empty.
% - If output_format=='afni': struct with field
% .BRICK_STATAUX.
% - If output_format=='bv': Nx1 cell with structs, each
% with fieldnames .name, .DF1 and .DF2.
% - If output_format=='nifti': struct with fieldnames
% .intent_code and .intent_p{1,2,3} if all stat codes
% are the same; empty otherwise (NIFTI does not support
% different stat codes for different samples).
%
% Notes:
% - this function is intended for fmri datasets
%
% # For CoSMoMVPA's copyright information and license terms, #
% # see the COPYING file distributed with CoSMoMVPA. #
if nargin<2, output_format=''; end
% see what kind of input we got
if iscellstr(ds)
% string representation of stats
stat_repr=ds(:);
elseif isstruct(ds) && isfield(ds,'samples')
% dataset struct
if isfield(ds,'sa') && isfield(ds.sa,'stats')
% get stat information
stat_repr=ds.sa.stats;
else
% no stat information; return empty output
stat_repr=[];
end
else
% assume header from AFNI, BV or NIFTI; try to convert it to string
% representation, or fail in hdr2strs if that's not the case
stat_repr=hdr2strs(ds);
end
% convert to (name, df) pairs
name_df=stat_strs2name_df(stat_repr);
if isempty(output_format)
% we're done
return
end
% convert to package-specific header
stat_repr=name_df2hdr(name_df, output_format);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% helper functions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function name_df=stat_strs2name_df(stat_strs)
% converts string representations to name and degrees of freedom
match_regexp_failed=@(res) isempty(res) || isempty(res.name);
n=numel(stat_strs);
name_df=cell(n,2);
for k=1:n
stat_str=stat_strs{k};
if isempty(stat_str)
name_df{k,1}='';
else
% try with parentheses: 'stat_name(df1,df2,df3)'
names=regexp(stat_str,'^(?<name>\w+)\((?<df>[,\d]*)\)$',...
'names');
if match_regexp_failed(names)
% try without parentheses: 'stat_name'
names=regexp(stat_str,'^(?<name>\w+)$','names');
df=[];
else
df_split=cosmo_strsplit(names.df,',');
if isequal(df_split,{''})
df=[];
else
df=cellfun(@str2num,df_split);
end
end
if match_regexp_failed(names)
error('Unable to parse stat string %s', stat_str);
end
name_df{k,1}=names.name;
name_df{k,2}=df;
end
end
function hdr=name_df2hdr(name_df, output_format)
% converts statistic name and df to package-specific header
%
% hdr=name_df2hdr(name_df, output_format)
%
% Inputs:
% name_df: Kx2 cell with name and degrees of freedom, for K
% samples
% output_format 'afni', 'bv', or 'nifti'
%
% Returns:
% hdr Part of the package-specific header containing the
% statistical information
codes=get_codes(output_format);
nsamples=size(name_df,1);
stat_idxs=zeros(nsamples,1);
names=codes(:,1);
for k=1:nsamples
name=name_df{k,1};
if isempty(name)
continue;
end
idx=find(cosmo_match(names, name));
if numel(idx)==1
% ignore if not found or multiple matches
stat_idxs(k)=idx;
elseif ~strcmp(name,'none')
cosmo_warning('Unrecognized stat name %s at sample %d',name,k);
end
end
switch output_format
case 'afni'
hdr=struct();
% space for each sample
auxs=cell(1,nsamples);
% build BRICK_STATAUX struct
for k=1:nsamples
stat_idx=stat_idxs(k);
if stat_idx==0
continue;
end
df=name_df{k,2};
auxs{k}=[(k-1) stat_idx numel(df) df];
end
% concatenate information for all samples
hdr.BRICK_STATAUX=[auxs{:}];
case 'bv'
maps=cell(nsamples,1);
for k=1:nsamples
map=struct();
stat_idx=stat_idxs(k);
map.Type=stat_idx;
if stat_idx~=0
df=name_df{k,2};
if numel(df)>=1, map.DF1=df(1); end
if numel(df)>=2, map.DF2=df(2); end
end
maps{k}=map;
end
hdr=maps; % return a cell with maps
case 'nifti'
hdr=struct();
if isempty(stat_idxs)
return
end
unq_stat_idx=unique(stat_idxs);
df=name_df(:,2);
n=numel(df);
ndf=cellfun(@numel,df);
max_ndf=max(ndf);
df_matrix=zeros(n, max_ndf);
for k=1:n
dfk=df{k};
n_dfk=numel(dfk);
if n_dfk>0
df_matrix(k,1:n_dfk)=dfk;
end
end
unq_df=unique(df_matrix,'rows');
if numel(unq_stat_idx)>1 || numel(unique(ndf))>1 || ...
size(unq_df,1)>1
cosmo_warning(['Multiple stat codes found, unsupported '...
'by nifti']);
return
end
hdr.intent_code=unq_stat_idx;
df=zeros(1,3);
if ~isempty(unq_df)
df(1:numel(unq_df))=unq_df;
end
hdr.intent_p1=df(1);
hdr.intent_p2=df(2);
hdr.intent_p3=df(3);
otherwise
assert(false,'should never come here');
end
function stat_strs=hdr2strs(hdr)
% Converts package-specific header to stat string representation
%
% stat_strs=hdr2strs(hdr)
%
% Inputs:
% hdr AFNI, BV vmp, or NIFTI header
%
% Returns:
% stat_strs Kx1
if isfield(hdr, 'Map') && (isfield(hdr,'VMRDimX') ||...
isfield(hdr,'NrOfVertices'))
% bv vmp or smp
codes=get_codes('bv');
nsamples=numel(hdr.Map);
stat_strs=cell(nsamples,1);
for k=1:nsamples
map=hdr.Map(k);
tp=map.Type;
df=[map.DF1 map.DF2];
stat_strs{k}=stat_code2str(codes, tp, df);
end
elseif isstruct(hdr) && isfield(hdr,'DATASET_RANK')
% afni
nsamples=hdr.DATASET_RANK(2);
stat_strs=repmat({''},nsamples,1);
if isfield(hdr,'BRICK_STATAUX')
codes=get_codes('afni');
aux=hdr.BRICK_STATAUX;
naux=numel(aux);
pos=0;
while (pos+2)<naux
pos=pos+1;
brik_idx=aux(pos)+1; % base0 => base1
stat_code=aux(pos+1); % stat code
n_df=aux(pos+2); % # of df
if naux<pos+2+n_df
break;
end
df=aux(pos+2+(1:n_df)); % get dfs
stat_strs{brik_idx}=stat_code2str(codes, stat_code, df);
pos=pos+2+n_df;
end
% sanity check
if pos~=naux
error('Not all elements were processed in STATAUX %s',...
sprintf('%d ', aux));
end
end
elseif isstruct(hdr) && isfield(hdr,'dime') && isfield(hdr.dime,'dim')
% nifti
codes=get_codes('nifti');
dime=hdr.dime;
% deal with potential data in >4th dimension
% (AFNI-NIFTI conversion syndrome)
nsamples=max(prod(dime.dim(5:end)),dime.dim(5));
% get stat codes
stat_code=dime.intent_code;
stat_df=[dime.intent_p1 dime.intent_p2 dime.intent_p3];
stat_str=stat_code2str(codes,stat_code,stat_df);
% repeat for as many samples as there are in the dataset
stat_strs=repmat({stat_str},nsamples,1);
else
error('Unsupported input');
end
function str=stat_code2str(codes, stat_code, df)
% give string representation of stat code with degrees of freedom
% (df can have more elements than required; superfluous ones are ignored.)
if stat_code==0
str='';
return
end
name=codes{stat_code,1};
if isempty(name)
error('Illegal stat_code: %d', stat_code);
end
% build string representation of df
df_count=codes{stat_code,2};
df_str=cellfun(@(x) sprintf('%d',x),num2cell(df(1:df_count)),...
'UniformOutput',false);
% join with stat name
str=sprintf('%s(%s)',name,cosmo_strjoin(df_str,','));
function codes=get_codes(package)
% get names and number of degrees of freedom for an analysis package
%
% codes=get_codes(package)
%
% Inputs:
% package 'nifti', 'afni' or 'bv'.
%
% Returns:
% codes Px2 cell, with codes{k,1} the name of the k-th stat
% and codes{k,2} the number of degrees of freedom
switch package
case 'nifti'
% http://nifti.nimh.nih.gov/nifti-1/documentation/...
% nifti_stats.pdf
codes={'',0;... % not defined
'Correl',1;... % first one starts at 2
'Ttest',1;...
'Ftest',2;...
'Zscore',0;...
'Chisq',1;...
'Beta',2;...
'Binom',2;...
'Gamma',2;...
'Poisson',1;...
'Normal',2,;...
'Ftest_Nonc',3;...
'Chisq_Nonc',2;...
'Logistic',2;...
'Laplace',2;...
'Uniform',2;...
'Ttest_Nonc',2;...
'Weibull',3;...
'Chi',1;...
'Invgauss',2;...
'Extval',2;...
'Pval',0;...
'Logpval',0;...
'Log10pval',0};
case 'afni'
% http://afni.nimh.nih.gov/pub/dist/doc/program_help/...
% README.attributes.html
% AFNI uses a subset of NIFTI.
nifti_codes=get_codes('nifti'); % recursive call
codes=nifti_codes(1:10,:);
case 'bv'
% http://support.brainvoyager.com/documents/Available_Tools/...
% Available_Plugins/niftiplugin_manual_v12.pdf
% why follow a standard (like NIFTI) if one can
% come up with something incompatible too?
codes={'Ttest',1;... % 1
'Correl',1;... % 2 correlation
'Correl',1;... % 3 "cross-correlation"
'Ftest',2;... % 4
'Zscore',0;... % 5
'',0;... % 6 [lots of undefined codes here]
'',0;... % 7
'',0;... % 8
'',0;... % 9
'',0;... % 10
'',0;... % 11 "percent signal change" ignored
'',0;... % 12 "ICA" ignored
'',0;... % 13
'Chisq',1;... % 14
'Beta',2;... % 15
'Pval',0}; % 16
otherwise
error('Unsupported analysis package %s', package);
end
