function clust_labels = cluster_data_em(data2clust, CantClusters, iter_times)

if nargin < 3 || isempty(iter_times)
    iter_times = 1;
end
if nargin < 2 || isempty(CantClusters)
    CantClusters = 5;
end

cant_ejemplos = getsize(data2clust,1);

std_train = std(+data2clust);
data2clust = setdata(data2clust, bsxfun(@rdivide, +data2clust, std_train ));

clustered_labels = repmat('0', cant_ejemplos, iter_times);
jj = 1;
%         for CantClusters = 2:5
for jj = 1:iter_times

    %EMclust
    ii = 0;
    bContinuar = true;
    while(bContinuar && ii < CantClusters)
        try
            dbclear if caught error            
            clust_labels = em_loop(data2clust, CantClusters-ii );
            dbstop if caught error            
            bContinuar = false;
        catch ME
            if(strcmpi(ME.message, 'Not possible to find desired number of components'))
                ii = ii + 1;
            else
                rethrow(ME)
            end
        end
    end
    clustered_labels(:,jj) = char(97+clust_labels);
%             jj = jj + 1;
end

bClustered = true;
while(bClustered)
    %analizo las coincidencias en las distintas iteraciones respecto data2clust
    %los etiquetados.
    [~, sort_idx] = sort( cellstr(clustered_labels) );
    [all_clusters, aux_location] = unique(clustered_labels(sort_idx,:), 'rows', 'first');

    aux_location = [colvec(aux_location); cant_ejemplos+1];
    cluster_sizes = diff(aux_location);
    [~, clust_sort_idx] = sort(cluster_sizes, 'descend');

    %agrupo todos los subclusters data2clust una dist maxima max_distance
    max_distance = round(0.2*iter_times); %distancia maxima para considerarse un cluster
    cant_clusters = size(all_clusters,1);
    cant_iter = size(all_clusters,2);
    aux_1 = repmat( all_clusters, cant_clusters, 1);
    aux_idx = colvec(repmat(1:cant_clusters,cant_clusters,1));
    aux_2 = all_clusters(aux_idx,:);
    distances = reshape( sum(aux_1 ~= aux_2,2), cant_clusters, cant_clusters );
    remaining_clusters = 1:cant_clusters;

    %clusterizo igual todo los clusters iguales hasta max_distance
    for ii = rowvec(clust_sort_idx)
        bClustered = false;
        cluster2fusion_idx = find(distances(:,ii) <= max_distance);
        cluster2fusion_idx = cluster2fusion_idx(cluster2fusion_idx < ii | cluster2fusion_idx > ii );
        for jj = 1:length(cluster2fusion_idx)
            aux_idx = find(strcmpi(cellstr(clustered_labels), cellstr(all_clusters(cluster2fusion_idx(jj),:))));
            if( ~isempty(aux_idx) )
                bClustered = true;
                remaining_clusters( remaining_clusters == cluster2fusion_idx(jj)) = [];
                clustered_labels( aux_idx ,:) = repmat(all_clusters(ii,:), length(aux_idx), 1 );
                all_clusters(cluster2fusion_idx(jj),:) = all_clusters(ii,:);
            end
        end
        if(bClustered)
            remaining_clusters( remaining_clusters == ii ) = [];
            %fuerzo el recalclo de distancias, para que no se junte todo y solo
            %los clusters grandes se coman data2clust los mas chicos.
            break
        end
    end
end

%REanalizo las coincidencias en las distintas iteraciones respecto data2clust
%los etiquetados.
[~, sort_idx] = sort( cellstr(clustered_labels) );
[all_clusters, aux_location] = unique(clustered_labels(sort_idx,:), 'rows', 'first');

aux_location = [colvec(aux_location); cant_ejemplos+1];
cluster_sizes = diff(aux_location);

cant_clusters = size(all_clusters,1);

clust_labels = nan(cant_ejemplos,1);
for ii = 1:cant_clusters
    clust_labels(sort_idx(aux_location(ii):(aux_location(ii+1)-1))) = ii;
end

if(any(isnan(clust_labels)))
   error() 
end



function new_lab = em_loop(data2clust, CantClusters )

n_ini		= 500;			% Maximum size of subset to use for initialisation.

cant_ejemplos = size(data2clust,1);

not_found = true;
iter = 0;
while(not_found)
    
    % try to find an initialisation with all class sizes > 1

    if (cant_ejemplos > n_ini)       						% ... data2clust random subset of data2clust.
        %sampleo uniformemente el feature-space
        a_ini = data2clust;
        a_ini = sum(abs(a_ini),2);
        [~, sort_idx ] = sort(a_ini);
        a_ini = data2clust(sort_idx(randsample(cant_ejemplos, n_ini)),:);
    else
        a_ini = data2clust;								% ... the entire set data2clust.
    end

    iter = iter + 1;
    if iter > 3
        error('Not possible to find desired number of components')
    end
    % 50 trials
    assign = zeros(1,CantClusters-1);
    std_a_ini = std(a_ini);
    while( length(unique(assign)) ~= CantClusters )
        % add some noise to data to avoid problems
        a_ini = a_ini.*(ones(size(a_ini))+ bsxfun(@times, max( [repmat(0.001,size(std_a_ini)); 0.01*std_a_ini] ), randn(size(a_ini))));
        %construct the distance matrix
        % The order of operations below is good for the accuracy.
        dist_mat = ones(n_ini,1)*sum(a_ini'.^2,1);
        dist_mat = dist_mat + sum(a_ini.^2,2)*ones(1,n_ini);
        dist_mat = dist_mat - 2 .* (a_ini)*(a_ini)';

        % Check for a numerical inaccuracy. 
        dist_mat(dist_mat<0) = zeros(size(J));          % dist_mat should be nonnegative.

        % take care of symmetric distance matrix
        dist_mat(1:n_ini+1:n_ini*n_ini) = zeros(1,n_ini);
        
        try
            assign  = kcentres(dist_mat,CantClusters,50);
        catch ME
            error('Not possible to find desired number of components')
        end
    end
    % Train initial classifier on labels generated by KCENTRES and find
    % initial hard labels. Use NMC instead of W_CLUST to make sure that we 
    % always have enough data to estimate the parameters.
    a_ini = dataset(a_ini,assign); 
    a_ini = setprior(a_ini,getprior(a_ini,0));
    dist_mat = data2clust*(a_ini*nmc);
    
    new_lab = dist_mat*labeld;
    cs_new_lab = classsizes(dataset(dist_mat,new_lab));
    if length(cs_new_lab) == CantClusters && all( cs_new_lab > 1)
        not_found = 0;
    end
end

lablist_org = [];

% Ready for the work.
iter = 0;
change = 1;
lab = ones(cant_ejemplos,1);

while (change ~= 0 && any(lab ~= new_lab)) % EM loop, run until labeling is stable.

    prwaitbar(100,100-100*exp(-iter/10));
    data2clust = setlabels(data2clust,new_lab); 		% 0. Set labels and store old labels.
    data2clust = setprior(data2clust,getprior(data2clust,0));%    Set priors to class frequencies
    lab = new_lab;								% 
    data2clust = remclass(data2clust,1);            %    demand class sizes > 2 objects
    iter = 0;

    while getsize(data2clust,3) < CantClusters        %    increase number of classes if necessary

        iter = iter + 1;

        if iter > 5
            error('Not possible to find desired number of components')
        end
        laba = getlablist(data2clust);
        labmax = max(laba);
        CantClusters = classsizes(data2clust);
        [Nmax,cmax] = max(CantClusters);        % find largest class
        aa = seldat(data2clust,cmax);         % select just that one

        std_aa = std(+aa);
        bContinue = true;
        while( bContinue )
            try
                new_lab_aa = kmeans(aa .* (ones(size(+aa))+ bsxfun(@times, max( [repmat(0.001,size(+aa)); 0.01*std_aa] ), randn(size(+aa)))) ,2);   % split it by kmeans
                bContinue = false;
            catch ME
            %                         if( strcmpi(ME.message, 'kcentres fails as some objects are identical: add some noise'))
            %                             std_aa = std_aa * 2;
            %                         else
            %                             rethrow(ME)
            %                         end
                error('Not possible to find desired number of components')
            end
        end

        N1 = sum(new_lab_aa == 1);   
        N2 = sum(new_lab_aa == 2);
        if (N1 > 1 & N2 > 1) % use it if both classes have more than one sample
            J = findlabels(data2clust,laba(cmax,:));
            data2clust = setlabels(data2clust,new_lab_aa + labmax,J);
        end
    end

    std_a = std(+data2clust);
    w_em = (data2clust .* (ones(size(+data2clust))+ bsxfun(@times, max( [repmat(0.001,size(+data2clust)); 0.01*std_a] ), randn(size(+data2clust))))) * w_clust;             % 1. Compute classifier, crisp output.
    b = data2clust*w_em;               		% 2. Classify training samples.
    new_lab = labeld(b);      		% 3. Insert classification as new labels.
    iter = iter+1;             %DXD Added also the iter for the crisp labels
    
    if iter > 50
        change = 0;
    end

    %para ver la velocidad de convergencia.
    %             disp(sum(lab ~= new_lab))

end


if ~isempty(lablist_org) % substitute original labels if desired
    new_lab = lablist_org(new_lab);
    wlab = getlabels(w_em);
    wlab = lablist_org(wlab);
    w_em = setlabels(w_em,wlab);
end


function labels = kcentres( dist_mat, CantClusters, itern)

dmax = max(max(dist_mat));
dopt = inf;

for tri = 1:itern

    M = randperm(cant_ejemplos); M = M(1:CantClusters);		% Random initializations

    J = zeros(1,CantClusters); 					    	  % Center objects to be found.

    % Iterate until J == M. See below.
    while true
        
        [dm,I] = min(dist_mat(M,:));

        % Find CantClusters centers. 
        for ii = 1:CantClusters                   
            JJ = find(I==ii); 
            if (isempty(JJ))
                %JJ can be empty if two or more objects are in the same position of 
                % feature space in dataset
                J(ii) = 0;									
            else
                % Find objects closest to the object M(ii)
                [dummy,j,dm] = kcentres(dist_mat(JJ,JJ),1,1);
                J(ii) = JJ(j);
            end
        end

        J(find(J==0)) = [];
        CantClusters = length(J);
        if CantClusters == 0
            error('kcentres fails as some objects are identical: add some noise')
        end
        if (length(M) == CantClusters) & (all(M == J))
            % CantClusters centers are found.
            [dmin,labs] = min(dist_mat(J,:));
            dmin = max(dmin);
            break;
        end
        M = J;
    end

    % Store the optimal centers in JOPT.

    if (dmin <= dopt)
        dopt = dmin;
        labels = labs';
        Jopt = J;
    end

end

prwaitbar(0)

% Determine the best centers over the N trials.
dm = zeros(1,CantClusters);   
for ii=1:CantClusters
    L = find(labels==ii);
    dm(ii) = max(dist_mat(Jopt(ii),L));
end