% Perceptron training function.
%
% X is a [n,f] matrix where n is the number of samples and f is a number of features.
% alpha is the learning rate
% D is a desired output [n,dim]; dim is a dimensionality of the output
% epochs limit is limit on the max number of epochs when the weights do not change
%
function W = perceptron(X, D, alpha, epochs_limit)
% get values for n, k, m
% Add bias x0 = 1 to input
[n,f] = size(X);
X = [ones(n,1) X];
[n,f] = size(X);
[N,dim] = size(D);
if (n ~=N)
    disp('Error: dimensionalities of X and D do not match.');
    return;
end

t = 0;% Initialize iteration number: t=0
W = rand(dim,f); % set weights to arbitrary values, f x dim dimensional vector
% Show some information about each iteration
disp(sprintf('Epoch#\tM.S.E.\t\tElapsed Time (s)\tIterations w/o change'));
% not_done = epochs_limit since we need certain number of iterations where
% weights don't change
not_done = epochs_limit; 
time = cputime; % get the current CPU time
totalmse = 0; % total the MSE.
maxrun = 0; % longest run without updating W

while (not_done > 0)
    i = mod(t,epochs_limit) + 1; % current iteration
    % Compute W(t)=ae(k)x(k), i.e. wi(t)=ae(k)x(ik) for i=0,...,n, where
    % e(k)=D(k)-y(k) and y(k)=f(i=0...nwi(t)x(ik));
   
    y = ((W * X(i,:)') >= 0);
    e = D(i,:)' - y;
    dW = alpha * e * X(i,:);
   
    % Increment iteration number t=t+1;
    t = t + 1;
    % Update W(t)=W(t-1)+~W(t-1), i.e. wi(t)=wi(t-1)+~wi(t-1) for i=0,...,n;
    W = W + dW;
    
    % If no weight changes occurred during last epoch or other stopping 
    % condition is true, then stop; otherwise, continue
    if (nnz(dW) == 0)% check the number of non zero elements 
        not_done = not_done - 1; % the weights did not change
    else % some weights changed, start counting epochs to the stoping condition
        not_done = epochs_limit;
    end
    
    % in case we were looking for solution for too long
    if (cputime - time) > 3600
        disp(sprintf('Giving up, over one hour of CPU time used.\n'));
        not_done = 0;
    end
    
    % Display data about the iteration
    thismse = mse(dW);
    totalmse = totalmse + thismse;
    
    if (maxrun < (epochs_limit - not_done))
        maxrun = (epochs_limit - not_done);
    end
    
    disp(sprintf('%d\t\t%f\t\t%f\t\t%d', t, thismse, (cputime - time), (epochs_limit - not_done)));
end

% Display average MSE.
disp(sprintf('Average M.S.E. is %f for %d epochs', (totalmse / t), t));
% Display time
disp(sprintf('Total CPU Time: %f (seconds).', (cputime - time)));