function [xk,supp] = INPMAT( A,y,mu,spar)
%---------------------------------------------------------------
% INPMAT implementation in MATLAB
%---------------------------------------------------------------
% Usage :
% [xk,supp] = INPMAT( A,y,mu,spar,maxit,tol)
%---------------------------------------------------------------
% INPUT :
% A     :  Measurement matrix.
% y     :  The measurement/observed vector.
%---------------------------------------------------------------
% Optionals:
% mu    :  The Tikhonov Regularization parameter
%        - Use 0 to set the default (least square solution).
% spar  : The sparsity number of the input signal.
%---------------------------------------------------------------
% OUTPUT :
% xk        : Recovered  sparse signal
% supp      : indices of nonzero elements
%---------------------------------------------------------------
% Date created: 27th August, 2016
% COPYRIGHT 2016 Asadi-Kangarshahi E, Esmaeili A, Marvasti F.
%---------------------------------------------------------------

    tol       =     1e-30;
% tol   :  Stop iterations if norm(y - Ax)^2 <= tol, where tol is the
% stopping tolerance.

    maxit     =     0.001;
    
% maxit :  The max. no. of iterations.

n=size(A,2);% length of the desired vector
if mu == 0
    B        =       pinv(A);
else
    B        =       pinv(mu*eye(n,n) + A'*A)*A';
end
    z        =       B*y; %initiation point derived from tikhonov regularized problem

threshold = max(abs(z));

%% parameters def: 
i       =       0; % iteration counter
xk = zeros(size(z)); % the vector holding updates of the desired signal
T = zeros(n,1); % support of the signal at each iteration
lastT = zeros(n,1);
prnumber=0; % convergence control counter
k = 0; % size of support
% while(k<=spar) % stopping criterion

Proj = B*A;

%% We can set two stopping criteria: 
% Case 1: We (approximately/precisely) know the sparsity level of the desired signal.
% Case 2: We do not know the sparsity level of the desired signal and use ...
... stopping criterion based on the squared residual,...
... i.e. (norm(y-A(:,Supp)*x(Supp)))^2>tol would be the stopping criterion.
%-------------------------------------------------------------------------------------
indis=[1:size(A,2)];
% while (norm(y-A(:,indis)*xk(indis)))>tol
while ( k <= spar )
    %---------------------------------------------------------------
    % Convergence control in scarce cases; the convergence could be speeded up
    % while(prnumber<=500)
    prnumber=prnumber+1;
    if prnumber == 100
        threshold = threshold * 0.9;
        prnumber = 0;
    elseif threshold<1e-14
        break
    end
    %     sum(T)
    %---------------------------------------------------------------
%% INPMAT Logistics:
    k;
    i = i+1;
    T = (abs(xk) >= threshold);
    k = sum(T);
    [~,indis] = sort(T,'descend');
    indis = indis(1:k);
    xkp = pinv(A(:,indis))*y;
    xk = zeros(n,1);
    xk(indis) = xkp;
    projected = Proj(:,indis)*xkp;
    xk = z + xk - projected;
    % threshold update:
    if lastT == T
        threshold = max(abs((1-T).*xk));
    end
    %----------------------------------------------------------------
    lastT = T;
end
spar = sum(T);
[vals,index]=sort(xk,'descend');
index = index(1:spar);
xkp = pinv(A(:,index))*y;
xk = zeros(n,1);
xk(index) = xkp;
supp = index;
% keeping the spars-largest components
end