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
