Make minimal changes to DenseNeuralNetwork class here to implement batch gradient descent (instead of SGD)
In batch gradient descent:-
- the whole training data is treated as a single batch, 
- loss is averaged for each of instance of the training data,
- one weight update occurs in the entire epoch based on that average loss.

class DenseNeuralNetwork:
    def __init__(self, layer_sizes, activations, loss_f):
        self.learning_rate = None      # Learning rate
        self.layer_sizes = layer_sizes # Array storing the layer sizes
        self.loss_f = loss_f           # Loss function

        # The array of weights, biases, net inputs, and activation functions have their first element as None to represent that
        # input layers don't have them, this also helps match the indexes of expressions here match with their mathematical formulas

        # To store weights(W) and biases(b)
        np.random.seed(42)
        self.W = [None] + [np.random.uniform(-1, 1, (layer_sizes[i], layer_sizes[i-1])) for i in range(1, len(layer_sizes))]
        self.b = [None] + [np.random.uniform(-1, 1, layer_sizes[i]) for i in range(1, len(layer_sizes))]

        # To store activations(a) and net inputs(z)
        self.a = [np.zeros(size) for size in layer_sizes] # This one does not start with None, because input layer activates with same value as the inputs
        self.z = [None] + [np.zeros(layer_sizes[i]) for i in range(1, len(layer_sizes))]

        # Set the activation functions for each layer
        self.activation_f = activations # The first None element is passed as input here

    def propagate_forward(self, x):
        # Set input layer activations
        self.a[0] = x

        for i in range(1, len(self.W)):
            self.z[i] = self.W[i] @ self.a[i-1] + self.b[i]  # Weighted sum
            self.a[i] = self.activation_f[i].function(self.z[i])  # Activation

    def propagate_backward(self, y_true):
        error_term = {}

        # Store the index of the last layer as L
        L = len(self.W) - 1

        # For output layer:-
        if self.activation_f[L] == Softmax:
            error_term[L] = self.activation_f[L].derivative(self.z[L]) @ self.loss_f.gradient(y_true, self.a[L])
            # print('LAST LAYER ERROR TERM (FROM MATRIX MULTIPLICATION): ', error_term[L]) # For debugging
        else:
            error_term[L] = self.activation_f[L].derivative(self.z[L]) * self.loss_f.gradient(y_true, self.a[L])

        # Update weights and biases for output layer  based on those error terms
        self.W[L] -= self.learning_rate * np.outer(error_term[L], self.a[L-1])
        self.b[L] -= self.learning_rate * error_term[L]

        # For hidden layers:-
        for i in reversed(range(1, L)):
            ## Propgating error terms backward
            if  self.activation_f[i] == Softmax:
                error_term[i] = self.activation_f[i].derivative(self.z[i]) @ (self.W[i+1].T @ error_term[i+1])
            else:
                error_term[i] = self.activation_f[i].derivative(self.z[i]) * (self.W[i+1].T @ error_term[i+1])

            ## Updating the weights and biases based on those error terms
            self.W[i] -= self.learning_rate * np.outer(error_term[i], self.a[i-1])
            self.b[i] -= self.learning_rate * error_term[i]

    # Training loop
    def train(self, X_train, y_train, X_val, y_val, epochs=100, learning_rate=0.01, verbose=True):
        training_losses, validation_losses, training_accuracies, validation_accuracies = [], [], [], []
        self.learning_rate = learning_rate

        for epoch in range(epochs):
            total_loss, correct_predictions = 0, 0

            for i in range(X_train.shape[0]):
                x, y_true = X_train[i], y_train[i]
                self.propagate_forward(x)
                total_loss += self.loss_f.function(y_true, self.a[-1])

                if self.is_correct_prediction(self.a[-1], y_true):
                    correct_predictions += 1

                self.propagate_backward(y_true)

            avg_loss = total_loss / X_train.shape[0]
            accuracy = correct_predictions / X_train.shape[0] * 100
            training_losses.append(avg_loss)
            training_accuracies.append(accuracy)

            # Validation phase
            val_loss, val_correct_predictions = 0, 0
            for i in range(X_val.shape[0]):
                x, y_true = X_val[i], y_val[i]
                self.propagate_forward(x)
                val_loss += self.loss_f.function(y_true, self.a[-1])

                # Use np.argmax for validation comparison as well
                if self.is_correct_prediction(self.a[-1], y_true):
                    val_correct_predictions += 1

            val_avg_loss = val_loss / X_val.shape[0]
            val_accuracy = val_correct_predictions / X_val.shape[0] * 100
            validation_losses.append(val_avg_loss)
            validation_accuracies.append(val_accuracy)

            if verbose:
                print(f"Epoch {epoch+1}")
                print(f"Training Loss = {avg_loss}, Validation Loss = {val_avg_loss}")
                print(f"Training Accuracy: {accuracy:.3f}%, Validation Accuracy: {val_accuracy:.3f}%")
                print("----------------------------------------------------------------------")

        return training_losses, validation_losses, training_accuracies, validation_accuracies

    # Checks if the prediction has been made correctly; used to evalutue accuracy
    # This function may need to be altered based on the type of problem being solved
    def is_correct_prediction(self, predicted, true):
        if len(predicted) == 1:
            # For binary classification, if neuron o/p > 0.5 we consider it postive class
            return ((predicted > 0.5).astype(int)) == true
        else:
            # For multiclass classification
            # Use np.argmax to convert predicted and true values to class labels for comparison
            return np.argmax(predicted) == np.argmax(true)

    # For making predctions after training the network
    def predict(self, x):
        # Forward propagate through all layers
        self.propagate_forward(x)

        # The output of the last layer is the prediction
        y_pred = self.a[-1]

        return y_pred

    # For evaluating test accuracy later on
    def get_test_accuracy(self, X_test, y_test):
        correct_predictions = 0
        total_loss = 0

        for i in range(X_test.shape[0]):
            x, y_true = X_test[i], y_test[i]
            self.propagate_forward(x)

            # Accumulate loss
            total_loss += self.loss_f.function(y_true, self.a[-1])

            # Check if predicted class matches the true class
            if self.is_correct_prediction(self.a[-1], y_true):
                correct_predictions += 1

        # Calculate accuracy
        accuracy = correct_predictions / X_test.shape[0] * 100

        # Average loss
        avg_loss = total_loss / X_test.shape[0]

        print(f"Test Loss: {avg_loss}")
        print(f"Test Accuracy: {accuracy:.3f}%")

        return avg_loss, accuracy


Only for reference - this is how I performed batch gradient descent for logistic regressison:
class BatchGradientDescent:
    def __init__(self, loss_function=BinaryCrossEntropy, accuracy=accuracy):
        self.weights = None
        self.learning_rate = None
        self.loss_f = loss_function
        self.accuracy = accuracy

    def fit(self, X_train, y_train, X_val, y_val, epochs=100, initial_weights=None, learning_rate = 0.003, 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:
                # Calculate predictions
                y_predicted = compute_predictions(X_train, self.weights)
                # Gradient of loss function w.r.t. weights
                delta_w = self.loss_f.gradient_wrt_weights(y_train, y_predicted, X_train)
                # Update the weights based on that gradient
                self.weights -= self.learning_rate * delta_w

            # Make predictions from the current weights
            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)