[varargout]=projOnClusters(x,W,varargin)

function [varargout]=projOnClusters(x,W,varargin)
% projOnClusters: A projection pursuit type algorithm that decomposes a
% signal by projecting the residual at each level onto a convex cone and
% the respective polar cone.
%
% The set giving cone k is described by a vector w = W(:,k) and scalar rho:
%
% C_k = { x :  x'w  >= rho*|| x || };
%
% The polar cone is:
%
% C*_k = { y: y'x <= 0, x element of C_k };
%
% where w is the reference vector for the set C. The projector for this set
% is described in:
% 
% Henry Stark and Youngyi Yang, "Vector Space Projections: A Numerical 
% Approach to Signal and Image Processing, Neural Nets, and Optics", Wiley
% Interscienc Publications, 1998.
%
% On pages 109-113, Section 3.5, in the discusssion on similarity
% constraints.
%
%--------------------------------------------------------------------------
% USAGE
%--------------------------------------------------------------------------
%
% [u]=projOnClusters(x,W);
%
%       ...
%
% [u,e,b,ind]=projOnClusters(x,W,rho,'plot');
%
% INPUT
% x: the signal to be analyzed. Must be a column vector or a matrix.
% W: A matrix whose columns are the cone representers (see [].pdf for a
%    description of these). W must have the same number of rows as x does.
% rho: A row vector of correlation coefficients corresponding to the columns
%      of W.
%
% OUTPUT
% u:    The projection onto the convex cones defined by W and rho.
% e:    The projections onto the polar cones correpsonding to the convex
%       cones
% b:    The norm of each projection onto a convex cone: b(k) = ||u(:,k)||
% ind:  The index of cluster that corresponds to u,e, and b.
%
% If no output is requested, a plot is made.  Alternatively, a plot is
% produced if an input argument 'plot' is given.
%
%% Examples:
%%Example 1: Using 5 signals to obtain a rough estimate of a step function
% t   = [linspace(-2*pi,2*pi,256)]';
% x   = zeros(size(t)); x(128-20:128+20)=4;
% y   = zeros(size(t)); y(128-20:128+20)=tukeywin(41,0.4);
% W   = [sin(2*t), y, sinc(2*t), t.*(t-1).*exp(-abs(t)/2),8*t.*[tukeywin(length(t))]];
% rho = 0.95*ones(1,size(W,2)); rho(2:3) = 0.94; %describes cone geometry
% 
% clf;
% [u,e,b,ind]=projOnClusters(x,W,rho,'plot');
% %Check energy preservation
% disp('Check Energy Preservation....');
% disp('Energy in reconstruction');
% recon_erg = sum(b.^2);
% disp('Energy in signal');
% orig_erg = norm(x).^2
% disp('Error : ')
% recon_erg - orig_erg
%--------------------------------------------------------------------------

%Get level of algorithm and length of vectors
J = size(W,2);
M = size(W,1);

%Initialize Outputs
b    = zeros(J+1,1);
u    = zeros(M,J+1);
e    = zeros(M,J+1);
ind  = zeros(J+1,1);
plotopt     = 'none';
rho         = 0.9*ones(1,J);

%% Default variable inputs
if(norm(W)<eps)
    error('W effectively zero');
end;

if( max(rho)>1||min(rho)<0 )
    error('rho must be between 0 and 1');
end;

if(nargin>2)
    for k=1:length(varargin)

        if(ischar(varargin{k}))
            plotopt = varargin{k};
        end;

        if(isnumeric(varargin{k}))
            rho     = varargin{k};
        end;
    end;
end;

%% Normalize Cluster Waveforms (without loops)
temp = W.^2;
temp = (sum(temp)).^(1/2);
s    = diag(temp.^(-1));
W    = W*s;

%make a copy of input W;
C   = W;

%Initialize the first error vector to be the signal
e(:,1)=x;

%% The Convex Set Projection Pursuit Algorithm
for k = 2:J+1

    %--------------------------------------------------------------------------
    %% Edit 21.Feb.2009: Optional to keep xcor-portion of code
%     temp    = [e(:,k-1),W];
%     [temp]  = xcorAlign(temp);
%     W       = temp(:,2:end);
    %--------------------------------------------------------------------------

    %Find the maximum projection of e_k-1 onto a cluster
    [m,i]   = max(rho.*(e(:,k-1)'*W));

    %Record the cluster index
    ind(k)  = i;

    %Extract the waveform for C_k, which is W(:,i)
    w       = W(:,i);

    %Find it's orthongonal complement
    wp      = e(:,k-1) - (e(:,k-1)'*w)*w;

    %Normalize it's orthongonal complement
    wp      = wp/norm(wp);

    %Test to see if e_k-1 is in the convex set defined by W(:,i)
    inSet   = (e(:,k-1)'*w >= rho(i)*norm(e(:,k-1)));

    %if e_k-1 is in C_k then u_k is e_k-1, and the algorithm stops.
    if(inSet)

        u(:,k)   = e(:,k-1);
        b(k)     = norm(u(:,k));
        continue;

    else
        %if e_k-1 is not in C_k then in may be in the C*_k, or have
        %a nonzero projection into both C_k and C*_k
        beta     = [e(:,k-1)'*[w,wp]]';
        r        = [rho(i);sqrt(1-rho(i)^2)];
        inPol    = (r'*beta <= 0);

        %if e_k-1 is in C*_k then it is equal to e_k
        if(inPol)
            e(:,k) = e(:,k-1);
        end;

        %if it is neither in C_k nor in C*_k then it has a projection onto
        %the boundary of C_k, and the boundary of C*_k
        if(~inPol)
            u(:,k) = [w,wp]*r*r'*beta;
            e(:,k) = [w,wp]*(eye(2,2)-r*r')*beta;
            b(k)   = norm(u(:,k));
        end;

    end;

    %Remove the i_th cluster so that it will not be tested again
    W(:,i)=NaN(M,1);

end;

u           = u(:,2:end);
e           = e(:,2:end);
ind         = ind(ind>0);
temp        = zeros(size(ind));
temp(ind)   = b(2:end);
b           = temp;

if(nargout>0)

    for k = 1:nargout

        if(nargout>=1), varargout{1}=u;     end;
        if(nargout>=2), varargout{2}=e;     end;
        if(nargout>=3), varargout{3}=b;     end;
        if(nargout>=4), varargout{4}=ind;   end;

    end;

end;


%% Plot Synthesis and Energy Distribution

if(nargout<1)
    plotopt='plot';
end;

ind = ind(ind>0);

if(strcmp(plotopt,'plot'))

    subplot(2,2,[1,2]);
    plot(x,'k','linewidth',2); hold on;
    %plot(sum(u,2)+e(:,end),'--r','linewidth',2);
    plot(sum(u,2),'--r','linewidth',2);
    title('Signal With Cluster Reconstruction','fontsize',16,'fontweight','bold');
    set(gca,'fontsize',16,'fontweight','bold');
    s=legend('Signal','Reconstruction');
    set(s,'fontsize',16);
    legend('boxoff');
    set(gca,'ytick',[]);
    set(gca,'Color', [0.85,0.9,1]);
    axis tight;

    for k=1:size(ind>0);

        if(b(k)>0.05*max(b))
            color='r';
        else
            color='k';
        end;

        subplot(2,2,3); hold on;

        plot(C(:,k)+k,char(color),'linewidth',1);

        subplot(2,2,4); hold on;

        bar(k,b(k)^2,char(color));

    end;

    subplot(2,2,3);

    set(gca,'xtick',[]);
    xlabel('Cluster Waveforms','fontsize',18,'fontweight','bold');
    set(gca,'fontsize',18,'fontweight','bold');
    axis tight;
    set(gca,'Color', [0.85,0.9,1]);

    subplot(2,2,4);

    %set(gca,'xtick',[]);
    xlabel('Energy in Cluster','fontsize',18,'fontweight','bold');
    set(gca,'fontsize',18,'fontweight','bold');
    axis tight;
    set(gca,'Color', [0.85,0.9,1]);

    set(gcf,'color','w');
end;