function perceptron_delta
%
% NSCI 850 - Computational Approaches to Neuroscience
%
% simple perceptron for face recognition
%
% G. Blohm (2010)
% adapted from S. Seung (MIT)
%-------------------------------------------------------------------------

clear

%% load training and test sets
load('ps9faces.mat'); % NOTE: the matlab and training set files must be in the current directory

%% define transfer function
g = inline('1./(1+exp(-x))');   % sigmoid transfer function of output node

%% initialize
[n,m] = size(train);    % number of pixels in each picture (n) and number of examples (m)
s = sqrt(n);            % pictures are square images
eps = 0.1;              % learning rate
etainit = 0.01;         % parameter controlling magnitude of initial conditions
w = etainit*randn(n,1); % initialization with Gaussian random values for weights
b = etainit*randn;      % initialization of bias input to network

tmax = 10000;           % max number of iterations to train the network
errsq = zeros(tmax,1);  % squared error on training set (initialization for memory opimization)

%% start training
for t = 1:tmax
    i = ceil(m*rand);               % pick a random image for online training
    x = double(train(:,i))/255;     % uint8 -> double, normalize max value to one
    y = double(trainlabels(i));     % this is the training goal associated with the random image
    actual = g(w'*x+b);             % current network performance
    error = y - actual;             % network error
    slope = actual*(1-actual);      % this simple expression results from the choice of the transfer function
    errsq(t) = error'*error;        % current squared error
    delta = slope*error;            % error signal scaled by derivative
    w = w + eps*delta*x;            % weight updates
    b = b + eps*delta;              % bias input updates
    if rem(t,100) == 0              % plot every 100 iterations
        subplot(2,1,1)
        imagesc(reshape(w,s,s))     % you can plot the actual images using the same function...
        title(sprintf('w, t=%d',t))
        axis image off
        colormap gray
        subplot(2,1,2)
        loglog(1:t,cumsum(errsq(1:t))./(1:t)') % cumulative mean squared error
        drawnow
    end
end

disp(sprintf('mean squared error: %.4f',mean(errsq)))
