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

% Yuan Yao, Jul-8-2009

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

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

numLevels = length(superLevelSet);
idx_sort = [];
len = [];
for i=1:numLevels,
	len = [len length(superLevelSet{i}.index)];
	idx_sort = [idx_sort,reshape(superLevelSet{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(superLevelSet)));
disp('');

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

tic; 
t0 = toc;

for level=1:length(superLevelSet),
%	if level>1, 
%		opt.v0 = X(superLevelSet{level}.index,superLevelSet{level-1}.index)  * specdata{level-1}.vecs(:,1);
%		disp(sprintf('size v0: %d ',length(opt.v0)));
%	end

	data = X(superLevelSet{level}.index,superLevelSet{level}.index);
	
	if isfield(superLevelSet{level},'cluster'),
		[ci{level},csize{level}, specdata{level}]=specClust(data,superLevelSet{level}.cluster,opt);
	else
		[ci{level},csize{level}, specdata{level}]=specClust(data);
	end
	
%	big_cluster{level} = find(csize{level}>opt.clusterSize*length(superLevelSet{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 = superLevelSet{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, superLevelSet{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);
				N_ab = length(intersect(nodeInfo{a}.set, nodeInfo{b}.set));
%                if( N_ab> 0)
                    adja(a, b) = N_ab;
                    adja(b, a) = N_ab;
%                end
            end
        end

    end % if(level > 1)
end

if numLevels == 1,
	adja = zeros(length(nodeInfo),length(nodeInfo));
end

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


