import numpy as np

# ── Illustrating the cost difference ──
def finite_difference_cost(n_params, forward_time_ms=1.0):
    """Naive gradient computation: 2 forward passes per parameter."""
    return 2 * n_params * forward_time_ms / 1000 / 3600 / 24  # days

def backprop_cost(forward_time_ms=1.0):
    """Backprop: roughly 2-3x a single forward pass, regardless of n_params."""
    return 2 * forward_time_ms / 1000  # seconds

# Small network: 100,000 parameters
n_small = 100_000
print(f"Small network ({n_small:,} params):")
print(f"  Finite differences: {finite_difference_cost(n_small):.2f} days")
print(f"  Backprop:           {backprop_cost():.3f} seconds")

# Large network: 100 million parameters (e.g., ResNet-50)
n_large = 100_000_000
print(f"\nLarge network ({n_large:,} params):")
print(f"  Finite differences: {finite_difference_cost(n_large):.0f} days")
print(f"  Backprop:           {backprop_cost():.3f} seconds")
print(f"\nBackprop speedup: {finite_difference_cost(n_large) * 86400 / backprop_cost():.2e}x faster!")

import numpy as np

# ── Chain rule example 1: y = sin(x²) ──
# Let u = x², so y = sin(u)
# dy/du = cos(u) = cos(x²)
# du/dx = 2x
# By chain rule: dy/dx = cos(x²) · 2x

x_val = 2.0
dy_dx_analytical = np.cos(x_val**2) * 2 * x_val
dy_dx_finite_diff = (np.sin((x_val+1e-5)**2) - np.sin((x_val-1e-5)**2)) / (2e-5)
print(f"dy/dx = cos(x²)·2x at x=2:")
print(f"  Analytical:       {dy_dx_analytical:.8f}")
print(f"  Finite difference: {dy_dx_finite_diff:.8f}")
print()

# ── Chain rule example 2: Multi-step chain ──
# y = loss(output(hidden(input)))
# loss = (ŷ - y_true)²  → dloss/dŷ = 2(ŷ - y_true)
# output = sigmoid(a)   → doutput/da = sigmoid(a) * (1 - sigmoid(a))
# hidden = w * x + b    → dhidden/dw = x

# Example numbers
x = 3.0
w = 0.5
b = 0.1
y_true = 1.0

def sigmoid(z): return 1 / (1 + np.exp(-z))

# Forward pass (computing and caching intermediate values)
a = w * x + b                  # pre-activation
yhat = sigmoid(a)              # output
loss = (yhat - y_true) ** 2    # loss

print(f"Forward pass:")
print(f"  a = w*x + b = {a:.4f}")
print(f"  ŷ = sigmoid(a) = {yhat:.4f}")
print(f"  loss = (ŷ - y_true)² = {loss:.4f}")

# Backward pass (chain rule applied step by step)
dloss_dyhat = 2 * (yhat - y_true)                     # ∂loss/∂ŷ
dyhat_da    = sigmoid(a) * (1 - sigmoid(a))           # ∂ŷ/∂a (sigmoid derivative)
da_dw       = x                                        # ∂a/∂w
da_db       = 1.0                                      # ∂a/∂b

# Chain rule: ∂loss/∂w = ∂loss/∂ŷ × ∂ŷ/∂a × ∂a/∂w
dloss_dw = dloss_dyhat * dyhat_da * da_dw
dloss_db = dloss_dyhat * dyhat_da * da_db

print(f"\nBackward pass:")
print(f"  ∂loss/∂ŷ = {dloss_dyhat:.4f}")
print(f"  ∂ŷ/∂a    = {dyhat_da:.4f}  (sigmoid derivative)")
print(f"  ∂a/∂w    = {da_dw:.4f}    (= x, the input)")
print(f"  ∂loss/∂w = {dloss_dw:.4f}  (full chain)")
print(f"  ∂loss/∂b = {dloss_db:.4f}  (full chain)")

import numpy as np

# ── Network architecture ──
# Input:  x  (2 inputs)  — shape (2,)
# Layer 1: W1 (2×2), b1 (2,) → pre-activation a1 = W1@x + b1 → activation z1 = relu(a1)
# Layer 2: W2 (1×2), b2 (1,) → pre-activation a2 = W2@z1 + b2 → output ŷ = a2 (linear)
# Loss: MSE = (ŷ - y_true)²

def relu(x):
    return np.maximum(0, x)

def relu_prime(x):
    """Derivative of ReLU: 1 if x > 0, else 0"""
    return (x > 0).astype(float)

# ── Small, specific example ──
np.random.seed(42)

# Network parameters (small random values)
W1 = np.array([[ 0.5, -0.3],
               [ 0.2,  0.8]])    # shape (2, 2)
b1 = np.array([0.1, -0.1])      # shape (2,)

W2 = np.array([[0.6, -0.4]])    # shape (1, 2)
b2 = np.array([0.0])            # shape (1,)

# Input and target
x = np.array([1.0, 2.0])        # input: shape (2,)
y_true = np.array([1.5])        # target: shape (1,)

# ── FORWARD PASS — compute and cache EVERYTHING ──
print("=" * 50)
print("FORWARD PASS")
print("=" * 50)

# Layer 1
a1 = W1 @ x + b1                # pre-activation: shape (2,)
z1 = relu(a1)                   # post-activation (hidden representation): shape (2,)

print(f"Layer 1 pre-activation a1 = W1@x + b1:")
print(f"  a1 = {a1}")
print(f"Layer 1 activation z1 = relu(a1):")
print(f"  z1 = {z1}")

# Layer 2
a2 = W2 @ z1 + b2              # output pre-activation: shape (1,)
yhat = a2                       # for regression, output = a2 directly

print(f"\nLayer 2 pre-activation a2 = W2@z1 + b2:")
print(f"  a2 = {a2}")
print(f"  ŷ  = {yhat}")

# Loss
loss = 0.5 * (yhat - y_true) ** 2   # 0.5 for clean gradient
print(f"\nMSE loss (×0.5) = (ŷ - y_true)² / 2 = {loss[0]:.6f}")

# Cache for backward pass
cache = {
    'x':  x,
    'a1': a1,
    'z1': z1,
    'a2': a2,
    'W1': W1, 'b1': b1,
    'W2': W2, 'b2': b2,
}
print(f"\nCached for backward: {list(cache.keys())}")

import numpy as np

# Continuing from Section 3 — all variables still defined

# ── BACKWARD PASS — derive all gradients using chain rule ──
print("=" * 50)
print("BACKWARD PASS")
print("=" * 50)

# ── Step 1: Gradient of loss w.r.t. output ŷ ──
# L = 0.5 * (ŷ - y_true)²
# ∂L/∂ŷ = (ŷ - y_true)
dL_dyhat = yhat - y_true                    # shape (1,)
print(f"∂L/∂ŷ = ŷ - y_true = {dL_dyhat}")

# ── Step 2: Gradients for Layer 2 (W2, b2) ──
# ŷ = W2 @ z1 + b2
# ∂ŷ/∂W2 = z1 (the input to this layer)
# ∂ŷ/∂b2 = 1
# ∂L/∂W2 = ∂L/∂ŷ · ∂ŷ/∂W2 = dL_dyhat · z1.T (outer product for matrix)
# ∂L/∂b2 = ∂L/∂ŷ

dL_dW2 = dL_dyhat.reshape(-1, 1) @ z1.reshape(1, -1)   # shape (1, 2)
dL_db2 = dL_dyhat                                        # shape (1,)
print(f"\n∂L/∂W2 = {dL_dW2}")
print(f"∂L/∂b2 = {dL_db2}")

# ── Step 3: Gradient flowing back into hidden layer (∂L/∂z1) ──
# ŷ = W2 @ z1 + b2
# ∂ŷ/∂z1 = W2.T
# ∂L/∂z1 = W2.T @ ∂L/∂ŷ  (gradient backpropagates through W2)
dL_dz1 = W2.T @ dL_dyhat                   # shape (2,)
print(f"\n∂L/∂z1 = W2.T @ ∂L/∂ŷ = {dL_dz1}")

# ── Step 4: Through ReLU activation (∂L/∂a1) ──
# z1 = relu(a1)
# ∂z1/∂a1 = relu'(a1) = 1 if a1 > 0, else 0
# ∂L/∂a1 = ∂L/∂z1 * relu'(a1)  (element-wise, not matrix multiply)
dL_da1 = dL_dz1 * relu_prime(a1)           # shape (2,)
print(f"\nrelu'(a1) = {relu_prime(a1)}")     # [0, 1] — first neuron was dead!
print(f"∂L/∂a1 = ∂L/∂z1 * relu'(a1) = {dL_da1}")
# Notice: gradient for dead neuron (a1[0] < 0) is 0 — it gets no update!

# ── Step 5: Gradients for Layer 1 (W1, b1) ──
# a1 = W1 @ x + b1
# ∂a1/∂W1 = x
# ∂L/∂W1 = ∂L/∂a1 · x.T (outer product)
# ∂L/∂b1 = ∂L/∂a1
dL_dW1 = dL_da1.reshape(-1, 1) @ x.reshape(1, -1)      # shape (2, 2)
dL_db1 = dL_da1                                          # shape (2,)
print(f"\n∂L/∂W1 = {dL_dW1}")
print(f"∂L/∂b1 = {dL_db1}")

print(f"\n{'='*50}")
print(f"Summary of gradients:")
print(f"  dL/dW1 shape: {dL_dW1.shape}")
print(f"  dL/db1 shape: {dL_db1.shape}")
print(f"  dL/dW2 shape: {dL_dW2.shape}")
print(f"  dL/db2 shape: {dL_db2.shape}")

import torch
import numpy as np

# Use same weights and data as Section 3/4
W1_np = np.array([[ 0.5, -0.3],
                   [ 0.2,  0.8]], dtype=np.float64)
b1_np = np.array([0.1, -0.1], dtype=np.float64)
W2_np = np.array([[0.6, -0.4]], dtype=np.float64)
b2_np = np.array([0.0], dtype=np.float64)
x_np  = np.array([1.0, 2.0], dtype=np.float64)
y_np  = np.array([1.5], dtype=np.float64)

# Convert to PyTorch tensors WITH requires_grad
W1 = torch.tensor(W1_np, requires_grad=True)
b1 = torch.tensor(b1_np, requires_grad=True)
W2 = torch.tensor(W2_np, requires_grad=True)
b2 = torch.tensor(b2_np, requires_grad=True)
x  = torch.tensor(x_np)
y_true = torch.tensor(y_np)

# Forward pass
a1   = W1 @ x + b1
z1   = torch.relu(a1)
a2   = W2 @ z1 + b2
yhat = a2
loss = 0.5 * ((yhat - y_true) ** 2).sum()

# Backward pass
loss.backward()

# Compare gradients
print("Gradient Comparison: Manual NumPy vs PyTorch Autograd")
print("=" * 55)

print(f"\ndL/dW2:")
print(f"  Manual:  {dL_dW2.round(6)}")
print(f"  PyTorch: {W2.grad.numpy().round(6)}")
match_W2 = np.allclose(dL_dW2, W2.grad.numpy(), atol=1e-10)
print(f"  Match: {'✓ YES' if match_W2 else '✗ NO'}")

print(f"\ndL/db2:")
print(f"  Manual:  {dL_db2.round(6)}")
print(f"  PyTorch: {b2.grad.numpy().round(6)}")
match_b2 = np.allclose(dL_db2, b2.grad.numpy(), atol=1e-10)
print(f"  Match: {'✓ YES' if match_b2 else '✗ NO'}")

print(f"\ndL/dW1:")
print(f"  Manual:\n{dL_dW1.round(6)}")
print(f"  PyTorch:\n{W1.grad.numpy().round(6)}")
match_W1 = np.allclose(dL_dW1, W1.grad.numpy(), atol=1e-10)
print(f"  Match: {'✓ YES' if match_W1 else '✗ NO'}")

print(f"\ndL/db1:")
print(f"  Manual:  {dL_db1.round(6)}")
print(f"  PyTorch: {b1.grad.numpy().round(6)}")
match_b1 = np.allclose(dL_db1, b1.grad.numpy(), atol=1e-10)
print(f"  Match: {'✓ YES' if match_b1 else '✗ NO'}")

all_match = match_W1 and match_b1 and match_W2 and match_b2
print(f"\nAll gradients verified: {'✓ PASS' if all_match else '✗ FAIL'}")

import numpy as np

def sigmoid_derivative_at_x(x):
    s = 1 / (1 + np.exp(-x))
    return s * (1 - s)  # max = 0.25

def relu_derivative_at_x(x):
    return 1.0 if x > 0 else 0.0

# Simulate gradient magnitude after n layers
def gradient_magnitude_after_layers(n_layers, activation='sigmoid', x_val=0.0):
    """How large is the gradient after passing through n layers?"""
    gradient = 1.0
    for _ in range(n_layers):
        if activation == 'sigmoid':
            gradient *= sigmoid_derivative_at_x(x_val)    # multiply by sigmoid'
        elif activation == 'relu':
            gradient *= relu_derivative_at_x(abs(x_val))  # 1 for + inputs
    return gradient

depths = [1, 5, 10, 20, 50]
print(f"{'Depth':>8} | {'Sigmoid grad':>15} | {'ReLU grad':>12}")
print("-" * 42)
for d in depths:
    sig_grad  = gradient_magnitude_after_layers(d, 'sigmoid', x_val=0.0)  # worst case (max)
    relu_grad = gradient_magnitude_after_layers(d, 'relu', x_val=1.0)     # positive inputs
    print(f"{d:>8} | {sig_grad:>15.2e} | {relu_grad:>12.1f}")

# The difference is catastrophic for deep networks

import numpy as np

# ── Illustrating residual connection gradient flow ──
def residual_block_gradient(x_val, W_scale=0.9):
    """
    Standard block: y = relu(W@x)  — gradient must pass through relu and W
    Residual block: y = relu(W@x) + x  — gradient has a direct path through +x
    """
    # Forward
    Wx = W_scale * x_val
    activated = max(0, Wx)     # relu
    y_standard = activated          # standard block output
    y_residual = activated + x_val  # residual block output

    # Gradients
    # Standard: dL/dx = dL/dy * dy/dx = dL/dy * W_scale * relu'(Wx)
    relu_prime = 1.0 if Wx > 0 else 0.0
    grad_standard = W_scale * relu_prime   # can be 0 if Wx < 0

    # Residual: dL/dx = dL/dy * (W_scale * relu'(Wx) + 1)
    # The "+1" comes from the skip connection derivative
    grad_residual = W_scale * relu_prime + 1.0   # always ≥ 1.0!

    return grad_standard, grad_residual

# Compare over 20 layers starting from gradient=1
layers = 20
grad_standard = 1.0
grad_residual = 1.0
print(f"{'Layer':>7} | {'Standard grad':>15} | {'Residual grad':>15}")
print("-" * 45)
for i in range(1, layers+1):
    gs, gr = residual_block_gradient(0.5)  # positive input
    grad_standard *= gs
    grad_residual *= gr
    if i in [1, 5, 10, 15, 20]:
        print(f"{i:>7} | {grad_standard:>15.4f} | {grad_residual:>15.1f}")

print("\nResidual gradients stay healthy through all 20 layers!")

import torch
import torch.nn as nn

# ── Exploring the computation graph manually ──
x = torch.tensor([1.0, 2.0], requires_grad=True)
W = torch.tensor([[0.5, -0.3], [0.2, 0.8]], requires_grad=True)
b = torch.tensor([0.1, -0.1], requires_grad=True)

# Layer 1 operations
a = W @ x + b          # linear transformation
z = torch.relu(a)      # activation

# Look at the graph
print("Computation graph nodes (grad_fn attributes):")
print(f"  x.grad_fn:     {x.grad_fn}")        # None — x is a leaf
print(f"  W.grad_fn:     {W.grad_fn}")        # None — W is a leaf
print(f"  a.grad_fn:     {a.grad_fn}")        # AddmmBackward (linear)
print(f"  z.grad_fn:     {z.grad_fn}")        # ReluBackward (relu)

# Each grad_fn knows its inputs
print(f"\n  a.grad_fn.next_functions:")
for fn, _ in a.grad_fn.next_functions:
    print(f"    {fn}")   # MmBackward (matmul) and AccumulateGrad (for b)

# The full graph: z → ReluBack → AddBack → MmBack → AccumulateGrad(x, W, b)

# ── Torchviz for visualization (optional, requires pip install torchviz) ──
# from torchviz import make_dot
# loss = (z.sum() - 1.0) ** 2
# dot = make_dot(loss, params={'W': W, 'b': b, 'x': x})
# dot.render("computation_graph", format="png")

# ── Counting parameters in a network ──
model = nn.Sequential(
    nn.Linear(784, 256),
    nn.ReLU(),
    nn.Linear(256, 128),
    nn.ReLU(),
    nn.Linear(128, 10)
)

total_params = 0
for name, param in model.named_parameters():
    print(f"  {name}: {param.shape}, {param.numel():,} params")
    total_params += param.numel()
print(f"\nTotal parameters: {total_params:,}")

import torch
import torch.nn as nn

# ── Bug 1: Forgot zero_grad — gradient accumulation ──
model = nn.Linear(3, 1)
x = torch.randn(5, 3)
y = torch.randn(5, 1)
loss_fn = nn.MSELoss()
optimizer = torch.optim.SGD(model.parameters(), lr=0.01)

print("Bug 1: Gradient accumulation (forgot zero_grad)")
for step in range(3):
    # BAD: no optimizer.zero_grad()
    loss = loss_fn(model(x), y)
    loss.backward()
    print(f"  Step {step+1} | grad norm: {model.weight.grad.norm().item():.4f}  ← grows!")

# Reset
optimizer.zero_grad()
print("\nCorrect: with zero_grad()")
for step in range(3):
    optimizer.zero_grad()   # ← CRITICAL
    loss = loss_fn(model(x), y)
    loss.backward()
    print(f"  Step {step+1} | grad norm: {model.weight.grad.norm().item():.4f}  ← stable")

# ── Bug 2: Calling backward on non-scalar ──
print("\nBug 2: backward on vector output")
x = torch.randn(3, requires_grad=True)
y = x ** 2   # vector output (shape (3,))
try:
    y.backward()  # RuntimeError!
except RuntimeError as e:
    print(f"  Error: {e}")

# Fix: reduce to scalar first
y.sum().backward()   # or y.mean().backward()
print(f"  Fix: y.sum().backward() → x.grad = {x.grad}")

# ── Bug 3: In-place operations breaking autograd ──
print("\nBug 3: In-place operations")
x = torch.tensor([1.0, 2.0], requires_grad=True)
z = x * 2
# BAD: in-place modification of a tensor in the graph
try:
    z += 1   # RuntimeError: in-place op on leaf that requires grad
    z.backward(torch.ones(2))
except RuntimeError as e:
    print(f"  Error: {e}")

# Fix: use out-of-place operations
x = torch.tensor([1.0, 2.0], requires_grad=True)
z = x * 2
z_new = z + 1    # ← out-of-place, creates new tensor
z_new.sum().backward()
print(f"  Fix: out-of-place op → x.grad = {x.grad}")

# ── Bug 4: Detached tensor doesn't participate in gradient ──
print("\nBug 4: Accidentally detaching from graph")
x = torch.tensor(2.0, requires_grad=True)
y = x * 3
y_detached = y.detach()   # y_detached has the value but not the graph
z = y_detached * 2        # this operation doesn't track x!
z.backward()
print(f"  x.grad after detach: {x.grad}")   # None — gradient can't flow through detach

# ── Gradient checking: verify your manual backprop ──
print("\nGradient checking (finite differences verification)")
def manual_loss(x_np):
    """A function we want to differentiate."""
    return float(np.sin(x_np[0]) * np.exp(-x_np[1]**2))

import numpy as np
x_check = np.array([1.0, 0.5])
eps = 1e-5

# Numerical gradient (finite differences)
grad_numerical = np.zeros(2)
for i in range(2):
    x_plus = x_check.copy(); x_plus[i] += eps
    x_minus = x_check.copy(); x_minus[i] -= eps
    grad_numerical[i] = (manual_loss(x_plus) - manual_loss(x_minus)) / (2 * eps)

# Analytical gradient: d/dx0 = cos(x0)*exp(-x1²), d/dx1 = sin(x0)*(-2x1)*exp(-x1²)
grad_analytical = np.array([
    np.cos(x_check[0]) * np.exp(-x_check[1]**2),
    np.sin(x_check[0]) * (-2*x_check[1]) * np.exp(-x_check[1]**2)
])

print(f"  Numerical:  {grad_numerical.round(8)}")
print(f"  Analytical: {grad_analytical.round(8)}")
print(f"  Difference: {np.abs(grad_numerical - grad_analytical).max():.2e}")