% T-61.5130 Machine Learning and Neural Networks
% Helsinki University of Technology
%
% Multilayer Perceptron (MLP) demo - generate a noisy
% data set and approximate the data using a neural network
%
% Matti Pöllä <matti.polla@tkk.fi> 2007-11-22


%% Generate noisy data
Ndat=100;
x=[0:8/(Ndat-1):8];
a=0.1;
y=sin(x)./(x+a);
t=y+0.01*randn(size(y));

xaxis = [0:.1:8];
plot(xaxis,sin(xaxis)./(xaxis+a),'-k',x,t,'b+')

%% Divide data to a training set and testing set
rind=randperm(Ndat);
x_train = x;
t_train = t;

%% Create a feed-forward MLP network
% net = newff(PR,[S1 S2...SNl],{TF1 TF2...TFNl},BTF,BLF,PF)
% PR - R x 2 matrix of min and max values for R input elements.
% Si - Size of ith layer, for Nl layers.
% TFi - Transfer function of ith layer, default = 'tansig'.
% BTF - Backpropagation network training function, default = 'traingdx'.
% BLF - Backpropagation weight/bias learning function, default = 'learngdm'.
% PF - Performance function, default = 'mse'.

% This is the number of hidden neurons, experiment with different values!
N=12;

net=newff([min(x),max(x)],[N 1],{'tansig' 'purelin'},'trainbfg','learngd','mse');

%% Initialize weights to a small value
% Hidden layer (input) weights
net.IW{1,1}
net.IW{1,1}=0.001*randn([N 1]);
net.IW{1,1}
% Output layer weights
net.LW{2,1}
net.LW{2,1}=0.001*randn([1 N]);
net.LW{2,1}
% Biases
net.b{1,1}=0.001*randn([N 1]);
net.b{2,1}=0.001*randn;

%% Train the network using the training data
%train(NET,P,T,Pi,Ai,VV,TV) takes,
%     net - Neural Network.
%     P - Network inputs.
%     T - Network targets, default = zeros.
%     Pi - Initial input delay conditions, default = zeros.
%     Ai - Initial layer delay conditions, default = zeros.
%     VV - Structure of validation vectors, default = [].
%     TV - Structure of test vectors, default = [].

% Maximum number of training rounds
net.trainParam.epochs = 200000;

% Maximum allowed error
net.trainParam.goal = 0.001;

[net,tr,Y,E]=train(net,x,t);


%% plot the results for training set and test set
xaxis = [0:.1:8];
plot(xaxis,sin(xaxis)./(xaxis+a),'-k',x,Y,'-r',x_train,t_train,'b+')
legend('Function y','Network output','Training samples')
