% This code implements a simple perceptron model to determine whether a valve should open or not based on three input features: Temperature, Pressure, and Concentration.
%% Input and Output Data

X = [
    0 0 1; 
    1 1 1; 
    1 0 1; 
    0 1 1
    ];                             % Inputs: [Temperature, Pressure, Concentration]
y = [0; 1; 1; 0];                  % Outputs: Valve position (1 = open, 0 = closed)

% X: This matrix contains the input data for training, where each row corresponds to a different observation, and each column represents a feature (Temperature, Pressure, Concentration).
% y: This vector contains the target output values corresponding to whether the valve should be open (1) or closed (0).
%% Initialization

rng(0);                % Set random seed for reproducibility
weights = rand(1, 3);  % Random initialization of weights
bias = rand(1);        % Random initialization of bias
learning_rate = 0.1;   % Learning rate
epochs = 1;            % Number of training iterations

% weights: Randomly initialized weights for each input feature.
% bias: A random bias term, added to the linear combination of inputs.
% learning_rate: A small factor that controls the step size during weight adjustment.
% epochs: Number of times the entire dataset is used to update the model’s weights.
%% Training the Perceptron

% Sigmoid function
sigmoid = @(x) 1 ./ (1 + exp(-x));

for epoch = 1:epochs
    % Forward pass
    weighted_sum = X * weights' + bias; 
    output = sigmoid(weighted_sum);

    % Calculate the error
    error = y - output;

    % Gradient descent
    % Calculate the derivative of the sigmoid function 
    % with respect to the weighted sum
    derivative = output .* (1 - output);

    % Update the weights and bias using the derivative
    % transpose back to row vector for consistent operations
    weights = weights + learning_rate * (error .* derivative)' * X;
    bias = bias + learning_rate * sum(error .* derivative);
end

% Sigmoid Function: Maps weighted inputs to a value between 0 and 1, representing the probability of belonging to a certain class.
% Epoch Loop: The outer loop runs for the specified number of epochs (iterations).
% Forward Pass: Calculates the weighted sum of inputs and biases, then applies the sigmoid function to get the model's output.
% Error Calculation: Measures how far the model's predictions are from the actual target values.
% Gradient Descent: Adjusts the weights and bias using the gradients of the error, helping the model learn and reduce the error over time.
% Weight and Bias Updates: These are iteratively refined across epochs to improve the model's predictions.
%% Testing the Perceptron

weighted_sum = X * weights' + bias; 
outputs = sigmoid(weighted_sum);
% outputs = sigmoid(X * weights') >= 0.5;

% The model is tested on the same dataset after training. The procedure is similar to the training process, but here, we only calculate the predictions without updating the weights and bias.
%% Results Display

disp('Predicted Valve Positions:');
disp(outputs);

disp('Final Weights:');
disp(weights);
disp('Final Bias:');
disp(bias);

%% With two hidden layers

% Input and Output Data
X = [
    0 0 1; 
    1 1 1; 
    1 0 1; 
    0 1 1
];                             % Inputs: [Temperature, Pressure, Concentration]
y = [0; 1; 1; 0];              % Outputs: Valve position (1 = open, 0 = closed)

% Initialization
rng(0);                        % Set random seed for reproducibility
input_layer_size = 3;          % Number of input features
hidden_layer1_size = 4;        % Number of neurons in the first hidden layer
hidden_layer2_size = 4;        % Number of neurons in the second hidden layer
output_layer_size = 1;         % Number of output neurons

% Random initialization of weights and biases
weights1 = rand(hidden_layer1_size, input_layer_size);   % Weights for input to first hidden layer
bias1 = rand(1, hidden_layer1_size);                     % Bias for first hidden layer

weights2 = rand(hidden_layer2_size, hidden_layer1_size); % Weights for first hidden layer to second hidden layer
bias2 = rand(1, hidden_layer2_size);                     % Bias for second hidden layer

weights3 = rand(output_layer_size, hidden_layer2_size);  % Weights for second hidden layer to output layer
bias3 = rand(1, output_layer_size);                      % Bias for output layer

learning_rate = 0.1;   % Learning rate
epochs = 20000;        % Number of training iterations

% Sigmoid function
sigmoid = @(x) 1 ./ (1 + exp(-x));

% Training the Neural Network
for epoch = 1:epochs
    % Forward pass
    z1 = X * weights1' + bias1;    % Input to first hidden layer
    a1 = sigmoid(z1);              % Output of first hidden layer
    
    z2 = a1 * weights2' + bias2;   % Input to second hidden layer
    a2 = sigmoid(z2);              % Output of second hidden layer
    
    z3 = a2 * weights3' + bias3;   % Input to output layer
    output = sigmoid(z3);          % Final output

    % Calculate the error
    error = y - output;

    % Backpropagation
    % Output layer
    delta3 = error .* (output .* (1 - output)); 
    
    % Second hidden layer
    delta2 = (delta3 * weights3) .* (a2 .* (1 - a2));
    
    % First hidden layer
    delta1 = (delta2 * weights2) .* (a1 .* (1 - a1));
    
    % Update the weights and biases using the derivatives
    weights3 = weights3 + learning_rate * delta3' * a2;
    bias3 = bias3 + learning_rate * sum(delta3);
    
    weights2 = weights2 + learning_rate * delta2' * a1;
    bias2 = bias2 + learning_rate * sum(delta2);
    
    weights1 = weights1 + learning_rate * delta1' * X;
    bias1 = bias1 + learning_rate * sum(delta1);
end

% Testing the Neural Network
z1 = X * weights1' + bias1;    
a1 = sigmoid(z1);              

z2 = a1 * weights2' + bias2;   
a2 = sigmoid(z2);              

z3 = a2 * weights3' + bias3;   
outputs = sigmoid(z3);         

% Results Display
disp('Predicted Valve Positions:');
disp(outputs);

% disp('Final Weights Layer 1:');
% disp(weights1);
% disp('Final Bias Layer 1:');
% disp(bias1);
% 
% disp('Final Weights Layer 2:');
% disp(weights2);
% disp('Final Bias Layer 2:');
% disp(bias2);
% 
% disp('Final Weights Layer 3:');
% disp(weights3);
% disp('Final Bias Layer 3:');
% disp(bias3);