% generate matrix with k clusters of size m(k)-by-n(k)
k=3;
m(1)=200; n(1)=100;
m(2)=100; n(2)=200;
m(3)=100; n(3)=100;

% treshold for sparsity
tresh=0.8;
% treshold for noise
treshn=0.2;
% put weigths 1=yes, 0=no
weights=0;

% incidence matrix
tmp=rand(m(1),n(1))<tresh;
mask=ones(m(1),n(1));
% for i=1:m(1), for j=1:n(1),
%     if (i-100)^2/10000+(j-50)^2/2500>1,
%        mask(i,j)=0;
%     end
% end
% end
tmp=tmp.*mask;
if weights==1, tmp=tmp.*rand(m(1),n(1)); end
a(1:m(1),1:n(1))=tmp;
for i=2:k  
   tmp=rand(m(i),n(i))<tresh;
   if weights==1, tmp=tmp.*rand(m(i),n(i)); end 
   a(sum(m(1:i-1))+1:sum(m(1:i)),sum(n(1:i-1))+1:sum(n(1:i)))=tmp;
end

% additional block
%   a(301:400,1:100)=rand(100,100)<0.5;


noise=rand(sum(m),sum(n))<treshn;
a=a+noise;

% figure(1);
[i,j,s]=find(a);
plot(j,i,'.')

% random premutation
[q,r,pm]=qr(rand(sum(m)));
[q,r,pn]=qr(rand(sum(n)));

ap=pm*a*pn;
% figure(2); 
[i,j,s]=find(ap);
plot(j,i,'.')

% clustering
[u,s,v]=svd(ap);
[usort,uindex]=sort(u(:,2));
[vsort,vindex]=sort(v(:,2));
apsort=ap(uindex,vindex);
% figure(3);
[i,j,s]=find(apsort);
plot(j,i,'.')

