Perceptron

"""The perceptron - mistake-driven updates, transcribed.

If the prediction is right, change nothing. If it is wrong, push the boundary
toward the answer by adding the example to the weights. That is the whole
algorithm, and it is also stochastic gradient descent on the criterion above -
decades before it was called that.
"""


def predict(weights, bias, x):
    """Sign of the score. NOT a probability - the magnitude means nothing."""
    score = bias
    for j in range(len(x)):
        score += weights[j] * x[j]
    return 1 if score > 0 else -1


def fit(X, y, learning_rate=1.0, epochs=20):
    """y must be in {-1, +1}: the update is built on the sign convention."""
    d = len(X[0])
    weights = [0.0] * d
    bias = 0.0

    for epoch in range(epochs):
        mistakes = 0

        for i in range(len(X)):
            # margin = y_i (w . x_i + b); at most zero means misclassified
            score = bias
            for j in range(d):
                score += weights[j] * X[i][j]
            margin = y[i] * score

            if margin <= 0.0:
                mistakes += 1
                for j in range(d):
                    weights[j] += learning_rate * y[i] * X[i][j]
                bias += learning_rate * y[i]

        # On separable data the mistakes stop and the algorithm terminates.
        # On non-separable data they never do - the weight vector cycles
        # forever, and nothing here reports that. The epoch budget is the
        # only thing that stops it.
        if mistakes == 0:
            return weights, bias, epoch + 1

    return weights, bias, epochs


def mistake_bound(X, margin):
    """(R / gamma)^2 - the theoretical bound, purely geometric.

    R is the largest example norm and gamma the margin of the best separating
    hyperplane. Note what is absent: the number of examples, and the number of
    features. Neither appears.
    """
    radius = 0.0
    for row in X:
        norm = sum(value * value for value in row) ** 0.5
        radius = max(radius, norm)
    return (radius / margin) ** 2

"""Averaged perceptron - typed, sparse, with the lazy averaging trick."""

from dataclasses import dataclass, field

import numpy as np
from numpy.typing import NDArray

Vector = NDArray[np.float64]


@dataclass
class AveragedPerceptron:
    """Returns the AVERAGE of every intermediate weight vector, not the last.

    The plain perceptron ends wherever its final update left it, which may
    have come from a single odd example. Averaging over the whole run closes
    most of the accuracy gap to far more expensive methods, and the naive
    implementation costs O(d) per example. The trick below makes it O(nnz)
    per MISTAKE instead, which is what makes averaging free enough to be the
    default rather than an option.
    """

    n_features: int
    weights: Vector = field(init=False)
    _totals: Vector = field(init=False)      # accumulated weights * time
    _stamps: NDArray = field(init=False)     # when each weight last changed
    _step: int = field(default=0, init=False)

    def __post_init__(self) -> None:
        if self.n_features < 1:
            raise ValueError(f"n_features must be at least 1, got {self.n_features}")
        self.weights = np.zeros(self.n_features, dtype=np.float64)
        self._totals = np.zeros(self.n_features, dtype=np.float64)
        self._stamps = np.zeros(self.n_features, dtype=np.int64)

    def score(self, indices: NDArray, values: Vector) -> float:
        """Sparse dot product: only present features are touched at all."""
        return float(self.weights[indices] @ values)

    def update(self, indices: NDArray, values: Vector, label: int) -> None:
        """One mistake-driven update, with averaging maintained lazily.

        Instead of adding the whole weight vector to a running total on every
        example, each coordinate carries the timestep it last changed. When it
        changes again, the intervening span is settled in one multiply - so
        the accumulated total stays exact while only touched coordinates are
        ever visited.
        """
        self._step += 1

        span = self._step - self._stamps[indices]
        self._totals[indices] += span * self.weights[indices]
        self._stamps[indices] = self._step

        self.weights[indices] += label * values

    def averaged(self) -> Vector:
        """Settle every coordinate to the current step, then average."""
        span = self._step - self._stamps
        return (self._totals + span * self.weights) / max(self._step, 1)


def fit(
    rows: list[tuple[NDArray, Vector]],
    labels: NDArray,
    n_features: int,
    epochs: int = 10,
    seed: int = 0,
) -> tuple[Vector, list[int]]:
    """Train on sparse rows. Raises ValueError on malformed input."""
    if len(rows) != labels.size:
        raise ValueError(f"{len(rows)} rows but {labels.size} labels")
    if not np.isin(labels, (-1, 1)).all():
        raise ValueError("labels must be -1 or +1; 0/1 silently optimizes something else")

    model = AveragedPerceptron(n_features=n_features)
    rng = np.random.default_rng(seed)
    order = np.arange(len(rows))
    history: list[int] = []

    for _ in range(epochs):
        # Order matters: the algorithm is order-dependent by construction, so
        # reshuffling each epoch is part of the method rather than hygiene.
        rng.shuffle(order)
        mistakes = 0

        for index in order:
            indices, values = rows[index]
            if labels[index] * model.score(indices, values) <= 0.0:
                model.update(indices, values, int(labels[index]))
                mistakes += 1

        history.append(mistakes)
        # A mistake count that stops falling means the data is not separable
        # under these features - which is information about the features.
        if mistakes == 0:
            break

    return model.averaged(), history

"""Perceptron - one contiguous CSR block, in-place scatter updates."""

import numpy as np
from numpy.typing import NDArray

Vector = NDArray[np.float64]


class CsrPerceptron:
    """Sparse rows in one flat block, not a list of arrays.

    The previous stage held each row as its own pair of arrays, which means a
    pointer chase and a bounds check per example. A CSR layout puts every
    index and value in one contiguous buffer, so a row is a slice and the
    whole dataset streams through cache in order.

    The algorithm itself stays sequential - each update depends on the
    boundary the previous updates produced, and no rewrite changes that. What
    is available is making every individual update as cheap as possible.
    """

    def __init__(self, n_features: int) -> None:
        # A single dtype throughout: float64 weights, int32 indices. Mixed
        # dtypes force a conversion on every gather.
        self._weights = np.zeros(n_features, dtype=np.float64)
        self._totals = np.zeros(n_features, dtype=np.float64)
        self._stamps = np.zeros(n_features, dtype=np.int64)
        self._step = 0
        # Scratch reused by every update rather than allocated per row.
        self._scratch = np.empty(n_features, dtype=np.float64)

    def fit(
        self,
        indptr: NDArray,
        indices: NDArray,
        values: Vector,
        labels: NDArray,
        epochs: int = 10,
        seed: int = 0,
    ) -> Vector:
        rng = np.random.default_rng(seed)
        order = np.arange(indptr.size - 1)

        for _ in range(epochs):
            rng.shuffle(order)
            mistakes = 0

            for row in order:
                start, stop = indptr[row], indptr[row + 1]
                cols = indices[start:stop]
                vals = values[start:stop]

                # Gather-and-dot over present features only. For a row with
                # 30 non-zeros out of a million features, this is 30 loads.
                if labels[row] * (self._weights[cols] @ vals) > 0.0:
                    continue

                self._step += 1
                mistakes += 1

                # Lazy averaging: settle only the coordinates being touched.
                span = self._step - self._stamps[cols]
                # In place, through a view rather than a temporary: at high
                # update rates these allocations dominate the arithmetic.
                np.multiply(span, self._weights[cols], out=self._scratch[: cols.size])
                np.add.at(self._totals, cols, self._scratch[: cols.size])
                self._stamps[cols] = self._step

                np.add.at(self._weights, cols, labels[row] * vals)

            if mistakes == 0:
                break

        span = self._step - self._stamps
        return (self._totals + span * self._weights) / max(self._step, 1)