function x=MCEJ(N,rho,alpha,k,y)

%rel=1e-9 ;

if nargin <5,
    fprintf('\nPlease supply parameters and data.\n')
    x=0;
    return;
end

	if k<2,
		 fprintf('Please set k>=2.\n')
	  x=0;return;
end

g=ceil(rho*N);

[n,d]=size(y);

p=ones(n,k)./k; % initial probabilities

ind=1;

dd=5;
st=rand(dd,1)+n*n;
dd2=5;
st2=rand(dd2,1)+n*n;
count=0;
injc=0;
inja=0.95;
inje=1;
while ind~=0;

    count=count+1;
    X=zeros(N,n);
    S=zeros(N,1);
    for i=1:N,
        
        % Assign points to clusters
        X(i,:)=genMCE(p);
        
        % Score points
        S(i)=scoreMCE(X(i,:),y,k);
        
    end
    
    % Order Scores
    [S,I]=sort(S);
    
    S(1:g);
    
    % Update p's based on elite samples
        fprintf('t = %d , g_hat_t = %g , S_star_t = %g\n',count,S(g),S(1))
        st(1:(dd-1))=st(2:dd);
        st(dd)=S(1);

w=zeros(n,k);
for i=1:n,
    for j=1:k,
    w(i,j)=w(i,j)+sum(X(I(1:g),i)==j);
end
end

    w=w./g;
    p=alpha*w + (1-alpha)*p;

   if (sum(st==S(1))==dd)&&(S(1)==S(g)), % require that we are truly "stuck"

%            break;
%        end

%      if (abs((st(end)./100-st(end-1)./100)/st(end-1)./100)<rel)&&(count>dd)
        st2(1:(dd2-1))=st2(2:dd2);
        st2(dd2)=S(1);

       if (sum(st2==S(1))==dd2),
            break;
        end

injc=injc+1;

[pmax,pmaxi]=max(p,[],2);
if k>2,
if (1-((inja)^(injc.^inje))*((k-2)/k))<1/(10*N),break;end
pmax=(1-((inja)^(injc.^inje))*((k-2)/k))*pmax;
fprintf('Injection: %g %% of largest probabilities redistributed.\n',(((inja)^(injc.^inje))*((k-2)/k))*100)
else
if (1-((inja)^(injc.^inje))*((k-1)/k))<1/(10*N),break;end
%pmax=pmax;
pmax=(1-((inja)^(injc.^inje))*((k-1)/k))*pmax;
fprintf('Injection2: %g %% of largest probabilities redistributed.\n',(((inja)^(injc.^inje))*((k-1)/k)))
end
pmin=(1-pmax)./(k-1);
p=pmin*ones(1,k);

for j=1:n,
p(j,pmaxi(j))=pmax(j);
end

p=p./(sum(p,2)*ones(1,k));

end
 
end

x=X(I(1),:);
