%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Demo Codes For:
%
%
%  "A New Algorithm for Sparsity Promoting Interpolation of Graph Signals" 
%           
%
%
%
% Written By: Maryam Fallah
%             Mahdi Boloursaz Mashhadi
%
% Affiliation:
%          Advanced Communications Research Institute (ACRI)
% Electrical Engineering Department, Sharif University of Technology
%                              Tehran, Iran
%
% For any problems, contact me at maryamfallah1994@gmail.com
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                             Initialization
close all
clear all
clc
%%                          Simulation Parameters
p = 0.45;                % Random (Missing) Sampling Rate
alpha = 0.03;            % Exponential Threshold Reduction Factor
IterNum = 200;           % Number of Algorithm Iterations
landa=1;                 % Relaxation parameter
%%                             Producing the random graph
N=100;                              % Number of nodes of the graph
a=round(rand(N,N));                 % Produces a random adjacent matrix
a(logical(eye(size(a))))=0;
A=floor((a+a')/2);                  % Generates a symmetrical adjacent matrix
B=zeros(N,1);                       % Finds the degree of each node to generate the degree matrix
for i=1:N
    
    bb=length(find(A(i,:))) ;
    B(i,1)=bb;
    
end

d=diag(B,0);
L=d-A;                              % Laplacian matrix
[U,D]=eig(L);                       % Generating the eigenvectors(V) and eigenvalues(D)
s=rand(1,N)<p;                      % Generates a random sampling vector
S=diag(s,0);                        % Diagnol matrix with ones for sampled vertices and zeros for not sampled vertices
fk=zeros(N,1);                      % Reconstructed signal
%%                            IMATGI Reconstruction
c=0;
for koverN=10:60                             % Reconstructs the signal for different sparsity factors (koverN)
    for j=1:100
        Supp=rand(1,N)<(koverN*.01);         % Produces a sparsity factor
        f=U*(Supp'.*(U'*rand(N,1)));         % Random sparse signal assigned to each vertex
        betta=p*max(U'*f);                   % Threshold coefficient
        measure=s'.*f;                       % Signifies available graph signals
        
        for k=1:IterNum                      % Iteration algorithm
            
            Th=betta*exp(-alpha*(k-1));      % Adaptive threshold
            Txs=U'*fk;                       % Transfer into GFT domain
            Tc =abs(Txs)>Th;                 % Thresholding in the GFT domain
            g=U*(Tc.*Txs);                   % Transfer back to the vertex domain
            fk=landa*measure+g-landa*S*g;    % Iterative update
            
            error(k)=norm(f-fk);        
            
        end
        error2=norm(f-fk);              
        psig=norm(f);                   
        SNR(j)=10*log10(psig./error2);  % Computes Reconstruction SNR
        
    end
    
    if mod(koverN,10)==0
        c=c+1;
        Progress_percentage=16.6*c
    end
    SNR2(koverN)=mean(SNR);
end
%%                             Simulation Results

plot(10:60,SNR2(10:end),'k-*');        % Plots SNR versus sparsity factor
xlabel('Sparsity Factor (%)')
ylabel('SNR (dB)')
