function ds_stacked=cosmo_stack(ds_cell,varargin)
% stacks multiple datasets to yield a single dataset
%
% ds_stacked=cosmo_stack(datasets[,dim,merge])
%
% Inputs
% datasets cell with input datasets. Each input dataset must be a
% struct and have a field .samples
% dim dimension over which data is stacked. This should be
% either 1 (stack .samples and .sa) or 2 (stack features
% and .fa).
% The default is 1.
% merge Merge strategy for the dimension other than which is
% stacked (i.e. .fa if dim==1 and .sa if dim==2), as well as
% the dataset attributes (.a). Must be one of:
% - 'drop' drop all elements
% - 'drop_nonunique' elements differing are dropped
% - 'unique' raise an exception if elements differ
% - I integer; use data from datasets{I}
% The default is 'unique'
% check Boolean indicating whether to check for the proper size of
% the input datasets. The default is true; setting this to
% false makes this function run faster but, if the input
% dataset do not have proper sizes, will result in less
% informative error messages (if any). Use this option only
% if you are sure that the inputs have proper dimensions.
%
% Ouput:
% ds_stacked Stacked dataset. If dim==1 [or dim==2] and the K-th
% dataset has .samples with size M_K x N_K, then it is
% required that all N_* [M_*] are the same, and the output
% has size sum(M_*) x N_1 [M_1 x sum(N_*)].
% It is also required that
% all input datasets have the same values across all
% the feature [sample] attributes.
%
% Example:
% % Build example data for split-half analysis
% % Normally, half1 and half2 could be from real brain data
% % that is imported using, e.g., cosmo_fmri_dataset
% ds=cosmo_synthetic_dataset('nchunks',2,'ntargets',3);
% half1=cosmo_slice(ds,ds.sa.chunks==1); % .samples is 3x6
% half2=cosmo_slice(ds,ds.sa.chunks==2); % .samples is 3x6
% %
% % join the data in the two halves, over samples.
% % stacking makes .samples a (3+3)x6 matrix
% merged=cosmo_stack({half1,half2});
% cosmo_disp(merged.samples)
% %|| [ 2.03 -0.892 -0.826 -1.08 3.4 -1.29
% %|| 0.584 1.84 1.17 -0.848 1.25 2.04
% %|| -3.68 -0.262 0.321 0.844 -1.37 1.73
% %|| 1.72 0.0975 0.441 1.86 0.479 0.0832
% %|| -1.05 2.04 -0.209 -0.486 -0.955 2.74
% %|| -1.33 0.482 2.39 0.502 1.17 -0.48 ]
% cosmo_disp(merged.sa)
% %|| .chunks
% %|| [ 1
% %|| 1
% %|| 1
% %|| 2
% %|| 2
% %|| 2 ]
% %|| .targets
% %|| [ 1
% %|| 2
% %|| 3
% %|| 1
% %|| 2
% %|| 3 ]
% %
% % data can also be merged over features, by using dim=2
% % here generated data simulates two (tiny) regions of interest
% roi1=cosmo_slice(ds,[1 3],2); % .samples is 6x2
% roi2=cosmo_slice(ds,[2 4 5],2); % .samples is 6x3
% % stacking makes .samples a 6x(2+3) matrix
% roi_both=cosmo_stack({roi1,roi2},2);
% cosmo_disp(roi_both.samples)
% %|| [ 2.03 -0.826 -0.892 -1.08 3.4
% %|| 0.584 1.17 1.84 -0.848 1.25
% %|| -3.68 0.321 -0.262 0.844 -1.37
% %|| 1.72 0.441 0.0975 1.86 0.479
% %|| -1.05 -0.209 2.04 -0.486 -0.955
% %|| -1.33 2.39 0.482 0.502 1.17 ]
% cosmo_disp(roi_both.fa)
% %|| .i
% %|| [ 1 3 2 1 2 ]
% %|| .j
% %|| [ 1 1 1 2 2 ]
% %|| .k
% %|| [ 1 1 1 1 1 ]
% %
% % stacking incompatible datasets gives an error
% cosmo_stack({roi1,half1})
% %|| error('value mismatch: .fa.i')
% cosmo_stack({roi1,half1},2)
% %|| error('value mismatch: .sa.targets')
%
% Note:
% - This function is like the inverse of cosmo_split, i.e. if
% ds_splits=cosmo_split(ds, split_by, dim),
% produces output (i.e., does not throw an error), then using
% ds_humpty_dumpty=cosmo_stack(ds_splits,dim)
% means that ds and ds_humpty_dumpty contain the same data, except that
% the order of the data (in the rows [or columns] of .samples, and
% values in .sa [.fa]) may be in different order if dim==1 [dim==2].
%
% See also: cosmo_split
%
% # For CoSMoMVPA's copyright information and license terms, #
% # see the COPYING file distributed with CoSMoMVPA. #
