# Things to note
- I prefer initialising weights at 0th epoch, and calculating loss at the end of each epoch, because that gives me the respective loss values of each 
weight (including the first), which ensures that I can plot the entire gradient descent process from first to last.
- During optimisation, I adjust gradient instead of learning rate, because Adam adjusts the learning rate differently in different directions, 
so it seems like dividing the gradient value makes sense because it has multiple components, while learning rate is a single scalar value. 
Unless this somehow results to a wrong value, I want it done this way.

The above two are my preferred design choices, please do not suggest me to change them (unless they’re giving wrong results, and the only 
solution is to not use them). Instead fix the code while maintaining this design choice.



# Introduction
I'm solving a simplified logical regression problem, that only uses weights w1 and w2 and no bias term. The data used here is generated synthetically 
to nicely fit into such model. (Not having a bias terms, reduces the parameters  to update during training to w1 and w2, making it possible for me to 
visualise the entire gradient descent process with only 3 dimension axes later on).

Theres two separate implementations of this here: 1st one does not use helper functions, and it is able to converge to 87% accuracy. While the second
implementation uses helper functions to reduce code repetititon, but something sneaky is going on in the second implementation, that is causing the 
accuracy to be stuck at 20% no matter how many epochs you train it for.



# Implementation 1: does not use helper functions, able to converge to 87% accuracy
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from sklearn.model_selection import train_test_split

"""# Set the data"""

# Generate Synthetic data
np.random.seed(42)
true_weights = np.array([-3, 3.5])
X = np.random.rand(500, len(true_weights)) * 10
y = (1 / (1 + np.exp(-np.dot(X, true_weights))) > 0.5).astype(int)
y[np.random.choice(500, int(0.1 * 500), replace=False)] ^= 1

# Split the data into training (60%), validation (20%), and test (20%) sets
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.4, random_state=42)
X_val, X_test, y_val, y_test = train_test_split(X_test, y_test, test_size=0.5, random_state=42)

"""# Gradient Descent Class"""

class GradientDescent:
    def __init__(self, optimizer=None):
        self.weights = None
        self.learning_rate = None
        self.optimizer = optimizer

    def fit(self, X_train, y_train, X_val, y_val, epochs=100, initial_weights=None, learning_rate=0.002, batch_size=32, verbose=True):
        self.learning_rate = learning_rate

        # Initialize lists to store metrics for each epoch
        train_losses, train_accuracies, val_losses, val_accuracies, weight_history = [], [], [], [], []

        for epoch in range(epochs):
            if epoch == 0:
                np.random.seed(42)
                self.weights = np.array(initial_weights) if initial_weights else np.random.uniform(-1, 1, num_features)
            else:
                for i in range(0, X_train.shape[0], batch_size):
                    X_batch = X_train[i:i + batch_size]
                    y_true_batch = y_train[i:i + batch_size]

                    # Compute linear combination and apply sigmoid
                    linear_combination = np.dot(X_batch, self.weights)
                    y_pred_batch = 1 / (1 + np.exp(-linear_combination))

                    # Compute gradient of binary-crossentropy loss w.r.t weights
                    grad_loss = np.dot(X_batch.T, (y_pred_batch - y_true_batch)) / X_batch.shape[0]

                    # Apply optimizer
                    if self.optimizer:
                        grad_loss = self.optimizer.update(grad_loss)

                    # Update weights
                    self.weights -= self.learning_rate * grad_loss

            # Training evaluation after the epoch
            train_loss, train_accuracy = self.evaluate(X_train, y_train)
            val_loss, val_accuracy = self.evaluate(X_val, y_val)

            # Store metrics & display the metrics for the current epoch
            weight_history.append(self.weights.copy()); train_losses.append(train_loss); train_accuracies.append(train_accuracy); val_losses.append(val_loss); val_accuracies.append(val_accuracy)
            if verbose:
                print(f"Epoch {epoch + 1}:\nWeights: {self.weights}\nTrain Loss: {train_loss:.8f}, Train Accuracy: {train_accuracy:.4f}\nValidation Loss: {val_loss:.8f}, Validation Accuracy: {val_accuracy:.4f}\n")

        return train_losses, train_accuracies, val_losses, val_accuracies, weight_history

    def evaluate(self, X, y_true):
        linear_combination = np.dot(X, self.weights)
        y_pred = 1 / (1 + np.exp(-linear_combination))
        loss = -np.mean(y_true * np.log(y_pred + 1e-8) + (1 - y_true) * np.log(1 - y_pred + 1e-8)) # Cross-entropy loss
        accuracy = np.mean((y_pred >= 0.5) == y_true)
        return loss, accuracy

"""# Define the optimizer"""

class Adam:
    def __init__(self, beta1=0.9, beta2=0.999, epsilon=1e-7):
        self.beta1 = beta1
        self.beta2 = beta2
        self.epsilon = epsilon
        self.m_t = None  # First moment estimate (momentum)
        self.v_t = None  # Second moment estimate (RMSProp-like)
        self.t = 0       # Timestep for bias correction

    def update(self, gradient):
        if self.m_t is None:
            self.m_t = np.zeros_like(gradient)
        if self.v_t is None:
            self.v_t = np.zeros_like(gradient)

        self.t += 1

        # Update biased first and second moment estimates
        self.m_t = self.beta1 * self.m_t + (1 - self.beta1) * gradient
        self.v_t = self.beta2 * self.v_t + (1 - self.beta2) * (gradient ** 2)

        # Bias-corrected moment estimates
        m_hat = self.m_t / (1 - self.beta1 ** self.t)
        v_hat = self.v_t / (1 - self.beta2 ** self.t)

        # Update weights
        adjusted_gradient = m_hat / (np.sqrt(v_hat) + self.epsilon)

        return adjusted_gradient

"""# Perform Gradient Descent
Here the accuracy slowly converges to around 87%
"""

gd = GradientDescent(optimizer=Adam())

# slowly converges at around 90th epoch (learning rate=0.03)
metrics = gd.fit(X_train, y_train, X_val, y_val, epochs=100, batch_size=32, initial_weights=[15.0, -20.0], learning_rate=0.03)

# slowly converges at around 130th epoch (learning rate=0.02)
metrics = gd.fit(X_train, y_train, X_val, y_val, epochs=200, batch_size=32, initial_weights=[15.0, -20.0], learning_rate=0.02)

# slowly converges at around 2540th epoch (learning rate=0.001)
metrics = gd.fit(X_train, y_train, X_val, y_val, epochs=5000, batch_size=32, initial_weights=[15.0, -20.0], learning_rate=0.001)



# Implementation 2: uses helper functions, accuracy stuck at 20%
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from sklearn.model_selection import train_test_split

"""# Set the data"""

# Generate Synthetic data
np.random.seed(42)
true_weights = np.array([-3, 3.5])
X = np.random.rand(500, len(true_weights)) * 10
y = (1 / (1 + np.exp(-np.dot(X, true_weights))) > 0.5).astype(int)
y[np.random.choice(500, int(0.1 * 500), replace=False)] ^= 1

# Split the data into training (60%), validation (20%), and test (20%) sets
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.4, random_state=42)
X_val, X_test, y_val, y_test = train_test_split(X_test, y_test, test_size=0.5, random_state=42)

"""# Define Helper Functions"""

# We will be using binary crossentrophy as the loss function to minimize
class BinaryCrossEntropy:
    @staticmethod
    def value(y_true, y_pred):
        # Average loss over all of the samples
        return -np.mean(y_true * np.log(y_pred + 1e-8) + (1 - y_true) * np.log(1 - y_pred + 1e-8))

    @staticmethod
    def gradient(y_true, y_pred):
        # Gradient of loss function w.r.t predictions
        return (y_pred - y_true) / (y_pred * (1 - y_pred) + 1e-8)

# Sigmoid function
def sigmoid(x):
    x = np.clip(x, -500, 500)
    return 1 / (1 + np.exp(-x))

# Given weights, computes the prediction of our model
def compute_predictions(X, weights):
    linear_combination = np.dot(X, weights)
    return sigmoid(linear_combination)

# Calculates accuracy as the proportion of correct predictions
def accuracy(y_true, y_pred):
    return np.mean((y_pred >= 0.5) == y_true)

"""# Gradient Descent Class"""

class GradientDescent:
    def __init__(self, loss_function=BinaryCrossEntropy, accuracy=accuracy, optimizer=None):
        self.weights = None
        self.learning_rate = None
        self.loss_f = loss_function
        self.accuracy = accuracy
        self.optimizer = optimizer

    def fit(self, X_train, y_train, X_val, y_val, epochs=100, initial_weights=None, learning_rate=0.003, batch_size=32, verbose=True):
        self.learning_rate = learning_rate
        num_samples, num_features = X_train.shape

        # Initialize lists to store metrics for each epoch
        train_losses, train_accuracies, val_losses, val_accuracies = [], [], [], []
        weight_history = []

        for epoch in range(epochs):
            # If it's the first epoch, initialize weights
            if epoch == 0:
                np.random.seed(42)
                self.weights = np.array(initial_weights) if initial_weights else np.random.uniform(-1, 1, num_features)
            # In the beginning of other epochs, we update the weights
            else:
                # Process data in mini-batches
                for i in range(0, num_samples, batch_size):
                    X_batch = X_train[i:i + batch_size]
                    y_true_batch = y_train[i:i + batch_size]

                    # Calculate predictions for the batch
                    y_predicted_batch = compute_predictions(X_batch, self.weights)

                    # Gradient of loss function w.r.t. predicted y
                    grad_loss_batch = self.loss_f.gradient(y_true_batch, y_predicted_batch)

                    # Use chain rule to get gradient of loss function w.r.t. weights
                    delta_w_batch = np.dot(X_batch.T, grad_loss_batch) / batch_size

                    # Get the adjusted gradient from the optimizer
                    if self.optimizer:
                        delta_w_batch = self.optimizer.update(delta_w_batch)

                    # Update the weights using that gradient
                    self.weights -= self.learning_rate * delta_w_batch

            # Make predictions for the whole training set after all updates in the epoch
            y_predicted = compute_predictions(X_train, self.weights)

            # Compute the binary cross-entropy loss and accuracy for training data
            train_loss = self.loss_f.value(y_train, y_predicted)
            train_accuracy = self.accuracy(y_train, y_predicted)

            # Validation
            val_y_predicted = compute_predictions(X_val, self.weights)
            val_loss = self.loss_f.value(y_val, val_y_predicted)
            val_accuracy = self.accuracy(y_val, val_y_predicted)

            # Store the weights and corresponding metrics
            weight_history.append(self.weights.copy())
            train_losses.append(train_loss)
            train_accuracies.append(train_accuracy)
            val_losses.append(val_loss)
            val_accuracies.append(val_accuracy)

            # Display the metrics for the current epoch
            if verbose:
                print(f"Epoch {epoch + 1}:")
                print(f"  Weights: {self.weights}")
                print(f"  Train Loss: {train_loss:.8f}, Train Accuracy: {train_accuracy:.4f}")
                print(f"  Validation Loss: {val_loss:.8f}, Validation Accuracy: {val_accuracy:.4f}")
                print("--------------------------")

        # Return all the training metrics
        return train_losses, train_accuracies, val_losses, val_accuracies, weight_history

    def predict(self, X):
        return compute_predictions(X, self.weights)

"""# Define the optimizer"""

class Adam:
    def __init__(self, learning_rate=0.001, beta1=0.9, beta2=0.999, epsilon=1e-8):
        self.learning_rate = learning_rate
        self.beta1 = beta1
        self.beta2 = beta2
        self.epsilon = epsilon
        self.m_t = None  # First moment estimate (momentum)
        self.v_t = None  # Second moment estimate (RMSProp-like)
        self.t = 0       # Timestep for bias correction

    def update(self, gradient):
        if self.m_t is None:
            self.m_t = np.zeros_like(gradient)
        if self.v_t is None:
            self.v_t = np.zeros_like(gradient)

        self.t += 1

        # Update biased first and second moment estimates
        self.m_t = self.beta1 * self.m_t + (1 - self.beta1) * gradient
        self.v_t = self.beta2 * self.v_t + (1 - self.beta2) * (gradient ** 2)

        # Bias-corrected moment estimates
        m_hat = self.m_t / (1 - self.beta1 ** self.t)
        v_hat = self.v_t / (1 - self.beta2 ** self.t)

        # Update weights
        adjusted_gradient = m_hat / (np.sqrt(v_hat) + self.epsilon)

        return self.learning_rate * adjusted_gradient

"""# Perform Gradient Descent
Here the accuracy never improves beyoned 20%
"""

gd = GradientDescent(optimizer=Adam())

# Does not converge even in 200th epoch
metrics = gd.fit(X_train, y_train, X_val, y_val, epochs=100, batch_size=32, initial_weights=[15.0, -20.0], learning_rate=0.03)

# Does not converge even in 200th epoch
metrics = gd.fit(X_train, y_train, X_val, y_val, epochs=200, batch_size=32, initial_weights=[15.0, -20.0], learning_rate=0.02)

# does not converge even in 500th epoch (learning rate=0.001)
metrics = gd.fit(X_train, y_train, X_val, y_val, epochs=5000, batch_size=32, initial_weights=[15.0, -20.0], learning_rate=0.001)