function [adja, nodeInfo, levelIdx, numCluster, specGap, spectrum, opt] = mapper_spec(X, filter, numLevels, overlap, opt)
% [adja, nodeInfo, levelIdx, opt] = mapper_spec(X, filter, numLevels, overlap, opt)
%
%	This implements mapper by spectral clustering. 
%
% INPUTS
%  X: n-by-p data matrix
%  filter: The filter function used to decompose the space
%  numLevels: Defines the length of each interval
%              length  = (max(filter) - min(filter))*resolution
%              Default : 1/5
%  overlap : The % overlap between successive intervals
%            Default : 50%
%  opt.magicFudge: The number of histograms to divide each homology plot into
%              Increase this number if you expect to see more clusters
%              and decrease it if you are seeing more clusters.
%              Default : 10
%  opt.metric: The way to compute pairwise distance of data points, which will
%				call pdist(X,'euclidean').
%				Default: 'euclidean'
%  opt.levelSize: the minimal size of each level set
%				Default: 100
%  opt.clustSize: the minimal size of each cluster
%				Default: 10
%  opt.histThreshold: the threshold to cut the histogram [0,1], or a vector
%               of length level numbers
%               Default: 0
%  opt.binNumber: choose the zero-bin which cut-off the minimal spanning tree. 
%               Default: 1
%  opt.linkage: ['single']|'average'|'complete'

% Yuan Yao 3/12/2008

%% Check input and set defaults
if nargin < 2
    error('Too few input arguments');
elseif nargin < 3
    % Set resolution to 1/5
    resolution = 1/5;
    % Set overlap to 50%
    overlap = 50;
    % Set magicFudge to 10
    opt.magicFudge = 10
	opt.metric = 'euclidean';
	opt.levelSize = 100;
	opt.clustSize = 10;
elseif nargin < 4
    % Set overlap to 50%
    overlap = 50;
    % Set magicFudge to 10
    opt.magicFudge = 10
	opt.metric = 'euclidean';
	opt.levelSize = 100;
	opt.clustSize = 10;
elseif nargin < 5
    % Set magicFudge to 10
    opt.magicFudge = 10
	opt.metric = 'euclidean';
	opt.levelSize = 100;
	opt.clustSize = 10;
end

opt = parse_opt(opt);

numPoints = size(X,1);
dim = size(X,2);

if (dim~=numPoints),
    error('X must be a symmetric transition count matrix!');
end

if iscell(filter),
    levelSet = filter;
else
    %h = filter;
    [ig,id_sort]=sort(filter,'descend');
    step = round(numPoints/numLevels);
    step_overlap = round(step*(1+overlap/100));
    level = 0;
    for i=1:step:numPoints,
        level = level + 1;
        end_overlap = min([i+step_overlap numPoints]);
        end_nonoverlap = min([i+step numPoints]);
        levelSet{level}.index = id_sort([i:end_nonoverlap]);
        levelSet{level}.overlap = id_sort([end_nonoverlap+1:end_overlap]);	
    end
end

numGraphNodes = 0;
nodeInfo = {};
levelIdx = {};
numCluster = {};
numComponent = {};
adja = [];
nodeInfo={};

if(overlap <=0 | overlap >= 100)
    error('Overlap must be between (0,100)');
end

disp(sprintf('Mapper : number of points %d', numPoints));
disp(sprintf('Mapper : Overlap : %2.1f%%', overlap));
disp('');

tic; 
t0 = toc;

for level=1:length(levelSet),
    
%    t0 = toc;
    
    idx_overlap = union(levelSet{level}.index,levelSet{level}.overlap);
    X1 = X(idx_overlap,idx_overlap);
    d = sum(X1,2);

    % remove zero row if any
    % if (sum(d==0)>0),
    %    idx_nz=find(d>0);
    %    idx_overlap=idx_overlap(idx_nz);
    %    X1 = X1(idx_nz,idx_nz);
    %    d = d(idx_nz);
    % end

    levelSize = length(idx_overlap);
    % perturb X1 to X1 + 0.0001*eye(levelSize) to avoid ill-conditioning
    if levelSet{level}.randsurf>0,
	X1 = X1 + levelSet{level}.randsurf*ones(size(X1))/length(X1);
	d = d + levelSet{level}.randsurf;
    end
    eps = 1e-4;
    [revec,reval]=eigs(X1+eps*eye(levelSize),sparse(1:levelSize,1:levelSize,sum(X1+1e-4*eye(levelSize),2)),min(20,levelSize),1+eps);
    spectrum{level}.vals=diag(reval);
    spectrum{level}.vecs=revec;
    specGap{level} = abs(diff(diag(reval)));
	numComponent{level} = sum(specGap{level}<=eps)+1;
	[ig,idd]=max(abs(specGap{level}));
	numCluster{level} = idd(1);
%	DISCONNECTED = 0;
    
    cluster_tree={};
    cluster_tree{1}.data=[1:levelSize];
    cluster_tree{1}.left = 0;
    cluster_tree{1}.right = 0;

	% Find those disconnected components.
	if numComponent{level}>1,
		cluster_tree{1}.left = 2;
		cluster_tree{1}.right = 3;
		for i=1:numComponent{level},
			mn = mean(abs(spectrum{level}.vecs(:,i)));
			cluster_tree{i+1}.data = find(abs(spectrum{level}.vecs(:,i))>mn);
			cluster_tree{i+1}.left = 0;
			cluster_tree{i+1}.right = 0;
		end		
	end

%    if ~DISCONNECTED,   
	if numCluster{level}>numComponent{level},
		for i = numComponent{level}+1:numCluster{level},
			tree_size = length(cluster_tree);
			variation = zeros(tree_size, 1);             
			for j = 1:tree_size,
				if (cluster_tree{j}.left == 0 & cluster_tree{j}.right == 0),
					tmp = revec(cluster_tree{j}.data, i);
					variation(j) = var( tmp );
				end;
			end;
			[mx ind] = max(variation);

			indices = cluster_tree{ind}.data;
%			if i==2,
%				vec_Fiedler = revec(:,i);
%			else
				X_split = X1(indices,indices);
				n_split = size(X_split,1);
				E = eps*eye(n_split);
				d_split = sum(X_split+E,2);
				[revec_split,reval_split]=eigs(X_split+E,sparse(1:n_split,1:n_split,d_split),min(3,n_split),1+eps);
				vec_Fiedler = revec_split(:,2);
%			end

			left_indices = find(vec_Fiedler >= 0);
			right_indices = find(vec_Fiedler < 0);

			% If neither of the bipartite sets is nonempty, then split the node to leaves 
			if min(length(left_indices),length(right_indices))>0,
				lind = tree_size + 1;
				rind = tree_size + 2;
				cluster_tree{ind}.left = lind;
				cluster_tree{ind}.right = rind;

				cluster_tree{lind}.data = indices(left_indices);
				cluster_tree{rind}.data = indices(right_indices);
				cluster_tree{lind}.left = 0;
				cluster_tree{lind}.right = 0;
				cluster_tree{rind}.left = 0;
				cluster_tree{rind}.right = 0;
			end
		end
	end
	
	leaves =[];
	leaveSize = [];  
	for i=1:length(cluster_tree),
		if (cluster_tree{i}.left==0 & cluster_tree{i}.right==0),
			leaves = [leaves, i];
			leaveSize = [leaveSize, length(setdiff(idx_overlap(cluster_tree{i}.data),levelSet{level}.overlap)) ];
		end
	end
	[leaveSize,id]=sort(leaveSize,'descend');
	leaves = leaves(id(find(leaveSize>opt.clustSize*length(levelSet{level}.index))));
	numCluster{level}=length(leaves);
	
	disp(sprintf('Level %d has %d clusters.',level,numCluster{level}));
	
	numGraphNodes = length(nodeInfo);
	levelIdx{level} = [];

	for i=1:length(leaves),
		new_node = i + numGraphNodes;
		levelIdx{level} = [levelIdx{level}, new_node];
		nodeInfo{new_node}.level=level;
		nodeInfo{new_node}.levelSize = length(levelSet{level}.index);
		nodeInfo{new_node}.set=setdiff(idx_overlap(cluster_tree{leaves(i)}.data),levelSet{level}.overlap);
		nodeInfo{new_node}.set_overlap=setdiff(idx_overlap(cluster_tree{leaves(i)}.data),nodeInfo{new_node}.set);
		if iscell(filter),
			if isfield(opt,'filter'),
				nodeInfo{new_node}.filter = mean(opt.filter(nodeInfo{new_node}.set));
			else
				nodeInfo{new_node}.filter = length(levelSet)-level;
			end
		else
			nodeInfo{new_node}.filter = mean(filter(nodeInfo{new_node}.set));
		end
	end
%    else 
%        [ci,csize] = components(sparse(X1>0));
%        for i=1:max(ci),
%            csize_nonoverlap(i) = length(setdiff(idx_overlap(find(ci==i)),levelSet{level}.overlap));
%        end
%        csize_nonoverlap = reshape(csize_nonoverlap,[length(csize_nonoverlap),1]);
%        large_clusters = find((csize>opt.clustSize*levelSize) .* (csize_nonoverlap>0));
%        numCluster{level}=length(large_clusters);
%
%        disp(sprintf('Level %d has %d clusters.',level,numCluster{level}));
%                
%        numGraphNodes = length(nodeInfo);
%        levelIdx{level} = [];
%
%        for i=1:numCluster{level},
%            new_node = i + numGraphNodes;
%            levelIdx{level} = [levelIdx{level}, new_node];
%            dataset = find(ci==large_clusters(i));
%            nodeInfo{new_node}.level=level;
%            nodeInfo{new_node}.levelSize = length(levelSet{level}.index);
%            nodeInfo{new_node}.set=setdiff(idx_overlap(dataset),levelSet{level}.overlap);
%            nodeInfo{new_node}.set_overlap=setdiff(idx_overlap(dataset),nodeInfo{new_node}.set);
%            if iscell(filter),
%                if isfield(opt,'filter'),
%                    nodeInfo{new_node}.filter = mean(opt.filter(nodeInfo{new_node}.set));
%                else
%                    nodeInfo{new_node}.filter = length(levelSet)-level;
%                end
%            else
%                nodeInfo{new_node}.filter = mean(filter(nodeInfo{new_node}.set));
%            end
%        end
%    end
    
    if(level > 1)
        prevLvlIdx  = levelIdx{level-1};
        thisLvlIdx  = levelIdx{level};
        for i1 = 1:length(prevLvlIdx)
            for i2 = 1:length(thisLvlIdx)
                a = prevLvlIdx(i1);
                b = thisLvlIdx(i2);
                if(length(intersect(union(nodeInfo{a}.set,nodeInfo{a}.set_overlap), union(nodeInfo{b}.set,nodeInfo{b}.set_overlap))) ~= 0)
                    adja(a, b) = 1;
                    adja(b, a) = 1;
                end
            end
        end

    end % if(level > 1)
end

t1 = toc;
dt = t1 - t0;
disp(sprintf('Elapsed time: %.1f seconds', dt));
    
disp('Mapper : Finished');

function option = parse_opt(opt)
% Parsing Option.

option = opt;
option_names = {'metric', 'clustSize', 'levelSize', 'clustSize', 'magicFudge'};

if ~isfield(option,'metric'),
	option = setfield(option,'metric','euclidean');
end

if ~isfield(option,'clustSize'),
	option = setfield(option,'clustSize',10);
end

if ~isfield(option,'levelSize'),
	option = setfield(option,'levelSize', 100);
end

if ~isfield(option,'magicFudge'),
	option = setfield(option,'magicFudge',10);
end

if ~isfield(option,'emptyBin'),
	option = setfield(option,'histThreshold',0);
end

if ~isfield(option,'binNumber'),
	option = setfield(option,'binNumber',1);
end

if ~isfield(option,'linkage'),
	option = setfield(option,'linkage','single');
end

if ~isfield(option,'clusterNumber'),
	option = setfield(option,'clusterNumber',1);
end
