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.           #


    [dim, merge_method, check]=process_parameters(ds_cell, varargin{:});


    n=numel(ds_cell);

    if n==0
        error('empty cell input');
    elseif n==1
        ds_stacked=ds_cell{1};
        return;
    end

    sample_values=get_struct_values(ds_cell, 'samples');
    if isempty(sample_values)
        error('missing field .samples');
    end
    if check
        sample_sizes=get_dimension_sizes(dim, sample_values);
    else
        sample_sizes=[];
    end
    ds_stacked.samples=stack_values(dim, sample_values, '.samples', check);

    % set the field names for the dimension to be stacked, and the other
    % one
    attrs_keys={'sa','fa'};
    stack_key=attrs_keys{dim};

    to_stack_attr=get_struct_values(ds_cell, stack_key);
    if ~isempty(to_stack_attr)

        ds_stacked.(stack_key)=stack_structs(dim,to_stack_attr,...
                                                sample_sizes,...
                                                ['.' stack_key]);
    end

    other_dim=3-dim; % the other dimension
    merge_key=attrs_keys{other_dim};

    to_merge_attr=get_struct_values(ds_cell, merge_key);
    if ~isempty(to_merge_attr)
        ds_stacked.(merge_key)=merge_structs(to_merge_attr,merge_method,...
                                    ['.' merge_key]);
    end



    to_merge_a=get_struct_values(ds_cell, 'a');
    if ~isempty(to_merge_a)
        ds_stacked.a=merge_structs(to_merge_a,merge_method,'.a');
    end



function [dim, merge, check]=process_parameters(ds_cell, varargin)
    narg=numel(varargin);
    if narg<1 || isempty(varargin{1})
        dim=1;
    else
        dim=varargin{1};
        if ~(isequal(dim,1) || isequal(dim,2))
            error('second argument must be 1 or 2');
        end
    end

    if narg<2 || isempty(varargin{2})
        merge='unique';
    else
        merge=varargin{2};
        if ~isnumeric(merge) && ...
                ~cosmo_match({merge},{'drop','drop_nonunique','unique'})
            error('illegal value for third argument ''merge''');
        end
    end

    if narg<3 || isempty(varargin)
        check=true;
    else
        check=varargin{3};

        if ~(islogical(check) && isscalar(check))
            error('fourth argument ''check'' must be a scalar boolean');
        end

    end

    if ~iscell(ds_cell)
        error('expected cell input');
    end




function values=get_struct_values(struct_cell, key)
    n=numel(struct_cell);
    values=cell(n,1);
    has_values=false;
    for k=1:n
        s=struct_cell{k};
        if isfield(s,key)
            has_values=true;
            values{k}=s.(key);
        end
    end

    if ~has_values
        values=[];
    end


function keys=get_struct_keys(struct_cell)
    n=numel(struct_cell);
    keys=cell(n,1);
    for k=1:n
        keys{k}=fieldnames(struct_cell{k});
    end


function s=merge_structs(vs, merge_method, where)
    if strcmp(merge_method,'drop')
        s=struct();
        return;
    elseif isnumeric(merge_method)
        if ~isscalar(merge_method) || merge_method<1 || ...
                merge_method>numel(vs) || round(merge_method)~=merge_method
            error(['''merge'' parameter, when an integer, must be in'...
                        'the range 1:%d'],numel(vs));
        end
        s=vs{merge_method};
        return;
    elseif ~cosmo_match({merge_method},{'drop_nonunique','unique'})
        error('illegal value for ''merge'' parameter');
    end


    all_keys_cell=get_struct_keys(vs);

    keys=unique(cat(1,all_keys_cell{:}));
    nkeys=numel(keys);

    s=struct();

    for k=1:nkeys
        key=keys{k};
        values=get_struct_values(vs, key);

        [has_unique_elem, unique_elem]=get_single_unique_element(values);

        if ~has_unique_elem
            switch merge_method
                case 'drop_nonunique'
                    % do not store values
                    continue;

                case 'unique'
                    error('non-unique elements in %s.%s',...
                            where, key);

                otherwise
                    error('illegal value for merge: %s', merge_method);
            end
        end

        s.(key)=unique_elem;
    end

function [has_unique_elem, unique_elem]=get_single_unique_element(vs)
    n=numel(vs);

    has_elem=false;
    unique_elem=[];

    for k=1:n
        v=vs{k};
        if isempty(v)
            continue;
        end

        if has_elem
            if ~cosmo_isequaln(v, unique_elem)
                has_unique_elem=false;
                return;
            end
        else
            has_elem=true;
            unique_elem=v;
        end
    end

    has_unique_elem=true;



function s=stack_structs(dim, structs, expected_sizes, where)
    check=~isempty(expected_sizes);
    n=numel(structs);

    for k=1:n
        s=structs{k};

        if ~isstruct(s)
            error('%d-th input is not a struct', k);
        end

        keys=sort(fieldnames(s));

        if k==1
            nkeys=numel(keys);

            stack_args=cell(n, nkeys);

            % keep track of the keys in the first input,
            % and how many values it has along the other dimension
            first_keys=keys;

        elseif ~isequal(keys, first_keys)
            delta=setxor(keys, first_keys);
            error(['key mismatch between 1st and %d-th input in %s ' ...
                        'for key ''%s'''], k, where, delta{1});
        end

        for j=1:nkeys
            key=keys{j};
            v=s.(key);

            if check && size(v,dim)~=expected_sizes(k)
                error(['size mismatch in %d-th input: size(%s.%s,%d)=%d'...
                        ', but size(.samples,%d)=%d'],...
                            k, where, key, dim, size(v,dim), ...
                            dim, expected_sizes(k));
            end

            stack_args{k,j}=v;
        end
    end

    % allocate space for output
    s=struct();

    for j=1:nkeys
        key=keys{j};
        s.(key)=stack_values(dim, stack_args(:,j), where, check);
    end

function sizes=get_dimension_sizes(dim, vs)
    n=numel(vs);
    sizes=zeros(n,1);
    for k=1:n
        sizes(k)=size(vs{k},dim);
    end

function ensure_same_size_along_dim(dim, vs, where)
    % throw an error if element sizes along vs is not the same

    sizes=get_dimension_sizes(dim, vs);
    if numel(sizes)>1
        m=sizes(1)~=sizes(2:end);
        if any(m)
            i=find(m,1)+1;
            error(['size mismatch along dimension %d between 1st '...
                        'and %-dth input for %s'], dim, i, where);
        end
    end


function c=stack_values(dim, vs, where, check)
    % stacks the contents of vs along dimension dim, or throw an error
    % if sizes are not compatible
    if check
        other_dim=3-dim;
        ensure_same_size_along_dim(other_dim, vs, where);
    end

    c=cat(dim, vs{:});

