% function HG = uniqueWardClust(data)
% uniqueWardClust performs Hierarchical Clustering using Ward's criterion
% on the dataset contained stored in 'data' and returns HG 
% (Hierarchical Groupings), which contains the significantly different
% hierarchical groupings (henceforth, dendrograms) that were found 
% by separately developing one dendrogram for each non-equivalent
% minimal-dissimilarity pair (refer to "A novel solution to the problem of
% non-uniqueness of the solution in Hierarchical Clustering" (Cattinelli,
% Valentini, Paulesu, Borghese) for more details about the algorithm and
% the specific terminology). 
%
% 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 HG = uniqueWardClust(data)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%% PREPARING DATA %%%%%%%%%%%%%%%%%%%%%%%%%%%
n = size(data, 1);

% compute squared Euclidean distances and store them to a vector
% (vector form of the dissimilarity matrix)
D = pdist(data, 'euclidean');
D = D .* D;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%% INITIALIZATIONS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% initialize (and pre-allocate) HG
% start assuming LARGE_NUMBER dendrograms will be created
LARGE_NUMBER = 5000;
dendro = zeros(n-1, 3); 
initHG = struct('dendrogram', dendro);
HG(1:LARGE_NUMBER) = deal(initHG);

% Other variables:
% maxCl refers to the total number of clusters that are created in the 
% clustering process:apart the initial n singleton clusters, we are 
% going to create a new cluster at each merging step (these are n-1)
maxCl = 2*n-1;                  
                                
% CARD counts the number of elements in each cluster
CARD = zeros(1, maxCl); 
CARD(1:n) = 1;                  % ... this is for our singleton clusters

% numExCl stores the number of existing clusters at each merging step
% (this will decrease by one at each step, starting from n)
numExCl = n;

% When a new cluster is created, it gets a new ID (n+step). Its 
% dissimilarity value for any other cluster c gets stored in the
% i-th row of dissimilarity matrix, where i is the smaller ID between those
% of the clusters just merged. The row (and column) corresponding to the
% other cluster, j, is removed, so that all clusters IDs from j+1 to n get
% moved back of one position. Because of all this renaming, one must keep
% track of changes in cluster IDs:
% idxMap maps each cluster ID into the appropriate index in the
% dissimilarity matrix
idxMap = 1:n;

% Dendrograms are developed in a depth-first fashion; therefore, we need
% to save the current state of the clustering process each time we need to
% start developing a new dendrogram. A stack structure is used.
initState = struct('dissField', D, 'cardField', CARD, ...
    'numExClField', numExCl, 'idxMapField', idxMap, 'stepField', 1, ...
    'parDendroField', 1, 'altPairsField', [], 'curMinField', -1);
stateStack(1:LARGE_NUMBER) = deal(initState);
stackSize = 1;

% Number of developed dendrograms
numDendro = 0;

% Service vectors, storing, respectively, non-critical and critical pairs 
% for each clustering step 
nonCritPairs = [];
critPairs = [];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%% MAIN CYCLE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% until all dendrograms have been developed...
while stackSize > 0

    % load a previously stored state
    currentState = stateStack(stackSize);
    stackSize = stackSize - 1;    
    [D, CARD, numExCl, idxMap, step, parDendro, critPairs, md] = ...
        deal(currentState.dissField, currentState.cardField, ...
        currentState.numExClField, currentState.idxMapField, ...
        currentState.stepField, currentState.parDendroField,  ...
        currentState.altPairsField, currentState.curMinField);
    
    % a new tree is created (print to assess progress)
    numDendro = numDendro + 1

    % the new dendrogram inherits previous steps from its parent tree;
    % the same goes for the collection of open equivalences
    if (step > 1)        
        HG(numDendro).dendrogram(1:(step-1), :) = ...
            HG(parDendro).dendrogram(1:(step-1), :);
    end
    
    % until the current dendrogram is completed... 
    while (step < n)

         % Compute MD pairs: this is done only once, i.e. these 
         % computations are not repeated at each step if the
         % minimum dissimilarity value is not changing - unless a critical
         % pair was selected for merging. In this case, the other critical 
         % pairs have to be re-assessed.
        if (isempty(critPairs) && isempty(nonCritPairs))    
           
           % Find the minimum dissimilarity value in the current matrix
           md = min(D);  
 
           % Find all MD pairs 
           MDpairs = find(D == md);
           numMD = size(MDpairs, 2);
           
           % if there is one MD pair only, it gets selected for merging
           if (numMD == 1)
                pairIdx = MDpairs(1);
                [p1 p2] = getMatrixIndexes(pairIdx, numExCl);         
           else

               % If more than one MD pair exists, we might have to develop
               % more than one dendrogram.
               %
               % We first distinguish between non-critical pairs (having no
               % extremes in common) and critical pairs

               % Preallocation
               extremes = zeros(numMD, 2);
               nonCritPairs = zeros(numMD, 2);
               critPairs = zeros(numMD, 2);

               % Find extremes for each MD pair
               for q = 1:numMD    
                   pairIdx = MDpairs(q);
                   [p1 p2] = getMatrixIndexes(pairIdx, numExCl);
                   extremes(q, :) = [p1 p2];
               end

               NCcount = 1;
               Ccount = 1;

               % this array will store only critical extremes: 
               % initially it containts all extremes (no repetitions)
               critExtremes = unique(reshape(extremes, 1, ...
                   size(extremes, 1)*size(extremes, 2)));   

               % Separate non-critical and critical pairs 
               for q = 1:numMD
                   p1 = extremes(q, 1);
                   p2 = extremes(q, 2);
                   countP1 = sum(sum(extremes == p1));
                   countP2 = sum(sum(extremes == p2)); 

                   % no extreme of the pair (p1, p2) is shared with
                   % other MD pairs: (p1, p2) is non-critical 
                   % (note the use of idxMap to take into account subsequent 
                   % remappings of indexes)
                   if (max(countP1, countP2) == 1) 
                       nonCritPairs(NCcount, :) = [idxMap(p1) idxMap(p2)];
                       NCcount = NCcount + 1;
                       % erase p1 and p2 from the list of critical extremes
                       critExtremes = setdiff(critExtremes, [p1 p2]);

                   else
                       % (p1, p2) is a critical pair
                       critPairs(Ccount, :) = [p1 p2];
                       Ccount = Ccount + 1;

                       % ... and we check both extremes for assessing 
                       % whether they are both critical
                       if (countP1 == 1)
                           critExtremes = setdiff(critExtremes, p1);
                       end
                       if (countP2 == 1)
                           critExtremes = setdiff(critExtremes, p2);
                       end                 

                   end
               end

               % resizing arrays
               nonCritPairs(NCcount:size(nonCritPairs, 1), :) = [];
               critPairs(Ccount:size(critPairs, 1), :) = [];

               % Then we look into critical pairs to find non-equivalent ones:
               % if two pairs are equivalent, only one is kept
               if (~isempty(critPairs))
                   
                    % special case: at the last-but-one step, two critical 
                    % pairs are necessarily equivalent, so we can simply 
                    % discard the second one
                    if (numExCl <= 3) 
                        critPairs = [];
                        pairIdx = MDpairs(1);
                        [p1 p2] = getMatrixIndexes(pairIdx, numExCl);
                    else
                        % initialization for non-equivalent pairs    
                        countNE = 1;
                        nonequiv = zeros(size(critPairs, 1), 2);

                        % for each critical extreme, we isolate equivalent classes;
                        % basically, we discard pairs that are found to be
                        % equivalent to the stored ones, so that only
                        % non-equivalent pairs are retained.
                        %
                        % Non-equivalent pairs are computed by considering
                        % the first critical extreme. Remaining critical
                        % extremes will be considered at subsequent
                        % iterations
                        sharedExtreme = critExtremes(1);

                        % find all pairs involving sharedExtreme...
                        [r, c] = find(critPairs == sharedExtreme);
                        % ... and extract from these all other (unshared) extremes
                        numUnshared = size(r, 1);
                        unsharedExtremes = zeros(1, numUnshared);
                        % c(p) points to the column position of sharedExtreme for each pair p; we
                        % need the unshared ones...
                        for p=1:numUnshared
                            if (c(p) == 1)
                                c(p) = 2;
                            else
                                c(p) = 1;
                            end
                            unsharedExtremes(p) = critPairs(r(p), c(p));
                            % save pairs into nonequiv
                            nonequiv(countNE, :) = [sharedExtreme unsharedExtremes(p)];
                            countNE = countNE + 1;
                        end
                        
                        % let us now compare pairs of those MD pairs that 
                        % share the same extreme (sharedExtreme)
                        for p=1:numUnshared-1    
                            chosenExtreme = unsharedExtremes(p);
                            
                            % order extremes
                            if (sharedExtreme < chosenExtreme)
                                i1 = sharedExtreme;
                                i2 = chosenExtreme;
                            else
                                i1 = chosenExtreme;
                                i2 = sharedExtreme;
                            end
                            % compute new dissimilarity matrix assuming
                            % that sharedExtreme and chosenExtreme are
                            % merged
                            newD1 = updateDissMatrixWard(D, i1, i2, ...
                                numExCl, md, idxMap, CARD);
                            cl1 = i1;

                            % find cluster(s) having minimum dissimilarity
                            % from cl1 = {sharedExtreme, chosenExtreme} 
                            % in newD1
                            [minFromCl1, nextFromCl1] = ...
                                lookAhead(newD1, cl1, numExCl-1, i2);

                            
                            for q=p+1:numUnshared
                                
                                otherExtreme = unsharedExtremes(q);
                                ot1 = otherExtreme;
                             
                                % find cluster(s) having minimum dissimilarity
                                % from ot1 in newD1 (remapping included)
                                if (ot1 > i2)
                                    [minFromOt1, nextFromOt1] = ...
                                        lookAhead(newD1, ot1-1, numExCl-1, i2);
                                else
                                    [minFromOt1, nextFromOt1] = ...
                                        lookAhead(newD1, ot1, numExCl-1, i2);                                   
                                end
                               

                                % compute new distance matrix assuming
                                % that sharedExtreme and otherExtreme are
                                % merged
                                
                                % order extremes
                                if (sharedExtreme < otherExtreme)
                                    j1 = sharedExtreme;
                                    j2 = otherExtreme;
                                else
                                    j1 = otherExtreme;
                                    j2 = sharedExtreme;
                                end
                                
                                newD2 = updateDissMatrixWard(D, j1, j2, ...
                                    numExCl, md, idxMap, CARD);
                                cl2 = j1;
                                ot2 = chosenExtreme;

                                % find cluster(s) having minimum dissimilarity
                                % from cl2 = {sharedExtreme, otherExtreme} 
                                % in newD2
                                [minFromCl2, nextFromCl2] = ...
                                    lookAhead(newD2, cl2, numExCl-1, j2);

                                % find cluster(s) having minimum dissimilarity
                                % from ot2 in newD2 (remapping included)
                                if (ot2 > j2)
                                    [minFromOt2, nextFromOt2] = ...
                                        lookAhead(newD2, ot2-1, numExCl-1, j2);
                                else
                                    [minFromOt2, nextFromOt2] = ...
                                        lookAhead(newD2, ot2, numExCl-1, j2);                                   
                                end
                               
                     
                                % EQUIVALENCE CHECK 
                                % if the closest point to cl1 is the
                                % excluded extreme ot1, and the closest 
                                % point to ot1 is cl1 - and the same goes 
                                % for cl2 and ot2, the two pairs 
                                % {sharedExtreme, chosenExtreme}, 
                                % {sharedExtreme, otherExtreme} are 
                                % equivalent - and the second one
                                % can be simply discarded
       
                                if (any(nextFromOt1  == cl1) && any(nextFromCl1 == ot1) ...
                                    && any(nextFromOt2 == cl2) && any(nextFromCl2 == ot2) ...
                                    && minFromOt1 == minFromOt2)
                               
                                    % the two pairs are equivalent: 
                                    % discard the second one
                                    [r, c] = find(nonequiv == otherExtreme);
                                    if (~isempty(r))
                                         nonequiv(r, :) = [];
                                         countNE = countNE - 1;
                                    end
 
                                 end
         
                             end

                        end
                               
                        % resizing and remapping (non-equivalent pairs are
                        % stored in the critPairs array
                        nonequiv(countNE:size(nonequiv, 1), :) = [];
                        critPairs = idxMap(nonequiv);
                        % just in case multiple copies are present...
                        critPairs = unique(critPairs, 'rows');       
                    end
                  
               end
               
           end
            
        end


       % if there are non-critical pairs, they are the first to be 
       % selected for merging: no new dendrogram is created, non-critical 
       % pairs are just merged in the (random) order they appear in
       if (~isempty(nonCritPairs)) 

           % select first non-critical pair
           % (re-mapping has to be taken into account)
           p1 = find(idxMap == nonCritPairs(1, 1));   
           p2 = find(idxMap == nonCritPairs(1, 2));
 
           % re-ordering
           if (p1 > p2)
               dummy = p1;
               p1 = p2;
               p2 = dummy;
           end
           
           % discard selected pair from the list of 
           % remaining non-critical pairs
           nonCritPairs(1, :) = [];
       
           
       % otherwise, we start merging critical pairs
       % (one for equivalence class, because of how we built 'critPairs')
       else
          if (~isempty(critPairs))
               
               % select first critical pair
               % (re-mapping has to be taken into account) 
               p1 = find(idxMap == critPairs(1, 1)); 
               p2 = find(idxMap == critPairs(1, 2));
                
               % re-ordering
               if (p1 > p2)
                   dummy = p1;
                   p1 = p2;
                   p2 = dummy;
               end
                
                % discard selected pair from the list of 
                % remaining critical (non-equivalent) pairs
                critPairs(1, :) = [];
                
                % if other critical pairs are left, new dendrograms will 
                % have to be created, one for each remaining pair. 
                % The current state is saved.
                if (~isempty(critPairs))
                   stackSize = stackSize + 1;
                   stateStack(stackSize) = ...
                       struct('dissField', D, 'cardField', CARD, ...
                       'numExClField', numExCl, 'idxMapField', idxMap, ...
                       'stepField', step, 'parDendroField', numDendro, ...
                       'altPairsField', critPairs, 'curMinField', md);
                end
                
                % the list of critical pairs is emptied - at the next step,
                % critical pairs will have to be re-assessed 
                % (the merging of a critical pair can cause 
                % other MD pairs to not disappear)
                critPairs = [];

           end
        end
 
        % store this merging step into the current dendrogram
        % representation
        HG(numDendro).dendrogram(step, :) = [idxMap(p1) idxMap(p2) md];

        % The dissimilarity matrix and the other variables
        % (idxMap, CARD, numExCl, step) are updated
        D = updateDissMatrixWard(D, p1, p2, numExCl, md, idxMap, CARD);
        
        % the newly created cluster gets ID 'n+step', and its cardinality
        % is computed as the sum of the cardinalities of the two
        % merged clusters
        CARD(n+step) = CARD(idxMap(p1)) + CARD(idxMap(p2));
        
        % in the new dissimilarity matrix, distances from the new cluster
        % (whose ID is n+step) are stored where p1 was previously stored.
        % Since p2 does not exist anymore, indexes are shifted to the left
        % starting from p2.
        idxMap(p1) = n+step;
        idxMap(p2:(n-1)) = idxMap((p2+1):n);
        idxMap(n) = 0;
        
        % the number of clusters decreases by one, 
        % and we go on to the next step
        numExCl = numExCl - 1;
        step = step + 1;

    end

end

% resizing 
HG(numDendro+1:LARGE_NUMBER) = [];

% unsquaring dissimilarity values and ordering cluster indexes
for i=1:numDendro
    HG(i).dendrogram(:, 3) = sqrt(HG(i).dendrogram(:, 3));
    HG(i).dendrogram(:, [1 2]) = sort(HG(i).dendrogram(:, [1 2]), 2);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%





