function [varargout]=matchingpursuit(data,dictionary,varargin)
% A matching pursuit algorithm.  It fits data from a set of dictionary
% elements by orthogonally projecting the data onto elements of the
% dictionary of vectors such that the energy in each projection is
% maximized. In constrast to G-S, the residual is orthogonalized, not
% the basis vectors.
%
%INPUTS
%data: A column vector containing a data vector (Hilbert-space vector)
%dictionary: A matrix whose columns contain dictionary elements.
%            Dictionary elements are Hilbert-space elements that are
%            being used to reconstruct the data.
%Optional Inputs:
%plotopt: If opt is NOT input, the string 'plot' will produce a plot.
%opt: A structure which contains fields that are options for plotting,
%     pre-aligning the data, pre-weighting the data, and weighting the
%     inner product. The fields are described below:
%opt fields:
%plot: string 'plot' produces a plot
%align: Two string options:
%       'max' aligns all dictionary elements in the center of the domain
%       by the maximum value so reconstruction begins element-by-element
%       in the middle of the data; 'corr' aligns all the data in the center
%       of the domain, by the maximum cross correlation value.  This option
%       is only recommended if the support of the data is well-known.
%window: A weight applied to the data.  Reconstruction is then applied to
%        product of the data and the window. window must be a vector the same
%        size as the data. Good for filtering.
%weight: A weight applied to each dictionary element. weight can be a
%        vector the same length as data or a matrix with the same number
%        of columns as the dictionary. If some dictionary elements are
%        down-weighted, then 'less of them' will be used in the data
%        reconstruction. Good for filtering.
%index: Three string options: 'pdf', 'width', 'erg'. This option orders the
%       dictionary elements by (respectively) noise content, signal width,
%       or energy localization (the max amplitude divided by the signal
%       width). If a plot is made, the projection coefficients are plotted
%       against the index option.
%
%!!
%  NOTE: There are some very simple functions that this function uses:
%  (1) sortwvfrms (2) maxalign (3) alignXmatrix (4) plotXmatrix.
%  All can be downloaded from the file exchange.
%!!
%--------------------------------------------------------------------%
%EXAMPLE: A signal with noise, w/out an opt structure:
% t=linspace(0,1,200);
% x=zeros(size(t));
% x(50:149)=triang(100)+(0.3)*(1-rand(100,1))-1/2;
% x=x';
% Q=zeros(200,100);
% for(k=1:200)
%  Q(:,k)=sin(k*pi*t);
% end;
% [coefficients,basis]=matchingpursuit(x,Q,'plot');
%
%EXAMPLE: Use the signal above, with opt that weights data to filter
% out corner effects:
% opt.plot='plot';
% opt.weight=triang(length(x));
% [coefficients,basis]=matchingpursuit(x,Q,opt);
%
%EXAMPLE: Use the signal above, with opt that weights the dictionary to
%reconstruct with/out corner effects
% opt.plot='plot';
% opt.weight=0.001+tukeywin(length(x));
% [coefficients,basis]=matchingpursuit(x,Q,opt);
%
%EXAMPLE: Use the signal above, with opt that windows the data to
%reconstruct just a portion of it:
% opt=rmfield(opt,'weight');
% opt.plot='plot';
% w=zeros(200,1); w(75:125)=1;
% opt.window=w;
% [coefficients,basis]=matchingpursuit(x,Q,opt);
%
%EXAMPLE: Include the weighting:
% opt.weight=w;
% [coefficients,basis]=matchingpursuit(x,Q,opt);
% Note that the data amplitude is larger: the transform is energy
% preserving.
%
%EXAMPLE: sort waveforms according to signal width:
% opt=rmfield(opt,'weight');
% opt=rmfield(opt,'window');
% opt.index='width';
% [coefficients,basis,Iout]=matchingpursuit(x,Q,opt);
%--------------------------------------------------------------------%

%Initializing Variables.
[M,N]=size(dictionary);
[m,n]=size(data);
colors={'k','r','b'};
indexopt=struct('');
proj=zeros(N,1);
phi=zeros(M,1);
coeffs=zeros(N,1);
%resid=zeros(N,1);
D=dictionary;
Idict=1:N; %dictionary index initialization
x=data;

if(nargin>0)
    plotopt='none';
end;
if(nargout==0)
    plotopt='plot';
end;

if(m~=M)
    error('The dictionary and data must have the same number of rows')
end;

if(n~=1)
    error('Data input must be column vector');
end;

if(nargout==3&&(~isstruct(varargin{1})))
    error('3 outputs require a structure input with the field index');
end;

if(nargout==3&&isstruct(varargin{1}))
    if(~isfield(varargin{1},'index'))
        error('3 outputs require a structure input with the field index');
    end;
end;
%% Check Input for 'Opt' Structure
if(nargin==3)
    if(isstruct(varargin{1}))
        opt=varargin{1};
        if(isfield(opt,'plot'))
            plotopt=opt.plot;
        end;
        if(isfield(opt,'align'))
            alignopt=opt.align;
            if(strncmpi(alignopt,'corr',4))
                D=alignXmatrix([x,D],'align');
                D(:,1)=[];
            end;
            if(strncmpi(alignopt,'max',3))
                D=maxalign(D);
                x=maxalign(x);
                data=maxalign(data);
            end;
        end;
        if(isfield(opt,'window'))
            wind=opt.window;
            x=wind.*x;
        end;
        if(isfield(opt,'weight'))
            wt=opt.weight;
            if(min(size(wt))==1)
                D=diag(wt)*D;
            end;
            if(~(min(size(wt))==1))
                D=wt*D;
            end;
        end;
    end;
    if(strncmpi(varargin{1},'plot',4))
        plotopt=varargin{1};
    end;
end;
% if(nargout==2)
%     error('Need opt structure to produce 3 outputs');
% end;
%% Forming the Orthogonal Projections
Iproj=zeros(N,1); %Indices of max. projection
sumind=1:N; %Indices to sum over
%Finds the maximum projection of data onto elements of the dictionary.
%Calls the maximum index "I".
for j=1:N;
    proj(1:N)=zeros(N,1); %reinitialize proj.
    for k=sumind
        if(j==1)
            D(:,k)=(D(:,k))/norm(D(:,k));
        end;
        proj(k)=x'*D(:,k);
    end;
    %Extract the dictionary element with the maximum projection
    [projmax,I]=max(abs(proj));
    phi(:,j)=D(:,I);
    coeffs(j)=proj(I);
    x=x-projmax*D(:,I);
    Iproj(j)=I;
    sumind=setdiff(sumind,I);
end;
%% Using Sortopt to Determine if to use some Indexing on Dictionary
if(nargin==3)
    if(isstruct(varargin{1})&&isstruct(opt))
        if(isfield(opt,'index'))
            indexopt=opt.index;
            if(ischar(indexopt))
                ID=sortwvfrms(phi,char(indexopt));
                Isort=ID(:,1);
                Imeas=ID(:,2);
                Imeas=Imeas(Isort);
                phi=phi(:,Isort);
                coeffs=coeffs(Isort);
                Idict=Imeas;
            end;
        end;
        if(~isfield(opt,'index')&&(nargout==3))
            error('Need index field to produce 3 outputs');
        end;
    end;
end;
%%
if(strcmp(plotopt,'plot'))
    Ndisp=20;
    figure;
    subplot(2,2,[1 3]);
    for k=1:Ndisp
        s=plot(Idict(k),abs(coeffs(k))/norm(abs(coeffs)),'o',...
            'MarkerFaceColor',...
            char(colors(mod(k-1,3)+1)), ...
            'MarkerEdgeColor',...
            char(colors(mod(k-1,3)+1)));
        set(s,'markersize',5);
        hold on;
        if(k==Ndisp)
            legend(sprintf('1st %i Projections',Ndisp));
        end;
    end;
    s=plot(Idict,abs(coeffs)/norm(abs(coeffs)),'b');
    set(s,'linewidth',2); axis tight;
    title('All Basis Coefficients','fontsize',14);
    set(gca,'Color', [0.85,0.9,0.95]);
    subplot(2,2,2);
    s=plot(data,'k');
    set(s,'linewidth',1.8); hold on;
    dataout=phi*coeffs;
    s=plot(dataout,'r');
    set(s,'linewidth',1.8);
    legend('Data','Reconstruction');
    title(sprintf('Reconstruction Using the Largest Projections'),'fontsize',14);
    subplot(2,2,4);
    plotXmatrix(phi(:,1:Ndisp));
    title(sprintf('First %i Dictionary Vectors',Ndisp),'fontsize',14);
    set(gcf,'Color', [1,1,1]);
end;

%% Assigning Outputs to Function
if(nargout==1)
    varargout{1}=coeffs;
end;
if(nargout==2)
    varargout{1}=coeffs;
    varargout{2}=phi;
end;
if(nargout==3&&nargin==3)
    if(isstruct(varargin{1}))
        varargout{1}=coeffs;
        varargout{2}=phi;
        varargout{3}=Imeas;
    end;
end;