%The data set
x1 = [0; 0; 1; 1];
x2 = [0; 1; 0; 1];
y =  [0; 0; 0; 1];

%Initialisation weights
W1 = stdnormal_rnd(2,3); 
W2 = stdnormal_rnd(3,1);
b1 = stdnormal_rnd(1,3); 
b2 = stdnormal_rnd(1,1);

%Initialisation variables
alpha = 1;
labda = 0.01;
max_iterations = 5000;
m = rows(x1);

%sigmoid function
function s = f(z) 
  s = 1./(1 + exp(-z));
end


for i = 1:maxi 
    
    %Forward propagation
    
    % layer 1 (input)
    a1 = [exx1, exx2];

    % layer 2
    z2 = a1 * W1 + b1;
    a2 = f(z2);
    
    % layer 3 (output)
    z3 = a2 * W2 + b2;
    h = f(z3);
    a3 = h;
    
    %Backward propagation
    
    % Delta values
    delta3 = -(exy - h)     .* (a3 .* (1 - a3)); % f'(z) = f(z)(1 − f(z)). See slide: http://i.imgur.com/K5QHTK9.png
    delta2 = (delta3 * W2') .* (a2 .* (1 - a2)); 

    %Update functions
    W2 = W2 - alpha * (a2' * delta3 * 1/m);  % + lambda * W2);
    W1 = W1 - alpha * (a1' * delta2 * 1/m);  % + lambda * W1);
    b2 = b2 - alpha * (sum(delta3) * 1/m);
    b1 = b1 - alpha * (sum(delta2) * 1/m);
    
end



%Forward propagation

% layer 1 (input)
a1 = [exx1, exx2];

% layer 2
z2 = a1 * W1 + b1;
a2 = f(z2);

% layer 3 (output)
z3 = a2 * W2 + b2;
h = f(z3);
a3 = h;

%Backward propagation

% Delta values
delta3 = -(exy - h)     .* (a3 .* (1 - a3)); % f'(z) = f(z)(1 − f(z)). See slide: http://i.imgur.com/K5QHTK9.png
delta2 = (delta3 * W2') .* (a2 .* (1 - a2)); 

%Update terms
deltaW2 = a2' * delta3;
deltaW1 = a1' * delta2;
deltab2 = sum(delta3);
deltab1 = sum(delta2);

%Update functions
W2 = W2 - alpha * (gradW2 * 1/m);  % + lambda * W2);
W1 = W1 - alpha * (gradW1 * 1/m);  % + lambda * W1);
b2 = b2 - alpha * (gradb2 * 1/m);
b1 = b1 - alpha * (gradb1 * 1/m);