function [BestLbs, BestCntrs, sumd, dist, iterArray,timeArray] = Kmeans2pp(Data,k,Replicates,distance,convTh,maxIter,tol,spamitmax)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% K-medians, K-means, K-spatialmedians clustering (Tommi K., Jun 12 2017)
% modified by Joonas H, 26.6.2017
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[N,n] = size(Data);
nrep = Replicates;
if(nargin<7)
    tol = 1.d-5;
    spamitmax = 100;
end

disp('Clustering started...')
Besterror = realmax;
iterArray = zeros(Replicates,1);
timeArray = zeros(Replicates,1);
I = zeros(k,1);
for rep=1:nrep
    tic
    %Generate initial prototypes using k-means++ like of method
    I(1) = randi(N,1); %Random index for first prototype
    Dist = inf(N,1);
    if(strcmp(distance,'seuclidean'))
        for i = 2:k
                Dist = min(Dist,pdist2(Data,Data(I(i-1),:)).^2);
                Distp = cumsum(Dist);
                if (Distp(end) > 0)
                    Distp = Distp/Distp(end);
                    ind = find(Distp > rand, 1);
                    I(i) = ind;
                else %Distances to all remaining observations equal to zero
                    disp('Not enough distinct values for initial prototypes')
                    I(i:k) = randi(N,k-i+1);
                end
        end
    else
        for i = 2:k
                Dist = min(Dist,pdist2(Data,Data(I(i-1),:),distance));
                Distp = cumsum(Dist);
                if (Distp(end) > 0)
                    Distp = Distp/Distp(end);
                    ind = find(Distp > rand, 1);
                    I(i) = ind;
                else %Distances to all remaining observations equal to zero
                    disp('Not enough distinct values for initial prototypes')
                    I(i:k) = randi(N,k-i+1);
                end
        end 
    end
    
    InitCntrs = Data(I,:);
    [L,C,iters] = kmeans3(Data,k,InitCntrs,distance,convTh,maxIter,tol,spamitmax);
    t = toc;
    if(iters < 0)
        timeArray(rep) = -1;
        iterArray(rep) = -1;
        Clusterror = realmax; 
    else
        timeArray(rep) = t;
        iterArray(rep) = iters;
        if(strcmp(distance,'seuclidean'))
            D = pdist2(Data,C).^2;
        elseif(strcmp(distance,'cityblock'))
            D = pdist2(Data,C,distance);
        elseif(strcmp(distance, 'euclidean'))
            D = pdist2(Data,C,distance);
        end
        Ds = min(D,[],2);
        for i=1:k
            SumD(i) = sum(Ds(L==i));
        end
        Clusterror = sum(SumD);
    end
    
    if Clusterror < Besterror
        Besterror = Clusterror;
        BestCntrs = C;
        BestLbs = L;
        sumd = SumD;
        dist = D;
    end
end