function y = commGSP_reconstruct(A, type, BW, S, x)
% Reconstructs a graph signal from the provided subsampled version (on
%   a subgraph)
% INPUT:
%   A - (symmetric) adjacency matrix
%   type - string defining type of the signal to reconstruct:
%          'smooth', 'non-smooth' (Laplacian-based),
%          'modular', 'anti-modular' (modularity-based)
%   BW - number of eigenvectors kept in the signal bandwidth
%   S - diagonal matrix defining subgraph on which the signal values
%       are used to reconstruct the signal: entry (i,i) equal
%       to 1 if i^th node belongs to the subgraph, 0 otherwise
%   x - signal to reconstruct (length equal to size of A). Unknown values
%       need to be set to 0.

% OUTPUT
%   y - reconstructed graph signal

N=size(A,1);
if strcmp(type, 'smooth') || strcmp(type, 'non-smooth')
    [U, ~] = commGSP_eig_implicit(A, 'laplacian');
else
    if strcmp(type, 'modular') || strcmp(type, 'anti-modular')
        [U, ~] = commGSP_eig_implicit(A, 'modularity');
    end
end

if strcmp(type, 'smooth') || strcmp(type, 'modular')
    w=1:BW;
else
    if strcmp(type, 'non-smooth') || strcmp(type, 'anti-modular')
        w=(N-BW+1):N;
    end
end
W = sparse(w,w,ones(length(BW),1),N,N);


B=sparse(U*W*U');
[V,PSI] = eigs(sparse(B*S*B),N); V=sparse(V);
eigval=diag(sparse(PSI));
eigval(find(abs(eigval)<1e-9))=0;
PSI=diag(eigval);

tmp=diag(PSI.*W); ind=find(tmp~=0);
rec=sparse(size(tmp,1),size(tmp,2));
rec(ind)=1./eigval(ind);
y=V*diag(rec)*V'*x;

y(find(diag(S)==1))=x(find(diag(S)==1));
end