function [nh_split,masks]=cosmo_neighborhood_split(nh, varargin)
% partitions a neighborhood in a cell with (smaller) neigborhoods
%
% cosmo_neighborhood_split(nh, ...)
%
% Inputs:
%   nh                  neighborhood struct with fields .neighbors and
%                       .origin
%   'divisions', d      Number of divisions along each feature dimension.
%                       For example, if d=4 for an fMRI dataset (with
%                       three feature dimensions corresponding to the three
%                       spatial dimensions), then the output will have
%                       at most 4^3=64 neighborhood structs and masks.
%
% Output:
%   nh_split            Cell with neighborhood structs. When all elements
%                       in nh_split are combined they should contain the
%                       same information in .neighbors and .fa as in the
%                       input nh, but with a re-indexing of the feature ids
%                       corresponding to the masks.
%   masks               Cell with mask datasets, corresponding to the
%                       features returned in neighborhoods in nh_split.
%                       When the dataset derived from the input
%                       neighborhood nh is masked by the k-th mask, then
%                       the indices in the resulting dataset match the
%                       indices in .neighbors
%
%
% Example:
%     % (This example cannot be documentation tested using Octave,
%     %  since Octave does not allow for-loops with evalc)
%     cosmo_skip_test_if_no_external('matlab');
%     %
%     % the following example shows how a neighborhood can be split in
%     % parts, a searchlight run on each part, and the different outputs
%     % joined together. This gives identical results as running the
%     % searchlight on the entire dataset
%     %
%     % generate a tiny dataset with 6 voxels
%     ds=cosmo_synthetic_dataset();
%     %
%     % define tiny neighborhoods with radii of 1 voxel
%     % (more typical values are ~3 voxels)
%     nh=cosmo_spherical_neighborhood(ds,'radius',1,'progress',false);
%     %
%     % Use a simple 'measure' that just averages data along the feature
%     % dimension. (Any other measure can be used as well)
%     averager=@(x,opt)cosmo_structjoin('samples',mean(x.samples,2));
%     %
%     % disable progress reporting for the searchlight
%     opt=struct();
%     opt.progress=false;
%     %
%     % run the searchlight
%     result=cosmo_searchlight(ds,nh,averager,opt);
%     %
%     % show output from searchlight
%     cosmo_disp({result.samples; result.fa});
%     %|| { [ 0.768     0.368        -1      1.45    0.0342     -0.32
%     %||     0.526      1.77     0.937      1.08      1.07      1.49
%     %||      0.46     -1.25    -0.152    0.0899     0.795    -0.522
%     %||      1.23     0.685     0.954     0.606      1.19     0.335
%     %||     0.914   -0.0442    0.0295     0.664     0.274    -0.221
%     %||     0.633      1.24     0.797     0.861      1.54      1.02 ]
%     %||   .nvoxels
%     %||     [ 3         4         3         3         4         3 ]
%     %||   .radius
%     %||     [ 1         1         1         1         1         1 ]
%     %||   .center_ids
%     %||     [ 1         2         3         4         5         6 ]
%     %||   .i
%     %||     [ 1         2         3         1         2         3 ]
%     %||   .j
%     %||     [ 1         1         1         2         2         2 ]
%     %||   .k
%     %||     [ 1         1         1         1         1         1 ]   }
%     %
%     [nh_splits,masks_split]=cosmo_neighborhood_split(nh);
%     n_split=numel(nh_splits);
%     %
%     % Allocate space for searchlight output from each part
%     res2_cell=cell(1,n_split);
%     %
%     for k=1:n_split
%         % get the k-th neighborhood
%         nh_sel=nh_splits{k};
%         %
%         % slice the dataset to select only the features indexed by
%         % the nh_sel neighborhood.
%         %
%         % (when using fmri datasets, it is also possible to use
%         %  cosmo_fmri_dataset with a filename and the mask from masks,
%         %  as in:
%         %
%         %     ds_sel=cosmo_fmri_dataset('data.nii','mask',masks_split{k});
%         %
%         %  Such an approach may result in significant reduction of memory
%         %  usage, if the file format supports loading partial files
%         %  (currently unzipped NIFTI (.nii), ANALYZE (.hdr) and AFNI)))
%         ds_sel=cosmo_slice(ds,masks_split{k}.samples,2);
%         %
%         % run the searchlight for this split, and store the result
%         res2_cell{k}=cosmo_searchlight(ds_sel,nh_sel,averager,opt);
%     end;
%     %
%     % join the results along the second (feature) dimension
%     result2=cosmo_stack(res2_cell,2);
%     %
%     % show the results; they are identical to the original output, modulo
%     % a possible permutation of the features.
%     cosmo_disp({result2.samples; result2.fa});
%     %|| { [ 0.768     0.368      1.45    0.0342        -1     -0.32
%     %||     0.526      1.77      1.08      1.07     0.937      1.49
%     %||      0.46     -1.25    0.0899     0.795    -0.152    -0.522
%     %||      1.23     0.685     0.606      1.19     0.954     0.335
%     %||     0.914   -0.0442     0.664     0.274    0.0295    -0.221
%     %||     0.633      1.24     0.861      1.54     0.797      1.02 ]
%     %||   .center_ids
%     %||     [ 1         2         3         4         1         2 ]
%     %||   .i
%     %||     [ 1         2         1         2         3         3 ]
%     %||   .j
%     %||     [ 1         1         2         2         1         2 ]
%     %||   .k
%     %||     [ 1         1         1         1         1         1 ]
%     %||   .nvoxels
%     %||     [ 3         4         3         4         3         3 ]
%     %||   .radius
%     %||     [ 1         1         1         1         1         1 ]   }
%
%
% Notes:
%   - typical use cases of this function are:
%     * searchlights on fMRI dataset with many samples and features
%       (i.e. too many to load all of them in memory). In this case, and
%       if the file format supports it (NIFTI, ANALYZE, AFNI), parts
%       of the data can be loaded with cosmo_fmri_datset, where the masks
%       are based on the second output from this function
%     * parallel searchlights on a cluster using multiple computing nodes.
%
% See also: cosmo_searchlight, cosmo_fmri_dataset
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #

    defaults=struct();
    defaults.divisions=4;

    opt=cosmo_structjoin(defaults,varargin);
    validate_input(nh, opt);

    % see how big each block is along each dimension
    [first,last,nfeatures]=get_dim_range(nh);

    min_first=min(first,[],2);
    max_last=max(last,[],2);
    max_size=max(last-first,[],2)+1;
    extent=max_last-min_first+1;

    % Examples of how max_size, step and extent are related:
    % max_size=1, step=5, extent=25, 5 div: 1-5, 6-10, 11-15, 16-20, 21-25
    % max_size=1, step=4, extent=20, 5 div: 1-4,  5-8,  9-12, 13-16, 17-20
    % max_size=2, step=5, extent=26, 5 div: 1-6, 6-11, 11-16, 16-21, 21-26
    % max_size=2, step=4, extent=21, 5 div: 1-5,  5-9,  9-13, 13-17, 17-21
    % max_size=3, step=4, extent=22, 5 div: 1-6, 5-10,  9-14, 13-18, 17-22

    % first element:                step+max_size-1     elements
    % all other (div-1) elements:   step                elements

    step=ceil((extent+1-max_size)/opt.divisions);
    assert(all(step>=1));

    % for each element in nh.neighbors, assign it to a block
    % indexed by ndim integers (where ndim is the number of dimensions)
    dim_pos=floor(bsxfun(@rdivide,first-1,step))+1;

    % get the indices for each block.
    % idxs{k}=[a_k1, ..., a_kN] means that nh.neighbors{a_k1}, ...,
    % nh.neighbors{a_kN} all refer to neighboring features
    %
    [idxs,unq]=cosmo_index_unique(dim_pos');

    nsplit=size(unq,1);

    nh_split=cell(nsplit,1);
    masks=cell(nsplit,1);

    keep_split_mask=false(nsplit,1);

    for split_id=1:nsplit
        msk=false(1,nfeatures);
        idx=idxs{split_id};
        n_nbrs=numel(idx);
        for j=1:n_nbrs
            nb=nh.neighbors{idx(j)};
            msk(1,nb)=true;
        end

        if ~any(msk)
            % skip empty neighbors
            continue;
        end

        keep_split_mask(split_id)=true;


        m_k=struct();
        m_k.samples=msk;
        m_k.a=nh.origin.a;
        m_k.fa=nh.origin.fa;

        nh_k=struct();
        nh_k.origin.a=m_k.a;
        nh_k.origin.fa=cosmo_slice(m_k.fa,msk,2,'struct');
        nh_k.a=nh.a;
        nh_k.fa=cosmo_slice(nh.fa, idx, 2, 'struct');

        all2some=zeros(1,nfeatures);
        all2some(msk)=1:sum(msk);
        neighbors=cell(n_nbrs,1);
        for j=1:n_nbrs
            some=all2some(nh.neighbors{idx(j)});
            assert(all(some~=0));
            neighbors{j}=some;
        end
        nh_k.neighbors=neighbors;

        nh_split{split_id}=nh_k;
        masks{split_id}=m_k;
    end

    % only keep non-empty masks
    nh_split=nh_split(keep_split_mask);
    masks=masks(keep_split_mask);


function [first,last,nfeatures]=get_dim_range(nh)
    % for each feature dimension, find the minimum and maximum value
    % in each neighborhood
    % If nh has N neighbors and M dimensions, then the output of first and
    % last in M x N.
    fdim=nh.origin.a.fdim;
    labels=fdim.labels;

    n_labels=numel(labels);
    n_neighbors=numel(nh.neighbors);

    fa=nh.origin.fa;

    % pre-store values for each feature attribute
    fa_idxs_cell=cell(n_labels,1);

    for k=1:n_labels
        label=labels{k};
        fa_idxs=fa.(label);
        fa_idxs_cell{k}=fa_idxs;
    end

    first=zeros(n_labels,n_neighbors);
    last=zeros(n_labels,n_neighbors);
    nfeatures=0;



    for j=1:n_neighbors
        nb=nh.neighbors{j};

        mx_nb=max(nb);
        if mx_nb>nfeatures
            nfeatures=mx_nb;
        end

        for k=1:n_labels
            % (this is a timing-critical part of the loop)
            nb_values=fa_idxs_cell{k}(nb);

            first(k,j)=min(nb_values);
            last(k,j)=max(nb_values);
        end
    end


function opt=validate_input(nh, opt)
    % some checks on the input
    cosmo_check_neighborhood(nh);

    allowed_fieldnames={'divisions'};
    delta=setdiff(fieldnames(opt), allowed_fieldnames);
    if ~isempty(delta)
        error('illegal field %s', delta{1});
    end

    divisions=opt.divisions;
    if ~(isnumeric(divisions) && ...
                isscalar(divisions) && ...
                isfinite(divisions) && ...
                divisions > 1)
        error('count must be finite scalar >1');
    end


    origin_labels={'origin.fa',...
                    'origin.a.fdim.labels',...
                    'origin.a.fdim.values'};
    if ~all(cosmo_isfield(nh,origin_labels))
        error('origin has missing labels');
    end

    is_numeric_vector=@(x)isnumeric(x) && isvector(x);
    if ~all(cellfun(is_numeric_vector, nh.a.fdim.values))
        error('elements in .a.fdim.values must be numeric vectors');
    end



