function crossed_nbrhood=cosmo_cross_neighborhood(ds, nbrhoods, varargin)
% cross neighborhoods along dataset dimensions
%
% crossed_nbrhood=cosmo_cross_neighborhood(ds,nbrhoods,...)
%
% Inputs:
%   ds            dataset struct
%   nbrhoods      1xK cell with neighborhood structs. Each element can be
%                 the output from cosmo_spherical_neighborhood,
%                 cosmo_meeg_chan_neighborhood,
%                 cosmo_surficial_neighborhood, or
%                 cosmo_interval_neighborhood.
%   'progress',p  if p is true, then progress is shown
%
% Returns:
%   crossed_nbrhood neighborhood struct with fields .[f]a and .neighbors,
%                   constructed by intersecting the neighborhoods from the
%                   input.
%
% Example:
%     % Illustrate neighborhood by crossing freq and time, with freq
%     % 5 bins wide and time 3 bins wide. Each neighborhood contains all
%     % the channels, repeated up to 5*3=15 times (fewer at the border)
%     ds=cosmo_synthetic_dataset('type','timefreq','size','big');
%     freq_nbrhood=cosmo_interval_neighborhood(ds,'freq','radius',3);
%     time_nbrhood=cosmo_interval_neighborhood(ds,'time','radius',5);
%     nbrhood=cosmo_cross_neighborhood(ds, {freq_nbrhood, time_nbrhood},...
%                                                    'progress',false);
%     cosmo_disp(nbrhood.a.fdim)
%     %|| .values
%     %||   { [ 2         4         6  ...  10        12        14 ]@1x7
%     %||     [ -0.2     -0.15      -0.1     -0.05         0 ]           }
%     %|| .labels
%     %||   { 'freq'
%     %||     'time' }
%     cosmo_disp(nbrhood.fa)
%     %|| .freq
%     %||   [ 1         2         3  ...  5         6         7 ]@1x35
%     %|| .time
%     %||   [ 1         1         1  ...  5         5         5 ]@1x35
%     cosmo_disp(nbrhood.neighbors)
%     %|| { [ 1   2   3  ...  9.79e+03  9.79e+03  9.79e+03 ]@1x6120
%     %||   [ 1   2   3  ...  1.01e+04  1.01e+04  1.01e+04 ]@1x7650
%     %||   [ 1   2   3  ...  1.04e+04  1.04e+04  1.04e+04 ]@1x9180
%     %||                                    :
%     %||   [ 307 308 309  ...  1.07e+04  1.07e+04  1.07e+04 ]@1x9180
%     %||   [ 613 614 615  ...  1.07e+04  1.07e+04  1.07e+04 ]@1x7650
%     %||   [ 919 920 921  ...  1.07e+04  1.07e+04  1.07e+04 ]@1x6120 }@35x1
%
% See also: cosmo_spherical_neighborhood, cosmo_meeg_chan_neighborhood,
%           cosmo_interval_neighborhood
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #

    cosmo_check_dataset(ds);
    check_nbrhoods(nbrhoods,ds);

    default.progress=1000;
    opt=cosmo_structjoin(default, varargin{:});

    ndim=numel(nbrhoods);

    dims=struct();
    dims.nbrs=cell(1,ndim);
    dims.labels=cell(1,ndim);
    dims.values=cell(1,ndim);
    dims.fa=cell(1,ndim);

    for k=1:ndim
        nbrhood=ensure_sorted(nbrhoods{k});

        dims.nbrs{k}=nbrhood.neighbors;
        dims.values{k}=nbrhood.a.fdim.values(:);
        dims.labels{k}=nbrhood.a.fdim.labels(:);
        dims.fa{k}=nbrhood.fa;
    end

    % merge labels and values
    dim_labels=cat(1,dims.labels{:});
    dim_values=cat(1,dims.values{:});

    ds_labels=ds.a.fdim.labels;

    % check labels
    check_labels(dim_labels,ds_labels);

    % optimization: compute conjunctions differently in Matlab and Octave
    is_matlab=cosmo_wtf('is_matlab');

    % compute conjunctions of neighborhoods
    [nbr_idxs, nbr_map_idxs]=conj_indices(dims.nbrs, ...
                                        opt.progress, is_matlab);

    % slice feature attributes
    fa_nbrs=cell(ndim,1);
    for k=1:ndim
        fa=dims.fa{k};
        fa_nbrs{k}=cosmo_slice(fa,nbr_map_idxs(k,:),2,'struct');
    end

    crossed_nbrhood=struct();
    crossed_nbrhood.neighbors=nbr_idxs;
    crossed_nbrhood.fa=cosmo_structjoin(fa_nbrs);
    crossed_nbrhood.a=ds.a;
    crossed_nbrhood.a.fdim=struct();
    crossed_nbrhood.a.fdim.values=dim_values;
    crossed_nbrhood.a.fdim.labels=dim_labels;

    origin=struct();
    origin.a=ds.a;
    origin.fa=ds.fa;
    crossed_nbrhood.origin=origin;


function check_nbrhoods(nbrhoods,ds)
    if ~iscell(nbrhoods)
        error(['second argument must a be cell of the form '...
                '{nbrhood1, nbrhood2, ...}, where each nbrhood* '...
                'is a neighborhood structure']);
    end

    for k=1:numel(nbrhoods)
        nbrhood=nbrhoods{k};
        cosmo_check_neighborhood(nbrhood,ds);
    end




function nbrhood=ensure_sorted(nbrhood)
    % ensure everything is sorted, as the helper function
    % 'conj_indices' requires that
    for j=1:numel(nbrhood.neighbors)
        if ~issorted(nbrhood.neighbors{j})
            nbrhood.neighbors{j}=sort(nbrhood.neighbors{j});
        end
    end

function check_labels(dim_labels,ds_labels)
    % ensure no duplicate or missing labels
    if ~isequal(sort(dim_labels), unique(dim_labels)) && ...
                    ~(isempty(dim_labels))
        error('Duplicate dimension labels in %s', ...
                    cosmo_strjoin(dim_labels,','));
    elseif ~all(cosmo_match(dim_labels, ds_labels))
        delta=setdiff(dim_labels, ds_labels);
        error('dimension label unknown in dataset: %s', delta{1});
    end

    show_warning_if_weird_dim_label_order(dim_labels);


function show_warning_if_weird_dim_label_order(dim_labels)
    expected_order={'chan','freq','time'};

    norder=numel(expected_order);

    pos_in_ds=@(needle)find(cosmo_match(dim_labels,{needle}));

    for i=1:(norder-1)
        for j=(i+1):norder
            pre_label=expected_order{i};
            pre=pos_in_ds(pre_label);

            post_label=expected_order{j};
            post=pos_in_ds(post_label);

            if ~isempty(pre) && ...
                    ~isempty(post) && ...
                    pre>post
                cosmo_warning(['dimension labels ''%s'' and ''%s'' are '...
                        'in an uncommon order, which may cause '...
                        'issues when exporting data to FieldTrip '...
                        'or EEGLAB. It is strongly adviced '...
                        'to change the order of the input '...
                        'neighborhoods so that dimension ''%s'''...
                        'preceeds dimension ''%s'''],...
                         pre_label,post_label,post_label,pre_label);

            end
        end
    end







function [flat_idxs, map_idxs]=conj_indices(dim_idxs, show_progress, use_fast)
    % computes conjunction indices
    %
    % Input:
    %   dim_idxs     A NDIMx1 cell, each with X_v cells with indices
    %                As used in this function, dim_idxs{dim}{j} are the
    %                sorted indices with feature attribute for the
    %                dim-th value equal to j.
    %
    % Outputs:
    %   flat_idxs    Nx1 cell values where N=prod(X_*), each of which
    %                has the linear indices of the neighbors of an output
    %                feature.
    %   map_idxs     N*ndim matrix with values in the dim-th column
    %                in the range 1..X_dim. This can be used to index the
    %                values in flat_idxs through sub-indices.

    ndim=numel(dim_idxs);

    % consider first dimension ('head')
    head=dim_idxs{1};
    nhead=numel(head);
    head_map=1:nhead;

    if ndim==1
        % done
        flat_idxs=head;
        map_idxs=head_map;
        return
    end

    % compute indices for remaining dimensions ('tail'), using recursion
    [tail, tail_map]=conj_indices(dim_idxs(2:end), false, use_fast);
    ntail=numel(tail);

    % allocate space for output
    n=nhead*ntail;
    flat_idxs=cell(n,1);
    map_idxs=zeros(ndim,n);

    if show_progress
        prev_msg='';
        clock_start=clock();
    end

    % combine the dimensions
    pos=0;

    for j=1:ntail
        for k=1:nhead
            pos=pos+1;
            headk=head{k};

            if use_fast
                flat_idxs{pos}=fast_intersect(headk, tail{j});
            else
                flat_idxs{pos}=intersect(headk, tail{j});
            end

            map_idxs(1,pos)=head_map(k);
            map_idxs(2:end,pos)=tail_map(:,j);
        end
        if show_progress
            msg=sprintf('crossing neighborhoods');
            prev_msg=cosmo_show_progress(clock_start,j/ntail,msg,prev_msg);
        end
    end


function xy=fast_intersect(x,y)
    % finds the intersection between two vectors
    %
    % xy=fast_intersect(x,y)
    %
    % Inputs:
    %   x     numeric vector with elements sorted.
    %   y     "                                  "
    %
    % Returns:
    %   xy    numeric vector containing elements present in both x and y,
    %         without duplicates and with elements sorted.
    %
    % Notes:
    %  - this function runs in O(n) compared to O(n*log(n)) with
    %    n=max(numel(x),numel(y)) for the built-in function 'intersect',
    %    as that function sorts the input data first.
    %  - in matlab it runs a factor 2 or 3 faster
    %  - in Octave it is very very slow for large inputs

    nx=numel(x);
    ny=numel(y);
    n=min(nx,ny); % maximum size possible for output

    xy=zeros(1,n); % allocate space for output

    pos=0; % last position where a value was stored in xy
    xi=1;  % position in x
    yi=1;  % position in y

    while xi<=nx && yi<=ny
        if x(xi)<y(yi)
            xi=xi+1;
        elseif x(xi)>y(yi)
            yi=yi+1;
        else % x(xi)==y(yi); keep the value
            pos=pos+1;
            xy(pos)=x(xi);

            xi=xi+1;
            yi=yi+1;
        end
    end

    xy=xy(1:pos); % only keep stored elements
