% backprop_sigmoid_xor.m
% Explain the backpropagation algorithm using a fully connected neural
% ntwork with one hidden layer. 
% Batch method is used.
% The activation function of each neuron is a sigmoid function
% h(z)=1/(1+exp(-z)), with the derivative being h'(z)=h(z).(1-h(z))
% The network can handle linearly non-separable cases, as in the XOR data.

% David Pan, UAH 

clear all;
close all;

alpha = 1;  % learning rate
max_iter = 1000;

% Linearly separable example (1,000 epochs are enough)
% Input data pattern
%X = [1 -1 -1 1; 1 -1 1 -1];
% Response
%R = [1  0  1 0; 0  1 0  1];

% Linearly non-separable example (with XOR pattern)
% Input data pattern
X = [1 -1 -1 1; 1 -1 1 -1];
% Response
R = [1  1  0 0; 0 0  1  1];

rng('default');
Std = 0.2;

% Initial weights and biases
W2 = Std*randn(2,2);
b2 = Std*randn(2,1);

W3 = Std*randn(2,2);
b3 = Std*randn(2,1);


mse = zeros(1, max_iter);

for epoch = 1: max_iter
    E = 0;
    W3_update = zeros(2,2);
    W2_update = zeros(2,2);
    b3_update = zeros(2,1);
    b2_update = zeros(2,1);
    
    for i = 1: 4
        A1 = X(:,i);
        % Z2 to replace A2 in 'backprop.m'
        Z2 = W2*A1 + b2;   % Z2 is the net input to the neuron
        % A2 is the output of the activation function on the input
        A2 = 1./(1+exp(-Z2));  
      
        Z3 = W3*A2 + b3;
        A3 = 1./(1+exp(-Z3));
        
        Deriv_A3 = A3 - R(:,i);
        E = E + 0.5*norm(Deriv_A3)^2;
        
        % Derivative of the activation function
        hd_Z3 = A3.*(1-A3);
        D3 = Deriv_A3 .* hd_Z3;
        
        % Backpropagation (now with sigmoid activation functions)
        hd_Z2 = A2.*(1-A2);
        D2 = W3'*D3.*hd_Z2;
        
        % Update the weights and biases
        W3_update = W3_update + alpha*D3*A2';
        W2_update = W2_update + alpha*D2*A1';
        
        b3_update  = b3_update + alpha*D3;
        b2_update  = b2_update + alpha*D2;    
    end
    
    % Update once after the whole batch of 4 pattern are processed
    W3 = W3 - W3_update;
    W2 = W2 - W2_update;
    b3 = b3 - b3_update;
    b2 = b2 - b2_update;
    
    mse(epoch) = E/4; 
    
end
mse(epoch)
plot(mse); grid
xlabel('Epochs');
ylabel('Mean Squared Error');

final_output = zeros(2, 4);

for i = 1: 4
    A1 = X(:,i);
        
    Z2 = W2*A1 + b2;
    A2 = 1./(1+exp(-Z2));

    Z3 = W3*A2 + b3;
    A3 = 1./(1+exp(-Z3));
    
    final_output(:,i) = A3;
end

final_output

% Compared to the desired response
R

% Scatter plot of the classified regions

figure;
hold on;
for x1 = -1.5:0.05:1.5
    for x2 = -1.5:0.05:1.5
        X_test = [x1; x2];
        A1 = X_test;
        Z2 = W2*A1 + b2;
        A2 = 1./(1+exp(-Z2));

        Z3 = W3*A2 + b3;
        A3 = 1./(1+exp(-Z3));
        
        if (A3(1)>=0.5)
            plot(x1, x2, 'r.');
        else
            plot(x1, x2, 'b.');
        end
    end
end
        
plot(X(1,1),X(2,1),'ro','MarkerSize',12, 'MarkerFaceColor','r');
plot(X(1,2),X(2,2),'ro','MarkerSize',12, 'MarkerFaceColor','r');
plot(X(1,3),X(2,3),'bo','MarkerSize',12, 'MarkerFaceColor','b');
plot(X(1,4),X(2,4),'bo','MarkerSize',12, 'MarkerFaceColor','b');


   