function [x,S] = omp(A,y,Parameters)

% Orthogonal Matching Pursuit (OMP) with various options
%
% A..........dictionary (normally, l2-normalized columns of dictionary is assumed)
% y..........signal
% x..........sparse solution found by OMP
%
% Parameters. (structure)
% StopMaxerror....maximum error - stopping criteria
% StopMaxiter.....maximum number of iterations - stopping criteria
% (Iterations stop if one of the criteria is met.
%  If one of them is not present, OMP is driven by the other one.)
%
% NormalizeDict...boolean; should the dictionary be normalized before OMP? 
% 
% % ompDisplay.update - display each .update iteration or disable (0)
% % ompDisplay.pause 
% %     - 'off'    - disable option 
% %     - seconds  - pause after each OMP iteration 
% %     - 'key'    - wait for key press
% % ompDisplay.originalSegment  -   reference segment
% % ompDisplay.originalDict     -   original dictionary
% % ompDisplay.missSel          -   assign the missing signal support
%
% Example:
%   A = randn(5,10);
%   y = randn(5,1);
%   param.StopMaxerror = 10e-3;
%   [x, S] = omp(normalize_dict(A),y,param)

% (c) 2012-2013 Pavel Rajmic, Jan Spirik, Marie Dankova, Vaclav Mach,
% Brno University of Technology

% Last revision 2013-05-02 Pavel Rajmic


%% Check and defaults
if nargin < 3
    error('Not enough input parameters')
end

if ~isstruct(Parameters)
    error('Parameters is not a field structure')
end

if ~isfield(Parameters,'StopMaxerror')
    Parameters.StopMaxerror = 0; %10e-5;  %default
end

if ~isfield(Parameters,'StopMaxiter')
    Parameters.StopMaxiter = NaN;  %default
end

if ~isfield(Parameters,'NormalizeDict')
    Parameters.NormalizeDict = false;  %default
end

n = size(A,2); %number of columns, i.e. atoms

%% Normalization of dictionary 
if Parameters.NormalizeDict == true
%     l2norms = sqrt(sum(A.^2)); %l2-norms of columns
%     if any( abs(l2norms-ones(1,n)) > 10e-5 ) % if norms are too far from 1
%         disp(' ')
        disp('Normalizing the columns of dictionary')
        [A, origl2norms] = normalize_dict(A);
%     end
%     clear l2norms
end

%% Initialization
x = zeros(n,1); %initial solution sloupcovy vektor plny nul o velikosti poctu sloupcu slovniku
r = y; %initial residual r je vstupni sign?l
S = []; %initial support prazdna mnozina
epsilon = norm(r); %initial error velikost y
k = 1; % iteration counter
%disp(' ')
AA=A;   %%helping dictionary for counting jo

%% Main loop
while (k <= Parameters.StopMaxiter) && (epsilon > Parameters.StopMaxerror)
%    disp(['Starting iteration ' num2str(k) ', previous residual norm ' num2str(epsilon)]);
    
    % Computing the errors for atoms
    e = abs(AA'*r); %%%%!!!zbytecne ze vsech sloupcu ackoliv u tech co jsou uz v S zde vyjde nula
    
    % Finding the error minimizer  
    jo = find(e == max(e),1,'first'); %najde prvni nejvesti?
    
    %Update support
    S = [S, jo]; %ulozi do s      

    %Update provisional solution (smallest l2-norm using support S)
    A_S = A(:,S);  %choose only atoms corresponding to S
    x(S) = pinv(A_S)*y;

    %Update residual    
    r = y - A_S*x(S);
    
    %Update helping dictionary
    AA(:,jo) = 0;

    %Norm of the residual
    epsilon = norm(r);
    
    %Iteration counter
    k = k+1;
    
%     %Display update
%     if nargin==4
%         if (ompDisplay.update~=0) && (mod(k,ompDisplay.update)==0)
%             ompMissSup = find(ompDisplay.missSel);  %support for displaying just missing samples
%             %ompMissSup = [1:1:length(ompDisplay.missSel)]; %support is of length of the segLen
%             figure(1)
%             plot(ompDisplay.originalSegment(ompMissSup), 'b');
%             hold on
%             x_tmp = ompDisplay.originalDict*x;
%             plot(x_tmp(ompMissSup), 'r');
%             legend('Original signal','Reconstructed signal');
%             title('Temporal reconstruction');
%             hold off
%                 if strcmp(ompDisplay.pause,'key');
%                     key=waitforbuttonpress;
%                 else
%                     pause(ompDisplay.pause);
%                 end
%         end
%     end
end

%Display which criterium was fulfilled
%disp('  ')
if (epsilon < Parameters.StopMaxerror)
 %   disp(['Stopping criterium ''StopMaxerror'' fulfilled, residual norm is ' num2str(epsilon)]);
%elseif  ((k-1) == Parameters.StopMaxiter)
%    disp(['Stopping criterium ''StopMaxiter'' fulfilled, residual norm is ' num2str(epsilon)]);
end

%% Sort the indexes of support
if nargout > 1
    S = sort(S)';
end