function [gene_weight, gene_group, sample_group, res] = ncis (X, N, ncg, ncs, d, T, eps);

[ng,ns]=size(X);

% train weight
expmad=(mad(X',1))';
expmm=max(expmad);
expmad=expmad/expmm;
s=sum(N);
index=find(s==0);
s(index)=1;
Deg=sparse(diag(1./s));
w=(sparse(eye(ng))-d*N*Deg)\expmad*(1-d);
W=sparse(diag(w));
res=[0];
clear Deg;

% randomly initialize clusters
% G=sparse(1:ng,c1(randperm(ng)),ones(1,ng),ng,ncg);
% F=sparse(1:ns,c2(randperm(ns)),ones(1,ns),ns,ncs);

% initialize the clusters using k-means
G=sparse(1:ng,kmeans(X,ncg)',ones(1,ng),ng,ncg);
F=sparse(1:ns,kmeans(X',ncs)',ones(1,ns),ns,ncs);

% co-clustering
for i=1:T
	S=(G'*W*G)\G'*W*X*F/(F'*F);

	A=X'*W*G*S;
	B=S'*G'*W*G*S;
	Aa=abs(A);
	Ba=abs(B);
	F=F.*sqrt((Aa+A+F*(Ba-B))./(Aa-A+F*(Ba+B)+eps));
    
	C=W*X*F*S';
	D=S*F'*F*S';
	Ca=abs(C);
	Da=abs(D);
	G=G.*sqrt((Ca+C+W*G*(Da-D))./(Ca-C+W*G*(Da+D)+eps));
   
	s=sum(F,2);
	iz=find(s<=0);
	ip=find(s>0);
	F(ip,:)=F(ip,:)./repmat(s(ip),1,ncs);
	F(iz,:)=1/ncs;
    
	s=sum(G,2);
	iz=find(s<=0);
	ip=find(s>0);
	G(ip,:)=G(ip,:)./repmat(s(ip),1,ncg);
	G(iz,:)=1/ncg;
    
	res1=trace(X'*W*X-2*X'*W*G*S*F'+F*S'*G'*W'*G*S*F');
	res=[res res1];
	if abs(res1-res(i))<eps
		break;
	end
end

gene_group=zeros(1,ng);
sample_group=zeros(1,ns);

for i=1:ng
	[m index]=max(G(i,:));
	gene_group(i)=index;
end

for i=1:ns
	[m index]=max(F(i,:));
	sample_group(i)=index;
end

gene_weight=w/max(w);
res=res(2:length(res));