Backpropagation &
Training Deep Networks

import numpy as np

# ==========================================
# COMPLETE NEURAL NETWORK WITH BACKPROPAGATION
# ==========================================

class DeepNeuralNetwork:
    """
    L-layer neural network with backpropagation.
    Supports: ReLU, Sigmoid, Tanh activations
    Optimizers: SGD, Momentum, RMSprop, Adam
    Regularization: L2, Dropout, Batch Normalization
    """

    def __init__(self, layer_dims, activations=None):
        """
        layer_dims: list of integers [n_x, n_h1, n_h2, ..., n_y]
        activations: list of strings ['relu', 'relu', ..., 'sigmoid']
        """
        self.L = len(layer_dims) - 1  # Number of layers
        self.layer_dims = layer_dims
        self.parameters = {}
        self.caches = {}
        self.grads = {}

        # Default activations: ReLU hidden, Sigmoid output
        if activations is None:
            self.activations = ['relu'] * (self.L - 1) + ['sigmoid']
        else:
            self.activations = activations

        self._initialize_parameters()

    def _initialize_parameters(self):
        """He initialization for ReLU, Xavier for Sigmoid/Tanh"""
        np.random.seed(42)
        for l in range(1, self.L + 1):
            if self.activations[l-1] == 'relu':
                # He initialization: scale by sqrt(2/n_prev)
                scale = np.sqrt(2.0 / self.layer_dims[l-1])
            else:
                # Xavier initialization: scale by sqrt(1/n_prev)
                scale = np.sqrt(1.0 / self.layer_dims[l-1])

            self.parameters[f'W{l}'] = np.random.randn(
                self.layer_dims[l], self.layer_dims[l-1]
            ) * scale
            self.parameters[f'b{l}'] = np.zeros((self.layer_dims[l], 1))

    # ---- Activation Functions ----

    def _sigmoid(self, Z):
        A = 1.0 / (1.0 + np.exp(-np.clip(Z, -500, 500)))
        return A

    def _relu(self, Z):
        return np.maximum(0, Z)

    def _tanh(self, Z):
        return np.tanh(Z)

    def _activation_forward(self, Z, activation):
        if activation == 'sigmoid':
            return self._sigmoid(Z)
        elif activation == 'relu':
            return self._relu(Z)
        elif activation == 'tanh':
            return self._tanh(Z)

    def _activation_backward(self, dA, Z, activation):
        if activation == 'sigmoid':
            s = self._sigmoid(Z)
            return dA * s * (1 - s)
        elif activation == 'relu':
            return dA * (Z > 0).astype(float)
        elif activation == 'tanh':
            return dA * (1 - np.tanh(Z) ** 2)

    # ---- Forward Pass ----

    def forward(self, X, keep_prob=1.0):
        """
        Full forward pass through L layers.
        Stores caches for backpropagation.
        """
        self.caches = {}
        A = X
        self.caches['A0'] = X

        for l in range(1, self.L + 1):
            A_prev = A
            W = self.parameters[f'W{l}']
            b = self.parameters[f'b{l}']

            # Linear step: Z = W·A_prev + b
            Z = np.dot(W, A_prev) + b
            self.caches[f'Z{l}'] = Z

            # Activation step
            A = self._activation_forward(Z, self.activations[l-1])

            # Dropout (not on output layer)
            if keep_prob < 1.0 and l < self.L:
                D = (np.random.rand(*A.shape) < keep_prob).astype(float)
                A = A * D / keep_prob  # Inverted dropout
                self.caches[f'D{l}'] = D

            self.caches[f'A{l}'] = A

        return A

    # ---- Loss Computation ----

    def compute_cost(self, AL, Y, lambd=0):
        """Binary cross-entropy with optional L2 regularization"""
        m = Y.shape[1]
        # Cross-entropy loss
        cost = -(1/m) * np.sum(
            Y * np.log(AL + 1e-8) + (1 - Y) * np.log(1 - AL + 1e-8)
        )

        # L2 regularization
        if lambd > 0:
            l2_cost = 0
            for l in range(1, self.L + 1):
                l2_cost += np.sum(np.square(self.parameters[f'W{l}']))
            cost += (lambd / (2 * m)) * l2_cost

        return np.squeeze(cost)

    # ---- Backward Pass ----

    def backward(self, Y, lambd=0, keep_prob=1.0):
        """
        Full backward pass — computes all gradients.
        This is the CORE of backpropagation.
        """
        m = Y.shape[1]
        AL = self.caches[f'A{self.L}']

        # Output layer: dZ[L] = A[L] - Y (sigmoid + BCE)
        dZ = AL - Y  # Shape: (n[L] × m)

        for l in range(self.L, 0, -1):
            A_prev = self.caches[f'A{l-1}']

            # Gradient of weights: dW = (1/m) · dZ · A_prev^T
            dW = (1/m) * np.dot(dZ, A_prev.T)

            # Add L2 regularization gradient
            if lambd > 0:
                dW += (lambd / m) * self.parameters[f'W{l}']

            # Gradient of biases: db = (1/m) · sum(dZ, axis=1)
            db = (1/m) * np.sum(dZ, axis=1, keepdims=True)

            # Store gradients
            self.grads[f'dW{l}'] = dW
            self.grads[f'db{l}'] = db

            # Propagate gradient to previous layer (if not at layer 1)
            if l > 1:
                # dA_prev = W^T · dZ
                dA_prev = np.dot(self.parameters[f'W{l}'].T, dZ)

                # Apply dropout mask
                if keep_prob < 1.0 and f'D{l-1}' in self.caches:
                    dA_prev = dA_prev * self.caches[f'D{l-1}'] / keep_prob

                # dZ_prev = dA_prev ⊙ g'(Z_prev)
                dZ = self._activation_backward(
                    dA_prev,
                    self.caches[f'Z{l-1}'],
                    self.activations[l-2]
                )

    # ---- Optimizers ----

    def update_parameters_sgd(self, learning_rate):
        """Vanilla SGD"""
        for l in range(1, self.L + 1):
            self.parameters[f'W{l}'] -= learning_rate * self.grads[f'dW{l}']
            self.parameters[f'b{l}'] -= learning_rate * self.grads[f'db{l}']

    def initialize_velocity(self):
        """Initialize velocity for Momentum"""
        self.v = {}
        for l in range(1, self.L + 1):
            self.v[f'dW{l}'] = np.zeros_like(self.parameters[f'W{l}'])
            self.v[f'db{l}'] = np.zeros_like(self.parameters[f'b{l}'])

    def update_parameters_momentum(self, learning_rate, beta=0.9):
        """SGD with Momentum"""
        for l in range(1, self.L + 1):
            # Update velocity
            self.v[f'dW{l}'] = beta * self.v[f'dW{l}'] + (1 - beta) * self.grads[f'dW{l}']
            self.v[f'db{l}'] = beta * self.v[f'db{l}'] + (1 - beta) * self.grads[f'db{l}']
            # Update parameters
            self.parameters[f'W{l}'] -= learning_rate * self.v[f'dW{l}']
            self.parameters[f'b{l}'] -= learning_rate * self.v[f'db{l}']

    def initialize_adam(self):
        """Initialize Adam moment estimates"""
        self.m = {}  # First moment
        self.v_adam = {}  # Second moment
        for l in range(1, self.L + 1):
            self.m[f'dW{l}'] = np.zeros_like(self.parameters[f'W{l}'])
            self.m[f'db{l}'] = np.zeros_like(self.parameters[f'b{l}'])
            self.v_adam[f'dW{l}'] = np.zeros_like(self.parameters[f'W{l}'])
            self.v_adam[f'db{l}'] = np.zeros_like(self.parameters[f'b{l}'])

    def update_parameters_adam(self, lr, t, beta1=0.9, beta2=0.999, eps=1e-8):
        """Adam optimizer with bias correction"""
        for l in range(1, self.L + 1):
            for param in ['dW', 'db']:
                key = f'{param}{l}'
                # Update first moment (mean of gradients)
                self.m[key] = beta1 * self.m[key] + (1-beta1) * self.grads[key]
                # Update second moment (mean of squared gradients)
                self.v_adam[key] = beta2 * self.v_adam[key] + (1-beta2) * (self.grads[key]**2)
                # Bias correction
                m_corrected = self.m[key] / (1 - beta1**t)
                v_corrected = self.v_adam[key] / (1 - beta2**t)
                # Update parameter
                p_key = f'W{l}' if 'W' in param else f'b{l}'
                self.parameters[p_key] -= lr * m_corrected / (np.sqrt(v_corrected) + eps)

    # ---- Gradient Checking ----

    def gradient_check(self, X, Y, epsilon=1e-7):
        """Verify backprop using numerical gradients"""
        # Compute analytic gradients
        self.forward(X)
        self.backward(Y)

        # Flatten all parameters and gradients
        params_values = []
        grads_values = []
        for l in range(1, self.L + 1):
            params_values.append(self.parameters[f'W{l}'].flatten())
            params_values.append(self.parameters[f'b{l}'].flatten())
            grads_values.append(self.grads[f'dW{l}'].flatten())
            grads_values.append(self.grads[f'db{l}'].flatten())

        theta = np.concatenate(params_values)
        d_theta = np.concatenate(grads_values)

        # Compute numerical gradients
        num_grads = np.zeros_like(theta)
        for i in range(len(theta)):
            # θ+
            theta_plus = theta.copy()
            theta_plus[i] += epsilon
            self._set_params_from_vector(theta_plus)
            AL_plus = self.forward(X)
            cost_plus = self.compute_cost(AL_plus, Y)

            # θ-
            theta_minus = theta.copy()
            theta_minus[i] -= epsilon
            self._set_params_from_vector(theta_minus)
            AL_minus = self.forward(X)
            cost_minus = self.compute_cost(AL_minus, Y)

            num_grads[i] = (cost_plus - cost_minus) / (2 * epsilon)

        # Restore original parameters
        self._set_params_from_vector(theta)

        # Compute difference
        diff = np.linalg.norm(d_theta - num_grads) / \
               (np.linalg.norm(d_theta) + np.linalg.norm(num_grads) + 1e-8)

        if diff < 1e-7:
            print(f"✅ Gradient check PASSED! Difference: {diff:.2e}")
        elif diff < 1e-5:
            print(f"⚠️  Gradient check SUSPICIOUS. Difference: {diff:.2e}")
        else:
            print(f"❌ Gradient check FAILED! Difference: {diff:.2e}")

        return diff

    def _set_params_from_vector(self, theta):
        """Helper: set parameters from flat vector"""
        idx = 0
        for l in range(1, self.L + 1):
            w_shape = self.parameters[f'W{l}'].shape
            w_size = np.prod(w_shape)
            self.parameters[f'W{l}'] = theta[idx:idx+w_size].reshape(w_shape)
            idx += w_size
            b_shape = self.parameters[f'b{l}'].shape
            b_size = np.prod(b_shape)
            self.parameters[f'b{l}'] = theta[idx:idx+b_size].reshape(b_shape)
            idx += b_size

    # ---- Full Training Loop ----

    def train(self, X, Y, epochs=1000, lr=0.01, optimizer='adam',
              batch_size=32, lambd=0, keep_prob=1.0,
              lr_decay=None, print_cost=True):
        """Complete training loop with mini-batches"""
        m = X.shape[1]
        costs = []
        t = 0  # Adam counter

        # Initialize optimizer state
        if optimizer == 'momentum':
            self.initialize_velocity()
        elif optimizer == 'adam':
            self.initialize_adam()

        for epoch in range(epochs):
            # Shuffle data
            permutation = np.random.permutation(m)
            X_shuffled = X[:, permutation]
            Y_shuffled = Y[:, permutation]

            # Create mini-batches
            num_batches = m // batch_size
            epoch_cost = 0

            for k in range(num_batches + (1 if m % batch_size != 0 else 0)):
                X_batch = X_shuffled[:, k*batch_size : (k+1)*batch_size]
                Y_batch = Y_shuffled[:, k*batch_size : (k+1)*batch_size]

                if X_batch.shape[1] == 0:
                    continue

                # Forward pass
                AL = self.forward(X_batch, keep_prob)

                # Compute cost
                batch_cost = self.compute_cost(AL, Y_batch, lambd)
                epoch_cost += batch_cost

                # Backward pass
                self.backward(Y_batch, lambd, keep_prob)

                # Update parameters
                t += 1
                current_lr = self._get_lr(lr, epoch, epochs, lr_decay)

                if optimizer == 'sgd':
                    self.update_parameters_sgd(current_lr)
                elif optimizer == 'momentum':
                    self.update_parameters_momentum(current_lr)
                elif optimizer == 'adam':
                    self.update_parameters_adam(current_lr, t)

            avg_cost = epoch_cost / max(num_batches, 1)
            costs.append(avg_cost)

            if print_cost and epoch % 100 == 0:
                print(f"Epoch {epoch:4d} | Cost: {avg_cost:.6f} | LR: {current_lr:.6f}")

        return costs

    def _get_lr(self, lr_init, epoch, total_epochs, schedule):
        """Learning rate scheduler"""
        if schedule is None:
            return lr_init
        elif schedule == 'step':
            return lr_init * (0.5 ** (epoch // 200))
        elif schedule == 'exponential':
            return lr_init * np.exp(-0.01 * epoch)
        elif schedule == 'cosine':
            return lr_init * 0.5 * (1 + np.cos(np.pi * epoch / total_epochs))
        elif schedule == 'warmup_cosine':
            warmup = total_epochs // 10
            if epoch < warmup:
                return lr_init * (epoch / warmup)
            return lr_init * 0.5 * (1 + np.cos(np.pi * (epoch - warmup) / (total_epochs - warmup)))

    def predict(self, X, threshold=0.5):
        """Make predictions"""
        AL = self.forward(X, keep_prob=1.0)
        return (AL > threshold).astype(int)

    def accuracy(self, X, Y):
        preds = self.predict(X)
        return np.mean(preds == Y) * 100


# ==========================================
# DEMO: Train on synthetic data
# ==========================================
if __name__ == "__main__":
    # Generate spiral dataset (non-linearly separable)
    np.random.seed(1)
    N = 200  # points per class
    D = 2   # dimensions
    K = 2   # classes

    X = np.zeros((N*K, D))
    Y = np.zeros((N*K, 1))
    for j in range(K):
        ix = range(N*j, N*(j+1))
        r = np.linspace(0.0, 1, N)
        t = np.linspace(j*4, (j+1)*4, N) + np.random.randn(N)*0.2
        X[ix] = np.c_[r*np.sin(t), r*np.cos(t)]
        Y[ix] = j

    X = X.T  # (2, 400)
    Y = Y.T  # (1, 400)

    # Create network: 2 inputs → 16 → 8 → 1 output
    nn = DeepNeuralNetwork([2, 16, 8, 1])

    # Verify backpropagation
    print("=== Gradient Check ===")
    nn.gradient_check(X[:, :5], Y[:, :5])

    # Train with Adam
    print("\n=== Training with Adam ===")
    costs = nn.train(X, Y, epochs=2000, lr=0.001,
                     optimizer='adam', batch_size=64,
                     lr_decay='cosine')

    print(f"\nFinal accuracy: {nn.accuracy(X, Y):.1f}%")

import numpy as np

class BatchNorm:
    """Batch Normalization layer — forward and backward"""

    def __init__(self, num_features, momentum=0.9, eps=1e-5):
        self.gamma = np.ones(num_features)       # Learnable scale
        self.beta = np.zeros(num_features)        # Learnable shift
        self.eps = eps
        self.momentum = momentum
        # Running stats for inference
        self.running_mean = np.zeros(num_features)
        self.running_var = np.ones(num_features)

    def forward(self, Z, training=True):
        """
        Forward pass of batch normalization.
        Z: (n, m) — n features, m examples
        """
        if training:
            # Step 1: Batch mean
            self.mu = np.mean(Z, axis=1, keepdims=True)

            # Step 2: Batch variance
            self.var = np.var(Z, axis=1, keepdims=True)

            # Step 3: Normalize
            self.Z_centered = Z - self.mu
            self.std_inv = 1.0 / np.sqrt(self.var + self.eps)
            self.Z_norm = self.Z_centered * self.std_inv

            # Step 4: Scale and shift
            out = self.gamma.reshape(-1, 1) * self.Z_norm + self.beta.reshape(-1, 1)

            # Update running statistics
            self.running_mean = (self.momentum * self.running_mean +
                                (1 - self.momentum) * self.mu.flatten())
            self.running_var = (self.momentum * self.running_var +
                               (1 - self.momentum) * self.var.flatten())
        else:
            # Inference: use running statistics
            Z_norm = (Z - self.running_mean.reshape(-1, 1)) / \
                     np.sqrt(self.running_var.reshape(-1, 1) + self.eps)
            out = self.gamma.reshape(-1, 1) * Z_norm + self.beta.reshape(-1, 1)

        return out

    def backward(self, dout):
        """
        Backward pass of batch normalization.
        dout: gradient from upstream, same shape as Z
        """
        m = dout.shape[1]
        gamma = self.gamma.reshape(-1, 1)

        # Gradient of gamma and beta
        self.dgamma = np.sum(dout * self.Z_norm, axis=1)
        self.dbeta = np.sum(dout, axis=1)

        # Gradient through normalization
        dZ_norm = dout * gamma
        dvar = np.sum(dZ_norm * self.Z_centered * (-0.5) *
                      (self.var + self.eps)**(-1.5), axis=1, keepdims=True)
        dmu = (np.sum(dZ_norm * (-self.std_inv), axis=1, keepdims=True) +
               dvar * np.mean(-2 * self.Z_centered, axis=1, keepdims=True))
        dZ = dZ_norm * self.std_inv + dvar * 2 * self.Z_centered / m + dmu / m

        return dZ


# Demo
bn = BatchNorm(num_features=3)
Z = np.array([[2.0, 4.0, 6.0, 8.0],
              [1.0, 3.0, 5.0, 7.0],
              [0.5, 1.5, 2.5, 3.5]])

out = bn.forward(Z, training=True)
print("BN output:\n", out)
print("Mean (should be ~0):", np.mean(out, axis=1))
print("Var (should be ~1):", np.var(out, axis=1))

import tensorflow as tf
import numpy as np

# ==========================================
# CUSTOM TRAINING LOOP WITH GRADIENTTAPE
# ==========================================

# 1. Build model
model = tf.keras.Sequential([
    tf.keras.layers.Dense(64, activation='relu', input_shape=(10,),
                          kernel_initializer='he_normal'),
    tf.keras.layers.BatchNormalization(),
    tf.keras.layers.Dropout(0.3),
    tf.keras.layers.Dense(32, activation='relu',
                          kernel_initializer='he_normal'),
    tf.keras.layers.BatchNormalization(),
    tf.keras.layers.Dropout(0.3),
    tf.keras.layers.Dense(1, activation='sigmoid')
])

# 2. Setup optimizer with learning rate schedule
lr_schedule = tf.keras.optimizers.schedules.CosineDecay(
    initial_learning_rate=0.001,
    decay_steps=1000,
    alpha=0.0001  # minimum learning rate
)
optimizer = tf.keras.optimizers.Adam(
    learning_rate=lr_schedule,
    beta_1=0.9,
    beta_2=0.999,
    epsilon=1e-8
)

loss_fn = tf.keras.losses.BinaryCrossentropy()

# 3. Generate synthetic data
np.random.seed(42)
X_train = np.random.randn(1000, 10).astype(np.float32)
Y_train = (np.sum(X_train[:, :3], axis=1) > 0).astype(np.float32).reshape(-1, 1)

train_dataset = tf.data.Dataset.from_tensor_slices((X_train, Y_train))
train_dataset = train_dataset.shuffle(1000).batch(32)

# 4. Custom training step
@tf.function
def train_step(x_batch, y_batch):
    """One training step with explicit gradient computation"""
    with tf.GradientTape() as tape:
        # Forward pass
        predictions = model(x_batch, training=True)
        loss = loss_fn(y_batch, predictions)

        # Add L2 regularization manually
        l2_loss = tf.add_n([tf.nn.l2_loss(v)
                           for v in model.trainable_variables
                           if 'kernel' in v.name])
        total_loss = loss + 1e-4 * l2_loss

    # Backward pass — THIS IS BACKPROPAGATION!
    gradients = tape.gradient(total_loss, model.trainable_variables)

    # Gradient clipping (prevent exploding gradients)
    gradients, global_norm = tf.clip_by_global_norm(gradients, 1.0)

    # Update weights (Adam optimizer handles momentum + RMSprop internally)
    optimizer.apply_gradients(zip(gradients, model.trainable_variables))

    return total_loss, predictions

# 5. Training loop
for epoch in range(50):
    epoch_loss = 0
    num_batches = 0

    for x_batch, y_batch in train_dataset:
        loss, preds = train_step(x_batch, y_batch)
        epoch_loss += loss.numpy()
        num_batches += 1

    if epoch % 10 == 0:
        avg_loss = epoch_loss / num_batches
        current_lr = optimizer.learning_rate(optimizer.iterations).numpy()
        print(f"Epoch {epoch:3d} | Loss: {avg_loss:.4f} | LR: {current_lr:.6f}")

# 6. Inspect gradients for diagnosis
print("\n=== Gradient Analysis ===")
with tf.GradientTape() as tape:
    preds = model(X_train[:32], training=False)
    loss = loss_fn(Y_train[:32], preds)

grads = tape.gradient(loss, model.trainable_variables)
for var, grad in zip(model.trainable_variables, grads):
    if grad is not None:
        print(f"{var.name:30s} | grad mean: {tf.reduce_mean(grad):.6f} | "
              f"grad std: {tf.math.reduce_std(grad):.6f}")

import tensorflow as tf
import numpy as np

def create_model():
    return tf.keras.Sequential([
        tf.keras.layers.Dense(32, activation='relu', input_shape=(10,)),
        tf.keras.layers.Dense(16, activation='relu'),
        tf.keras.layers.Dense(1, activation='sigmoid')
    ])

optimizers = {
    'SGD':          tf.keras.optimizers.SGD(learning_rate=0.01),
    'SGD+Momentum': tf.keras.optimizers.SGD(learning_rate=0.01, momentum=0.9),
    'RMSprop':      tf.keras.optimizers.RMSprop(learning_rate=0.001),
    'Adam':         tf.keras.optimizers.Adam(learning_rate=0.001),
}

results = {}
for name, opt in optimizers.items():
    model = create_model()
    model.compile(optimizer=opt, loss='binary_crossentropy', metrics=['accuracy'])
    history = model.fit(X_train, Y_train, epochs=50, batch_size=32, verbose=0)
    results[name] = history.history['loss']
    print(f"{name:15s} | Final loss: {history.history['loss'][-1]:.4f} | "
          f"Final acc: {history.history['accuracy'][-1]:.4f}")

from sklearn.neural_network import MLPClassifier
from sklearn.datasets import make_moons
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
import numpy as np

# Generate dataset
X, y = make_moons(n_samples=1000, noise=0.2, random_state=42)
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)

scaler = StandardScaler()
X_train = scaler.fit_transform(X_train)
X_test = scaler.transform(X_test)

# Compare different solvers (optimizers) in scikit-learn
solvers = {
    'sgd': {'solver': 'sgd', 'learning_rate_init': 0.01,
            'momentum': 0.9, 'nesterovs_momentum': True},
    'adam': {'solver': 'adam', 'learning_rate_init': 0.001,
            'beta_1': 0.9, 'beta_2': 0.999},
    'lbfgs': {'solver': 'lbfgs'},  # Quasi-Newton method
}

print("Scikit-Learn MLPClassifier — Optimizer Comparison")
print("=" * 55)

for name, params in solvers.items():
    mlp = MLPClassifier(
        hidden_layer_sizes=(64, 32),  # Two hidden layers
        activation='relu',
        max_iter=500,
        random_state=42,
        alpha=0.0001,       # L2 regularization
        batch_size=32,
        **params
    )
    mlp.fit(X_train, y_train)
    train_acc = mlp.score(X_train, y_train) * 100
    test_acc = mlp.score(X_test, y_test) * 100
    print(f"{name:8s} | Train: {train_acc:.1f}% | Test: {test_acc:.1f}% | "
          f"Iters: {mlp.n_iter_}")

# Access the loss curve
print(f"\nAdam loss curve (first 10 epochs): {mlp.loss_curve_[:10]}")
print(f"Final loss: {mlp.loss_curve_[-1]:.6f}")