function [X , S , r]= OMP(A,b,s,ep)
%OMP Summary of this function goes here
%   Detailed explanation goes here

if(nargin<4)
    ep = 1e-12;
end

[m,N] = size(A);
itermax = s;
%%

r = b;
X = zeros(N,1);

S = [];

for iter=1:itermax
    %step1
    mu = abs(A'*r);
   
    [q,lambda] = max(mu);
    
    %step2
    S = [S,lambda]; % display(Lda); %Lda = union(Lda,lambda);
    A_S = A(:,S);
    
    
    %step3, pseudo inverse
    x = lscov(A_S,b);
    
    %step4, residual
    r = b - A_S*x;
    
    %dump
    res = norm(r,2);
    %display(sprintf('iter=%d, res=%g', iter, res))
    
    X(S) = x;
%    if(res<ep) % || length(S)==s)
%        return;
%    end
end

S = S';

