function clusters=cosmo_clusterize(sample,nbrhood_mat)
% fast depth-first clustering based on equal values of neighbors
%
% clusters=cosmo_clusterize(ds,nbrhood,sizes)
%
% Inputs:
%   ds        A vector of size 1xN, or a dataset struct with a field
%             .samples of size 1xN.
%   nbrhood   PxN neighborhood definition, if each feature has P
%             neighbors at most. Values of zero can be used as
%             fillers if a feature has less than P neighbors.
%
% Output:
%   clusters  1xQ cell with cluster indices, if Q clusters are found.
%             C{k} denotes the indices in ds that belong to the k-th
%             cluster. Each value in C{k} is in the range 1:N.
%
% Example:
%     % generate tiny fmri dataset with 6 voxels
%     ds=cosmo_synthetic_dataset('ntargets',1,'nchunks',1);
%     % compute neighborhood over voxels, so that two voxels are neighbors
%     % if they share a vertex
%     nh=cosmo_cluster_neighborhood(ds,'progress',false);
%     %
%     % set the data to be, when unflatten, a 2D matrix:
%     % [-1  0 -1]
%     %  3  0 -1]
%     % with three clusters:
%     % - linear index 1         (matrix element (1,1)) with element -1
%     % - linear index 4         (matrix element (2,1)) with element 3
%     % - linear indices 3 and 6 (matrix elments (:,3)) with element -1
%     % Note that elements with indices:
%     % - 1 and 4 are not neighbors because they have different values
%     % - 1 and 3 are not neighbors because they are not connected, as the
%     % neighborhood definition
%     %
%     ds.samples=[2 0 -1 3 0 -1];
%     %
%     % convert neighborhood definition to matrix form
%     nh_mat=cosmo_convert_neighborhood(nh,'matrix');
%     %
%     % clusterize the data
%     cl=cosmo_clusterize(ds.samples,nh_mat);
%     %
%     % show cluster indices
%     cosmo_disp(cl)
%     %|| { [ 1 ]  [ 3    [ 4 ]
%     %||            6 ]        }
%
% Notes:
%   - Two features j and k are in the same cluster if:
%     * their values ds(j) and ds(k) are the same (ds(j)==ds(k)),
%     * features j and k are connected.
%     Features j and k are connected if:
%     * j==k, or
%     % j and k are neigbors, or
%     * there is a feature m, so that j and m are connected, and m and k
%       are connected
%   - Values of zero and NaN in ds are ignored for clusters.
%
% See also: cosmo_cluster_neighborhood, cosmo_convert_neighborhood
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #

    check_input_sizes(sample,nbrhood_mat);

    nfeatures=size(sample,2);
    [maxnbrs,nfeatures_matrix]=size(nbrhood_mat);

    assert(nfeatures==nfeatures_matrix);

    queue=zeros(nfeatures,1); % queue of nodes to add to current cluster
    queue_end=0; % last position where a value in q was stored
    queue_next=1; % first free position

    cl_start=zeros(nfeatures+1,1); % start index of i-th cluster.
                                  % one bigger for ease in setting output
    cl_count=0; % total number of clusters found so far
    cl_idxs=zeros(nfeatures,1); % label for each cluster
    in_queue=false(nfeatures,1); % indicates if an element is in the queue

    data_pos=1; % first non-visited node

    % go over all nodes from 1 to n. xpos is the candidate node for a new
    % cluster, and we add its neighbours in a breadth-first search
    while data_pos<=nfeatures
        if cl_idxs(data_pos)>0 % is already part of a cluster, continue
            data_pos=data_pos+1;
            continue;
        end

        data_val=sample(data_pos); % value of candidate cluster

        if data_val~=0 && ~isnan(data_val)
            % found an element for a new cluster

            queue(queue_next)=data_pos; % first element in the queue
            queue_end=queue_next; % last element of the queue
            cl_count=cl_count+1; % number of clusters
            cl_start(cl_count)=queue_next; % first index of this cluster
            cl_idx=cl_count; % index for current clustter

            in_queue(data_pos)=true;

            % now visit all neighbours of q(qnext),
            % iteratively, until the queue is exhausted
            while true
                qpos=queue(queue_next); % position in queue
                cl_idxs(qpos)=cl_idx; % assign cluster index

                for j=1:maxnbrs
                    nbr=nbrhood_mat(j,qpos);
                    % neighbour should have positive index, not queued yet,
                    % and have the same value as the first element of the
                    % current cluster
                    if nbr>0 && ~in_queue(nbr) && sample(nbr)==data_val
                        % add to queue of nodes to be visited
                        queue_end=queue_end+1;
                        queue(queue_end)=nbr; %

                        % mark it as queued, so it will not be added
                        % to the queue multiple times
                        in_queue(nbr)=true;
                    end
                end

                % get ready for next element in queue
                queue_next=queue_next+1;

                % if queue is exhausted, break out and find next cluster
                if queue_next>queue_end
                    break;
                end
            end
        end
        data_pos=data_pos+1;
    end

    % allocate space for output
    clusters=cell(1,cl_count); % cluster index

    % set extra boundary for last cluster
    cl_start(cl_count+1)=queue_end+1;

    % assign cluster index, size, and value
    for k=1:cl_count
        clusters{k}=queue(cl_start(k):(cl_start(k+1)-1));
    end



function check_input_sizes(sample,nbrhood_mat)
    if size(sample,1)~=1
        error('sample input must be a row vector');
    end

    if ~isnumeric(sample) && ~islogical(sample)
        error('sample input must be numeric or logical');
    end

    if numel(size(nbrhood_mat))~=2
        error('neighborhood matrix must be a matrix');
    end

    if ~isnumeric(nbrhood_mat)
        error('neighborhood matrix must be numeric');
    end

    % ensure they have matching size
    if numel(sample)~=size(nbrhood_mat,2)
        error(['size mismatch between data (%d elements) '...
                'and neighborhood (%d elements)'],...
                    numel(sample),size(nbrhood_mat,2));
    end
