Example: A house with 3 features
  x = [size_in_sqft, number_of_bedrooms, age_in_years]
  x = [1500, 3, 10]

We write: x ∈ ℝ³  (x is in 3-dimensional real space)

Addition:    [1, 2] + [3, 4] = [4, 6]      (element-wise)
Scalar mult: 2 × [1, 2, 3]  = [2, 4, 6]     (multiply each element)
Dot product: [1, 2] · [3, 4] = 1×3 + 2×4 = 11   (scalar result)

import numpy as np

# Vectors in code
x = np.array([1500, 3, 10])  # a house
y = np.array([2000, 4, 5])   # another house

# Dot product = how similar are they?
similarity = np.dot(x, y)  # 1500×2000 + 3×4 + 10×5 = 3,000,062

Matrix A (2×3):    Matrix B (3×2):    A × B (2×2):
[1, 2, 3]          [7, 8]            [1×7+2×9+3×11,   1×8+2×10+3×12]
[4, 5, 6]          [9, 10]           [4×7+5×9+6×11,   4×8+5×10+6×12]
                   [11, 12]          
                                     = [58, 64]
                                       [139, 154]

import numpy as np

# Matrix multiplication
A = np.array([[1, 2, 3], [4, 5, 6]])
B = np.array([[7, 8], [9, 10], [11, 12]])
C = A @ B  # or np.matmul(A, B)

print(C)  # [[58, 64], [139, 154]]

# A neural network layer: h = X @ W + b
X = np.random.randn(32, 784)    # batch of 32 images, 784 pixels each
W = np.random.randn(784, 256)   # weight matrix
b = np.random.randn(256)        # bias

h = X @ W + b  # shape: (32, 256)
# This is called a "linear layer" or "fully-connected layer"

Property            Definition                         ML Use
─────────           ──────────                         ──────
Transpose (Aᵀ)      Swap rows and columns              Backpropagation, attention
Inverse (A⁻¹)       A × A⁻¹ = I                       Solving normal equations
Identity (I)        1 on diagonal, 0 elsewhere          The "do nothing" matrix
Trace (tr[A])       Sum of diagonal elements            Regularization, optimization
Determinant (|A|)   "Volume" of transformation          Normalizing flows, change of var.
Eigenvalues (λ)     Av = λv                            PCA, graph theory, stability
Singular Values     Diagonal of Σ in A = UΣVᵀ         SVD, compression, PCA

# SVD in one line
U, S, Vt = np.linalg.svd(A, full_matrices=False)

# Low-rank approximation: keep k largest singular values
k = 50
A_approx = U[:, :k] @ np.diag(S[:k]) @ Vt[:k, :]
# This is the foundation of PCA, recommendation systems, and compression

Scalar:   shape []        or ()
Vector:   shape [3]       or (3,)
Matrix:   shape [3, 4]    or (3, 4)
3-tensor: shape [3, 4, 5] or (3, 4, 5)

# In ML, everything is tensors:
batch_images = torch.randn(32, 3, 256, 256) # batch, channels, height, width
batch_text   = torch.randn(32, 128, 768)     # batch, seq_len, embedding_dim

L1 norm:   ||x||₁ = |x₁| + |x₂| + ... + |xₙ|    (sparsity)
L2 norm:   ||x||₂ = √(x₁² + x₂² + ... + xₙ²)     (Euclidean distance)
L∞ norm:   ||x||∞ = max(|x₁|, |x₂|, ..., |xₙ|)   (worst-case)
Frobenius: ||A||_F = √(ΣΣ Aᵢⱼ²)                   (matrix L2)

x = np.array([3, -4])
l1 = np.sum(np.abs(x))           # 7
l2 = np.sqrt(np.sum(x**2))       # 5
l2 = np.linalg.norm(x)           # same

Geometric: slope of the tangent line at point x
Intuitive: if I nudge x up a tiny bit, does f go up or down?

# Numerical derivative (approximate)
def derivative(f, x, h=1e-7):
    return (f(x + h) - f(x - h)) / (2 * h)

# Example: f(x) = x²
f = lambda x: x**2
print(derivative(f, 3))  # ≈ 6.0 (since derivative of x² is 2x)

Function       Derivative          Used In
────────       ──────────          ───────
xⁿ             n·xⁿ⁻¹              Polynomials
eˣ             eˣ                  Softmax, sigmoid
ln(x)          1/x                 Cross-entropy loss
sin(x)         cos(x)              Positional encoding
σ(x)=1/(1+e⁻ˣ) σ(x)(1-σ(x))       Sigmoid activation
tanh(x)        1-tanh²(x)          Tanh activation
ReLU(x)=max(0,x) 0 if x<0, 1 if x>0  ReLU activation

# f(w, b) = (w*x + b - y)²  (squared error for one data point)
def f(w, b, x=2, y=5):
    return (w*x + b - y)**2

# Partial derivative w.r.t w: 2*(w*x + b - y)*x
# Partial derivative w.r.t b: 2*(w*x + b - y)*1

# Gradient descent visualized
θ = random_initial_weights()
for step in range(steps):
    gradient = compute_gradient(loss_function, θ, data)
    θ = θ - learning_rate * gradient

P(rain tomorrow) = 0.3    means 30% chance
P(sunny) = 0.7            means 70% chance
                        total must = 1

Posterior = Likelihood × Prior / Evidence

θ = model parameters, D = data
P(θ|D): what we believe about parameters AFTER seeing data
P(D|θ): how likely the data is given parameters (likelihood)
P(θ): what we believed BEFORE seeing data (prior)
P(D): probability of the data (normalization)

import numpy as np
import matplotlib.pyplot as plt

# Generate samples from a Gaussian
mu, sigma = 0, 1
samples = np.random.normal(mu, sigma, 10000)

# Estimate parameters from data
mu_hat = np.mean(samples)     # ≈ 0
sigma_hat = np.std(samples)   # ≈ 1

# Why Gaussian? Central Limit Theorem:
# The sum of many independent random variables
# approaches a Gaussian, regardless of their individual distributions

# MSE = negative log-likelihood of Gaussian
# For y ~ N(ŷ, σ²):
#   log P(y|ŷ) = -½(y-ŷ)²/σ² - ½log(2πσ²)
# Minimizing MSE = maximizing log likelihood
# (assuming fixed σ²)

loss = np.mean((y_pred - y_true) ** 2)  # this IS MLE under Gaussian noise

Example: A fair coin
  P(heads) = 0.5, P(tails) = 0.5
  H = -0.5·log₂(0.5) - 0.5·log₂(0.5)
  H = -0.5·(-1) - 0.5·(-1) = 0.5 + 0.5 = 1 bit

Example: A always-heads coin
  P(heads) = 1, P(tails) = 0
  H = -1·log₂(1) - 0·log₂(0) = 0 bits (no uncertainty)

# Cross-entropy loss for classification
# True label: cat (index 2 in 10 classes)
# Model predicts: [0.01, 0.02, 0.93, 0.01, ..., 0.01]

true_class = 2
predicted_probs = [0.01, 0.02, 0.93, 0.01, ...]

loss = -np.log(predicted_probs[true_class])  # -log(0.93) = 0.073

def gradient_descent(gradient_fn, init_params, lr=0.01, steps=1000):
    """
    The simplest optimization algorithm.

    gradient_fn: function that returns ∇L at current parameters
    init_params: starting point θ₀
    lr: step size η (learning rate)
    steps: how many iterations
    """
    params = init_params.copy()
    history = [params.copy()]

    for i in range(steps):
        grad = gradient_fn(params)
        params = params - lr * grad
        history.append(params.copy())

    return params, history

def sgd_step(model, batch_X, batch_y, lr):
    """
    One step of SGD. 
    Key insight: gradient on mini-batch ≈ gradient on full dataset
    but ~1000x cheaper to compute.
    """
    predictions = model(batch_X)
    loss = compute_loss(predictions, batch_y)
    loss.backward()  # compute gradients

    with torch.no_grad():
        for param in model.parameters():
            param -= lr * param.grad
            param.grad.zero_()

A linear model with 2 features: h(x) = w₁x₁ + w₂x₂ + b
Hypothesis space: all possible (w₁, w₂, b) ∈ ℝ³

A decision tree of depth 3: all possible tree structures of depth ≤ 3
Hypothesis space: all decision trees fitting this constraint

import numpy as np

def linear_regression_closed_form(X, y):
    """
    Closed-form solution for linear regression.
    X: (n_samples, n_features)
    y: (n_samples,)
    """
    # Add bias term (column of 1s)
    X_with_bias = np.c_[np.ones(X.shape[0]), X]

    # Normal equation: w = (X^T X)^(-1) X^T y
    w = np.linalg.inv(X_with_bias.T @ X_with_bias) @ X_with_bias.T @ y

    return w[0], w[1:]  # bias, weights

def linear_regression_gradient(X, y, lr=0.01, epochs=1000):
    """
    Gradient descent solution for linear regression.
    More scalable than closed-form for large datasets.
    """
    n, d = X.shape
    w = np.zeros(d)
    b = 0.0

    for epoch in range(epochs):
        y_pred = X @ w + b
        error = y_pred - y

        # Gradients
        dw = (2/n) * X.T @ error
        db = (2/n) * np.sum(error)

        # Update
        w -= lr * dw
        b -= lr * db

        if epoch % 100 == 0:
            loss = np.mean(error ** 2)
            print(f"Epoch {epoch}: loss = {loss:.4f}")

    return w, b

def sigmoid(z):
    """Numerically stable sigmoid."""
    return 1 / (1 + np.exp(-np.clip(z, -100, 100)))

def logistic_regression_gradient(X, y, lr=0.01, epochs=1000):
    """
    Logistic regression via gradient descent.
    """
    n, d = X.shape
    w = np.zeros(d)
    b = 0.0

    for epoch in range(epochs):
        # Forward pass
        z = X @ w + b
        y_pred = sigmoid(z)

        # Binary cross-entropy loss
        loss = -np.mean(y * np.log(y_pred + 1e-8) + (1-y) * np.log(1 - y_pred + 1e-8))

        # Gradients
        error = y_pred - y
        dw = (1/n) * X.T @ error
        db = (1/n) * np.sum(error)

        # Update
        w -= lr * dw
        b -= lr * db

        if epoch % 200 == 0:
            print(f"Epoch {epoch}: loss = {loss:.4f}")

    return w, b

# Decision boundary: P(y=1|x) > 0.5
def predict(X, w, b):
    return (sigmoid(X @ w + b) >= 0.5).astype(int)

def softmax(z):
    """Numerically stable softmax."""
    shifted = z - np.max(z, axis=-1, keepdims=True)
    exp = np.exp(shifted)
    return exp / np.sum(exp, axis=-1, keepdims=True)

# Cross-entropy for K classes:
# L = -Σ y_k · log(ŷ_k)
# Where y is one-hot encoded (e.g., [0, 0, 1, 0])

# Common kernels
linear:     K(x, y) = x · y
polynomial: K(x, y) = (x · y + c)^d
RBF/Gaussian: K(x, y) = exp(-γ||x-y||²)
sigmoid:    K(x, y) = tanh(αx·y + c)

from sklearn.svm import SVC

# RBF SVM — the "I don't know what to use" champion
svm = SVC(kernel='rbf', C=1.0, gamma='scale')
svm.fit(X_train, y_train)

# C controls the margin-violation penalty
# Large C → hard margin (overfit potential)
# Small C → soft margin (more tolerant of errors)

class DecisionTreeStub:
    """A single-level decision tree (decision stump)."""
    def fit(self, X, y):
        self.best_feature = None
        self.best_threshold = None
        self.best_gini = 1.0

        n, d = X.shape
        for feature in range(d):
            thresholds = np.percentile(X[:, feature], np.linspace(10, 90, 20))
            for threshold in thresholds:
                left = y[X[:, feature] <= threshold]
                right = y[X[:, feature] > threshold]

                gini = self._weighted_gini(left, right)
                if gini < self.best_gini:
                    self.best_gini = gini
                    self.best_feature = feature
                    self.best_threshold = threshold

        self.left_pred = np.mean(y[X[:, self.best_feature] <= self.best_threshold])
        self.right_pred = np.mean(y[X[:, self.best_feature] > self.best_threshold])

    def predict(self, X):
        left_mask = X[:, self.best_feature] <= self.best_threshold
        preds = np.zeros(X.shape[0])
        preds[left_mask] = self.left_pred
        preds[~left_mask] = self.right_pred
        return preds

    def _weighted_gini(self, left, right):
        n_left, n_right = len(left), len(right)
        n_total = n_left + n_right
        gini_left = 1 - sum((np.mean(left == k))**2 for k in np.unique(left)) if n_left > 0 else 0
        gini_right = 1 - sum((np.mean(right == k))**2 for k in np.unique(right)) if n_right > 0 else 0
        return (n_left/n_total) * gini_left + (n_right/n_total) * gini_right

class RandomForest:
    """Simplified random forest for classification."""
    def __init__(self, n_trees=100, max_depth=5, max_features='sqrt'):
        self.n_trees = n_trees
        self.max_depth = max_depth
        self.max_features = max_features
        self.trees = []

    def fit(self, X, y):
        n, d = X.shape
        n_features = int(np.sqrt(d)) if self.max_features == 'sqrt' else d

        for _ in range(self.n_trees):
            # Bootstrap sample
            indices = np.random.choice(n, n, replace=True)
            X_boot, y_boot = X[indices], y[indices]

            # Random feature subset
            feature_subset = np.random.choice(d, n_features, replace=False)
            X_subset = X_boot[:, feature_subset]

            # Train tree
            tree = DecisionTree(max_depth=self.max_depth)
            tree.fit(X_subset, y_boot)
            self.trees.append((tree, feature_subset))

    def predict(self, X):
        predictions = np.zeros((X.shape[0], len(self.trees)))
        for i, (tree, features) in enumerate(self.trees):
            predictions[:, i] = tree.predict(X[:, features])
        # Majority vote
        return np.round(np.mean(predictions, axis=1))

# Simplified gradient boosting for regression
def gradient_boosting(X, y, n_estimators=100, lr=0.1, max_depth=3):
    """Train a gradient boosting ensemble."""
    # Start with mean prediction
    F = np.full(y.shape, np.mean(y))
    trees = []

    for m in range(n_estimators):
        # Compute pseudo-residuals (negative gradient)
        residuals = y - F  # For MSE: -dL/dF = y - F

        # Train tree on residuals
        tree = DecisionTreeRegressor(max_depth=max_depth)
        tree.fit(X, residuals)

        # Update ensemble
        F += lr * tree.predict(X)
        trees.append(tree)

    return trees

def predict_boosting(X, trees, lr=0.1):
    F = np.zeros(X.shape[0])
    for tree in trees:
        F += lr * tree.predict(X)
    return F

def kmeans(X, k, max_iters=100, tol=1e-4):
    """
    K-means clustering from scratch.

    Steps:
    1. Initialize k random centroids
    2. Assign each point to nearest centroid
    3. Update centroids to mean of assigned points
    4. Repeat 2-3 until convergence
    """
    n, d = X.shape

    # Step 1: Initialize centroids (k-means++ style)
    centroids = [X[np.random.choice(n)]]
    for _ in range(k - 1):
        distances = np.min([np.linalg.norm(X - c, axis=1) for c in centroids], axis=0)
        probs = distances / distances.sum()
        centroids.append(X[np.random.choice(n, p=probs)])
    centroids = np.array(centroids)

    for iteration in range(max_iters):
        # Step 2: Assign each point to nearest centroid
        distances = np.array([np.linalg.norm(X - c, axis=1) for c in centroids])
        labels = np.argmin(distances, axis=0)

        # Step 3: Update centroids
        new_centroids = np.array([X[labels == j].mean(axis=0) for j in range(k)])

        # Check convergence
        if np.all(np.linalg.norm(new_centroids - centroids, axis=1) < tol):
            break
        centroids = new_centroids

    return labels, centroids

inertias = []
for k in range(1, 11):
    labels, centroids = kmeans(X, k)
    inertia = sum(np.min([np.linalg.norm(X - c, axis=1) for c in centroids], axis=0)**2)
    inertias.append(inertia)

# Plot: look for the "elbow" where adding more clusters gives diminishing returns

def pca_from_scratch(X, n_components=2):
    """
    PCA via eigendecomposition of covariance matrix.
    """
    # Center the data
    X_centered = X - X.mean(axis=0)

    # Covariance matrix
    cov = X_centered.T @ X_centered / (X_centered.shape[0] - 1)

    # Eigendecomposition
    eigenvalues, eigenvectors = np.linalg.eigh(cov)

    # Sort by eigenvalue (descending)
    idx = np.argsort(eigenvalues)[::-1]
    eigenvalues = eigenvalues[idx]
    eigenvectors = eigenvectors[:, idx]

    # Select top n_components
    components = eigenvectors[:, :n_components]

    # Project data
    X_reduced = X_centered @ components

    # Explained variance ratio
    explained_variance_ratio = eigenvalues[:n_components] / eigenvalues.sum()

    return X_reduced, components, explained_variance_ratio

Example: Linear classifiers in 2D
- Can shatter 3 points? YES (any labeling of 3 points can be separated)
- Can shatter 4 points? NO (for any 4 points, at least one labeling is impossible)
- VC dimension = 3

Example: Decision trees of depth D
- VC dimension = O(2^D) (exponential in depth)

# ALWAYS split before anything else
from sklearn.model_selection import train_test_split

# Standard: 80/20 or 70/15/15
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

# For hyperparameter tuning: three-way split
X_train, X_temp, y_train, y_temp = train_test_split(X, y, test_size=0.3)
X_val, X_test, y_val, y_test = train_test_split(X_temp, y_temp, test_size=0.5)

from sklearn.model_selection import KFold

kf = KFold(n_splits=5, shuffle=True, random_state=42)
scores = []

for train_idx, val_idx in kf.split(X_train):
    X_fold_train, X_fold_val = X_train[train_idx], X_train[val_idx]
    y_fold_train, y_fold_val = y_train[train_idx], y_train[val_idx]

    model = train(X_fold_train, y_fold_train)
    score = model.score(X_fold_val, y_fold_val)
    scores.append(score)

print(f"Mean CV score: {np.mean(scores):.3f} ± {np.std(scores):.3f}")

              Predicted Positive    Predicted Negative
Actual Pos    TP (True Positive)    FN (False Negative)
Actual Neg    FP (False Positive)   TN (True Negative)

Accuracy  = (TP + TN) / (TP + TN + FP + FN)  [overall correctness]
Precision = TP / (TP + FP)                    [when we say yes, how often right?]
Recall    = TP / (TP + FN)                    [how many real positives do we catch?]
F1        = 2·P·R / (P + R)                  [harmonic mean of P and R]

def forward_neuron(x, w, b, activation='relu'):
    """
    A single neuron forward pass.
    x: input vector
    w: weight vector
    b: bias scalar
    """
    z = np.dot(x, w) + b
    if activation == 'relu':
        return np.maximum(0, z)
    elif activation == 'sigmoid':
        return 1 / (1 + np.exp(-np.clip(z, -100, 100)))
    elif activation == 'tanh':
        return np.tanh(z)
    elif activation == 'linear':
        return z

Layer 1: h = ReLU(x · W₁ + b₁)     x: (batch, D_in), W₁: (D_in, D_hidden)
Layer 2: ŷ = h · W₂ + b₂            W₂: (D_hidden, D_out)

z₁ = x · W₁ + b₁        [pre-activation, layer 1]
a₁ = σ(z₁)               [activation, layer 1]
z₂ = a₁ · W₂ + b₂       [pre-activation, layer 2]
a₂ = σ(z₂)               [activation, layer 2]
z₃ = a₂ · W₃ + b₃       [pre-activation, output layer]
ŷ = z₃                   [linear output for regression]
L = ½(ŷ - y)²            [MSE loss]

def backprop_full(X, y, W1, b1, W2, b2, W3, b3):
    """Full forward + backward pass for 3-layer network."""
    # Forward pass
    z1 = X @ W1 + b1
    a1 = np.maximum(0, z1)  # ReLU
    z2 = a1 @ W2 + b2
    a2 = np.maximum(0, z2)
    z3 = a2 @ W3 + b3
    y_pred = z3

    # Loss
    loss = 0.5 * np.mean((y_pred - y) ** 2)

    # Backward pass
    m = X.shape[0]  # batch size

    # Output layer
    dz3 = (y_pred - y) / m  # δ₃
    dW3 = a2.T @ dz3
    db3 = np.sum(dz3, axis=0)

    # Layer 2
    da2 = dz3 @ W3.T
    dz2 = da2 * (z2 > 0)  # ReLU derivative
    dW2 = a1.T @ dz2
    db2 = np.sum(dz2, axis=0)

    # Layer 1
    da1 = dz2 @ W2.T
    dz1 = da1 * (z1 > 0)  # ReLU derivative
    dW1 = X.T @ dz1
    db1 = np.sum(dz1, axis=0)

    return loss, (dW1, db1, dW2, db2, dW3, db3)

def conv2d_forward(image, kernel, stride=1, padding=0):
    """
    2D convolution from scratch.
    image: (H, W, C_in) — input image (height, width, channels)
    kernel: (K, K, C_in, C_out) — conv filters
    """
    H, W, C_in = image.shape
    K, _, _, C_out = kernel.shape

    # Apply padding
    if padding > 0:
        image = np.pad(image, ((padding, padding), (padding, padding), (0, 0)), mode='constant')

    # Output dimensions
    H_out = (H + 2*padding - K) // stride + 1
    W_out = (W + 2*padding - K) // stride + 1

    output = np.zeros((H_out, W_out, C_out))

    for h in range(H_out):
        for w in range(W_out):
            for c in range(C_out):
                # Extract region
                h_start = h * stride
                w_start = w * stride
                region = image[h_start:h_start+K, w_start:w_start+K, :]

                # Convolve
                output[h, w, c] = np.sum(region * kernel[:, :, :, c])

    return output

class ResidualBlock(nn.Module):
    def __init__(self, channels):
        super().__init__()
        self.conv1 = nn.Conv2d(channels, channels, 3, padding=1)
        self.bn1 = nn.BatchNorm2d(channels)
        self.conv2 = nn.Conv2d(channels, channels, 3, padding=1)
        self.bn2 = nn.BatchNorm2d(channels)

    def forward(self, x):
        residual = x
        out = F.relu(self.bn1(self.conv1(x)))
        out = self.bn2(self.conv2(out))
        out = out + residual  # skip connection
        return F.relu(out)

def rnn_step(x_t, h_prev, W_h, W_x, b):
    """
    One timestep of an RNN.
    x_t: input at time t
    h_prev: hidden state from t-1
    """
    return np.tanh(W_h @ h_prev + W_x @ x_t + b)

Forget gate:    f_t = σ(W_f · [h_{t-1}, x_t] + b_f)
Input gate:     i_t = σ(W_i · [h_{t-1}, x_t] + b_i)
Candidate:      Ċ_t = tanh(W_c · [h_{t-1}, x_t] + b_c)
Cell update:    c_t = f_t ⊙ c_{t-1} + i_t ⊙ Ċ_t
Output gate:    o_t = σ(W_o · [h_{t-1}, x_t] + b_o)
Hidden state:   h_t = o_t ⊙ tanh(c_t)

class LSTMCell:
    def __init__(self, input_dim, hidden_dim):
        # Concatenated weights for efficiency
        self.W = np.random.randn(input_dim + hidden_dim, 4 * hidden_dim) * 0.01
        self.b = np.zeros(4 * hidden_dim)
        self.hidden_dim = hidden_dim

    def forward(self, x, h_prev, c_prev):
        """
        x: (input_dim,)
        h_prev: (hidden_dim,) — previous hidden state
        c_prev: (hidden_dim,) — previous cell state
        """
        # Concatenate input and previous hidden
        combined = np.concatenate([x, h_prev])

        # Compute all gates at once
        gates = self.W.T @ combined + self.b

        # Split into forget, input, candidate, output gates
        f = sigmoid(gates[:self.hidden_dim])
        i = sigmoid(gates[self.hidden_dim:2*self.hidden_dim])
        c_tilde = np.tanh(gates[2*self.hidden_dim:3*self.hidden_dim])
        o = sigmoid(gates[3*self.hidden_dim:])

        # Update cell state and hidden state
        c = f * c_prev + i * c_tilde
        h = o * np.tanh(c)

        return h, c

# In practice: use weight_decay parameter in optimizer
optimizer = torch.optim.AdamW(model.parameters(), lr=1e-3, weight_decay=0.01)

def dropout(x, p=0.5, training=True):
    """
    Dropout: randomly zero out p fraction of neurons.
    During testing: scale by p (or use all neurons).
    """
    if not training:
        return x

    mask = np.random.binomial(1, 1-p, x.shape) / (1-p)
    return x * mask

class BatchNorm:
    def __init__(self, num_features, eps=1e-5, momentum=0.1):
        self.gamma = np.ones(num_features)
        self.beta = np.zeros(num_features)
        self.running_mean = np.zeros(num_features)
        self.running_var = np.ones(num_features)
        self.eps = eps
        self.momentum = momentum

    def forward(self, x, training=True):
        if training:
            batch_mean = x.mean(axis=0)
            batch_var = x.var(axis=0)

            self.running_mean = (1 - self.momentum) * self.running_mean + self.momentum * batch_mean
            self.running_var = (1 - self.momentum) * self.running_var + self.momentum * batch_var

            x_norm = (x - batch_mean) / np.sqrt(batch_var + self.eps)
        else:
            x_norm = (x - self.running_mean) / np.sqrt(self.running_var + self.eps)

        return self.gamma * x_norm + self.beta

# Vanilla SGD
θ = θ - lr * ∇L(θ)

v = momentum * v - lr * grad
θ = θ + v

def adam(params, grads, m, v, t, lr=0.001, beta1=0.9, beta2=0.999, eps=1e-8):
    """
    One step of Adam optimizer.
    """
    t += 1

    for i, (param, grad) in enumerate(zip(params, grads)):
        # Biased first moment estimate (momentum)
        m[i] = beta1 * m[i] + (1 - beta1) * grad

        # Biased second moment estimate (adaptive lr)
        v[i] = beta2 * v[i] + (1 - beta2) * grad**2

        # Bias correction
        m_hat = m[i] / (1 - beta1**t)
        v_hat = v[i] / (1 - beta2**t)

        # Update
        param[i] -= lr * m_hat / (np.sqrt(v_hat) + eps)

    return params, m, v, t

def scaled_dot_product_attention(Q, K, V, mask=None):
    """
    Q: (batch, seq_len_q, d_k) — queries
    K: (batch, seq_len_k, d_k) — keys  
    V: (batch, seq_len_v, d_v) — values
                (seq_len_k must equal seq_len_v)
    """
    # Compute attention scores
    scores = Q @ K.transpose(-2, -1)  # (batch, q_len, k_len)

    # Scale to prevent softmax saturation
    scores = scores / np.sqrt(K.shape[-1])

    # Optional mask (for causal attention)
    if mask is not None:
        scores = scores + mask  # mask has -inf for forbidden positions

    # Softmax over key dimension
    weights = np.exp(scores - scores.max(axis=-1, keepdims=True))
    weights = weights / weights.sum(axis=-1, keepdims=True)

    # Weighted sum of values
    output = weights @ V  # (batch, q_len, d_v)

    return output, weights

class MultiHeadAttention(nn.Module):
    def __init__(self, d_model, n_heads):
        super().__init__()
        self.n_heads = n_heads
        self.d_head = d_model // n_heads

        # All projections
        self.W_q = nn.Linear(d_model, d_model)
        self.W_k = nn.Linear(d_model, d_model)
        self.W_v = nn.Linear(d_model, d_model)
        self.W_o = nn.Linear(d_model, d_model)

    def forward(self, x, mask=None):
        B, L, D = x.shape

        # Project and reshape: (B, L, D) → (B, n_heads, L, d_head)
        Q = self.W_q(x).view(B, L, self.n_heads, self.d_head).transpose(1, 2)
        K = self.W_k(x).view(B, L, self.n_heads, self.d_head).transpose(1, 2)
        V = self.W_v(x).view(B, L, self.n_heads, self.d_head).transpose(1, 2)

        # Compute attention
        scores = Q @ K.transpose(-2, -1) / np.sqrt(self.d_head)

        if mask is not None:
            scores = scores + mask

        weights = F.softmax(scores, dim=-1)
        output = weights @ V  # (B, n_heads, L, d_head)

        # Concatenate heads: (B, n_heads, L, d_head) → (B, L, D)
        output = output.transpose(1, 2).contiguous().view(B, L, D)

        return self.W_o(output)

def sinusoidal_positional_encoding(seq_len, d_model):
    """Create sinusoidal position encodings."""
    pe = np.zeros((seq_len, d_model))
    position = np.arange(0, seq_len)[:, np.newaxis]
    div_term = np.exp(np.arange(0, d_model, 2) * -(np.log(10000.0) / d_model))

    pe[:, 0::2] = np.sin(position * div_term)
    pe[:, 1::2] = np.cos(position * div_term)
    return pe

def create_causal_mask(seq_len):
    """
    Causal mask: each token attends to itself + previous tokens only.
    Upper triangle = -inf (can't attend to future).
    """
    mask = np.triu(np.ones((seq_len, seq_len)) * -1e9, k=1)
    return mask

# With causal mask:
# Position 0 can attend to: [0]
# Position 1 can attend to: [0, 1]
# Position 2 can attend to: [0, 1, 2]
# ...

class Transformer(nn.Module):
    """
    Full encoder-decoder transformer (original Vaswani et al. 2017).
    """
    def __init__(self, d_model=512, n_heads=8, n_layers=6, d_ff=2048, vocab_size=30000):
        super().__init__()
        self.encoder = nn.TransformerEncoder(
            nn.TransformerEncoderLayer(d_model, n_heads, d_ff),
            num_layers=n_layers
        )
        self.decoder = nn.TransformerDecoder(
            nn.TransformerDecoderLayer(d_model, n_heads, d_ff),
            num_layers=n_layers
        )
        self.embedding = nn.Embedding(vocab_size, d_model)
        self.output = nn.Linear(d_model, vocab_size)

    def forward(self, src_tokens, tgt_tokens):
        src_emb = self.embedding(src_tokens)
        tgt_emb = self.embedding(tgt_tokens)

        # Encoder: process input bidirectionally
        memory = self.encoder(src_emb)

        # Decoder: generate output with cross-attention to encoder
        output = self.decoder(tgt_emb, memory)

        return self.output(output)

Input:  "The [MASK] chased the [MASK]"
Target: "The cat chased the mouse"

class BERTEmbedding(nn.Module):
    def __init__(self, vocab_size, d_model=768, max_len=512):
        super().__init__()
        self.token_embedding = nn.Embedding(vocab_size, d_model)
        self.position_embedding = nn.Embedding(max_len, d_model)
        self.segment_embedding = nn.Embedding(2, d_model)  # sentence A/B

    def forward(self, token_ids, segment_ids):
        positions = torch.arange(token_ids.shape[1], device=token_ids.device)

        return (self.token_embedding(token_ids) + 
                self.position_embedding(positions) + 
                self.segment_embedding(segment_ids))

class GroupedQueryAttention(nn.Module):
    def __init__(self, d_model, n_heads, n_kv_heads):
        super().__init__()
        self.n_heads = n_heads           # total query heads
        self.n_kv_heads = n_kv_heads     # key/value heads (n_heads // group_size)
        self.d_head = d_model // n_heads

        self.W_q = nn.Linear(d_model, n_heads * self.d_head, bias=False)
        self.W_k = nn.Linear(d_model, n_kv_heads * self.d_head, bias=False)
        self.W_v = nn.Linear(d_model, n_kv_heads * self.d_head, bias=False)
        self.W_o = nn.Linear(n_heads * self.d_head, d_model, bias=False)

    def forward(self, x):
        B, L, _ = x.shape

        Q = self.W_q(x).view(B, L, self.n_heads, self.d_head).transpose(1, 2)
        K = self.W_k(x).view(B, L, self.n_kv_heads, self.d_head).transpose(1, 2)
        V = self.W_v(x).view(B, L, self.n_kv_heads, self.d_head).transpose(1, 2)

        # Repeat K, V for each query head in the group
        n_groups = self.n_heads // self.n_kv_heads
        K = K.repeat_interleave(n_groups, dim=1)
        V = V.repeat_interleave(n_groups, dim=1)

        # Now all shapes match: (B, n_heads, L, d_head)
        scores = Q @ K.transpose(-2, -1) / (self.d_head ** 0.5)
        weights = F.softmax(scores, dim=-1)
        output = weights @ V

        output = output.transpose(1, 2).contiguous().view(B, L, -1)
        return self.W_o(output)

def generate_with_kv_cache(model, input_ids, max_new_tokens=100):
    """
    Efficient generation using KV-cache.
    Instead of recomputing K,V for all previous tokens each step,
    we cache them and only compute for the new token.
    """
    kv_cache = {}

    for i in range(max_new_tokens):
        if i == 0:
            # First step: process full input sequence
            logits, kv_cache = model(input_ids, use_cache=True)
        else:
            # Subsequent steps: only process the last token
            logits, kv_cache = model(input_ids[:, -1:], kv_cache=kv_cache, use_cache=True)

        next_token = sample(logits[:, -1, :])
        input_ids = torch.cat([input_ids, next_token], dim=-1)

        if next_token.item() == tokenizer.eos_token_id:
            break

    return input_ids

n = 1,000:  1M operations
n = 10,000: 100M operations
n = 100,000: 10B operations
n = 1,000,000: 1T operations ← infeasible

# Standard attention: writes N×N matrix to HBM (slow)
S = Q @ K.T     # read from HBM
P = softmax(S)  # write to HBM
O = P @ V       # read from HBM

# Flash attention: process in tiles, stay in SRAM
# Never materialize full N×N matrix
# Hardware-aware tiling

# Simplified Mamba block processing one sequence
def mamba_block(x, dt_linear, A, B_proj, C_proj):
    """
    x: (batch, seq_len, d_model)
    Returns: (batch, seq_len, d_model)
    """
    B, L, D = x.shape

    # Input-dependent parameters
    delta = F.softplus(dt_linear(x))  # (B, L, D) — step size, per token

    B = B_proj(x)  # (B, L, D) — input-dependent B
    C = C_proj(x)  # (B, L, D) — input-dependent C

    # Discretize A
    A_bar = torch.exp(delta.unsqueeze(-1) * A)  # (B, L, D, N)
    B_bar = delta * B  # simplified

    # Scan: h_t = A_bar_t * h_{t-1} + B_bar_t * x_t
    h = torch.zeros(B, D, device=x.device)
    ys = []

    for t in range(L):
        h = A_bar[:, t] * h + B_bar[:, t] * x[:, t]
        y = (C[:, t] * h).sum(-1)
        ys.append(y)

    return torch.stack(ys, dim=1)

# In practice: use parallel associative scan (O(n log n) instead of O(n))

class GCNLayer(nn.Module):
    """
    Graph Convolutional Network layer.
    AGGREGATE = mean of neighbors.
    """
    def __init__(self, in_dim, out_dim):
        super().__init__()
        self.W = nn.Linear(in_dim, out_dim)

    def forward(self, x, adj_matrix):
        """
        x: (num_nodes, in_dim)
        adj_matrix: (num_nodes, num_nodes) — normalized adjacency
        """
        # Message passing: neighbor features aggregated by adjacency
        neighbor_aggregated = adj_matrix @ x
        return torch.relu(self.W(neighbor_aggregated))

class VAE(nn.Module):
    def __init__(self, input_dim=784, latent_dim=20):
        super().__init__()
        # Encoder: x → z (Gaussian parameters)
        self.encoder = nn.Sequential(
            nn.Linear(input_dim, 256),
            nn.ReLU(),
        )
        self.mu_layer = nn.Linear(256, latent_dim)
        self.logvar_layer = nn.Linear(256, latent_dim)

        # Decoder: z → x
        self.decoder = nn.Sequential(
            nn.Linear(latent_dim, 256),
            nn.ReLU(),
            nn.Linear(256, input_dim),
            nn.Sigmoid(),
        )

    def encode(self, x):
        h = self.encoder(x)
        return self.mu_layer(h), self.logvar_layer(h)

    def reparameterize(self, mu, logvar):
        """Reparameterization trick: z = μ + σ * ε"""
        std = torch.exp(0.5 * logvar)
        eps = torch.randn_like(std)
        return mu + eps * std

    def decode(self, z):
        return self.decoder(z)

    def forward(self, x):
        mu, logvar = self.encode(x)
        z = self.reparameterize(mu, logvar)
        x_recon = self.decode(z)
        return x_recon, mu, logvar, z

def vae_loss(x, x_recon, mu, logvar):
    """Reconstruction + KL divergence loss."""
    # Reconstruction: binary cross-entropy
    recon_loss = F.binary_cross_entropy(x_recon, x, reduction='sum')

    # KL divergence: q(z|x) || N(0, I)
    # Closed form for Gaussian: -½Σ(1 + log(σ²) - μ² - σ²)
    kl_loss = -0.5 * torch.sum(1 + logvar - mu.pow(2) - logvar.exp())

    return recon_loss + kl_loss

class GAN(nn.Module):
    def __init__(self, latent_dim=100):
        super().__init__()
        self.generator = nn.Sequential(
            nn.Linear(latent_dim, 256),
            nn.ReLU(),
            nn.Linear(256, 512),
            nn.ReLU(),
            nn.Linear(512, 784),
            nn.Tanh(),  # images [-1, 1]
        )

        self.discriminator = nn.Sequential(
            nn.Linear(784, 512),
            nn.LeakyReLU(0.2),
            nn.Linear(512, 256),
            nn.LeakyReLU(0.2),
            nn.Linear(256, 1),
            nn.Sigmoid(),
        )

def train_gan_step(gan, real_images, optimizer_G, optimizer_D, latent_dim=100):
    batch_size = real_images.size(0)

    # --- Train Discriminator ---
    optimizer_D.zero_grad()

    # Real images: D should output 1
    real_pred = gan.discriminator(real_images)
    real_loss = F.binary_cross_entropy(real_pred, torch.ones_like(real_pred))

    # Fake images: D should output 0
    z = torch.randn(batch_size, latent_dim)
    fake_images = gan.generator(z)
    fake_pred = gan.discriminator(fake_images.detach())  # stop gradient to G
    fake_loss = F.binary_cross_entropy(fake_pred, torch.zeros_like(fake_pred))

    d_loss = (real_loss + fake_loss) / 2
    d_loss.backward()
    optimizer_D.step()

    # --- Train Generator ---
    optimizer_G.zero_grad()

    z = torch.randn(batch_size, latent_dim)
    fake_images = gan.generator(z)
    fake_pred = gan.discriminator(fake_images)

    # G wants D to output 1 on its fakes
    g_loss = F.binary_cross_entropy(fake_pred, torch.ones_like(fake_pred))
    g_loss.backward()
    optimizer_G.step()

    return d_loss.item(), g_loss.item()