perform_omp(D,Y,options)

function X = perform_omp(D,Y,options)

% perform_omp - perform orthogonal matching pursuit
%
%   X = perform_omp(D,Y,options);
%
%   D is the dictionary of size (n,p) of p atoms
%   Y are the m vectors to decompose of size (n,m)
%   X are the m coefficients of the decomposition of size (p,m).
%
%   Orthogonal matching pursuit is a greedy procedure to minimise
%       |x|_0   under the constraint that D*x=y.
%   (maximum sparsity).
%   You can stop the iteration after |x|_0=k  by setting options.nbr_max_atoms=k
%   You can set relative tolerance options.tol so that iterations stops when
%       |Y-D*X| < tol*|Y|
%
%	This code calls the fast mex version of Antoine Grolleau
%	when possible.
%
%   Copyright (c) 2006 Gabriel Peyre


if isfield(options, 'sparse_coder') && strcmp(options.sparse_coder, 'mp')
    X = perform_mp(D,Y,options);
    return;
end

options.null = 0;
nbr_max_atoms = getoptions(options, 'nbr_max_atoms', 50);

tol = getoptions(options, 'tol', 1e-3);
use_slow_code = getoptions(options, 'use_slow_code', 0);
verb = getoptions(options, 'verb', 1);

% P : number of signals, n dimension
[n,P]=size(Y);
% K : size of the dictionary
[n,K]=size(D);


if isfield(options, 'use_mex') && options.use_mex==1 && exist('perform_omp_mex')==3
    % use fast mex interface
    X = zeros(K,P);
    for k=1:P
        X(:,k) = perform_omp_mex(D,Y(:,k),nbr_max_atoms);
    end
    return;
end

if use_slow_code
    X = perform_omp_other(Y,D,tol,nbr_max_atoms);
    return;
end

for k=1:1:P,
    if P>20 && verb==1
        progressbar(k,P);
    end
    a=[];
    x = Y(:,k);
    r = x; % residual
    I = zeros(nbr_max_atoms,1); % index of the chosen atoms
    % norm of the original signal
    e0 = sqrt( sum( r.^2 ) );
    
    for j=1:1:nbr_max_atoms,
        proj = D'*r;
        pos = find(abs(proj)==max(abs(proj)));
        pos = pos(1);
        
        I(j) = pos;        
        a=pinv(D(:,I(1:j)))*x;
        r=x-D(:,I(1:j))*a;

        % compute error
        e = sqrt( sum( r.^2 ) );
        if e/e0 < tol
            break;
        end

    end;
    
    
    temp=zeros(K,1);
    temp(I(1:length(a)))=a;
    X(:,k)=sparse(temp);
end;
return;


function [D,Di,Q,beta,c] = perform_omp_other(f,D,tol,No,ind)

% perform_omp - Optimized Orthogonal Matching Pursuit
%
% It creates an atomic decomposition of a signal using OOMP method. You can choose a
% tolerance, the number of atoms to take in or an initial subspace to influence the OOMP
% algorithm. Non-selected atoms subtracted by their component in the chosen space are also
% available.
%
% Usage: 
%        x = perform_omp_other(f,D,tol,No,ind);
%
% Inputs:
%   f    analyzing signal of size (n,s) where s is the number of signals.
%   D    dictionary of normalized atoms of size (n,d) where d is the number
%       of atoms.
%   tol  desired distance between f and its approximation
%        the routine will stop if norm(f'-Dsub*(f*beta)')*sqrt(delta)<tol
%        where delta=1/nbr_max_atoms, nbr_max_atoms is number of points in a sample
%   No   (optional) maximal number of atoms to choose, if the number of chosen
%        atoms equals to No, routine will stop (default No=size(D,2)
%   ind  (optional) indices determining  the initial subspace,
%
% Outputs:
%   x    coefficients of the decomposition such that f \approx D*x of size (p,s)
%
% References:
%    nbr_max_atoms. Rebollo-Neira and D. Lowe, "Optimized Orthogonal Matching Pursuit Approach",
%      IEEE Signal Processing Letters, Vol(9,4), 137-140, (2002).
%
% See also OMPF.

%        [Dnew,Di]=oompf(f,D,tol);
%        [Dnew,Di,Q,beta,c]=oompf(f,D,tol,No,ind);

%   D    the dictionary D rearranged according to the selection process
%        D(:,1:k) contains the atoms chosen into the atomic decomposition
%   Di   indices of atoms in new D written w.r.t the original D
%   Q    Q(:,1:k) contains orthonormal functions spanning new D(:,1:k)
%        Q(:,k+1:N) contains new D(:,k+1:N) subtracted by the projection
%        onto the space generated by  Q(:,1:k)  (resp. D(:,1:k))
%   beta 'k' biorthogonal functions corresponding to new D(:,1:k)
%   c    'k' coefficients of the atomic decomposition

% See http://www.ncrg.aston.ac.uk/Projects/BiOrthog/ for more details

% verbosity
verb = 0;

name='OOMPF';  %name of routine
% get inputs and setup parameters
[nbr_max_atoms,N]=size(D);
beta=[];
Di=1:N;     % index vector of full-dictionary
Q=D;
%delta=1/nbr_max_atoms;  %uncomment in the case of analytical norm
delta=1;   %uncomment in the case of discrete norm
if nargin<5 ind=[];end
if nargin<4 No=N;end
if nargin<3 tol=0.01;end;


if size(f,2)>1
    % code several vectors
    s = size(f,2);
    x = zeros(size(D,2),s);
    for i=1:s
        if s>20
            progressbar(i,s);
        end
        x(:,i) = perform_omp_other(f(:,i),D,tol,No,ind);
    end
    D = x;
    return;
end

% the code use transposed data
f = f';

numind=length(ind);

%atoms having smaller norm than tol1 are supposed be zero ones
tol1=1e-7; %1e-5
%threshold for coefficients
tol2=1e-10;   %0.0001  %1e-5

if verb
tic;
fprintf('\n%s is running for tol=%g, tol1=%g and tol2=%g.\n',name,tol,tol1,tol2);
fprintf('tol  -> required precision, distance ||f-f_approx||\n')
fprintf('tol1 -> atoms having smaller norm are supposed to be zero ones.\n');
fprintf('tol2 -> algorithm stops if max(|<f,q>|/||q||)<= tol2.\n');
end

%===============================
% Main algorithm: at kth iteration
%===============================
H=min(No,N); %maximal number of function in sub-dictionary
for k=1:H
    %finding of maximal coefficient
    c=zeros(1,N);cc=zeros(1,N);
    if k<=numind
        [test,q]=ismember(ind(k),Di);
        if test~=1 error('Demanded index %d is out of dictionary',ind(k));end
        c(k)=norm(Q(:,q));
        if c(k)<tol1
            error('Demanded atom with index %d is dependent on the others.',ind(k));
        end;
    else
        for p=k:N, c(p)=norm(Q(:,p));
            if c(p)>tol1 cc(p)=abs(f*Q(:,p))/c(p);end;
        end
        [max_c,q]=max(cc);
        %stopping criterion (coefficient)
        if max_c<tol2
            k=k-1;
            if verb
            fprintf('%s stopped, max(|<f,q>|/||q||)<= tol2=%g.\n',name,tol2);
            end
            break;
        end
    end
    if q~=k
        Q(:,[k q])=Q(:,[q k]); % swap k-th and q-th columns
        D(:,[k q])=D(:,[q k]);
        Di([k q])=Di([q k]);
    end
    %re-orthogonalization of Q(:,k)  w.r.t Q(:,1),..., Q(:,k-1)
    if k>1
        for p=1:k-1
            Q(:,k)=Q(:,k)-(Q(:,p)'*Q(:,k))*Q(:,p);
        end
    end
    nork=norm(Q(:,k));
    Q(:,k)=Q(:,k)/nork; %normalization
    % compute biorthogonal functions beta from 1 to k-1
    if k>1
        beta=beta-Q(:,k)*(D(:,k)'*beta)/nork;
    end
    beta(:,k)=Q(:,k)/nork;          % kth biorthogonal function

    %orthogonalization of Q(:,k+1:n) wrt Q(:,k)
    if k==N, break; end
    for p=k+1:N
        % orthogonalization of Q(:,p) wrt Q(:,k)
        Q(:,p)=Q(:,p)-(Q(:,k)'*Q(:,p))*Q(:,k);
    end
    %stopping criterion (distance)
    if (tol~= 0) & (norm(f'-D(:,1:k)*(f*beta)')*sqrt(delta) < tol) break;end;
end

c=f*beta;
normf=norm(f'-D(:,1:k)*c')*sqrt(delta);
if verb
fprintf('From %d atoms in this dictionary %s has chosen %d atoms.\n',N,name,k);
fprintf('\n%s lasted %g seconds\n',name,toc);
fprintf('\nNorm ||f-f_approx|| is %g. \n',normf);
end
% error testing
% orthogonality
errortest(Q(:,1:k));
% biorthogonality
errortest(D(:,1:k),beta);

if nargout==1
    x = zeros( size(D,2),1 );
    Di = Di(1:length(c));
    x(Di(:)) = c(:);
    D = x;
end


%Copyright (C) 2006 Miroslav ANDRLE and Laura REBOLLO-NEIRA
%
%This program is free software; you can redistribute it and/or modify it under the terms
%of the GNU General Public License as published by the Free Software Foundation; either
%version 2 of the License, or (at your option) any later version.
%
%This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
%without even the implied warranty of MERCHANTABILITY or FITNESS FOR X PARTICULAR PURPOSE.
%See the GNU General Public License for more details.
%
%You should have received a copy of the GNU General Public License along with this program;
%if not, write to the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
%Boston, MA  02110-1301, USA.


function f=errortest(D,beta);
% ERRORTEST tests orthogonality of a sequence or biorthogonality of two sequences.
%
% Usage: errortest(D,beta);
%        errortest(Q);
%
% Inputs:
%   D    sequence of vectors
%   beta (optional) biorthogonal sequence
%
% Output:
%   f    orthogonality of D or biorthogonality D w.r.t beta

% See http://www.ncrg.aston.ac.uk/Projects/BiOrthog/ for more details

name='ERRORTEST';

verb = 0;

if nargin==1
    % orthogonality
    PP=D'*D;f=norm(eye(size(PP))-PP);
    if verb
        fprintf('Orthogonality:   %g \n',f);
    end
elseif nargin==2
    % biorthogonality
    PP=beta'*D;f=norm(eye(size(PP))-PP);
    if verb
        fprintf('Biorthogonality: %g \n',f);
    end
else
    if verb
        fprintf('%s: Wrong number of arguments',name);
    end
end


%Copyright (C) 2006 Miroslav ANDRLE and Laura REBOLLO-NEIRA
%
%This program is free software; you can redistribute it and/or modify it under the terms
%of the GNU General Public License as published by the Free Software Foundation; either
%version 2 of the License, or (at your option) any later version.
%
%This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
%without even the implied warranty of MERCHANTABILITY or FITNESS FOR X PARTICULAR PURPOSE.
%See the GNU General Public License for more details.
%
%You should have received a copy of the GNU General Public License along with this program;
%if not, write to the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
%Boston, MA  02110-1301, USA.