cds2.ml

Classical machine learning models and metrics.

Classical machine-learning models, utilities and metrics.

StandardScaler

Zero-mean unit-variance feature scaler.

Source code in src/cds2/ml.py

class StandardScaler:
    """Zero-mean unit-variance feature scaler."""

    def __init__(self) -> None:
        self.mean_: FloatArray | None = None
        self.scale_: FloatArray | None = None

    def fit(self, x: object) -> StandardScaler:
        arr = np.asarray(x, dtype=float)
        if arr.ndim == 1:
            arr = arr.reshape(-1, 1)
        self.mean_ = np.asarray(arr.mean(axis=0))
        std = arr.std(axis=0)
        self.scale_ = np.where(std == 0, 1.0, std)
        return self

    def transform(self, x: object) -> FloatArray:
        if self.mean_ is None or self.scale_ is None:
            msg = "scaler is not fitted"
            raise RuntimeError(msg)
        arr = np.asarray(x, dtype=float)
        if arr.ndim == 1:
            arr = arr.reshape(-1, 1)
        return np.asarray((arr - self.mean_) / self.scale_)

    def fit_transform(self, x: object) -> FloatArray:
        return self.fit(x).transform(x)

LinearRegression

Ordinary least-squares linear regression with an intercept.

Source code in src/cds2/ml.py

class LinearRegression:
    """Ordinary least-squares linear regression with an intercept."""

    def __init__(self) -> None:
        self.coef_: FloatArray | None = None
        self.intercept_: float = 0.0

    def fit(self, x: object, y: object) -> LinearRegression:
        design = self._design(x)
        target = np.asarray(y, dtype=float).ravel()
        solution, *_rest = np.linalg.lstsq(design, target, rcond=None)
        self.coef_ = solution[:-1]
        self.intercept_ = float(solution[-1])
        return self

    def predict(self, x: object) -> FloatArray:
        if self.coef_ is None:
            msg = "model is not fitted"
            raise RuntimeError(msg)
        design = self._design(x)
        return np.asarray(design @ np.concatenate([self.coef_, [self.intercept_]]))

    def score(self, x: object, y: object) -> float:
        return r2_score(np.asarray(y, dtype=float).ravel(), self.predict(x))

    @staticmethod
    def _design(x: object) -> FloatArray:
        arr = np.asarray(x, dtype=float)
        if arr.ndim == 1:
            arr = arr.reshape(-1, 1)
        return np.hstack([arr, np.ones((arr.shape[0], 1))])

LogisticRegression

Binary logistic regression trained by gradient descent on BCE loss.

Source code in src/cds2/ml.py

class LogisticRegression:
    """Binary logistic regression trained by gradient descent on BCE loss."""

    def __init__(
        self,
        learning_rate: float = 0.1,
        max_iter: int = 2000,
        l2: float = 0.0,
    ) -> None:
        self.learning_rate = learning_rate
        self.max_iter = max_iter
        self.l2 = l2
        self.coef_: FloatArray | None = None
        self.intercept_: float = 0.0

    def fit(self, x: object, y: object) -> LogisticRegression:
        if not np.isfinite(self.learning_rate) or self.learning_rate <= 0.0:
            raise ValueError("learning_rate must be a positive finite number")
        if self.max_iter < 1:
            raise ValueError("max_iter must be at least 1")
        if not np.isfinite(self.l2) or self.l2 < 0.0:
            raise ValueError("l2 must be a non-negative finite number")
        features = np.asarray(x, dtype=float)
        if features.ndim == 1:
            features = features.reshape(-1, 1)
        if features.ndim != 2 or features.shape[0] == 0 or features.shape[1] == 0:
            raise ValueError("x must contain a non-empty 2-D feature matrix")
        if not bool(np.all(np.isfinite(features))):
            raise ValueError("x must contain only finite values")
        target = np.asarray(y, dtype=float)
        if target.ndim != 1 or target.shape[0] != features.shape[0]:
            raise ValueError("y must be 1-D and match the number of samples")
        if not bool(np.all(np.isfinite(target))) or not bool(np.all(np.isin(target, [0.0, 1.0]))):
            raise ValueError("y must contain only binary labels 0 and 1")
        n_samples, n_features = features.shape
        weights = np.zeros(n_features)
        bias = 0.0
        for _ in range(self.max_iter):
            logits = features @ weights + bias
            probs = _sigmoid(logits)
            error = probs - target
            grad_w = (features.T @ error) / n_samples + self.l2 * weights
            grad_b = float(np.mean(error))
            weights -= self.learning_rate * grad_w
            bias -= self.learning_rate * grad_b
            if not bool(np.all(np.isfinite(weights))) or not np.isfinite(bias):
                raise FloatingPointError("gradient descent diverged to non-finite parameters")
        self.coef_ = np.asarray(weights, dtype=np.float64)
        self.intercept_ = bias
        return self

    def predict_proba(self, x: object) -> FloatArray:
        if self.coef_ is None:
            raise RuntimeError("model is not fitted")
        features = np.asarray(x, dtype=float)
        if features.ndim == 1:
            features = features.reshape(-1, 1)
        if features.ndim != 2 or features.shape[1] != self.coef_.size:
            raise ValueError("x has a different number of features than the fitted data")
        if not bool(np.all(np.isfinite(features))):
            raise ValueError("x must contain only finite values")
        return _sigmoid(features @ self.coef_ + self.intercept_)

    def predict(self, x: object, threshold: float = 0.5) -> NDArray[np.int64]:
        if not 0.0 <= threshold <= 1.0:
            raise ValueError("threshold must be between 0 and 1")
        return (self.predict_proba(x) >= threshold).astype(np.int64)

    def score(self, x: object, y: object) -> float:
        return accuracy_score(y, self.predict(x))

KMeans

Lloyd's k-means with k-means++ seeding.

Source code in src/cds2/ml.py

class KMeans:
    """Lloyd's k-means with k-means++ seeding."""

    def __init__(
        self, n_clusters: int = 8, max_iter: int = 300, tol: float = 1e-8, seed: int | None = None
    ) -> None:
        self.n_clusters = n_clusters
        self.max_iter = max_iter
        self.tol = tol
        self.seed = seed
        self.cluster_centers_: FloatArray | None = None
        self.labels_: NDArray[np.int64] | None = None
        self.inertia_: float | None = None

    def fit(self, x: object) -> KMeans:
        points = np.asarray(x, dtype=float)
        if points.ndim == 1:
            points = points.reshape(-1, 1)
        if points.ndim != 2 or points.shape[0] == 0 or points.shape[1] == 0:
            raise ValueError("x must contain a non-empty 2-D feature matrix")
        if not bool(np.all(np.isfinite(points))):
            raise ValueError("x must contain only finite values")
        n_samples = points.shape[0]
        if not 1 <= self.n_clusters <= n_samples:
            raise ValueError("n_clusters must be between 1 and number of samples")
        if self.max_iter < 1:
            raise ValueError("max_iter must be at least 1")
        if not np.isfinite(self.tol) or self.tol < 0.0:
            raise ValueError("tol must be a non-negative finite number")
        rng = np.random.default_rng(self.seed)
        centers = self._kmeans_pp_init(points, self.n_clusters, rng)
        if _HAS_C_KERNEL and _c_kernel is not None:
            _labels, centers = self._run_c_lloyd(points, centers)
        else:
            _labels, centers = self._run_numpy_lloyd(points, centers)
        distances = cdist(points, centers, "sqeuclidean")
        labels = distances.argmin(axis=1).astype(np.int64)
        inertia = float(distances[np.arange(n_samples), labels].sum())
        self.cluster_centers_ = np.asarray(centers, dtype=np.float64)
        self.labels_ = labels
        self.inertia_ = inertia
        return self

    def _run_c_lloyd(
        self, points: FloatArray, centers: FloatArray
    ) -> tuple[NDArray[np.int64], FloatArray]:
        points_contiguous = np.ascontiguousarray(points, dtype=np.float64)
        centers_contiguous = np.ascontiguousarray(centers, dtype=np.float64)
        labels_buffer, centers_buffer, _iterations = _c_kernel.lloyd(
            points_contiguous,
            centers_contiguous,
            self.max_iter,
            self.tol,
        )
        labels = np.frombuffer(labels_buffer, dtype=np.int64).copy()
        fitted_centers = (
            np.frombuffer(centers_buffer, dtype=np.float64)
            .copy()
            .reshape(self.n_clusters, points.shape[1])
        )
        return labels, fitted_centers

    def _run_numpy_lloyd(
        self, points: FloatArray, centers: FloatArray
    ) -> tuple[NDArray[np.int64], FloatArray]:
        n_samples = points.shape[0]
        sample_indices = np.arange(n_samples)
        for _ in range(self.max_iter):
            distances = cdist(points, centers, "sqeuclidean")
            new_labels = distances.argmin(axis=1).astype(np.int64)
            counts = np.bincount(new_labels, minlength=self.n_clusters).astype(float)
            assignment = sp_sparse.csr_matrix(
                (np.ones(n_samples), (new_labels, sample_indices)),
                shape=(self.n_clusters, n_samples),
            )
            sums = np.asarray(assignment @ points)
            if (counts == 0).any():
                min_distances = distances[np.arange(n_samples), new_labels].copy()
                for cluster_index in np.flatnonzero(counts == 0):
                    candidates = np.flatnonzero(counts[new_labels] > 1.0)
                    farthest = int(candidates[np.argmax(min_distances[candidates])])
                    old_cluster = int(new_labels[farthest])
                    sums[old_cluster] -= points[farthest]
                    counts[old_cluster] -= 1.0
                    sums[cluster_index] = points[farthest]
                    counts[cluster_index] = 1.0
                    new_labels[farthest] = cluster_index
                    min_distances[farthest] = -1.0
            new_centers = sums / counts[:, None]
            shift = float(np.abs(new_centers - centers).max())
            centers = np.asarray(new_centers, dtype=np.float64)
            if shift < self.tol:
                break
        final_distances = cdist(points, centers, "sqeuclidean")
        labels = final_distances.argmin(axis=1).astype(np.int64)
        return labels, centers

    def predict(self, x: object) -> NDArray[np.int64]:
        if self.cluster_centers_ is None:
            raise RuntimeError("model is not fitted")
        points = np.asarray(x, dtype=float)
        if points.ndim == 1:
            points = points.reshape(1, -1)
        if points.ndim != 2 or points.shape[1] != self.cluster_centers_.shape[1]:
            raise ValueError("x has a different number of features than the fitted data")
        if not bool(np.all(np.isfinite(points))):
            raise ValueError("x must contain only finite values")
        distances = cdist(points, self.cluster_centers_, "sqeuclidean")
        return np.asarray(distances.argmin(axis=1), dtype=np.int64)

    @staticmethod
    def _kmeans_pp_init(
        points: FloatArray, n_clusters: int, rng: np.random.Generator
    ) -> FloatArray:
        n_samples = points.shape[0]
        chosen = [int(rng.integers(n_samples))]
        closest = cdist(points, points[chosen], "sqeuclidean").ravel()
        while len(chosen) < n_clusters:
            total = float(closest.sum())
            if total == 0.0:
                taken = set(chosen)
                candidates = [i for i in range(n_samples) if i not in taken]
                next_index = int(rng.choice(candidates)) if candidates else 0
            else:
                next_index = int(rng.choice(n_samples, p=closest / total))
            chosen.append(next_index)
            closest = np.minimum(
                closest, cdist(points, points[next_index : next_index + 1], "sqeuclidean").ravel()
            )
        return np.asarray(points[chosen], dtype=float)

PCA

Principal component analysis via SVD of the centered data.

Source code in src/cds2/ml.py

class PCA:
    """Principal component analysis via SVD of the centered data."""

    def __init__(self, n_components: int = 2) -> None:
        self.n_components = n_components
        self.components_: FloatArray | None = None
        self.explained_variance_ratio_: FloatArray | None = None
        self.mean_: FloatArray | None = None

    def fit(self, x: object) -> PCA:
        arr = np.asarray(x, dtype=float)
        if arr.ndim != 2 or arr.shape[0] == 0 or arr.shape[1] == 0:
            raise ValueError("PCA expects a non-empty 2-D array of shape (n_samples, n_features)")
        if not bool(np.all(np.isfinite(arr))):
            raise ValueError("x must contain only finite values")
        k = self.n_components
        if not 1 <= k <= min(arr.shape):
            raise ValueError("n_components must be between 1 and min(n_samples, n_features)")
        self.mean_ = np.asarray(arr.mean(axis=0), dtype=np.float64)
        centered = arr - self.mean_
        _u, s, vt = np.linalg.svd(centered, full_matrices=False)
        variance = s**2 / max(arr.shape[0] - 1, 1)
        total_variance = float(variance.sum())
        self.components_ = np.asarray(vt[:k], dtype=np.float64)
        self.explained_variance_ratio_ = (
            np.zeros(k, dtype=np.float64)
            if total_variance == 0.0
            else np.asarray(variance[:k] / total_variance, dtype=np.float64)
        )
        return self

    def transform(self, x: object) -> FloatArray:
        if self.components_ is None or self.mean_ is None:
            raise RuntimeError("model is not fitted")
        arr = np.asarray(x, dtype=float)
        if arr.ndim != 2 or arr.shape[1] != self.mean_.size:
            raise ValueError("x has a different number of features than the fitted data")
        if not bool(np.all(np.isfinite(arr))):
            raise ValueError("x must contain only finite values")
        return np.asarray((arr - self.mean_) @ self.components_.T, dtype=np.float64)

    def fit_transform(self, x: object) -> FloatArray:
        return self.fit(x).transform(x)

KNeighborsClassifier

k-nearest-neighbour classifier backed by a KD-tree.

Source code in src/cds2/ml.py

class KNeighborsClassifier:
    """k-nearest-neighbour classifier backed by a KD-tree."""

    def __init__(self, n_neighbors: int = 5) -> None:
        self.n_neighbors = n_neighbors
        self._tree: cKDTree | None = None
        self._targets: NDArray[np.int64] | None = None
        self._n_features: int | None = None

    def fit(self, x: object, y: object) -> KNeighborsClassifier:
        points = np.asarray(x, dtype=float)
        if points.ndim == 1:
            points = points.reshape(-1, 1)
        if points.ndim != 2 or points.shape[0] == 0 or points.shape[1] == 0:
            raise ValueError("x must contain a non-empty 2-D feature matrix")
        if not bool(np.all(np.isfinite(points))):
            raise ValueError("x must contain only finite values")
        targets = np.asarray(y, dtype=np.int64)
        if targets.ndim != 1 or targets.shape[0] != points.shape[0]:
            raise ValueError("y must be 1-D and match the number of samples")
        if not 1 <= self.n_neighbors <= points.shape[0]:
            raise ValueError("n_neighbors must be between 1 and number of samples")
        self._tree = cKDTree(points)
        self._targets = targets.copy()
        self._n_features = points.shape[1]
        return self

    def predict(self, x: object) -> NDArray[np.int64]:
        if self._tree is None or self._targets is None or self._n_features is None:
            raise RuntimeError("model is not fitted")
        queries = np.asarray(x, dtype=float)
        if queries.ndim == 1:
            queries = queries.reshape(1, -1)
        if queries.ndim != 2 or queries.shape[1] != self._n_features:
            raise ValueError("x has a different number of features than the fitted data")
        if not bool(np.all(np.isfinite(queries))):
            raise ValueError("x must contain only finite values")
        _distances, indices = self._tree.query(queries, k=self.n_neighbors)
        if self.n_neighbors == 1:
            indices = np.asarray(indices, dtype=np.int64).reshape(-1, 1)
        else:
            indices = np.asarray(indices, dtype=np.int64)
        predictions = np.zeros(queries.shape[0], dtype=np.int64)
        for row_index, neighbor_indices in enumerate(indices):
            votes = self._targets[neighbor_indices]
            values, counts = np.unique(votes, return_counts=True)
            predictions[row_index] = values[counts.argmax()]
        return predictions

train_test_split(*arrays, test_size=0.25, shuffle=True, seed=None)

Split arrays into non-empty train/test partitions along axis 0.

Source code in src/cds2/ml.py

def train_test_split(
    *arrays: object,
    test_size: float = 0.25,
    shuffle: bool = True,
    seed: int | None = None,
) -> tuple[NDArray[np.float64], ...]:
    """Split arrays into non-empty train/test partitions along axis 0."""
    if not arrays:
        raise ValueError("at least one array is required")
    prepared = [np.asarray(a) for a in arrays]
    if any(a.ndim == 0 for a in prepared):
        raise ValueError("arrays must have at least one dimension")
    n = prepared[0].shape[0]
    for a in prepared:
        if a.shape[0] != n:
            raise ValueError("all arrays must share the same first dimension")
    if n < 2:
        raise ValueError("at least two samples are required")
    if not 0 < test_size < 1:
        raise ValueError("test_size must be in (0, 1)")
    indices = np.arange(n)
    rng = np.random.default_rng(seed)
    if shuffle:
        rng.shuffle(indices)
    split = int(np.floor(n * (1 - test_size)))
    split = min(max(split, 1), n - 1)
    train_idx, test_idx = indices[:split], indices[split:]
    parts: list[NDArray[np.float64]] = []
    for a in prepared:
        parts.append(a[train_idx])
        parts.append(a[test_idx])
    return tuple(parts)

make_regression_data(n=100, n_features=2, noise=10.0, seed=None)

Synthetic linear-regression dataset y = X w + b + noise.

Source code in src/cds2/ml.py

def make_regression_data(
    n: int = 100,
    n_features: int = 2,
    noise: float = 10.0,
    seed: int | None = None,
) -> tuple[FloatArray, FloatArray]:
    """Synthetic linear-regression dataset y = X w + b + noise."""
    rng = np.random.default_rng(seed)
    x = rng.normal(size=(n, n_features))
    weights = rng.uniform(-4.0, 4.0, size=n_features)
    y = x @ weights + rng.normal(scale=noise, size=n)
    return x, np.asarray(y)

make_blobs(n_samples=300, centers=3, cluster_std=1.0, seed=None)

Gaussian blobs around random centers with integer labels.

Source code in src/cds2/ml.py

def make_blobs(
    n_samples: int = 300,
    centers: int = 3,
    cluster_std: float = 1.0,
    seed: int | None = None,
) -> tuple[FloatArray, NDArray[np.int64]]:
    """Gaussian blobs around random centers with integer labels."""
    rng = np.random.default_rng(seed)
    center_points = rng.uniform(-10.0, 10.0, size=(centers, 2))
    counts = np.full(centers, n_samples // centers, dtype=int)
    counts[: n_samples % centers] += 1
    xs: list[FloatArray] = []
    ys: list[NDArray[np.int64]] = []
    for label, (count, center) in enumerate(zip(counts, center_points, strict=True)):
        xs.append(rng.normal(loc=center, scale=cluster_std, size=(count, 2)))
        ys.append(np.full(count, label, dtype=np.int64))
    return np.vstack(xs), np.concatenate(ys)