% fap2D.m    
%         1 March 2017 <= 13 March 2005 <= 25 March 1999
% Example of approximation of two 2-dimensional functions
% using backpropagation in a batch mode
%
% Andrew P. Paplinski, Monash University, AUSTRALIA 

clear

D  = 2 ; % Dimensionality of input signals
Dh = 24; % Dimensionality of hidden signals 
M = 2 ;  % Dimensionality of outputs

na = 16 ;  N = na^2; nn = 0:na-1;  % Number of training cases

% Generation of the training cases as coordinates of points
% from two 2-D surfaces
% Specification of the sampling grid
X1 = nn*4/na - 2; % na  points from  -2 step (4/na)=.25
                  %            to (2 - 4/na)=1.75
[X1, X2] = meshgrid(X1);  % coordinates of the grid vertices
                          % X1 and X2  are  na by na

R  = (X1.^2 + X2.^2 +1e-5);  % R is a  na by na  matrix of
   % squares of distances of the grid vertices from the origin

Y1 = X1 .* exp(-R);  Y = (Y1(:))';
                  % Y1 is na by na, Y  is   1 by N

Y2 = 0.25*sin(2*R)./R ;  Y = [Y ; (Y2(:))'];
           % Y2 is na by na
           % Y  is  2 by N  matrix of target 2D output vectors

Z = zeros(size(Y)) ;  % memory allocation for the output matrix
Xh = [ X1(:)'; X2(:)'; ones(1,N) ];  % Appending the bias inputs
                                    % X  is  D+1 by N
figure(1), clf reset, hold off
myfigA(1,11,10)
surfc([X1-2 X1+2], [X2 X2], [Y1 Y2]), % Two 2-D target functions
                                      % are ploted side by side
grid on ,  title('Two 2-D target functions') 
text(-1,2,0.25, 'y_1')
text( 4,0,0.25, 'y_2')
print -dpdf -f1 dFap2Dy

% Initialization of the weight matrices
% Hidden layer weight matrix
 Wh = randn(Dh, D+1)/(D+1);     % Wh  is Dh by D+1
% Output layer weight matrix
 Wz = randn(M, Dh+1)/(Dh+1) ;   % Wz  is  M by Dh+1

nepchs = 2000;           % maximum number of training epochs
Er = zeros(1, nepchs);   % memory allocation for the error function
gamma = [0.0010 0.035];  % Training gains

% **********************  main loop *********************
figure(2), clf reset, hold off
tic
for ep = 1:nepchs  % the epoch loop

 % The forward pass
 Uh = Wh*Xh;  % hidden post synaptic
 % Hidden signals  (Dh by N)  
 Zh = ones(Dh, N)./(1+exp(-Uh));	
 dZh = Zh.*(1-Zh);       % Derivatives of hidden signals
 Zhh = [Zh ; ones(1,N)] ; % appending the bias input
 Z = tanh(Wz*Zhh)  ;      % Output signals        (M by N)
 dZ = 1 - Z.^2 ;	      % Derivatives of output signals
                                    
 % The backward pass
 Ey = Z - Y;		     % The output errors      (M by N)
 % the performance function: total mse
 Er(ep) = sum(sum(Ey.^2)/2)/N ;  % mse after each epoch
 
 delY = Ey.*dZ;		% Output delta signal    (M by N)
 dWz = -delY*Zhh';	% Update of the output matrix M by Dh+1                       
 Eh = Wz(:,1:end-1)'*delY; % The backpropagated hidden error Dh by N                      
 delH = Eh.*dZh ;   % Hidden delta signals  (Dh by N)
 dWh = -delH*Xh';	% Update of the hidden matrix                                              
 
 % The batch update of the weights:
 Wz = Wz + gamma(1)*dWz ;  Wh = Wh + gamma(2)*dWh ;
 Y1(:)=Z(1,:)'; Y2(:)=Z(2,:)';

 % Two 2-D aproximated functions are ploted side by side
 surfc([X1-2 X1+2], [X2 X2], [Y1 Y2]), grid on 
 title(['epoch: ', num2str(ep), ', error: ', num2str(Er(ep),'%0.4f  '), ...
       ', gamma: ', num2str(gamma,'%0.4f  ')]), drawnow
end  % of the epoch loop
toc

figure(3), clf reset
plot(Er'), grid on
title(['The approximation error. Final: ', num2str(Er(ep),'%0.4f  ')] )
xlabel('number of training epochs')

print -dpdf -f2 dFap2Dz
print -dpdf -f3 dFap2De

