function [hartK] = hartigan_k_index(X,cindsK,cindsKP1)
% HARTIGAN_K_INDEX - Implements the Hartigan cluster index for various clusterings
% 
% Inputs: 
% 
% Written by:
% -- 
% John L. Weatherwax                2008-11-05
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----

[n,p] = size(X); 

%--
% Compute the within class scatter matrix S_W under the clustering with "K" clusters:  
%--

cinds       = cindsK; 
uniq_labels = unique(cinds); 
K           = length(uniq_labels)        % <- the number of clusters 

S_W = zeros( p, p ); 
for jj=1:K,
  tlab = uniq_labels(jj); 
  indj = find( cinds==tlab );               % <- the indices of the samples from class j
  xbarj = mean( X(indj,:) );                % <- the mean of the elements in class j
  nj    = length( indj );                   % <- the number of samples in class j
  for ii=1:nj
    S_W = S_W + ( X(indj(ii),:).' - xbarj(:) ) * ( ( X(indj(ii),:).' - xbarj(:) ).' ); 
  end
end
S_W = (1/n)*S_W; clear ii jj tlab indj xbarj nj; 

% compute the trace of S_W:
tsk = trace(S_W);

%--
% Do the SAME THING but for the clustering with "K+1" clusters:
%--

cinds       = cindsKP1; 
uniq_labels = unique(cinds); 
KP1         = length(uniq_labels);       % <- the number of clusters 

S_W = zeros( p, p ); 
for jj=1:KP1,
  tlab = uniq_labels(jj); 
  indj = find( cinds==tlab );               % <- the indices of the samples from class j
  xbarj = mean( X(indj,:) );                % <- the mean of the elements in class j
  nj    = length( indj );                   % <- the number of samples in class j
  for ii=1:nj
    S_W = S_W + ( X(indj(ii),:).' - xbarj(:) ) * ( ( X(indj(ii),:).' - xbarj(:) ).' ); 
  end
end
S_W = (1/n)*S_W; 

% compute the trace of S_W:
tskp1 = trace(S_W);

% compute the Hartigan index: 
hartK = ( ( tsk / tskp1 ) - 1 ) * ( n - K - 1 );