function [U mu] = compute_subspace_basis(X,cn,th)
% This function computes an orthonormal basis for the affine subspace spanned by the columns of
% the data matrix X based on the chosen norm.
% Inputs
%   X   -  dxn data matrix
%   cn  -  selected norm 'l2' or 'l1'
%   th  -  threshold between 0 and 1 for determining the subspace size
% Ouput
%   U   - computed dxl orthormal basis matrix
%   mu  - mean of the samples in the projected subspace
%
% Written by Hakan Cevikalp, 11/18/2008, Eskisehir Osmangazi University

tol=10^-2; % stopping citeria (used in the iterative case involving 'l1' norm)
if nargin<2,
    cn='l2';
    th=1;
elseif nargin<3,
    th=1;
end
%keyboard
[d n]=size(X);
if d>n, % in case of high dimensinality we need some tricks
    mu=mean(X')';
    XX=X-repmat(mu, [1 n]);
    clear X
    LL=XX'*XX;
    [U,D,V]=svd(LL);
    [val,ind]=sort(diag(D),'descend');
    U=U(:,ind);
    te=sum(val); % total energy
    kk=cumsum(val)/te;
    ll=find(kk>=th); % finding the index where the ratio is larger than or equal to the selected threshold
    Uo=XX*U(:,1:ll(1)); % Basis for the subspace  
    l=size(Uo,2);
    % Obtaining an orthonormal basis set
    for i=1:l,
        Uo(:,i)=Uo(:,i)/norm(Uo(:,i));
    end
    if strcmp('l2',cn),
        % Basis based on L2 norm
        U=Uo;
        mu_projected=U*(U'*mu); % projection of the mean onto the subspace
    else
        % Basis based on L1 norm - R1-PCA  Algorithm
        % Intialization
        ss=diag(XX'*XX)-diag((Uo'*XX)'*(Uo'*XX));
        cc=sqrt(median(ss))+10^-5;
        % Computing weights
        Ut=Uo;
        loss_vec=sqrt(diag(XX'*XX)-diag((Ut'*XX)'*(Ut'*XX)));
        val_obj(1)=sum(loss_vec);
        for it=1:20;
            for i=1:n,
                dd=norm(XX(:,i)-Ut*(Ut'*XX(:,i)));
                if dd<=cc,
                    w(i)=1;
                else
                    w(i)=cc/dd;
                end
            end
            W=sqrt(diag(w));
            Unew=XX*W*((W*XX')*Ut);
            %sc=sqrt(diag(Unew'*Unew));
            Ut=orth(Unew);
            loss_vec=sqrt(diag(XX'*XX)-diag((Ut'*XX)'*(Ut'*XX)));
            val_obj(it+1)=sum(loss_vec);
            if abs(val_obj(it+1)-val_obj(it))<=tol;
                break;
            end
            %[it val_obj(it+1) val_obj(it)]
        end
        U=Ut;
        mu_projected=U*(U'*mu);
        Lambda=U'*XX*W*((W*XX')*Ut);
    end
else % if n>=d
    mu=mean(X')';
    XX=X-repmat(mu, [1 n]);
    S=cov(X');
    [U D V]=svd(S);
    [val,ind]=sort(diag(D),'descend');
    U=U(:,ind);
    te=sum(val); % total energy
    kk=cumsum(val)/te;
    ll=find(kk>=th); % finding where the ratio is larger than or equal to the selected threshold
    Uo=U(:,1:ll(1)); % Basis for the subspace  
    l=size(Uo,2);
    if strcmp('l2',cn),
        % Basis based on L2 norm
        U=Uo;
        mu_projected=U*(U'*mu); % projection of the mean onto the subspace
    else
        % Basis based on L1 norm - R1-PCA  Algorithm
        % Intialization
        ss=diag(XX'*XX)-diag((Uo'*XX)'*(Uo'*XX));
        cc=sqrt(median(ss))+eps;
        % Computing weights
        Ut=Uo;
        val_obj(1)=norm(Ut'*XX);
        for it=1:50;
            for i=1:n,
                dd=norm(XX(:,i)-(Ut*Ut')*XX(:,i));
                if dd<=cc,
                    w(i)=1;
                else
                    w(i)=cc/dd;
                end
            end
            W=sqrt(diag(w));
            Unew=(XX*W*W*XX')*Ut;
            %sc=sqrt(diag(Unew'*Unew));
            Ut=orth(Unew);
            val_obj(it+1)=norm(Ut'*XX);
            if abs(val_obj(it+1)-val_obj(it))<=tol;
                break;
            end
            [it abs(val_obj(it+1)-val_obj(it))]
        end
        U=Ut;
        mu_projected=(U*U')*mu;
        Lambda=U'*XX*W*((W*XX')*Ut);
    end
end







