function [x,count,sz]=NCE(N,g,alpha,k,y,modif,drplot,c,sig0)
 Betta = 0.95; 
 q = 6; 
 
if nargin <5,
    fprintf('\nPlease supply parameters and data.\n')
    x=0;
    return;
end

%g=ceil(rho*N); % use this is want rho instead of g input

[n,d]=size(y);

if nargin<8, % starting means if none supplied
c =   ones(k,1)*min(y)  + (diag(max(y) - min(y))*rand(d,k))';
end

if nargin<9,
sig0=max(max(y) - min(y))*1.0;
end

if nargin<6,modif=0;end
if nargin<7,drplot=0;end

sig=ones(k,d).*sig0;

epsilon=10^-4;

ind=1;

dd=10;
st=rand(dd,1)+n*n;
count=0;

if drplot~=0,
plot(y(:,1),y(:,2),'k.',c(:,1),c(:,2),'ro')
drawnow;
end

while max(max(sig))>epsilon,
    count=count+1;
    B_Mod = Betta - Betta*(1-1/count)^q;
    C=zeros(N,k,d);
    S=zeros(N,1);
    for i=1:N,
        % Assign points to clusters
        ccc=genNCE(c,sig);
        X(i,:,:)=ccc;
        % Score points
        S(i)=scoreNCE(ccc,y,k);
    end
    
    % Order Scores
    [S,I]=sort(S);
    % Update p's based on elite samples
    
        fprintf('t = %d , g_hat_t = %g , S_star_t = %g, maximum sigma = %g\n',count,S(g),S(1),max(max(sig)))
        st(1:(dd-1))=st(2:dd);
        st(dd)=floor(S(1)*100);
       if (sum(st==floor(S(1)*100))==dd),
            break;
        end
for i=1:k,
    mu(i,:)=mean(X(I(1:g),i,:),1);
end

co=zeros(k,d);
for i=1:k,
    for j=1:d,
        co(i,j)=std(X(I(1:g),i,j),1);
    end
end
		  if modif==0,      
        c=alpha*mu+(1-alpha)*c;
        sig=alpha*co+(1-alpha)*sig;
else
        c=alpha*mu+(1-alpha)*c;      
        sig=B_Mod*co+(1-B_Mod)*sig;
end

		  if drplot~=0,
plot(y(:,1),y(:,2),'k.',c(:,1),c(:,2),'ro')
drawnow;
end

end

for i=1:k,
    x(i,:)=X(I(1),i,:);
end
sz=S(1);
