function [adja, nodeInfo, levelIdx, opt] = mapper_link(X, filter, resolution, overlap, opt)
% [adja, nodeInfo, levelIdx, opt] = mapper_link(X, filter, resolution, overlap, opt)
%
%	This implements mapper by Linkage clustering. 
%
% INPUTS
%  X: n-by-p data matrix
%  filter: The filter function used to decompose the space
%  resolution: 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);

if(resolution >= 1 | resolution <=0)
    error('Resolution must be between (0,1)');
end

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

if(opt.magicFudge < 0)
    error('magicFudge must be >0');
end

%h = filter;
[ig,id_ascend]=sort(filter,'descend');
numPoints = length(filter);
%step = [0: fix(length(filter)/(1-overlap/100)/(resolution-1)): length(filter)];
%step = [step numPoints];
step = round(numPoints*resolution);
step_overlap = round(step*(1+overlap/100));
numGraphNodes = 0;
level = 1;
nodeInfo = cell(0);
levelIdx = cell(0);
adja = [];

disp(sprintf('Mapper : Filter Range [%0.2f-%0.2f]', min(filter), max(filter)));
disp(sprintf('Mapper : number of points %d', numPoints));
disp(sprintf('Mapper : resolutions %2.1f%%', resolution*100));
disp(sprintf('Mapper : Interval size %d', step));
disp(sprintf('Mapper : Overlap : %2.1f%%', overlap));
disp(sprintf('Mapper : magicFudge : %d', opt.magicFudge));
disp('');

t0=tic; 
dim = size(X,2);
for i=1:step:numPoints,
	t0 = toc;
%	idx = id_ascend(step(i-1)+1:step(i+1));
    %% Select points in this index range
	id0 = i;
%	disp(sprintf(' i+step_overlap: %d; numPoints: %d; min %d',i+step_overlap,numPoints,min(i+step_overlap,numPoints)));
	id1_overlap = min([i+step_overlap numPoints]);
	id1_nonoverlap = min([i+step numPoints]);
	idx = id_ascend([id0:id1_nonoverlap]);	 
	idx_overlap = id_ascend([id0:id1_overlap]);
    disp(sprintf('Mapper: Level %d: filter values from range [%0.2f-%0.2f]', level, min(filter(idx)),max(filter(idx))));
	disp(sprintf('	Level %d has %d points.',level,length(idx))); 

    %% If we have enough points, grow the graph.
    if(max(size(idx)) >= opt.levelSize),
		di = dim*pdist(X(idx_overlap,:),opt.metric);
		ZZ = linkage(di,opt.linkage);
        disp('ZZ size:'); disp(size(ZZ));
        lens = ZZ(:,3);
        
        % histogram of linkage edges with opt.magicFudge bins
        [numBin,cenBin] = hist(lens, opt.magicFudge);
            
        disp('numBin: ');    disp(numBin);
        disp('cenBin: ');    disp(cenBin);

        % opt.emptyBin*100 percentage of the maximal height is the
        % cut-off value 
        cutHist = round(opt.emptyBin*max(numBin));
        disp('Empty Bin at: '); disp(cutHist);
        
        % compute the cumulative sum of numBin, to decide the top edges 
        numBinAccum = cumsum(numBin);    
        
        %
        % Choose the opt.topEdges*100 percent edges
        nz = find(numBin > cutHist & numBinAccum > (1-opt.topEdges)*max(numBinAccum));
        numBinpart = numBin(nz(1):end);
        z = find(numBinpart <= cutHist) + (nz(1) - 1);
        
		if (length(z)==0),
			numClusts = 1;
			c_siz_o = length(idx_overlap);
			c_siz = length(idx);
			c_ind_o = ones(c_siz_o,1);
			c_ind = ones(c_siz,1);
			c_large_id = 1;
        else
            binNumber = 1;

            % Choose the first empty bin center as the cut-off if the
            % accumulated sum is less than 0.25 then the remainder;
            % otherwise choose the second empty bin center as the cut-off.
            if (sum(numBin(1:z(binNumber)))<0.25*sum(numBin(z(binNumber)+1:length(numBin)))),
                cutLink = cenBin(z(binNumber+1));              
            else
                cutLink = cenBin(z(binNumber));
            end

            disp('MST Threshold: '); disp(cutLink);
            c_ind_o = cluster(ZZ,'Cutoff',cutLink,'Criterion','distance');
            c_siz_o = zeros(max(c_ind_o),1);
            for i=1:max(c_ind_o),
                c_siz_o(i) = sum(c_ind_o == i);
            end
			c_ind = c_ind_o(1:length(idx));
            c_siz = c_siz_o;
			for i=1:length(c_siz_o),
                c_siz(i) = sum(c_ind == i);	
            end
            % opt.clustSize is the percentage of the overlapped level size
            % 
            % 			%%c_large_id = find(c_siz>=opt.clustSize);
            c_large_id = find(c_siz>=opt.clustSize*max(size(idx)));
 			numClusts = length(c_large_id);

        end
		disp(sprintf('	Level %d has %d clusters.', level,numClusts));

%         if length(opt.clusterNumber)>1,
%             disp('clusterNumber'); disp(opt.clusterNumber(level));
%      		c_ind_o = cluster(ZZ,'maxclust',opt.clusterNumber(level),'Criterion','distance');
%             c_siz_o = zeros(max(c_ind_o),1);
%             for i=1:max(c_ind_o),
%                 c_siz_o(i) = sum(c_ind_o == i);
%             end
% 			c_ind = c_ind_o(1:length(idx));
%             c_siz = c_siz_o;
% 			for i=1:length(c_siz_o),
% 				c_siz(i) = sum(c_ind == i);
% 				end
% 			%% c_large_id = find(c_siz>=opt.clustSize);
%             c_large_id = find(c_siz>=opt.clustSize*max(size(idx)));
% 			numClusts = length(c_large_id);       
%         else
%             
%     %%  z = find(numBin == 0);
%         if (length(opt.histThreshold)>1),
%             cenBin = opt.histThreshold(level);
%             disp('cenBin: ');    disp(cenBin);
%             z = 1;
%             else
%             [numBin,cenBin] = hist(lens, opt.magicFudge);
%             disp('numBin: ');    disp(numBin);
%             disp('cenBin: ');    disp(cenBin);
%             cutHist = round(opt.histThreshold*max(numBin));
%             disp('histThreshold: '); disp(cutHist);
%             numBinAccum = numBin;
%             for i=1:length(numBin),
%                 numBinAccum(i) = sum(numBin(1:i));
%             end
%             nz = find(numBin > cutHist & numBinAccum > 0.90*max(numBinAccum));
%             numBinpart = numBin(nz(1):end);
%             z = find(numBinpart <= cutHist) + (nz(1) - 1);
%             end
%         
% 		if (length(z)==0),
% 			numClusts = 1;
% 			c_siz_o = length(idx_overlap);
% 			c_siz = length(idx);
% 			c_ind_o = ones(c_siz_o,1);
% 			c_ind = ones(c_siz,1);
% 			c_large_id = 1;
%         else
% 			%binNumber = min(length(z),opt.binNumber);
%             binNumber = 1;
% 
%  %           if (sum(numBin(1:z(binNumber)))<0.25*sum(numBin(z(binNumber)+1:length(numBin)))),
%  %               cutLINK = cenBin(z(binNumber+1));
%  %           else
%                 cutLINK = cenBin(z(binNumber));
%  %           end
%             
% %             if level==6,
% %                 cutMST = 2.1;
% %             end
% 
%             disp('Cut-off value: '); disp(cutLINK);
% 			c_ind_o = cluster(ZZ,'Cutoff',cutLINK,'Criterion','distance');
%             c_siz_o = zeros(max(c_ind_o),1);
%             for i=1:max(c_ind_o),
%                 c_siz_o(i) = sum(c_ind_o == i);
%             end
% 			c_ind = c_ind_o(1:length(idx));
%             c_siz = c_siz_o;
% 			for i=1:length(c_siz_o),
% 				c_siz(i) = sum(c_ind == i);
% 				end
% 			%% c_large_id = find(c_siz>=opt.clustSize);
%             c_large_id = find(c_siz>=opt.clustSize*max(size(idx)));
% 			numClusts = length(c_large_id);
%         end
%         end
% 		disp(sprintf('	Level %d has %d clusters.', level,numClusts));
		
        %% Color the graph
		if numClusts>0,
			numGraphNodes = max(size(nodeInfo));
			levelIdx{level} = [];
			for j=1:numClusts
				levelIdx{level} = [levelIdx{level} (j+numGraphNodes)];
				nodeInfo{j+numGraphNodes}.level = level;
				nodeInfo{j+numGraphNodes}.fnval = [min(filter(idx)),max(filter(idx))];
				idx_tmp = find(c_ind_o(length(idx)+1:length(idx_overlap)) == c_large_id(j));
				if length(idx_tmp) == 0,
					nodeInfo{j+numGraphNodes}.set_overlap = [];
					else
					nodeInfo{j+numGraphNodes}.set_overlap = idx_overlap(length(idx)+idx_tmp); 
					end
				nodeInfo{j+numGraphNodes}.set = idx(find(c_ind == c_large_id(j)));
				nodeInfo{j+numGraphNodes}.filter = mean(filter(nodeInfo{j+numGraphNodes}.set));
                nodeInfo{j+numGraphNodes}.levelSize = length(idx);
				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)
				
				level = level + 1;
        end
        t1 = toc;
        dt = t1 - t0;
        info = sprintf('Elapsed time: %.1f seconds', dt);
        disp(info);
    end
end

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