function CVInd = inds(ClassLabs, Data, Centres, sumd, Method, p, q,tol,spamitmax)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Implementation of ESANN2017-papers' best internal cluster validation indices
%   (c) Tommi K., Nov 27, 2016
%Generalized for median, June 2017, Susanne Jauhiainen
%Modified, 15.6.2017, Joonas H
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if(nargin<8)
    tol = 1.d-5;
    spamitmax = 100;
end

%Define the dist function and center of whole data
if (p==2 && q == 2)
    distFun = @(X,Y) pdist2(X,Y).^2; %Squared Euclidean distances between X and Y
    s = mean(Data);
elseif (p==1 && q ==1)
    distFun = @(X,Y) pdist2(X,Y,'cityblock');
    s = median(Data);
elseif (p==2 && q ==1)
    distFun = @(X,Y) pdist2(X,Y,'euclidean');
    [s,~] = SpatMed3(Data,tol,spamitmax);
end

%Get the size of data and K
[N,n] = size(Data);
K = max(ClassLabs);

%Calculate the clustering error
FinalDists = min(distFun(Data,Centres),[],2);
ClustErr = sum(FinalDists); %sama kuin sum(sumd)

% The new K x Clustering Error CVAI
if strcmp(Method,'kCE') %k times ClusteError
    CVInd = K*ClustErr;
    
%WB index
elseif strcmp(Method,'WB')
    BGSS = 0;
    for i=1:K
        I = find(ClassLabs == i);
        BGSS = BGSS + length(I)*distFun(Centres(i,:),s);
    end
    CVInd = (K*ClustErr)/BGSS;
    
%Calinski-Harabasz
elseif strcmp(Method,'CH') 
    BGSS = 0;
    for i=1:K
        I = find(ClassLabs == i);
        BGSS = BGSS + length(I)*distFun(Centres(i,:),s);
    end
    intra = (K-1)*ClustErr;
    inter = (N-K)*BGSS;
    CVInd = inter/intra;


%Davies-Bouldin
elseif strcmp(Method,'DB')
    for i=1:K
        ni = length(find(ClassLabs == i));
        Svec(i) = (1/ni)*sumd(i);
    end
    for i=1:K
       for j=i+1:K
            Mscal = distFun(Centres(i,:),Centres(j,:));
            Ratios(i,j) = (Svec(i) + Svec(j))/Mscal;
            Ratios(j,i) = Ratios(i,j);
       end
    end  
    Rvec = max(Ratios); 
    CVInd = sum(Rvec)/K;


%Pakhira-Bandyopadhyay-Maulik
elseif strcmp(Method,'PBM') 
    D = sum(distFun(Data,s));
    if (p == 2 && q == 2)
        inter = max(pdist(Centres).^2)*D;
    elseif (p==1 && q == 1) 
        inter = max(pdist(Centres, 'cityblock'))*D;
    elseif (p==2 && q == 1)
        inter = max(pdist(Centres, 'euclidean'))*D;
    end
    intra = K*ClustErr;
    CVInd = (inter/intra)^2;

%Ray Turi
elseif strcmp(Method,'RT') 
    intra = (1/N)*ClustErr;
    if (p == 2 && q == 2)
        inter = min(pdist(Centres).^2);
    elseif (p==1 && q == 1) 
        inter = min(pdist(Centres, 'cityblock'));
    elseif (p==2 && q == 1)
        inter = min(pdist(Centres, 'euclidean'));
    end
    CVInd = intra/inter;
    

% Wemmert-Gancarski
elseif strcmp(Method,'WG')
    CVInd = 0;
    for i=1:K
        I = find(ClassLabs == i);
        Centres2 = Centres;
        Centres2(i,:) = realmax/(10^6)*ones(1,n); %Enforce not to be closest
        OthClustdists = min(distFun(Data(I,:),Centres2),[],2);
        Ratio = FinalDists(I)./OthClustdists;
        CVInd = CVInd + max(0,length(I) - sum(Ratio));
    end
    CVInd = CVInd/N; 

else
    disp('Requested Internal Cluster Validation Index not supported...')

end
end



