% sftMx2LyrPt.m
%             23 Mar 2017  <= 13 Jan 2017
% data inspired by   http://cs231n.github.io/neural-networks-case-study/
% Andrew P Paplinski, Monash University
% softmax 2-Layer nnet classifier 
% Pattern learning

clear
disp('Pattern Learning with a Softmax 2-Layer Nnet Classifier')
tic
load Sprl ;   % X, Y  noisy spiral data
[D, N] = size(X) ;     % dimensionality of X and the number of points
M = length(unique(Y)); % number of classes, dimensionality of the output Z

Xh = [X ; ones(1,N)] ; % adding bias input
% converting the expected class labels Y into outputs YY
YY = zeros(M,N) ;
for n = 1:N
  YY(Y(n), n) = 1;  
end

% NNet specification:
% Dh hidden neurons with ReLU, M output neurons with Softmax
Dh = 100 ; % number of hidden neurons
% Weights initialization
Wh = 0.01*(rand(Dh,D+1)-0.5) ; % 1st weight matrix is Dh by D+1
Wz  = 0.01*(rand(M,Dh+1)-0.5) ; % 2nd weight matrix is M by Dh+1

% learning parameters
etah = 5/N ;    % learning rate (hidden)
etaz = 1/N ;    % learning rate (output)
reg = 1e-2/N ;  % lambda in reqularization
nepchs = 3000 ; % number of learning epochs  

L = zeros(nepchs,1) ;     % total loss
LData = zeros(nepchs,1) ; % data loss
logpk = zeros (N,1) ;     % pattern log(pk_n)

% learning
for epch=1:nepchs   % the epoch loop 
  for n=1:N % pattern loop, for each pattern  
    % evaluate all signals in the 2-layer Neural Network
    % first layer
    uh = Wh*Xh(:,n) ;   % postsynaptic activities in the hidden layer
    zh = [max(0,uh); 1];% ReLU nonlinearity and bias
    % second layer
    % Calculating the postsynaptic activities (scores) M by 1
    uz = Wz*zh ; % each column (stimulus) gives the class scores  
                 % corresponding to the M classes 
    % the Softmax classifier -- the cross-entropy loss function
    ewx = exp(uz) ;     %  is M by 1
    s_ewx = sum(ewx) ;  %  is 1 by 1
    Pr = ewx./s_ewx ;   % probabilities, M by 1 = outputs z_n
           % selecting the correct probability for each class
    pk = sum(Pr.*YY(:,n)) ;  % is 1 by N , class probabilities
    logpk(n) = log(pk) ;
    
    % gradient of LData for the n-th pattern
    gL = Pr - YY(:,n) ;
    
    % The weight matrix update (including the regularization)
    dWz = gL*zh' + reg*Wz ; 
     Wz = Wz - etaz*dWz ;

    % backpropagation from 2nd layer to the 1st layer
    delh = Wz(:,1:end-1)'*gL ;  % backprop error
    delh(uh<=0) = 0 ;           % backprop through ReLU
    % The weight matrix update (including the regularization)
    dWh = delh*Xh(:,n)' + reg*Wh  ; 
    Wh = Wh - etah*dWh ;
  end   % pattern loop
  
  % data loss is an average cross entropy
  LData(epch) = -mean(logpk) ; % mean of -log(probs(pk))
  % regularization loss
  Lreg = 0.5*reg*(Wh(:)'*Wh(:) + Wz(:)'*Wz(:)) ;
  L(epch) = LData(epch) + Lreg  ;  % total loss
    
end     %  the epoch loop  

% Validation = training set accuracy
Uh = Wh*Xh ; % postsynaptic activities in the hidden layer
Zh = [max(0,Uh); ones(1,N)];   % ReLU nonlinearity and bias
Uz = Wz*Zh ; 
% evaluate learning set accuracy
[~, Umx] = max(Uz) ; % predicted classes
tsacc = sum(Umx == Y)/N ;

figure(1), clf
plot(1:nepchs, [L LData]), grid on
legend('L', 'Ldata')
text(20,0.9, num2str(tsacc, 'Learning accuracy: %1.3f'))
title('Loss function')

klr = hsv(M) ;  % colours for clusters
% calculating the boundries
% ax = [0 1.1 0 1.1] ;  % axis of the plot
ax = [-1 1 -1 1] ;  % axis of the plot
k = 50; % K = (2*k+1)^2 ;
kk = (ax(2)-ax(1))*((0:k)/k)+ax(1) ;
[x1, x2] = meshgrid(kk) ; nx = numel(x1) ;
XT = [x1(:) x2(:) ones(nx,1)]' ;

Uz = Wz*[max(0, Wh*XT); ones(1, nx)] ;
[~, kl] = max(Uz) ;

% plotting the data points and the class boundries
figure(2), clf 
axis(ax), grid on, hold on
for m = 1:M
  plot(Xh(1,Y==m),  Xh(2,Y==m),  'o', 'Color', klr(m,:))
  plot(XT(1,kl==m), XT(2,kl==m), '.', 'Color', klr(m,:))
end
xlabel('x_2'), ylabel('x_1')
hold off
title('Data points and calculated class boundries')
toc

% disp('Linear Classifier fails to learn the spiral dataset')
% print('-f1', '-dpdf', 'sftMxLLss')
% print('-f2', '-dpdf', 'sftMxLsep')
% print('-f1', '-dpdf', 'sftMxLspr1')
% print('-f2', '-dpdf', 'sftMxLspr2')
% print('-f1', '-dpdf', 'sftMx2Lyr1')
% print('-f2', '-dpdf', 'sftMx2Lyr2')
% print('-f1', '-dpdf', 'sftMx2LrP1')
% print('-f2', '-dpdf', 'sftMx2LrP2')


