% Multilayer feedforward network for function approximation

clear all
figure(1)
clf
rand('state',9)
% Create training data
% input data
x = -4:0.01:4;
% Output data
y = cos(2*pi*x).*exp(-4*x.^2);

%Display training data
plot(x,y)
xlabel('Input x')
ylabel('Output y')
axis tight
title('Unknown Function to be Approximated'); 
pause

% Building the Network 

% Range of input data
PR = [min(x) max(x)]; 

% No. of nodes in each layer
S1 = 8;  %layer 1
S2 = 1;  %layer 2

%Activation function of each layer
TF1 = 'logsig';
TF2 = 'purelin';

%Back propagation network learning function
BTF = 'trainlm';

%Back propagation weight/bias learning function
BLF = 'learngd';

%Cost function or performance function or error function
PF = 'mse';

rand('state',0);
% Create the network
net = newff(PR,[S1 S2],{TF1 TF2},BTF,BLF,PF);

% Before training
figure(2)
clf

w1 = net.IW{1}; % weights of w1 input to hidden layer
b  = net.b{1};  % bias b of hidden layer
w2 = net.LW{2}; % weights of w2 hidden to output layer

% color of graph
c = [1 0 0;
   0 1 0;
   0 0 1;
   1 0 1;
   0 1 1;
   0 0 0];

% Style of graph
g = {'-',':','--','-.'};
txt = cell([1 size(w1,1)]);
for i=1:size(w1,1)
   colorindex = rem(i-1,size(c,1))+1;
   styleindex = rem(i-1,size(g,2))+1;

   nodeout(i,:) = eval([TF1 '(w1(i)*x+b(i))*w2(i)']);
   
   t = plot(x,nodeout(i,:),g{styleindex});
   set(t,'color',c(colorindex,:));
   
   txt{i} = sprintf('Node %d',i);
   hold on
end
legend(txt)
title('Weighted Activation Functions of 1st Layer Nodes Before Training'); 
pause
clf

%Set no. of training rounds
net.trainParam.epochs = 2000;

%Set minimum error criterion
net.trainParam.goal = 0.0002;

% Training the network
net = train(net,x,y);
pause

% After training
w1 = net.IW{1}; % weights of w1 input to hidden layer
b  = net.b{1};  % bias b of hidden layer
w2 = net.LW{2}; % weights of w2 hidden to output layer

clf
% color of graph
c = [1 0 0;
   0 1 0;
   0 0 1;
   1 0 1;
   0 1 1;
   0 0 0];

% Style of graph
g = {'-',':','--','-.'};
txt = cell([1 size(w1,1)]);
for i=1:size(w1,1)
   colorindex = rem(i-1,size(c,1))+1;
   styleindex = rem(i-1,size(g,2))+1;

   nodeout(i,:) = eval(['w2(i)*' TF1 '(w1(i)*x+b(i))']);
   
   t = plot(x,nodeout(i,:),g{styleindex});
   set(t,'color',c(colorindex,:));
   
   txt{i} = sprintf('Node %d',i);
   hold on
end
legend(txt)
title('Weighted Activation Functions of 1st Layer Nodes After Training'); 
hold off
pause

nodesum = sum(nodeout,1);
plot(x,nodesum)
title('Sum of weighted output from all first layer nodes');
pause

%compute the network output after training
netout = sim(net,x);

% plot approximated function
plot(x,y,x,netout)
legend('Desired output','Network output')
title('Original Graph and Approximated Function Using the 2-Layer Network');
pause

