Logistic Regression

import numpy as np
import pandas as pd
from typing import Literal, Tuple, Optional

class LogisticRegression:
    """
    Logistic Regression with multiple optimization variants.

    Supports:
    - Optimizers: BGD, SGD, Mini-Batch SGD, SGD+Momentum, Adam
    - Regularization: None, L1, L2, Elastic Net

    Example:
        model = LogisticRegression(optimizer='adam', reg_type='l2')
        model.fit(X_train, y_train)
        y_pred = model.predict(X_test)
    """

    def __init__(
        self,
        learning_rate: float = 0.01,
        n_iterations: int = 1000,
        optimizer: Literal['bgd', 'sgd', 'minibatch', 'momentum', 'adam'] = 'bgd',
        batch_size: int = 32,
        reg_type: Literal['none', 'l1', 'l2', 'elasticnet'] = 'none',
        reg_lambda: float = 0.01,
        l1_ratio: float = 0.5,
        momentum_beta: float = 0.9,
        adam_beta1: float = 0.9,
        adam_beta2: float = 0.999,
        adam_epsilon: float = 1e-8,
        verbose: bool = False
    ):
        self.learning_rate = learning_rate
        self.n_iterations = n_iterations
        self.optimizer = optimizer
        self.batch_size = batch_size
        self.reg_type = reg_type
        self.reg_lambda = reg_lambda
        self.l1_ratio = l1_ratio
        self.momentum_beta = momentum_beta
        self.adam_beta1 = adam_beta1
        self.adam_beta2 = adam_beta2
        self.adam_epsilon = adam_epsilon
        self.verbose = verbose

        # Initialize weights and biases
        self.w = None
        self.b = 0.0
        self.loss_history = []

        # For momentum and Adam
        self.velocity = None  # for momentum
        self.m = None  # for Adam (1st moment)
        self.v = None  # for Adam (2nd moment)

    def _sigmoid(self, z: np.ndarray) -> np.ndarray:
        """Sigmoid activation: 1 / (1 + exp(-z))"""
        # Clip to prevent overflow
        z = np.clip(z, -500, 500)
        return 1 / (1 + np.exp(-z))

    def _compute_loss(self, X: np.ndarray, y: np.ndarray) -> float:
        """
        Compute binary cross-entropy loss + regularization penalty.

        J(w) = -(1/N) * Σ[y*log(p̂) + (1-y)*log(1-p̂)] + reg_penalty
        """
        m = X.shape[0]
        z = X @ self.w + self.b
        p_hat = self._sigmoid(z)

        # Clip to prevent log(0)
        p_hat = np.clip(p_hat, 1e-15, 1 - 1e-15)

        # BCE loss
        bce_loss = -(1/m) * np.sum(
            y * np.log(p_hat) + (1 - y) * np.log(1 - p_hat)
        )

        # Regularization penalty
        if self.reg_type == 'l1':
            reg_penalty = self.reg_lambda * np.sum(np.abs(self.w))
        elif self.reg_type == 'l2':
            reg_penalty = self.reg_lambda * np.sum(self.w ** 2)
        elif self.reg_type == 'elasticnet':
            reg_penalty = self.reg_lambda * (
                self.l1_ratio * np.sum(np.abs(self.w)) +
                (1 - self.l1_ratio) * np.sum(self.w ** 2)
            )
        else:
            reg_penalty = 0.0

        return bce_loss + reg_penalty

    def _compute_gradients(self, X: np.ndarray, y: np.ndarray) -> Tuple[np.ndarray, float]:
        """
        Compute gradient of loss w.r.t. weights and bias.

        ∂J/∂w_j = (1/N) * Σ(p̂_i - y_i) * x_{i,j} + reg_gradient
        ∂J/∂b = (1/N) * Σ(p̂_i - y_i)
        """
        m = X.shape[0]
        z = X @ self.w + self.b
        p_hat = self._sigmoid(z)

        error = p_hat - y
        grad_w = (1/m) * (X.T @ error)
        grad_b = (1/m) * np.sum(error)

        return grad_w, grad_b

    def _apply_regularization_gradient(self, grad_w: np.ndarray) -> np.ndarray:
        """Add regularization term to gradient."""
        if self.reg_type == 'l1':
            grad_w += self.reg_lambda * np.sign(self.w)
        elif self.reg_type == 'l2':
            grad_w += 2 * self.reg_lambda * self.w
        elif self.reg_type == 'elasticnet':
            grad_w += self.reg_lambda * (
                self.l1_ratio * np.sign(self.w) +
                2 * (1 - self.l1_ratio) * self.w
            )
        return grad_w

    def _update_bgd(self, X: np.ndarray, y: np.ndarray) -> None:
        """Batch gradient descent: update on entire dataset."""
        grad_w, grad_b = self._compute_gradients(X, y)
        grad_w = self._apply_regularization_gradient(grad_w)

        self.w -= self.learning_rate * grad_w
        self.b -= self.learning_rate * grad_b

    def _update_sgd(self, X: np.ndarray, y: np.ndarray) -> None:
        """Stochastic gradient descent: update on one sample at a time."""
        indices = np.random.permutation(X.shape[0])
        for i in indices:
            x_i = X[i:i+1]
            y_i = y[i:i+1]

            grad_w, grad_b = self._compute_gradients(x_i, y_i)
            grad_w = self._apply_regularization_gradient(grad_w)

            self.w -= self.learning_rate * grad_w
            self.b -= self.learning_rate * grad_b

    def _update_minibatch(self, X: np.ndarray, y: np.ndarray) -> None:
        """Mini-batch SGD: update on batch_size samples."""
        indices = np.random.permutation(X.shape[0])
        for start_idx in range(0, X.shape[0], self.batch_size):
            end_idx = min(start_idx + self.batch_size, X.shape[0])
            batch_indices = indices[start_idx:end_idx]

            X_batch = X[batch_indices]
            y_batch = y[batch_indices]

            grad_w, grad_b = self._compute_gradients(X_batch, y_batch)
            grad_w = self._apply_regularization_gradient(grad_w)

            self.w -= self.learning_rate * grad_w
            self.b -= self.learning_rate * grad_b

    def _update_momentum(self, X: np.ndarray, y: np.ndarray) -> None:
        """SGD with momentum: accumulate velocity direction."""
        if self.velocity is None:
            self.velocity = np.zeros_like(self.w)

        indices = np.random.permutation(X.shape[0])
        for i in indices:
            x_i = X[i:i+1]
            y_i = y[i:i+1]

            grad_w, grad_b = self._compute_gradients(x_i, y_i)
            grad_w = self._apply_regularization_gradient(grad_w)

            # Accumulate velocity
            self.velocity = self.momentum_beta * self.velocity - self.learning_rate * grad_w
            self.w += self.velocity
            self.b -= self.learning_rate * grad_b

    def _update_adam(self, X: np.ndarray, y: np.ndarray) -> None:
        """Adam optimizer: adaptive per-parameter learning rates."""
        if self.m is None:
            self.m = np.zeros_like(self.w)
            self.v = np.zeros_like(self.w)

        indices = np.random.permutation(X.shape[0])
        for t, i in enumerate(indices, 1):
            x_i = X[i:i+1]
            y_i = y[i:i+1]

            grad_w, grad_b = self._compute_gradients(x_i, y_i)
            grad_w = self._apply_regularization_gradient(grad_w)

            # Update biased first moment estimate (mean)
            self.m = self.adam_beta1 * self.m + (1 - self.adam_beta1) * grad_w

            # Update biased second moment estimate (variance)
            self.v = self.adam_beta2 * self.v + (1 - self.adam_beta2) * (grad_w ** 2)

            # Bias correction
            m_hat = self.m / (1 - self.adam_beta1 ** t)
            v_hat = self.v / (1 - self.adam_beta2 ** t)

            # Update parameters
            self.w -= self.learning_rate * m_hat / (np.sqrt(v_hat) + self.adam_epsilon)
            self.b -= self.learning_rate * grad_b

    def fit(self, X: np.ndarray, y: np.ndarray) -> None:
        """
        Train the logistic regression model.

        Args:
            X: (n_samples, n_features) feature matrix
            y: (n_samples,) binary labels {0, 1}
        """
        # Initialize weights
        if self.w is None:
            self.w = np.zeros(X.shape[1])

        # Training loop
        for iteration in range(self.n_iterations):
            # Update weights based on chosen optimizer
            if self.optimizer == 'bgd':
                self._update_bgd(X, y)
            elif self.optimizer == 'sgd':
                self._update_sgd(X, y)
            elif self.optimizer == 'minibatch':
                self._update_minibatch(X, y)
            elif self.optimizer == 'momentum':
                self._update_momentum(X, y)
            elif self.optimizer == 'adam':
                self._update_adam(X, y)

            # Record loss
            loss = self._compute_loss(X, y)
            self.loss_history.append(loss)

            if self.verbose and (iteration + 1) % max(1, self.n_iterations // 10) == 0:
                print(f"Iteration {iteration + 1}/{self.n_iterations}, Loss: {loss:.6f}")

    def predict_proba(self, X: np.ndarray) -> np.ndarray:
        """Return predicted probabilities P(y=1 | X)."""
        z = X @ self.w + self.b
        return self._sigmoid(z)

    def predict(self, X: np.ndarray, threshold: float = 0.5) -> np.ndarray:
        """Return binary predictions {0, 1}."""
        return (self.predict_proba(X) >= threshold).astype(int)


# ============ ADTECH EXAMPLE: PREDICTING AD CLICKS ============

# Generate synthetic AdTech click data
np.random.seed(42)
n_samples = 10000

# Features: user_age, hour_of_day, is_mobile, ad_relevance, user_past_ctr, page_score, ad_position
X = pd.DataFrame({
    'user_age_bucket': np.random.randint(1, 8, n_samples),
    'hour_of_day': np.random.randint(0, 24, n_samples),
    'is_mobile': np.random.binomial(1, 0.4, n_samples),
    'ad_relevance': np.random.rand(n_samples),
    'user_past_ctr': np.random.rand(n_samples) * 0.05,
    'page_context_score': np.random.rand(n_samples),
    'ad_position': np.random.randint(1, 6, n_samples),
})

# Normalize features
X = (X - X.mean()) / X.std()
X_array = X.values

# Generate labels with signal from features
true_weights = np.array([0.3, 0.1, -0.2, 0.5, 0.4, 0.2, -0.15])
logits = X_array @ true_weights + np.random.normal(0, 0.3, n_samples)
click_probs = 1 / (1 + np.exp(-logits))
y = (click_probs > np.random.rand(n_samples)).astype(int)

print(f"Dataset: {n_samples} samples, {X.shape[1]} features")
print(f"Click rate: {y.mean():.2%}")
print()

# Split into train/test
split = int(0.8 * n_samples)
X_train, X_test = X_array[:split], X_array[split:]
y_train, y_test = y[:split], y[split:]

# Compare all optimizers with L2 regularization
optimizers = ['bgd', 'sgd', 'minibatch', 'momentum', 'adam']
results = {}

print("=" * 70)
print("COMPARING OPTIMIZERS (L2 Regularization, λ=0.01)")
print("=" * 70)

for opt in optimizers:
    model = LogisticRegression(
        optimizer=opt,
        learning_rate=0.01,
        n_iterations=500,
        batch_size=64,
        reg_type='l2',
        reg_lambda=0.01,
        verbose=False
    )

    model.fit(X_train, y_train)

    # Compute metrics
    train_proba = model.predict_proba(X_train)
    test_proba = model.predict_proba(X_test)

    train_pred = (train_proba >= 0.5).astype(int)
    test_pred = (test_proba >= 0.5).astype(int)

    train_acc = (train_pred == y_train).mean()
    test_acc = (test_pred == y_test).mean()
    final_loss = model.loss_history[-1]

    results[opt] = {
        'train_acc': train_acc,
        'test_acc': test_acc,
        'final_loss': final_loss,
        'convergence': min(model.loss_history),
        'model': model
    }

    print(f"{opt.upper():12} | Loss: {final_loss:.6f} | "
          f"Train Acc: {train_acc:.4f} | Test Acc: {test_acc:.4f}")

print()
print("=" * 70)
print("COMPARING REGULARIZATION TYPES (Mini-Batch SGD)")
print("=" * 70)

reg_types = ['none', 'l1', 'l2', 'elasticnet']
for reg in reg_types:
    model = LogisticRegression(
        optimizer='minibatch',
        learning_rate=0.01,
        n_iterations=500,
        batch_size=64,
        reg_type=reg,
        reg_lambda=0.01 if reg != 'none' else 0,
        l1_ratio=0.5,  # for elastic net
        verbose=False
    )

    model.fit(X_train, y_train)

    test_pred = model.predict(X_test)
    test_acc = (test_pred == y_test).mean()
    final_loss = model.loss_history[-1]
    sparsity = (model.w == 0).mean()

    print(f"{reg.upper():12} | Loss: {final_loss:.6f} | "
          f"Test Acc: {test_acc:.4f} | Sparsity: {sparsity:.2%}")

print()
print("=" * 70)
print("BEST MODEL WEIGHTS (Adam, L2)")
print("=" * 70)
best_model = results['adam']['model']
feature_names = list(X.columns)
for fname, weight in zip(feature_names, best_model.w):
    print(f"{fname:20} : {weight:+.6f}")