% function [clustID, clusterData] = findOptimizedSolution(data,
% cutThreshold)
% findOptimizedSolution.m calls uniqueWardClust to compute all
% significantly different dendrograms obtained by hierarchical clustering
% on 'data'. These dendrograms are then cut according to the
% cutThreshold parameter to extract the corresponding clustering solutions.
% The best one (according to some criterion, here the maximum 
% between-cluster error sum of squares) among them is returned. ClustID
% lists for each input data point the corresponding cluster ID, and
% clusterData summarizes information for each cluster (mean, standard
% deviation, cardinality).
%
% 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, clusterData] = findOptimizedSolution(data, cutThreshold)

% running uniqueWardClust function
HG = uniqueWardClust(data);

% number of non-equivalent trees
numTrees = size(HG, 2);
% number of data points
numData = size(data, 1);

% creating a structure for storing clustering solutions
clID = zeros(numData, 1); 
initSol = struct('clusterID', clID, 'clusterData', []);
solutions(1:numTrees) = deal(initSol);

% cut dendrograms
solutions(1).clusterID = cutDendrogram(HG(1).dendrogram, data, cutThreshold, 0);
startingLevel = max(solutions(1).clusterID) + 5;
for i=2:numTrees
    solutions(i).clusterID = cutDendrogram(HG(i).dendrogram, data, cutThreshold, ...
        startingLevel);
end

% computer cluster data for each stored solution
for i=1:numTrees
    solutions(i).clusterData = computeClusterInfo(data, solutions(i).clusterID);
end


% if we have just one dendrogram, we are done... 
if (numTrees == 1)
    disp('***One dendrogram only!***')
    clustID = solutions(1).clusterID;
    clusterData = solutions(1). clusterData;
   
else     
    % find unique solutions
    uniqueSol = findUniqueSolutions(solutions, numTrees);
    solutions = [];   
    
    % compute the best solution according to the user-defined criterion
    % (here, the solution maximizing the b-ESS is chosen; see function
    % computeBESS.m for more details)
    % (notice that you could plug in a different function here)
    bestSol = maximizeBESS(uniqueSol, data);
    
    clustID = bestSol.clusterID;
    clusterData = bestSol.clusterData;
   
end

