% This is just a MLP with 2 inputs, one output and 3 hidden nodes in its
% middle hidden layer
clear;clc;

samples=100;
bias=-1;

for n=1:samples
   x(n)=n/10;
   u(n)=x(n)+n/10;
   y(n)=u(n)+x(n);
end

% My own samples
x = [3 5 10];
u = [5 1 2];
y = [75 82 93];
% Save it for later
x_ = x;
u_ = u;
y_ = y;

% Feedforward part
% Random weights between the input and hidden layer
% First row corresponds to the input to the first hidden node
weights1 = rand(3,3);

% Random weights between the hidden layer and the output
weights2 = rand(1,4);
eta = -0.2;
learningTimes=500;

% One way to normalize
%x = (x-mean(x))/var(x);
%u = (u-mean(u))/var(u);
%y = (y-mean(y))/var(y);

% Another way to normalize
x = (x-min(x))/(max(x)-min(x));
u = (u-min(u))/(max(u)-min(u));
y = (y-min(y))/(max(y)-min(y));

for t=1:learningTimes

    biasColumn = bias*ones(length(x),1);
    input = [biasColumn x' u'];

    z2 = input*weights1';

    % Get the activation in each element of the array
    a2=arrayfun(@(x) activation(x,'sigmoid'),z2);
    % Add the bias to the current activation
    a2=[bias*ones(length(x),1) a2];

    z3=a2*weights2';

    y_hat=arrayfun(@(x) activation(x,'linear'), z3);

    J = 0.5*sum((y'-y_hat).^2)

    % Now we calculate the derivative of the squared error (J) respect to
    % weights2:

    d3 = -1*(y'-y_hat).*arrayfun(@(x) activation(x,'Dsigmoid'), z3);
    dJ2 = a2'*d3;

    % Let's remove the bias from the hidden layer as well
    dJ1_tmp = d3*weights2(:,2:end);
    dJ1 = input'*(dJ1_tmp.*arrayfun(@(x) activation(x,'Dsigmoid'), z2));

    weights1 = weights1 + eta*dJ1;
    weights2 = weights2 + (eta*dJ2)';
    
end

% Check if it works
% Let's say that we have normalization:
x1 = 3; % Input 1
x1_normalized = (x1-min(x_))/(max(x_)-min(x_));
x2 = 5; % Input 2
x2_normalized = (x2-min(u_))/(max(u_)-min(u_));
y1 = 75; % Output
y1_normalized = (y1-min(y_))/(max(y_)-min(y_));
input = [bias x1_normalized x2_normalized]
z2 = input*weights1';
a2=arrayfun(@(x) activation(x,'sigmoid'),z2);
a2=[bias a2]
z3=a2*weights2';
n_y_hat=arrayfun(@(x) activation(x,'linear'), z3)

result=n_y_hat*(max(y_)-min(y_))+min(y_);
disp('Final result');
disp(result);
J = 0.5*sum((y1_normalized-n_y_hat).^2)
