% function [updatedD] = updateDissMatrixWard(D, i1, i2, num, d12, idxMap,
% CARD)
% updateDissMatrixWard updates the current dissimilarity matrix D 
% to reflect the merging of elements i1 and i2 (whose dissimilairty 
% value is d12).
% Ward's method is used for computing inter-cluster dissimilarity values.
% The function also requires to specify the number num of current clusters,
% their remapped indexes, idxMap, and their cardinality CARD.
%
% 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 [updatedD] = updateDissMatrixWard(D, i1, i2, num, d12, idxMap, CARD)


updatedD = D;

% computing indexes to be updated: all indexes except for i1 and i2
A = 1:(i1-1);
B = (i1+1):(i2-1);
C = (i2+1):num;
OTHER = [A B C];

% positions of dissimilarity values to be used in the 
% Lance-Williams update formula
dA_i1 = A.*(num - (A+1)/2) - num + i1;
di1_B = i1*(num - (i1+1)/2) - num + B;
di1_C = i1*(num - (i1+1)/2) - num + C;
dA_i2 = A.*(num - (A+1)/2) - num + i2;
dB_i2 = B.*(num - (B+1)/2) - num + i2;
di2_C = i2*(num - (i2+1)/2) - num + C;

I1 = [dA_i1 di1_B di1_C];
I2 = [dA_i2 dB_i2 di2_C];


% update formula (new values are stored in updatedD considering that i1 is
% now the index for the newly created cluster). Ward's criterion is used.
updatedD(I1) = ((CARD(idxMap(OTHER)) + CARD(idxMap(i1))).* D(I1) + ...
    (CARD(idxMap(OTHER)) + CARD(idxMap(i2))).* D(I2) - ...
    (CARD(idxMap(OTHER))* d12))./(CARD(idxMap(OTHER))+ ...
    CARD(idxMap(i1)) + CARD(idxMap(i2)));

% discard values corresponding to cluster i2 (as it does not exist anymore)
di1_di2 = i1 * (num - (i1+1)/2) - num + i2;
I2 = [I2 di1_di2];
updatedD(I2) = [];



