clear;

load('challenge2.mat');
[Ur,Sr,Vr]=svd(Ar);
[Uc,Sc,Vc]=svd(Ac);

% right-hand side B
subplot(1,1,1), imagesc(B), axis image, title('Right-hand side B'); colormap('gray'); 
print right-hand_side_B -dpsc 
fprintf('press any key to continue ...\n'); pause;

% S - singluar values of matrix A=(Ar kron Ac)
sc=diag(Sc);
sr=diag(Sr);
S=sc*sr';

% we are looking for V_i-basis, images, that we need to construct solution 
%
% A = Ar kron Ac
% A = USV', U = Ur kron Uc, S = Sr kron Sc, V = Vr kron Vc
% x = sum [phi_i] * (u_i^T * b / sigma_i) * v_i
% V_i = reshape (v_i)
%
% singluar valuses in variable S are disordered !!!

sp_Sc=sparse(Sc);
sp_Sr=sparse(Sr);
sp_S=kron(sp_Sr,sp_Sc);
diag_S=diag(sp_S);
%%subplot(1,1,1); semilogy(diag_S); title('Singular values are disordered'); 
%%fprintf('press any key to continue ...\n');
%%pause;    
%%
%%fprintf('disordered basis (please wait)\n');
%%
%%for k=[1:1:9,10:10:260,261:1:269,270:10:520,521:260:4941,4941:1:4949,4950:20:5190,5201:3128:77145,78000]
%%    i=k;
%%
%%    c_idx=mod(i-1,260)+1;
%%    r_idx=floor((i-1)/260)+1;
%% 
%%    vr_i=Vr(:,r_idx);
%%    vc_i=Vc(:,c_idx);
%%    v_i=kron(vr_i,vc_i);
%%    V_i=reshape(v_i,260,300);
%%    
%%    subplot(1,2,1); imagesc(V_i); colormap('gray'); title(sprintf('V_{%d}',k)); axis square;
%%    subplot(1,2,2); semilogy(1:78000,diag_S,'b-',k,diag_S(k),'ro'); title(sprintf('{%s}_{%d}','\sigma',k)); axis square;
%%    H=text(k,diag_S(k),sprintf('    %s_{%d} = %f','\sigma',k,diag_S(k))); set(H,'FontSize',8);
%%    
%%    if k<20, pause(0.5); else, pause(0.2); end;
%%end;
%%fprintf('press any key to continue ...\n');
%%pause;    

% sorting singlular values
% variable p contain sorting-permutation 

[diag_S,p]=sort(-diag_S); diag_S=-diag_S;
subplot(1,1,1); semilogy(diag_S); title('Singular values of A');
fprintf('press any key to continue ...\n');
pause;    

fprintf('singular images (please wait)\n');

for k=[1:1:9,10:10:90,100:100:900,1000:1000:9000,10000:10000:70000,78000]
    i=p(k);

    c_idx=mod(i-1,260)+1;
    r_idx=floor((i-1)/260)+1;
 
    vr_i=Vr(:,r_idx);
    vc_i=Vc(:,c_idx);
    v_i=kron(vr_i,vc_i);
    V_i=reshape(v_i,260,300);
    
    subplot(1,2,1); imagesc(V_i); colormap('gray'); title(sprintf('Singular image V_{%d}',k)); axis square;
    subplot(1,2,2); semilogy(1:78000,diag_S,'b-',k,diag_S(k),'ro'); title(sprintf('{%s}_{%d}','\sigma',k)); axis square;
    H=text(k,diag_S(k),sprintf('    %s_{%d} = %f','\sigma',k,diag_S(k))); set(H,'FontSize',8);
    
    if k<20, pause(0.5); else, pause(0.2); end;
end;
fprintf('press any key to continue ...\n');
pause;    

% image reconstruction, TSVD

UTB=Uc'*B*Ur;
uTb=UTB(:); uTb=uTb(p);
subplot(1,1,1), semilogy(1:78000,abs(uTb),'r.',1:78000,diag_S,'b-',1:78000,diag_S(2983)*ones(1,78000),'k-'), title('Singular values of A and projections u_i^Tb'); 
legend('projections u_i^Tb','singlular values \sigma_i','noise level');
fprintf('press any key to continue ...\n');
pause;    

fprintf('summation (please wait)\n');

X_TSVD=zeros(size(B));
for k=[1:49,50:5:195,200:10:290,300:100:2600,2700:10:3190,3200:100:3900,4000:1000:9000,10000:10000:70000,78000]
    i=p(k);

    c_idx=mod(i-1,260)+1;
    r_idx=floor((i-1)/260)+1;
 
    vr_i=Vr(:,r_idx);
    vc_i=Vc(:,c_idx);
    v_i=kron(vr_i,vc_i);
    V_i=reshape(v_i,260,300);
    
    phi=zeros(size(B)); phi(:,:)=S(:,:)>=diag_S(k); Sfilt=phi./S;
    X_TSVD=Vc*((Uc'*B*Ur).*Sfilt)*Vr';

    subplot(1,2,1); imagesc(X_TSVD); colormap('gray'); title(sprintf('\\Sigma_{j=1}^{%d} ( u_j^Tb / \\sigma_j ) V_j',k)); axis square;
    hor_X=[1,78000]; hor_Y=[1,1]*diag_S(k);
    ver_X=[1,1]*k; ver_Y=[1e-12,1];
    subplot(1,2,2); semilogy(1:78000,diag_S,'b-',1:k,diag_S(1:k),'r-',hor_X,hor_Y,'k--',ver_X,ver_Y,'k--',1:78000,diag_S(2983)*ones(1,78000),'k-'); 
    title(sprintf('j = %d',k)); axis square;
    
    pause(0.1);
end;
fprintf('press any key to continue ...\n'); pause;

% naive solution x=inv(A)*b
Xnaive=Vc*((Uc'*B*Ur)./S)*Vr';
subplot(1,1,1), imagesc(Xnaive), axis image, title('Naive solution'); colormap('gray');
print naive_solution -dpsc
fprintf('press any key to end ...\n');
pause;
