% function [clustID] = cutDendrogramBioInf(den, data, threshold, startingLevel)
% cutDendrogramBioInf cuts dendrogram den based on the threshold parameter, 
% to extract the corresponding clustering solution. Here a method based 
% on the average intra-cluster standard deviation is used: this value is 
% computed at each level of the dendrogram starting from the leaves, 
% averaging over all existent clusters for that level. When the computed
% value exceeds the set threshold, the process stops, and the 
% clustering solution at the previous level of the dendrogram is returned.
% The function also requires the original data that have been clustered,
% and the level (startingLevel) at which to start climbing the dendrogram
% (= 0 when starting from the leaves).
% clustID lists, for each data point, the ID of the cluster containing that
% point.
%
% This file is part of Cluster4MetaAnalysis.
%     Cluster4MetaAnalysis is free software: you can redistribute it 
%     and/or modify it under the terms of the GNU General Public License 
%     as published by the Free Software Foundation, either version 3 of 
%     the License, or (at your option) any later version.
% 
%     Cluster4MetaAnalysis is distributed in the hope that it will be 
%     useful, but WITHOUT ANY WARRANTY; without even the implied warranty 
%     of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
% 
%     You should have received a copy of the GNU General Public License
%     along with Cluster4MetaAnalysis.  If not, see
%     <http://www.gnu.org/licenses/>.


function [clustID] = cutDendrogramBioInf(den, data, threshold, startingLevel)

% Try all the cuts until the stopping criterion is met...
% If startingLevel = 0, or it is greater than the total number of data 
% points, let us just start from the lowest level. 
% Otherwise, let us take it as a starting level for cutting (for speed sake)

% number of data points
N = size(data, 1);
% dimensionality of the data
numDim = size(data, 2);

if (startingLevel == 0 || startingLevel >= N)
    level = N - 1;          % starting from N - 1 clusters
else
    level = startingLevel;
end

% preparing the data
merge = zeros(N, numDim+1);
merge(:, 1:numDim) = data;
avgStd = zeros(1, numDim);

% climbing the dendrogram to compute the norm of the average 
% standard deviation at each level, until the stopping criterion is met
while (sqrt(sum(avgStd.^2)) <= threshold)
    
    if (level == 0)
        disp('Warning: the threshold you set might be too high! Returning one cluster.');
        clustID = cluster(den, 'MaxClust', 1);
        return;
    else
        avgStd = zeros(1, numDim);
        % find clusters
        cl = cluster(den, 'MaxClust', level);
        merge(:, numDim+1) = cl;

        % Compute the average standard deviations over 
        % all clusters at this cut-level
        for i=1:level 
            idx = find(merge(:, numDim+1) == i);
            if (size(idx, 1) > 1)
                avgStd = avgStd + std(merge(idx, 1:numDim));
            end
        end
        avgStd = avgStd./level;

        % next level
        level = level - 1;
    end

end
level = level + 1;

% last cycle execution made us go over threshold... let us take the
% clustering solution at the previous level.
level = level + 1;
clustID = cluster(den, 'MaxClust', level);















