NSC(Train, Test, cs, repeat, type)

function [result2, relative_error2] = NSC(Train, Test, cs, repeat, type)
% 
% Nearest Subspace Classifier
% [1] A. Majumdar and R. K. Ward, Nearest Subspace Classifier submitted
% to International Conference on Image Processing (ICIP09).
% [2] A. Majumdar and R. K. Ward, Compressive Classification 
% in preparation.
% Angshul Majumdar (c) 2009

if nargin < 3 cs = 1; end
if nargin < 4 repeat = 1; end
if nargin < 5 type = 'general'; end

% --- INPUTS
% Train.X - training samples
% Train.y - training labels
% Test.X - testing samples
% Test.y - testing labels
% cs - number of compressive samples to be takes (<1)
% repeat - number of times to be repeated
% type - 'general' Random Undersampling
%        'gaussian' Gaussian Undersampling
% --- OUTPUTS
% result - recognition result for entire test data
% relative_error - relative recognition error

% y - normalised test sample
% A - normalised training samples
% x - sparse classification vector

d = size(Train.X,1);    % dimension of the feature vectors

for r = 1:repeat
        % finding restricted isometry
        j=1;
        R(j) = ceil(round(1/cs)*j*rand);
        while j ~= round(d*cs)
            temp = ceil(round(1/cs)*(j+1)*rand);
            if temp > R(j)
                j=j+1;
                R(j)=temp;
            end
        end
        I = eye(d);
        for j = 1:round(d*cs)
            IR(j,:) = I(R(j),:);    % IR restricted isometry matrix
        end
%     end

    G = randn(d,d);
    GR = IR*G;      % GR restricted gaussian matrix

    switch lower(type)
        case {'general'}
            trn.X = IR*Train.X;
            trn.y = Train.y;
            tst.X = IR*Test.X;
            tst.y = Test.y;
        case {'gaussian'}
            trn.X = GR*Train.X;
            trn.y = Train.y;
            tst.X = GR*Test.X;
            tst.y = Test.y;
    end

    % Normalising the training set
    for i=1:size(trn.X,2)
        Atilde(:,i)=trn.X(:,i)/norm(trn.X(:,i),2);
    end

    n = size(trn.X,2); % total number of samples in the training database
    k = max(trn.y); % number of classes
    onesvecs = full(ind2vec(trn.y)); % preparing classification target
    n_test = size(tst.X,2); % number of test data to be classified

    % start classifcation
    for iter = 1:n_test % looping through all the test samples
        ytilde = tst.X(:,iter);
        ytilde = ytilde/norm(ytilde);
        % find the error for each class
        for i = 1:k % number of classes
            ind = find(trn.y == i);
            V = Atilde(:,ind);
            alpha = V\ytilde;
%             alpha = inv(V'*V)*V'*ytilde;
            xp = zeros(n,1);    xp(ind) = alpha;
            deltavec(:,i) = onesvecs(i,:)'.* xp;
            residual2(r,iter,i) = norm(ytilde-Atilde*deltavec(:,i));
        end
    end
end

res2 = sum(residual2,1);

for iter = 1:n_test
    [B, Ind] = sort(res2(1,iter,:));
    result2(iter) = Ind(1);
%     result2(iter) = find(res2(1,iter,:) == min(res2(1,iter,:))); % storing the classification results
end

% finding the relative classification error
relative_error2 = size(find(abs(Test.y-result2)),2)/size(Test.y,2);