% generate matrix with clusters of size m1-by-n1 and m2-by-n2
m1=200; n1=100;
m2=100; n2=200;
% treshold for sparsity
tresh=0.6;
% treshold for noise
treshn=0.4;

% only incidence matrix
% a11=rand(m1,n1)<tresh;
% a22=rand(m2,n2)<tresh;
% noise=rand(m1+m2,n1+n2)<treshn;

% incidence matrix + weights
a11=(rand(m1,n1)<tresh).*rand(m1,n1);
mask=ones(m1,n1);
for i=1:m1, for j=1:n1,
    if (i-100)^2/10000+(j-50)^2/2500>1,
       mask(i,j)=0;
    end
end
end
a11=a11.*mask;
a22=(rand(m2,n2)<tresh).*rand(m2,n2);
noise=(rand(m1+m2,n1+n2)<treshn).*rand(m1+m2,n1+n2);

a=[a11 zeros(m1,n2); zeros(m2,n1) a22];
a=a+noise;

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

% random premutation
[q,r,pm]=qr(rand(m1+m2));
[q,r,pn]=qr(rand(n1+n2));

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,'.')

% Check Lanczos

% Compute U on A*A'
[Yu, ThetaU, Qu,Tu]=Lanczos(ap*ap',3);
[minu,iminu]=sort(diag(ThetaU),'descend');
[usortL,uindexL]=sort(Yu(:,iminu(2)));

% Compute V on A'*A
[Yv, ThetaV, Qv,Tv]=Lanczos(ap'*ap,3);
[minv,iminv]=sort(diag(ThetaV),'descend');
[vsortL,vindexL]=sort(Yv(:,iminv(2)));

apsortL=ap(uindexL,vindexL);
figure(4);
[i,j,s]=find(apsortL); plot(j,i,'.')
