function [arr, dim_labels, dim_values]=cosmo_unflatten(ds, dim, varargin)
% unflattens a dataset from 2 to (1+K) dimensions.
%
% [arr, dim_labels, dim_values]=cosmo_unflatten(ds, [dim, ][,...])
%
% Inputs:
% ds dataset structure, with fields:
% .samples PxQ for P samples and Q features.
% .a.Xdim.labels 1xK cell with string labels for each dimension,
% with X='s' for samples (dim=1) or X='f' for features
% (dim=2).
% .a.Xdim.values 1xK cell, with S_J values (J in 1:K) corresponding
% to the labels in each of the K dimensions.
% .Xa.(label) for each label in a.Xdim.labels it contains the
% sub-indices for the K dimensions. It is required
% that for every dimension J in 1:K, all values in
% ds.fa.(a.fdim.labels{J}) are in the range 1:S_K, and
% that every combination across labels is unique.
% dim dimension to be unflattened, either 1 (for samples)
% or 2 (for features; default)
% 'set_missing_to',s value to set missing values to (default: 0)
% 'matrix_labels',m Allow labels in the cell string m to be matrices
% rather than vectors. Currently the only use case is
% the 'pos' attribute for MEEG source space data.
%
% Returns:
% arr S_1 x ... x S_K x Q array if (dim==1), or
% P x S_1 x ... x S_K array if (dim==2), where
% Q=prod(S_*) if dim==1 and P=prod(S_*) if dim==2
% dim_labels the value of .a.Xdim.labels
% dim_values the value of .a.Xdim.values
%
% Example:
% % ds is an FMRI dataset with 6 samples, volumes are 3 x 2 x 5 voxels
% ds=cosmo_synthetic_dataset('size','normal','type','fmri');
% size(ds.samples)
% %|| [ 6 30 ]
% cosmo_disp(ds.a.fdim)
% %|| .labels
% %|| { 'i' 'j' 'k' }
% %|| .values
% %|| { [ 1 2 3 ] [ 1 2 ] [ 1 2 3 4 5 ] }
% %
% % flatten the dataset
% [unfl,labels,values]=cosmo_unflatten(ds);
% %
% % the unflattened dataset is of size 6 x 3 x 2 x 5
% size(unfl)
% %|| [ 6 3 2 5 ]
% cosmo_disp(labels)
% %|| { 'i' 'j' 'k' }
% cosmo_disp(values)
% %|| { [ 1 2 3 ] [ 1 2 ] [ 1 2 3 4 5 ] }
%
% % ds is a small dataset with 2 classes
% ds=cosmo_synthetic_dataset();
% %
% % compute all (2x2) split-half correlation values
% res=cosmo_correlation_measure(ds,'output','raw',...
% 'post_corr_func',[]);
% cosmo_disp(res)
% %|| .samples
% %|| [ 0.363
% %|| -0.404
% %|| -0.447
% %|| 0.606 ]
% %|| .sa
% %|| .half1
% %|| [ 1
% %|| 2
% %|| 1
% %|| 2 ]
% %|| .half2
% %|| [ 1
% %|| 1
% %|| 2
% %|| 2 ]
% %|| .a
% %|| .sdim
% %|| .labels
% %|| { 'half1' 'half2' }
% %|| .values
% %|| { [ 1 [ 1
% %|| 2 ] 2 ] }
% %
% % reshape the correlations into a square matrix
% [unfl,labels,values]=cosmo_unflatten(res,1);
% %
% % yields a 2x2x1 matrix (matlab omits the last, singleton dimension)
% cosmo_disp(unfl)
% %|| [ 0.363 -0.447
% %|| -0.404 0.606 ]
% %
% cosmo_disp(labels)
% %|| { 'half1' 'half2' }
% %
% cosmo_disp(values)
% %|| { [ 1 [ 1
% %|| 2 ] 2 ] }
%
%
% Notes:
% - A typical use case is mapping an fMRI or MEEG dataset struct
% back to a 3D or 4D array.
% - This function is the inverse of cosmo_flatten.
%
% See also: cosmo_flatten, cosmo_map2fmri, cosmo_map2meeg
%
% # For CoSMoMVPA's copyright information and license terms, #
% # see the COPYING file distributed with CoSMoMVPA. #