function [adja, nodeInfo, levelIdx, ci, csize, specGap] = mapper_level(X,levelSets)
% [adja, nodeInfo, levelIdx, ci, csize, specGap] = mapper_level(X,levelSets,opt)
%
%	Mapper graph by spectral clustering on super-levelsets 
%
% INPUTS
%  X: n-by-p data matrix
%  levelSets: the levelsets index
%  opt.clustSize: the minimal size of each cluster
%				Default: 0

% Yuan Yao 9/6/2008

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

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

numLevels = length(levelSets);
idx_sort = [];
len = [];
for i=1:numLevels,
	len = [len length(levelSets{i}.index)];
	idx_sort = [idx_sort,reshape(levelSets{i}.index,[1,len(i)])];
end

Y = X(idx_sort,idx_sort);

disp(sprintf('Mapper : number of points %d', numPoints));
disp(sprintf('		number of levels %d', length(levelSets)));
disp('');

numGraphNodes = 0;
nodeInfo = {};
levelIdx = {};
numCluster = {};
numComponent = {};
adja = [];
nodeInfo={};
opt.disp = 0;
opt.clusterSizeThreshold = 0;

tic; 
t0 = toc;

for level=1:length(levelSets),

	data = X(levelSets{level}.index,levelSets{level}.index);
	
	if isfield(levelSets{level},'cluster'),
		[ci{level},csize{level}, specGap{level}]=specClust(data,levelSets{level}.cluster,opt);
	else
		[ci{level},csize{level}, specGap{level}]=specClust(data);
	end
	
%	big_cluster{level} = find(csize{level}>opt.clusterSize*length(levelSets{level}.index));

	numCluster{level} = length(csize{level});	
	disp(sprintf('Level %d has %d macrostates out of %d microstates.',level,numCluster{level},length(data)));
	
	numGraphNodes = length(nodeInfo);
	levelIdx{level} = [];

	for i=1:numCluster{level},
		new_node = i + numGraphNodes;
		levelIdx{level} = [levelIdx{level}, new_node];
		nodeInfo{new_node}.level=level;
		nodeInfo{new_node}.levelSize = length(data);
		thisNodeIndex = find(ci{level}==i);
		nodeInfo{new_node}.set = levelSets{level}.index(thisNodeIndex);
%		MetaNodeIndex = thisNodeIndex(find(thisNodeIndex<=numGraphNodes));
%		nodeInfo{new_node}.set = [];
%		for j=1:length(MetaNodeIndex),
%			nodeInfo{new_node}.set = union(nodeInfo{new_node}.set, nodeInfo{MetaNodeIndex(j)}.set);
%		end
%		nodeInfo{new_node}.index = union(nodeInfo{new_node}.index, levelSets{level}.index(setdiff(thisNodeIndex,MetaNodeIndex)));
		nodeInfo{new_node}.filter = level;
	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);
				cut_ab = sum(sum(X(nodeInfo{a}.set, nodeInfo{b}.set)));
				N_ab = cut_ab;
%                if( N_ab> 0)
                    adja(a, b) = N_ab;
                    adja(b, a) = N_ab;
%                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
