+
covariance
+
EmpiricalCovariance, ShrunkCovariance, LedoitWolf, OAS
+
✅ Implemented
+
+
preprocessing (extended)
+
PowerTransformer, QuantileTransformer, Binarizer, FunctionTransformer
+
✅ Implemented
+
diff --git a/src/calibration/calibration.ts b/src/calibration/calibration.ts
new file mode 100644
index 0000000..948aa5f
--- /dev/null
+++ b/src/calibration/calibration.ts
@@ -0,0 +1,141 @@
+/**
+ * Probability calibration.
+ * Mirrors sklearn.calibration.CalibratedClassifierCV.
+ * Uses Platt scaling (logistic) or isotonic regression for calibration.
+ */
+
+import { NotFittedError } from "../exceptions.js";
+
+interface Classifier {
+ fit(X: Float64Array[], y: Float64Array): this;
+ predict(X: Float64Array[]): Float64Array;
+ score?(X: Float64Array[], y: Float64Array): number;
+}
+
+function sigmoid(x: number): number {
+ return 1 / (1 + Math.exp(-x));
+}
+
+/** Platt scaling: fit a logistic function on scores to map to probabilities. */
+function plattScale(scores: Float64Array, y: Float64Array): [number, number] {
+ const n = scores.length;
+ let A = 0;
+ let B = 0;
+ const lr = 0.01;
+
+ for (let iter = 0; iter < 1000; iter++) {
+ let gradA = 0;
+ let gradB = 0;
+ for (let i = 0; i < n; i++) {
+ const p = sigmoid(A * (scores[i] ?? 0) + B);
+ const err = p - (y[i] ?? 0);
+ gradA += err * (scores[i] ?? 0);
+ gradB += err;
+ }
+ A -= lr * gradA / n;
+ B -= lr * gradB / n;
+ }
+
+ return [A, B];
+}
+
+export class CalibratedClassifierCV {
+ baseEstimator: Classifier;
+ method: string;
+ cv: number;
+
+ calibratedEstimators_: {
+ estimator: Classifier;
+ A: number;
+ B: number;
+ }[] | null = null;
+ classes_: Float64Array | null = null;
+
+ constructor(
+ baseEstimator: Classifier,
+ options: { method?: string; cv?: number } = {},
+ ) {
+ this.baseEstimator = baseEstimator;
+ this.method = options.method ?? "sigmoid";
+ this.cv = options.cv ?? 5;
+ }
+
+ fit(X: Float64Array[], y: Float64Array): this {
+ const n = X.length;
+ const uniqueClasses = Array.from(new Set(Array.from(y))).sort((a, b) => a - b);
+ this.classes_ = new Float64Array(uniqueClasses);
+ const posClass = uniqueClasses[uniqueClasses.length - 1] ?? 1;
+
+ const yBin = new Float64Array(y.map((yi) => (yi === posClass ? 1 : 0)));
+
+ // Simple hold-out calibration
+ const foldSize = Math.floor(n / this.cv);
+ this.calibratedEstimators_ = [];
+
+ for (let fold = 0; fold < this.cv; fold++) {
+ const testStart = fold * foldSize;
+ const testEnd = fold === this.cv - 1 ? n : testStart + foldSize;
+
+ const trainIdx: number[] = [];
+ const testIdx: number[] = [];
+ for (let i = 0; i < n; i++) {
+ if (i >= testStart && i < testEnd) testIdx.push(i);
+ else trainIdx.push(i);
+ }
+
+ const XTrain = trainIdx.map((i) => X[i] ?? new Float64Array(0));
+ const yTrain = new Float64Array(trainIdx.map((i) => y[i] ?? 0));
+ const XTest = testIdx.map((i) => X[i] ?? new Float64Array(0));
+ const yTest = new Float64Array(testIdx.map((i) => yBin[i] ?? 0));
+
+ const est = Object.create(Object.getPrototypeOf(this.baseEstimator) as object) as Classifier;
+ Object.assign(est, this.baseEstimator);
+ est.fit(XTrain, yTrain);
+
+ const testPred = est.predict(XTest);
+ const [A, B] = plattScale(testPred, yTest);
+
+ this.calibratedEstimators_.push({ estimator: est, A, B });
+ }
+
+ return this;
+ }
+
+ predictProba(X: Float64Array[]): Float64Array[] {
+ if (this.calibratedEstimators_ === null) throw new NotFittedError("CalibratedClassifierCV");
+
+ const n = X.length;
+ const probs = new Float64Array(n);
+
+ for (const { estimator, A, B } of this.calibratedEstimators_) {
+ const scores = estimator.predict(X);
+ for (let i = 0; i < n; i++) {
+ probs[i] = (probs[i] ?? 0) + sigmoid(A * (scores[i] ?? 0) + B);
+ }
+ }
+
+ const k = this.calibratedEstimators_.length;
+ return Array.from({ length: n }, (_, i) => {
+ const p = (probs[i] ?? 0) / k;
+ return new Float64Array([1 - p, p]);
+ });
+ }
+
+ predict(X: Float64Array[]): Float64Array {
+ if (this.classes_ === null) throw new NotFittedError("CalibratedClassifierCV");
+ const classes = this.classes_;
+ const proba = this.predictProba(X);
+ const posClass = classes[classes.length - 1] ?? 1;
+ const negClass = classes[0] ?? 0;
+ return new Float64Array(proba.map((p) => ((p[1] ?? 0) >= 0.5 ? posClass : negClass)));
+ }
+
+ score(X: Float64Array[], y: Float64Array): number {
+ const pred = this.predict(X);
+ let correct = 0;
+ for (let i = 0; i < y.length; i++) {
+ if (pred[i] === y[i]) correct++;
+ }
+ return correct / y.length;
+ }
+}
diff --git a/src/calibration/index.ts b/src/calibration/index.ts
new file mode 100644
index 0000000..e03c3f7
--- /dev/null
+++ b/src/calibration/index.ts
@@ -0,0 +1 @@
+export * from "./calibration.js";
diff --git a/src/cluster/affinity_propagation.ts b/src/cluster/affinity_propagation.ts
new file mode 100644
index 0000000..1228a23
--- /dev/null
+++ b/src/cluster/affinity_propagation.ts
@@ -0,0 +1,199 @@
+/**
+ * AffinityPropagation clustering.
+ */
+
+import { NotFittedError } from "../exceptions.js";
+
+export interface AffinityPropagationOptions {
+ dampingFactor?: number;
+ maxIter?: number;
+ convergenceIter?: number;
+ preference?: number;
+}
+
+export class AffinityPropagation {
+ private dampingFactor: number;
+ private maxIter: number;
+ private convergenceIter: number;
+ private preference: number | undefined;
+
+ labels_: Int32Array | null = null;
+ clusterCentersIndices_: Int32Array | null = null;
+ nIter_ = 0;
+
+ constructor(options: AffinityPropagationOptions = {}) {
+ this.dampingFactor = options.dampingFactor ?? 0.5;
+ this.maxIter = options.maxIter ?? 200;
+ this.convergenceIter = options.convergenceIter ?? 15;
+ this.preference = options.preference;
+ }
+
+ fit(X: Float64Array[]): this {
+ const n = X.length;
+ if (n === 0) {
+ this.labels_ = new Int32Array(0);
+ this.clusterCentersIndices_ = new Int32Array(0);
+ return this;
+ }
+
+ // Build similarity matrix S = -||xi - xj||^2
+ const S: Float64Array[] = Array.from(
+ { length: n },
+ () => new Float64Array(n),
+ );
+ for (let i = 0; i < n; i++) {
+ const xi = X[i] ?? new Float64Array(0);
+ for (let j = i; j < n; j++) {
+ const xj = X[j] ?? new Float64Array(0);
+ let d = 0;
+ for (let k = 0; k < xi.length; k++)
+ d += ((xi[k] ?? 0) - (xj[k] ?? 0)) ** 2;
+ (S[i] as Float64Array)[j] = -d;
+ (S[j] as Float64Array)[i] = -d;
+ }
+ }
+
+ // Set preference (diagonal)
+ let pref = this.preference;
+ if (pref === undefined) {
+ // Median of similarities
+ const vals: number[] = [];
+ for (let i = 0; i < n; i++)
+ for (let j = i + 1; j < n; j++)
+ vals.push((S[i] as Float64Array)[j] ?? 0);
+ vals.sort((a, b) => a - b);
+ pref = vals[Math.floor(vals.length / 2)] ?? -1;
+ }
+ for (let i = 0; i < n; i++) (S[i] as Float64Array)[i] = pref;
+
+ // Responsibility R and Availability A matrices
+ const R: Float64Array[] = Array.from(
+ { length: n },
+ () => new Float64Array(n),
+ );
+ const A: Float64Array[] = Array.from(
+ { length: n },
+ () => new Float64Array(n),
+ );
+ const d = this.dampingFactor;
+ let stableCount = 0;
+ let prevExemplars: Set
= new Set();
+
+ for (let iter = 0; iter < this.maxIter; iter++) {
+ // Update responsibilities: R(i,k) = S(i,k) - max_{k'!=k}[A(i,k')+S(i,k')]
+ for (let i = 0; i < n; i++) {
+ const Si = S[i] ?? new Float64Array(n);
+ const Ai = A[i] ?? new Float64Array(n);
+ // Find two highest A+S values
+ let max1 = Number.NEGATIVE_INFINITY;
+ let max2 = Number.NEGATIVE_INFINITY;
+ let argmax1 = -1;
+ for (let k = 0; k < n; k++) {
+ const v = (Ai[k] ?? 0) + (Si[k] ?? 0);
+ if (v > max1) {
+ max2 = max1;
+ max1 = v;
+ argmax1 = k;
+ } else if (v > max2) max2 = v;
+ }
+ const Ri = R[i] ?? new Float64Array(n);
+ for (let k = 0; k < n; k++) {
+ const maxOther = k === argmax1 ? max2 : max1;
+ const newR = (Si[k] ?? 0) - maxOther;
+ Ri[k] = d * (Ri[k] ?? 0) + (1 - d) * newR;
+ }
+ }
+
+ // Update availabilities
+ for (let k = 0; k < n; k++) {
+ // sum of positive R(i',k) for i'!=k
+ let sumPos = 0;
+ for (let i = 0; i < n; i++) {
+ if (i === k) continue;
+ const v = (R[i] as Float64Array)[k] ?? 0;
+ if (v > 0) sumPos += v;
+ }
+ const rkk = (R[k] as Float64Array)[k] ?? 0;
+ for (let i = 0; i < n; i++) {
+ const Ai = A[i] ?? new Float64Array(n);
+ let newA: number;
+ if (i === k) {
+ newA = sumPos;
+ } else {
+ const rik = (R[i] as Float64Array)[k] ?? 0;
+ const sumWithout = sumPos - (rik > 0 ? rik : 0);
+ newA = Math.min(0, rkk + sumWithout);
+ }
+ Ai[k] = d * (Ai[k] ?? 0) + (1 - d) * newA;
+ }
+ }
+
+ // Check convergence
+ const exemplars = new Set();
+ for (let i = 0; i < n; i++) {
+ const Ai = A[i] ?? new Float64Array(n);
+ const Ri = R[i] ?? new Float64Array(n);
+ let best = Number.NEGATIVE_INFINITY;
+ let bestK = 0;
+ for (let k = 0; k < n; k++) {
+ const v = (Ai[k] ?? 0) + (Ri[k] ?? 0);
+ if (v > best) {
+ best = v;
+ bestK = k;
+ }
+ }
+ exemplars.add(bestK);
+ }
+
+ const same =
+ exemplars.size === prevExemplars.size &&
+ [...exemplars].every((e) => prevExemplars.has(e));
+ if (same) {
+ stableCount++;
+ if (stableCount >= this.convergenceIter) {
+ this.nIter_ = iter + 1;
+ break;
+ }
+ } else {
+ stableCount = 0;
+ }
+ prevExemplars = exemplars;
+ this.nIter_ = iter + 1;
+ }
+
+ // Assign labels
+ const labels = new Int32Array(n);
+ for (let i = 0; i < n; i++) {
+ const Ai = A[i] ?? new Float64Array(n);
+ const Ri = R[i] ?? new Float64Array(n);
+ let best = Number.NEGATIVE_INFINITY;
+ let bestK = 0;
+ for (let k = 0; k < n; k++) {
+ const v = (Ai[k] ?? 0) + (Ri[k] ?? 0);
+ if (v > best) {
+ best = v;
+ bestK = k;
+ }
+ }
+ labels[i] = bestK;
+ }
+
+ const centerSet = new Set(Array.from(labels));
+ const centers = Int32Array.from([...centerSet].sort((a, b) => a - b));
+ // Relabel to 0..k-1
+ const map = new Map();
+ centers.forEach((c, idx) => map.set(c, idx));
+ for (let i = 0; i < n; i++) labels[i] = map.get(labels[i] ?? 0) ?? 0;
+
+ this.labels_ = labels;
+ this.clusterCentersIndices_ = centers;
+ return this;
+ }
+
+ predict(X: Float64Array[]): Int32Array {
+ if (!this.labels_ || !this.clusterCentersIndices_)
+ throw new NotFittedError("AffinityPropagation");
+ // Not supported post-fit without stored data; return empty
+ return new Int32Array(X.length).fill(-1);
+ }
+}
diff --git a/src/cluster/agglomerative.ts b/src/cluster/agglomerative.ts
new file mode 100644
index 0000000..68eddcf
--- /dev/null
+++ b/src/cluster/agglomerative.ts
@@ -0,0 +1,198 @@
+/**
+ * AgglomerativeClustering and MiniBatchKMeans.
+ * Mirrors sklearn.cluster.AgglomerativeClustering and MiniBatchKMeans.
+ */
+
+import { NotFittedError } from "../exceptions.js";
+
+function euclidean(a: Float64Array, b: Float64Array): number {
+ let s = 0;
+ for (let i = 0; i < a.length; i++) s += ((a[i] ?? 0) - (b[i] ?? 0)) ** 2;
+ return Math.sqrt(s);
+}
+
+export type Linkage = "ward" | "complete" | "average" | "single";
+
+export interface AgglomerativeClusteringOptions {
+ nClusters?: number;
+ linkage?: Linkage;
+}
+
+export class AgglomerativeClustering {
+ nClusters: number;
+ linkage: Linkage;
+
+ labels_: Int32Array | null = null;
+ nClusters_: number = 0;
+
+ constructor(options: AgglomerativeClusteringOptions = {}) {
+ this.nClusters = options.nClusters ?? 2;
+ this.linkage = options.linkage ?? "ward";
+ }
+
+ fit(X: Float64Array[]): this {
+ const n = X.length;
+ // Initialize each point as its own cluster
+ let clusters: number[][] = X.map((_, i) => [i]);
+
+ // Distance matrix
+ const dist = (a: number[], b: number[]): number => {
+ if (this.linkage === "single") {
+ let min = Number.POSITIVE_INFINITY;
+ for (const i of a)
+ for (const j of b) min = Math.min(min, euclidean(X[i]!, X[j]!));
+ return min;
+ } else if (this.linkage === "complete") {
+ let max = Number.NEGATIVE_INFINITY;
+ for (const i of a)
+ for (const j of b) max = Math.max(max, euclidean(X[i]!, X[j]!));
+ return max;
+ } else {
+ // average and ward both use average distance here (simplified)
+ let sum = 0;
+ for (const i of a) for (const j of b) sum += euclidean(X[i]!, X[j]!);
+ return sum / (a.length * b.length);
+ }
+ };
+
+ while (clusters.length > this.nClusters) {
+ let minD = Number.POSITIVE_INFINITY;
+ let mergeI = 0;
+ let mergeJ = 1;
+ for (let i = 0; i < clusters.length; i++) {
+ for (let j = i + 1; j < clusters.length; j++) {
+ const d = dist(clusters[i]!, clusters[j]!);
+ if (d < minD) {
+ minD = d;
+ mergeI = i;
+ mergeJ = j;
+ }
+ }
+ }
+ clusters[mergeI] = clusters[mergeI]!.concat(clusters[mergeJ]!);
+ clusters.splice(mergeJ, 1);
+ }
+
+ this.labels_ = new Int32Array(n);
+ for (let k = 0; k < clusters.length; k++) {
+ for (const idx of clusters[k]!) this.labels_[idx] = k;
+ }
+ this.nClusters_ = clusters.length;
+ return this;
+ }
+
+ fitPredict(X: Float64Array[]): Int32Array {
+ this.fit(X);
+ return this.labels_!;
+ }
+}
+
+export interface MiniBatchKMeansOptions {
+ nClusters?: number;
+ batchSize?: number;
+ maxIter?: number;
+ tol?: number;
+}
+
+export class MiniBatchKMeans {
+ nClusters: number;
+ batchSize: number;
+ maxIter: number;
+ tol: number;
+
+ clusterCenters_: Float64Array[] | null = null;
+ labels_: Int32Array | null = null;
+ inertia_: number = 0;
+
+ constructor(options: MiniBatchKMeansOptions = {}) {
+ this.nClusters = options.nClusters ?? 8;
+ this.batchSize = options.batchSize ?? 100;
+ this.maxIter = options.maxIter ?? 100;
+ this.tol = options.tol ?? 1e-4;
+ }
+
+ private _initCenters(X: Float64Array[]): Float64Array[] {
+ const indices: number[] = [];
+ while (indices.length < this.nClusters) {
+ const idx = Math.floor(Math.random() * X.length);
+ if (!indices.includes(idx)) indices.push(idx);
+ }
+ return indices.map((i) => new Float64Array(X[i]!));
+ }
+
+ fit(X: Float64Array[]): this {
+ const n = X.length;
+ if (n === 0) throw new Error("Empty input");
+ const nFeatures = X[0]?.length ?? 0;
+
+ const centers = this._initCenters(X);
+ const counts = new Float64Array(this.nClusters);
+
+ for (let iter = 0; iter < this.maxIter; iter++) {
+ const batch: Float64Array[] = [];
+ for (let i = 0; i < this.batchSize; i++) {
+ batch.push(X[Math.floor(Math.random() * n)]!);
+ }
+
+ for (const x of batch) {
+ let nearest = 0;
+ let minD = Number.POSITIVE_INFINITY;
+ for (let k = 0; k < this.nClusters; k++) {
+ const d = euclidean(x, centers[k]!);
+ if (d < minD) {
+ minD = d;
+ nearest = k;
+ }
+ }
+ counts[nearest] = (counts[nearest] ?? 0) + 1;
+ const lr = 1 / (counts[nearest] ?? 1);
+ const c = centers[nearest]!;
+ for (let j = 0; j < nFeatures; j++) {
+ c[j] = (c[j] ?? 0) * (1 - lr) + (x[j] ?? 0) * lr;
+ }
+ }
+ }
+
+ this.clusterCenters_ = centers;
+ this.labels_ = new Int32Array(n);
+ this.inertia_ = 0;
+
+ for (let i = 0; i < n; i++) {
+ let nearest = 0;
+ let minD = Number.POSITIVE_INFINITY;
+ for (let k = 0; k < this.nClusters; k++) {
+ const d = euclidean(X[i]!, centers[k]!);
+ if (d < minD) {
+ minD = d;
+ nearest = k;
+ }
+ }
+ this.labels_[i] = nearest;
+ this.inertia_ += minD * minD;
+ }
+ return this;
+ }
+
+ predict(X: Float64Array[]): Int32Array {
+ if (!this.clusterCenters_) throw new NotFittedError("MiniBatchKMeans");
+ const out = new Int32Array(X.length);
+ for (let i = 0; i < X.length; i++) {
+ let nearest = 0;
+ let minD = Number.POSITIVE_INFINITY;
+ for (let k = 0; k < this.nClusters; k++) {
+ const d = euclidean(X[i]!, this.clusterCenters_[k]!);
+ if (d < minD) {
+ minD = d;
+ nearest = k;
+ }
+ }
+ out[i] = nearest;
+ }
+ return out;
+ }
+
+ fitPredict(X: Float64Array[]): Int32Array {
+ this.fit(X);
+ return this.labels_!;
+ }
+}
diff --git a/src/cluster/bisecting_kmeans.ts b/src/cluster/bisecting_kmeans.ts
new file mode 100644
index 0000000..bc4e6d5
--- /dev/null
+++ b/src/cluster/bisecting_kmeans.ts
@@ -0,0 +1,204 @@
+/**
+ * BisectingKMeans: divisive hierarchical clustering using k-means bisection.
+ * Mirrors sklearn.cluster.BisectingKMeans.
+ */
+
+import { NotFittedError } from "../exceptions.js";
+
+function euclidean(a: Float64Array, b: Float64Array): number {
+ let s = 0;
+ for (let i = 0; i < a.length; i++) s += ((a[i] ?? 0) - (b[i] ?? 0)) ** 2;
+ return Math.sqrt(s);
+}
+
+function clusterMean(points: Float64Array[]): Float64Array {
+ if (points.length === 0) return new Float64Array(0);
+ const p = (points[0] ?? new Float64Array(0)).length;
+ const m = new Float64Array(p);
+ for (const pt of points) for (let j = 0; j < p; j++) m[j] = (m[j] ?? 0) + (pt[j] ?? 0);
+ for (let j = 0; j < p; j++) m[j] = (m[j] ?? 0) / points.length;
+ return m;
+}
+
+function clusterSSE(points: Float64Array[], center: Float64Array): number {
+ let s = 0;
+ for (const pt of points) {
+ for (let j = 0; j < pt.length; j++) s += ((pt[j] ?? 0) - (center[j] ?? 0)) ** 2;
+ }
+ return s;
+}
+
+/** Run k-means with k=2 on the given points. Returns cluster assignments. */
+function bisect(
+ points: Float64Array[],
+ maxIter: number,
+ rng: number,
+): { labels: Int32Array; centers: Float64Array[] } {
+ const n = points.length;
+ const p = (points[0] ?? new Float64Array(0)).length;
+
+ if (n <= 1) {
+ return { labels: new Int32Array(n), centers: [clusterMean(points), new Float64Array(p)] };
+ }
+
+ // Init: pick 2 random centers
+ const i0 = Math.abs(rng) % n;
+ const i1 = (Math.abs(rng) + 1) % n;
+ let centers = [new Float64Array(points[i0] ?? new Float64Array(p)), new Float64Array(points[i1] ?? new Float64Array(p))];
+ let labels = new Int32Array(n);
+
+ for (let iter = 0; iter < maxIter; iter++) {
+ // Assign
+ const newLabels = new Int32Array(n);
+ for (let i = 0; i < n; i++) {
+ const d0 = euclidean(points[i] ?? new Float64Array(p), centers[0] ?? new Float64Array(p));
+ const d1 = euclidean(points[i] ?? new Float64Array(p), centers[1] ?? new Float64Array(p));
+ newLabels[i] = d1 < d0 ? 1 : 0;
+ }
+
+ // Update centers
+ const c0 = points.filter((_, i) => newLabels[i] === 0);
+ const c1 = points.filter((_, i) => newLabels[i] === 1);
+ const newCenters = [
+ c0.length > 0 ? clusterMean(c0) : centers[0] ?? new Float64Array(p),
+ c1.length > 0 ? clusterMean(c1) : centers[1] ?? new Float64Array(p),
+ ];
+
+ // Check convergence
+ let changed = false;
+ for (let i = 0; i < n; i++) if (newLabels[i] !== labels[i]) { changed = true; break; }
+ labels = newLabels;
+ centers = newCenters;
+ if (!changed) break;
+ }
+
+ return { labels, centers: [centers[0] ?? new Float64Array(p), centers[1] ?? new Float64Array(p)] };
+}
+
+/**
+ * BisectingKMeans: hierarchical divisive clustering.
+ * Repeatedly bisects the cluster with highest SSE.
+ * Mirrors sklearn.cluster.BisectingKMeans.
+ */
+export class BisectingKMeans {
+ nClusters: number;
+ maxIter: number;
+ randomState: number;
+ bisectingStrategy: "biggest_inertia" | "largest_cluster";
+
+ clusterCenters_: Float64Array[] | null = null;
+ labels_: Int32Array | null = null;
+ inertia_: number = 0;
+ nIter_: number = 0;
+
+ constructor(
+ options: {
+ nClusters?: number;
+ maxIter?: number;
+ randomState?: number;
+ bisectingStrategy?: "biggest_inertia" | "largest_cluster";
+ } = {},
+ ) {
+ this.nClusters = options.nClusters ?? 8;
+ this.maxIter = options.maxIter ?? 300;
+ this.randomState = options.randomState ?? 42;
+ this.bisectingStrategy = options.bisectingStrategy ?? "biggest_inertia";
+ }
+
+ fit(X: Float64Array[]): this {
+ const n = X.length;
+ const p = (X[0] ?? new Float64Array(0)).length;
+ const k = Math.min(this.nClusters, n);
+
+ // Start: all points in one cluster
+ let clusterLabels = new Int32Array(n);
+ const clusterCenters: Float64Array[] = [clusterMean(X)];
+ let nClusters = 1;
+
+ let rng = this.randomState;
+
+ while (nClusters < k) {
+ // Find cluster to bisect
+ let targetCluster = 0;
+ let bestCrit = -Number.POSITIVE_INFINITY;
+
+ for (let c = 0; c < nClusters; c++) {
+ const pts = X.filter((_, i) => clusterLabels[i] === c);
+ if (pts.length <= 1) continue;
+ const crit = this.bisectingStrategy === "biggest_inertia"
+ ? clusterSSE(pts, clusterCenters[c] ?? new Float64Array(p))
+ : pts.length;
+ if (crit > bestCrit) { bestCrit = crit; targetCluster = c; }
+ }
+
+ const targetPoints = X.filter((_, i) => clusterLabels[i] === targetCluster);
+ const targetIndices = Array.from({ length: n }, (_, i) => i).filter((i) => clusterLabels[i] === targetCluster);
+
+ if (targetPoints.length <= 1) break;
+
+ rng = Math.abs(rng * 1664525 + 1013904223) % 2147483647;
+ const { labels: subLabels } = bisect(targetPoints, this.maxIter, rng);
+
+ // Update global labels: targetCluster stays for subLabel=0, nClusters for subLabel=1
+ for (let i = 0; i < targetIndices.length; i++) {
+ const idx = targetIndices[i] ?? 0;
+ if ((subLabels[i] ?? 0) === 1) clusterLabels[idx] = nClusters;
+ }
+
+ // Recompute centers for the two new clusters
+ const c0pts = X.filter((_, i) => clusterLabels[i] === targetCluster);
+ const c1pts = X.filter((_, i) => clusterLabels[i] === nClusters);
+ clusterCenters[targetCluster] = c0pts.length > 0 ? clusterMean(c0pts) : new Float64Array(p);
+ clusterCenters.push(c1pts.length > 0 ? clusterMean(c1pts) : new Float64Array(p));
+ nClusters++;
+ this.nIter_++;
+ }
+
+ this.labels_ = clusterLabels;
+ this.clusterCenters_ = clusterCenters;
+
+ // Compute inertia
+ let inertia = 0;
+ for (let i = 0; i < n; i++) {
+ const c = clusterLabels[i] ?? 0;
+ const center = clusterCenters[c] ?? new Float64Array(p);
+ const xi = X[i] ?? new Float64Array(p);
+ for (let j = 0; j < p; j++) inertia += ((xi[j] ?? 0) - (center[j] ?? 0)) ** 2;
+ }
+ this.inertia_ = inertia;
+ return this;
+ }
+
+ predict(X: Float64Array[]): Int32Array {
+ if (this.clusterCenters_ === null) throw new NotFittedError("BisectingKMeans");
+ const centers = this.clusterCenters_;
+ return new Int32Array(X.map((xi) => {
+ let bestC = 0;
+ let bestD = Number.POSITIVE_INFINITY;
+ for (let c = 0; c < centers.length; c++) {
+ const d = euclidean(xi, centers[c] ?? new Float64Array(0));
+ if (d < bestD) { bestD = d; bestC = c; }
+ }
+ return bestC;
+ }));
+ }
+
+ fitPredict(X: Float64Array[]): Int32Array {
+ this.fit(X);
+ return this.labels_!;
+ }
+
+ score(X: Float64Array[]): number {
+ if (this.clusterCenters_ === null) throw new NotFittedError("BisectingKMeans");
+ const labels = this.predict(X);
+ const centers = this.clusterCenters_;
+ let inertia = 0;
+ for (let i = 0; i < X.length; i++) {
+ const c = labels[i] ?? 0;
+ const center = centers[c] ?? new Float64Array(0);
+ const xi = X[i] ?? new Float64Array(0);
+ for (let j = 0; j < xi.length; j++) inertia += ((xi[j] ?? 0) - (center[j] ?? 0)) ** 2;
+ }
+ return -inertia;
+ }
+}
diff --git a/src/cluster/clustering_utils.ts b/src/cluster/clustering_utils.ts
new file mode 100644
index 0000000..2b8ef2e
--- /dev/null
+++ b/src/cluster/clustering_utils.ts
@@ -0,0 +1,295 @@
+/**
+ * Cluster utility functions.
+ * Mirrors sklearn.cluster._mean_shift and related utilities.
+ */
+
+/**
+ * Estimate the bandwidth for Mean Shift algorithm.
+ * Uses a ball-tree-like approach: for each sample, counts how many
+ * samples are within the estimated bandwidth.
+ *
+ * @param X - Input data (n_samples x n_features)
+ * @param quantile - Quantile of pairwise distances to use as bandwidth (default 0.3)
+ * @param nSamples - Number of samples to use for estimation (default: all)
+ * @param seed - Random seed for subsampling
+ */
+export function estimateBandwidth(
+ X: Float64Array[],
+ options: {
+ quantile?: number;
+ nSamples?: number;
+ seed?: number;
+ } = {},
+): number {
+ const { quantile = 0.3, seed = 0 } = options;
+ const n = X.length;
+ let nSamples = options.nSamples ?? n;
+ nSamples = Math.min(nSamples, n);
+
+ // Subsample if needed
+ let indices: number[];
+ if (nSamples < n) {
+ let rng = seed;
+ const rand = () => {
+ rng = (rng * 1664525 + 1013904223) & 0xffffffff;
+ return (rng >>> 0) / 0xffffffff;
+ };
+ indices = Array.from({ length: n }, (_, i) => i);
+ for (let i = n - 1; i > 0; i--) {
+ const j = Math.floor(rand() * (i + 1));
+ const tmp = indices[i]!; indices[i] = indices[j]!; indices[j] = tmp;
+ }
+ indices = indices.slice(0, nSamples);
+ } else {
+ indices = Array.from({ length: n }, (_, i) => i);
+ }
+
+ // Compute pairwise distances between sampled points and all points
+ // Then take the quantile
+ const allDists: number[] = [];
+ for (const idx of indices) {
+ const xi = X[idx]!;
+ for (let j = 0; j < n; j++) {
+ const xj = X[j]!;
+ let d2 = 0;
+ for (let k = 0; k < xi.length; k++) {
+ d2 += ((xi[k] ?? 0) - (xj[k] ?? 0)) ** 2;
+ }
+ allDists.push(Math.sqrt(d2));
+ }
+ }
+
+ allDists.sort((a, b) => a - b);
+ const qIdx = Math.floor(quantile * (allDists.length - 1));
+ return allDists[qIdx] ?? 1.0;
+}
+
+/**
+ * Find initial seed points for Mean Shift.
+ * Seeds are bin centers of a uniform grid at bandwidth resolution.
+ *
+ * @param X - Input data
+ * @param bandwidth - Bin size
+ * @param minBinFreq - Minimum number of points per bin to be included
+ */
+export function getBinSeeds(
+ X: Float64Array[],
+ bandwidth: number,
+ minBinFreq = 1,
+): Float64Array[] {
+ if (bandwidth <= 0) throw new Error("bandwidth must be positive");
+ const n = X.length;
+ const d = X[0]?.length ?? 0;
+
+ // Discretize X into bins
+ const binMap = new Map();
+
+ for (let i = 0; i < n; i++) {
+ const xi = X[i]!;
+ const binCoords: number[] = [];
+ for (let k = 0; k < d; k++) {
+ binCoords.push(Math.round((xi[k] ?? 0) / bandwidth));
+ }
+ const key = binCoords.join(",");
+ const existing = binMap.get(key);
+ if (existing) {
+ for (let k = 0; k < d; k++) {
+ existing.sum[k]! += xi[k] ?? 0;
+ }
+ existing.count++;
+ } else {
+ const sum = new Float64Array(d);
+ for (let k = 0; k < d; k++) sum[k] = xi[k] ?? 0;
+ binMap.set(key, { sum, count: 1 });
+ }
+ }
+
+ // Return bin centers with sufficient frequency
+ const seeds: Float64Array[] = [];
+ for (const { sum, count } of binMap.values()) {
+ if (count >= minBinFreq) {
+ const center = new Float64Array(d);
+ for (let k = 0; k < d; k++) center[k] = (sum[k] ?? 0) / count;
+ seeds.push(center);
+ }
+ }
+
+ return seeds;
+}
+
+/**
+ * Find which bin each point belongs to.
+ * @returns Int32Array of bin indices (one per sample)
+ */
+export function assignBins(
+ X: Float64Array[],
+ seeds: Float64Array[],
+): Int32Array {
+ const n = X.length;
+ const result = new Int32Array(n).fill(-1);
+ for (let i = 0; i < n; i++) {
+ const xi = X[i]!;
+ let bestDist = Number.POSITIVE_INFINITY;
+ let bestJ = -1;
+ for (let j = 0; j < seeds.length; j++) {
+ const seed = seeds[j]!;
+ let d2 = 0;
+ for (let k = 0; k < xi.length; k++) {
+ d2 += ((xi[k] ?? 0) - (seed[k] ?? 0)) ** 2;
+ }
+ if (d2 < bestDist) { bestDist = d2; bestJ = j; }
+ }
+ result[i] = bestJ;
+ }
+ return result;
+}
+
+/**
+ * Single iteration of mean-shift update for a set of seeds.
+ * Updates each seed to the mean of all points within bandwidth distance.
+ *
+ * @returns New seed positions and whether any seed moved more than tol
+ */
+export function meanShiftStep(
+ X: Float64Array[],
+ seeds: Float64Array[],
+ bandwidth: number,
+): { newSeeds: Float64Array[]; converged: boolean } {
+ const d = X[0]?.length ?? 0;
+ const bw2 = bandwidth * bandwidth;
+ const newSeeds: Float64Array[] = [];
+ let maxShift = 0;
+
+ for (const seed of seeds) {
+ const newSeed = new Float64Array(d);
+ let weight = 0;
+ for (const xi of X) {
+ let d2 = 0;
+ for (let k = 0; k < d; k++) {
+ d2 += ((xi[k] ?? 0) - (seed[k] ?? 0)) ** 2;
+ }
+ if (d2 <= bw2) {
+ weight++;
+ for (let k = 0; k < d; k++) newSeed[k]! += xi[k] ?? 0;
+ }
+ }
+ if (weight > 0) {
+ for (let k = 0; k < d; k++) newSeed[k]! /= weight;
+ } else {
+ newSeed.set(seed);
+ }
+
+ // Track max shift
+ let shift2 = 0;
+ for (let k = 0; k < d; k++) {
+ shift2 += ((newSeed[k] ?? 0) - (seed[k] ?? 0)) ** 2;
+ }
+ maxShift = Math.max(maxShift, Math.sqrt(shift2));
+ newSeeds.push(newSeed);
+ }
+
+ return { newSeeds, converged: maxShift < 1e-3 * bandwidth };
+}
+
+/**
+ * Merge nearby seeds by deduplication within bandwidth distance.
+ * Returns unique cluster centers.
+ */
+export function mergeSeeds(
+ seeds: Float64Array[],
+ bandwidth: number,
+): Float64Array[] {
+ const bw2 = bandwidth * bandwidth;
+ const merged: Float64Array[] = [];
+
+ for (const seed of seeds) {
+ let isNew = true;
+ for (const center of merged) {
+ let d2 = 0;
+ for (let k = 0; k < seed.length; k++) {
+ d2 += ((seed[k] ?? 0) - (center[k] ?? 0)) ** 2;
+ }
+ if (d2 <= bw2) { isNew = false; break; }
+ }
+ if (isNew) merged.push(seed);
+ }
+
+ return merged;
+}
+
+/**
+ * Compute cluster labels for X given cluster centers.
+ * Each point is assigned to its nearest center.
+ */
+export function clusterLabels(
+ X: Float64Array[],
+ centers: Float64Array[],
+): Int32Array {
+ const labels = new Int32Array(X.length);
+ for (let i = 0; i < X.length; i++) {
+ const xi = X[i]!;
+ let best = -1;
+ let bestDist = Number.POSITIVE_INFINITY;
+ for (let j = 0; j < centers.length; j++) {
+ const c = centers[j]!;
+ let d2 = 0;
+ for (let k = 0; k < xi.length; k++) {
+ d2 += ((xi[k] ?? 0) - (c[k] ?? 0)) ** 2;
+ }
+ if (d2 < bestDist) { bestDist = d2; best = j; }
+ }
+ labels[i] = best;
+ }
+ return labels;
+}
+
+/**
+ * Compute inertia (within-cluster sum of squared distances to centers).
+ */
+export function computeInertia(
+ X: Float64Array[],
+ centers: Float64Array[],
+ labels: Int32Array,
+): number {
+ let inertia = 0;
+ for (let i = 0; i < X.length; i++) {
+ const xi = X[i]!;
+ const c = centers[labels[i]!]!;
+ let d2 = 0;
+ for (let k = 0; k < xi.length; k++) {
+ d2 += ((xi[k] ?? 0) - (c[k] ?? 0)) ** 2;
+ }
+ inertia += d2;
+ }
+ return inertia;
+}
+
+/**
+ * Compute cluster centers from assignments.
+ */
+export function computeCenters(
+ X: Float64Array[],
+ labels: Int32Array,
+ nClusters: number,
+): Float64Array[] {
+ const d = X[0]?.length ?? 0;
+ const sums: Float64Array[] = Array.from({ length: nClusters }, () => new Float64Array(d));
+ const counts = new Int32Array(nClusters);
+
+ for (let i = 0; i < X.length; i++) {
+ const xi = X[i]!;
+ const lbl = labels[i] ?? 0;
+ if (lbl >= 0 && lbl < nClusters) {
+ counts[lbl]!++;
+ for (let k = 0; k < d; k++) sums[lbl]![k]! += xi[k] ?? 0;
+ }
+ }
+
+ return sums.map((s, j) => {
+ const cnt = counts[j] ?? 1;
+ if (cnt === 0) return s;
+ const c = new Float64Array(d);
+ for (let k = 0; k < d; k++) c[k] = (s[k] ?? 0) / cnt;
+ return c;
+ });
+}
diff --git a/src/cluster/feature_agglomeration.ts b/src/cluster/feature_agglomeration.ts
new file mode 100644
index 0000000..0a0ca57
--- /dev/null
+++ b/src/cluster/feature_agglomeration.ts
@@ -0,0 +1,169 @@
+/**
+ * FeatureAgglomeration — hierarchical clustering applied to features (columns).
+ * Each sample's features are grouped; the representative value (mean/median/max)
+ * of each group becomes the transformed feature.
+ *
+ * Ports: FeatureAgglomeration
+ */
+
+import { BaseEstimator } from "../base.js";
+
+export interface FeatureAgglomerationOptions {
+ nClusters?: number;
+ poolingFunc?: "mean" | "median" | "max" | "min";
+ linkage?: "ward" | "complete" | "average" | "single";
+}
+
+function columnMean(X: Float64Array[], col: number): number {
+ let s = 0;
+ for (const row of X) s += row[col] ?? 0;
+ return s / X.length;
+}
+
+function colDist(X: Float64Array[], a: number, b: number): number {
+ const ma = columnMean(X, a);
+ const mb = columnMean(X, b);
+ return Math.abs(ma - mb);
+}
+
+/**
+ * Agglomerative (bottom-up) clustering on columns using average-column-value distance.
+ * Returns an array mapping each column → cluster index (0-based).
+ */
+function agglomerateCols(
+ X: Float64Array[],
+ nClusters: number,
+ _linkage: string,
+): Int32Array {
+ const nFeatures = X[0]?.length ?? 0;
+ if (nClusters >= nFeatures) {
+ return Int32Array.from({ length: nFeatures }, (_, i) => i);
+ }
+ // Start: each feature is its own cluster
+ const assignments = Int32Array.from({ length: nFeatures }, (_, i) => i);
+ let nActive = nFeatures;
+ // Track which features belong to each cluster
+ const clusters: number[][] = Array.from({ length: nFeatures }, (_, i) => [i]);
+
+ while (nActive > nClusters) {
+ // Find two closest clusters (by mean column distance)
+ let minDist = Number.POSITIVE_INFINITY;
+ let mergeA = -1;
+ let mergeB = -1;
+ const activeIds = [...new Set(Array.from(assignments))].sort((a, b) => a - b);
+ for (let ai = 0; ai < activeIds.length; ai++) {
+ for (let bi = ai + 1; bi < activeIds.length; bi++) {
+ const ca = activeIds[ai] ?? 0;
+ const cb = activeIds[bi] ?? 0;
+ const colsA = clusters[ca] ?? [];
+ const colsB = clusters[cb] ?? [];
+ // average linkage between column groups
+ let d = 0;
+ let count = 0;
+ for (const fa of colsA) {
+ for (const fb of colsB) {
+ d += colDist(X, fa, fb);
+ count++;
+ }
+ }
+ d = count > 0 ? d / count : Number.POSITIVE_INFINITY;
+ if (d < minDist) {
+ minDist = d;
+ mergeA = ca;
+ mergeB = cb;
+ }
+ }
+ }
+ if (mergeA < 0 || mergeB < 0) break;
+ // Merge mergeB into mergeA
+ const colsB = clusters[mergeB] ?? [];
+ for (const col of colsB) {
+ assignments[col] = mergeA;
+ }
+ clusters[mergeA] = [...(clusters[mergeA] ?? []), ...colsB];
+ clusters[mergeB] = [];
+ nActive--;
+ }
+ // Remap cluster IDs to 0..nClusters-1
+ const idMap = new Map();
+ let nextId = 0;
+ for (let i = 0; i < assignments.length; i++) {
+ const a = assignments[i] ?? 0;
+ if (!idMap.has(a)) idMap.set(a, nextId++);
+ assignments[i] = idMap.get(a) ?? 0;
+ }
+ return assignments;
+}
+
+/**
+ * Cluster features using hierarchical clustering and pool each group.
+ */
+export class FeatureAgglomeration extends BaseEstimator {
+ nClusters: number;
+ poolingFunc: "mean" | "median" | "max" | "min";
+ linkage: "ward" | "complete" | "average" | "single";
+
+ labels_!: Int32Array;
+ nClusters_!: number;
+
+ constructor(options: FeatureAgglomerationOptions = {}) {
+ super();
+ this.nClusters = options.nClusters ?? 2;
+ this.poolingFunc = options.poolingFunc ?? "mean";
+ this.linkage = options.linkage ?? "ward";
+ }
+
+ fit(X: Float64Array[]): this {
+ this.labels_ = agglomerateCols(X, this.nClusters, this.linkage);
+ this.nClusters_ = new Set(Array.from(this.labels_)).size;
+ return this;
+ }
+
+ transform(X: Float64Array[]): Float64Array[] {
+ if (this.labels_ === undefined) throw new Error("Not fitted");
+ const k = this.nClusters_;
+ return X.map((row) => {
+ const groups: number[][] = Array.from({ length: k }, () => []);
+ for (let j = 0; j < row.length; j++) {
+ const cid = this.labels_[j] ?? 0;
+ (groups[cid] ?? []).push(row[j] ?? 0);
+ }
+ const out = new Float64Array(k);
+ for (let c = 0; c < k; c++) {
+ const vals = groups[c] ?? [];
+ if (vals.length === 0) { out[c] = 0; continue; }
+ if (this.poolingFunc === "mean") {
+ out[c] = vals.reduce((a, b) => a + b, 0) / vals.length;
+ } else if (this.poolingFunc === "median") {
+ const s = [...vals].sort((a, b) => a - b);
+ const m = Math.floor(s.length / 2);
+ out[c] = s.length % 2 === 0
+ ? ((s[m - 1] ?? 0) + (s[m] ?? 0)) / 2
+ : (s[m] ?? 0);
+ } else if (this.poolingFunc === "max") {
+ out[c] = Math.max(...vals);
+ } else {
+ out[c] = Math.min(...vals);
+ }
+ }
+ return out;
+ });
+ }
+
+ fitTransform(X: Float64Array[]): Float64Array[] {
+ return this.fit(X).transform(X);
+ }
+
+ /** Reconstruct original shape from reduced representation. */
+ inverseTransform(Xred: Float64Array[]): Float64Array[] {
+ if (this.labels_ === undefined) throw new Error("Not fitted");
+ const nFeatures = this.labels_.length;
+ return Xred.map((row) => {
+ const out = new Float64Array(nFeatures);
+ for (let j = 0; j < nFeatures; j++) {
+ out[j] = row[this.labels_[j] ?? 0] ?? 0;
+ }
+ return out;
+ });
+ }
+}
diff --git a/src/cluster/hdbscan.ts b/src/cluster/hdbscan.ts
new file mode 100644
index 0000000..2a1f489
--- /dev/null
+++ b/src/cluster/hdbscan.ts
@@ -0,0 +1,189 @@
+/**
+ * HDBSCAN — Hierarchical Density-Based Spatial Clustering of Applications with Noise.
+ * Mirrors sklearn.cluster.HDBSCAN.
+ */
+
+import { NotFittedError } from "../exceptions.js";
+
+export interface HDBSCANOptions {
+ minClusterSize?: number;
+ minSamples?: number | null;
+ clusterSelectionEpsilon?: number;
+ maxClusterSize?: number | null;
+ alpha?: number;
+ clusterSelectionMethod?: "eom" | "leaf";
+ allowSingleCluster?: boolean;
+ metric?: "euclidean" | "manhattan" | "chebyshev";
+}
+
+/**
+ * HDBSCAN clustering algorithm.
+ * Extends DBSCAN by converting it into a hierarchical clustering then using a stability
+ * criterion to extract a flat clustering.
+ */
+export class HDBSCAN {
+ minClusterSize: number;
+ minSamples: number;
+ clusterSelectionEpsilon: number;
+ alpha: number;
+ clusterSelectionMethod: "eom" | "leaf";
+ allowSingleCluster: boolean;
+ metric: "euclidean" | "manhattan" | "chebyshev";
+
+ labels_: Int32Array | null = null;
+ probabilities_: Float64Array | null = null;
+ clusterPersistence_: Float64Array | null = null;
+ nFeatures_: number = 0;
+
+ constructor(options: HDBSCANOptions = {}) {
+ this.minClusterSize = options.minClusterSize ?? 5;
+ this.minSamples = options.minSamples ?? 5;
+ this.clusterSelectionEpsilon = options.clusterSelectionEpsilon ?? 0;
+ this.alpha = options.alpha ?? 1.0;
+ this.clusterSelectionMethod = options.clusterSelectionMethod ?? "eom";
+ this.allowSingleCluster = options.allowSingleCluster ?? false;
+ this.metric = options.metric ?? "euclidean";
+ }
+
+ private _dist(a: Float64Array, b: Float64Array): number {
+ const p = a.length;
+ if (this.metric === "manhattan") {
+ let s = 0;
+ for (let j = 0; j < p; j++) s += Math.abs((a[j] ?? 0) - (b[j] ?? 0));
+ return s;
+ }
+ if (this.metric === "chebyshev") {
+ let s = 0;
+ for (let j = 0; j < p; j++) s = Math.max(s, Math.abs((a[j] ?? 0) - (b[j] ?? 0)));
+ return s;
+ }
+ let s = 0;
+ for (let j = 0; j < p; j++) s += ((a[j] ?? 0) - (b[j] ?? 0)) ** 2;
+ return Math.sqrt(s);
+ }
+
+ fit(X: Float64Array[]): this {
+ const n = X.length;
+ this.nFeatures_ = X[0]?.length ?? 0;
+
+ // Compute pairwise distances
+ const dists: Float64Array[] = Array.from({ length: n }, () => new Float64Array(n));
+ for (let i = 0; i < n; i++) {
+ for (let j = i + 1; j < n; j++) {
+ const d = this._dist(X[i]!, X[j]!);
+ dists[i]![j]! = d;
+ dists[j]![i]! = d;
+ }
+ }
+
+ // Core distances (kth nearest neighbor distance)
+ const k = Math.min(this.minSamples, n - 1);
+ const coreDists = new Float64Array(n);
+ for (let i = 0; i < n; i++) {
+ const sorted = Array.from(dists[i]!).filter((_, j) => j !== i).sort((a, b) => a - b);
+ coreDists[i]! = sorted[k - 1] ?? 0;
+ }
+
+ // Mutual reachability distances
+ const mrd: Float64Array[] = Array.from({ length: n }, () => new Float64Array(n));
+ for (let i = 0; i < n; i++) {
+ for (let j = 0; j < n; j++) {
+ if (i === j) continue;
+ mrd[i]![j]! = Math.max(coreDists[i]!, coreDists[j]!, dists[i]![j]!);
+ }
+ }
+
+ // Build MST (Prim's algorithm)
+ const inMST = new Uint8Array(n);
+ const minEdge = new Float64Array(n).fill(Number.POSITIVE_INFINITY);
+ const parent = new Int32Array(n).fill(-1);
+ minEdge[0]! = 0;
+
+ const edges: Array<[number, number, number]> = [];
+ for (let step = 0; step < n; step++) {
+ let u = -1;
+ for (let i = 0; i < n; i++) {
+ if (!inMST[i] && (u < 0 || (minEdge[i] ?? 0) < (minEdge[u] ?? 0))) u = i;
+ }
+ if (u < 0) break;
+ inMST[u]! = 1;
+ if (parent[u]! >= 0) edges.push([parent[u]!, u, mrd[parent[u]!]![u]!]);
+ for (let v = 0; v < n; v++) {
+ if (!inMST[v] && (mrd[u]![v]! < (minEdge[v] ?? Number.POSITIVE_INFINITY))) {
+ minEdge[v]! = mrd[u]![v]!;
+ parent[v]! = u;
+ }
+ }
+ }
+
+ // Sort MST edges by weight
+ edges.sort((a, b) => (a[2] ?? 0) - (b[2] ?? 0));
+
+ // Build hierarchy via single-linkage (union-find)
+ const uf = Array.from({ length: n }, (_, i) => i);
+ const find = (x: number): number => {
+ while (uf[x] !== x) {
+ uf[x]! = uf[uf[x]!]!;
+ x = uf[x]!;
+ }
+ return x;
+ };
+ const clusterSizes = new Int32Array(n).fill(1);
+ const labels = new Int32Array(n).fill(-1);
+
+ // Simplified flat clustering: use density-based approach
+ // Group points where edge weight <= threshold
+ const threshold = this.clusterSelectionEpsilon > 0
+ ? this.clusterSelectionEpsilon
+ : (edges[Math.floor(edges.length * 0.5)]?.[2] ?? 0);
+
+ for (const [u, v, w] of edges) {
+ if (w <= threshold) {
+ const pu = find(u);
+ const pv = find(v);
+ if (pu !== pv) {
+ const newSize = (clusterSizes[pu] ?? 1) + (clusterSizes[pv] ?? 1);
+ if ((clusterSizes[pu] ?? 1) >= (clusterSizes[pv] ?? 1)) {
+ uf[pv]! = pu;
+ clusterSizes[pu]! = newSize;
+ } else {
+ uf[pu]! = pv;
+ clusterSizes[pv]! = newSize;
+ }
+ }
+ }
+ }
+
+ // Assign cluster labels
+ const rootToCluster = new Map();
+ let nextCluster = 0;
+ for (let i = 0; i < n; i++) {
+ const root = find(i);
+ const sz = clusterSizes[root] ?? 1;
+ if (sz >= this.minClusterSize) {
+ if (!rootToCluster.has(root)) rootToCluster.set(root, nextCluster++);
+ labels[i]! = rootToCluster.get(root)!;
+ }
+ }
+
+ this.labels_ = labels;
+ this.probabilities_ = new Float64Array(n).fill(1.0);
+ // Mark noise points
+ for (let i = 0; i < n; i++) {
+ if (labels[i] === -1) this.probabilities_[i]! = 0;
+ }
+ this.clusterPersistence_ = new Float64Array(nextCluster).fill(1.0);
+ return this;
+ }
+
+ fitPredict(X: Float64Array[]): Int32Array {
+ this.fit(X);
+ if (!this.labels_) throw new NotFittedError("HDBSCAN is not fitted");
+ return this.labels_;
+ }
+
+ get nClusters_(): number {
+ if (!this.labels_) return 0;
+ return Math.max(...Array.from(this.labels_)) + 1;
+ }
+}
diff --git a/src/cluster/index.ts b/src/cluster/index.ts
new file mode 100644
index 0000000..13b3dc2
--- /dev/null
+++ b/src/cluster/index.ts
@@ -0,0 +1,9 @@
+export * from "./kmeans.js";
+export * from "./agglomerative.js";
+export * from "./spectral.js";
+export * from "./hdbscan.js";
+export * from "./bisecting_kmeans.js";
+export * from "./affinity_propagation.js";
+export * from "./feature_agglomeration.js";
+export * from "./ward.js";
+export * from "./clustering_utils.js";
diff --git a/src/cluster/kmeans.ts b/src/cluster/kmeans.ts
new file mode 100644
index 0000000..3e043d0
--- /dev/null
+++ b/src/cluster/kmeans.ts
@@ -0,0 +1,301 @@
+/**
+ * KMeans and DBSCAN clustering.
+ * Mirrors sklearn.cluster.KMeans and DBSCAN.
+ */
+
+import { NotFittedError } from "../exceptions.js";
+
+function euclideanSq(a: Float64Array, b: Float64Array): number {
+ let s = 0;
+ for (let i = 0; i < a.length; i++) {
+ s += ((a[i] ?? 0) - (b[i] ?? 0)) ** 2;
+ }
+ return s;
+}
+
+function euclidean(a: Float64Array, b: Float64Array): number {
+ return Math.sqrt(euclideanSq(a, b));
+}
+
+export class KMeans {
+ nClusters: number;
+ maxIter: number;
+ tol: number;
+ nInit: number;
+
+ clusterCenters_: Float64Array[] | null = null;
+ labels_: Int32Array | null = null;
+ inertia_: number = 0;
+
+ constructor(
+ options: {
+ nClusters?: number;
+ maxIter?: number;
+ tol?: number;
+ nInit?: number;
+ } = {},
+ ) {
+ this.nClusters = options.nClusters ?? 8;
+ this.maxIter = options.maxIter ?? 300;
+ this.tol = options.tol ?? 1e-4;
+ this.nInit = options.nInit ?? 10;
+ }
+
+ private _kmeanspp(X: Float64Array[], k: number): Float64Array[] {
+ const n = X.length;
+ const p = (X[0] ?? new Float64Array(0)).length;
+ const centers: Float64Array[] = [];
+
+ // Pick first center randomly
+ centers.push(new Float64Array(X[Math.floor(Math.random() * n)] ?? new Float64Array(p)));
+
+ for (let c = 1; c < k; c++) {
+ const dists = X.map((xi) => {
+ let minD = Number.POSITIVE_INFINITY;
+ for (const center of centers) {
+ const d = euclideanSq(xi, center);
+ if (d < minD) minD = d;
+ }
+ return minD;
+ });
+ const totalDist = dists.reduce((a, b) => a + b, 0);
+ let rand = Math.random() * totalDist;
+ let selected = 0;
+ for (let i = 0; i < n; i++) {
+ rand -= dists[i] ?? 0;
+ if (rand <= 0) {
+ selected = i;
+ break;
+ }
+ }
+ centers.push(new Float64Array(X[selected] ?? new Float64Array(p)));
+ }
+ return centers;
+ }
+
+ private _run(
+ X: Float64Array[],
+ k: number,
+ ): { centers: Float64Array[]; labels: Int32Array; inertia: number } {
+ const n = X.length;
+ const p = (X[0] ?? new Float64Array(0)).length;
+ let centers = this._kmeanspp(X, k);
+ const labels = new Int32Array(n);
+
+ for (let iter = 0; iter < this.maxIter; iter++) {
+ // Assignment step
+ for (let i = 0; i < n; i++) {
+ let minDist = Number.POSITIVE_INFINITY;
+ let minIdx = 0;
+ for (let c = 0; c < centers.length; c++) {
+ const d = euclideanSq(X[i] ?? new Float64Array(p), centers[c] ?? new Float64Array(p));
+ if (d < minDist) {
+ minDist = d;
+ minIdx = c;
+ }
+ }
+ labels[i] = minIdx;
+ }
+
+ // Update step
+ const newCenters: Float64Array[] = Array.from({ length: k }, () => new Float64Array(p));
+ const counts = new Int32Array(k);
+ for (let i = 0; i < n; i++) {
+ const c = labels[i] ?? 0;
+ counts[c] = (counts[c] ?? 0) + 1;
+ const xi = X[i] ?? new Float64Array(p);
+ const center = newCenters[c] ?? new Float64Array(p);
+ for (let j = 0; j < p; j++) {
+ center[j] = (center[j] ?? 0) + (xi[j] ?? 0);
+ }
+ }
+
+ let maxShift = 0;
+ for (let c = 0; c < k; c++) {
+ const cnt = counts[c] ?? 0;
+ const center = newCenters[c] ?? new Float64Array(p);
+ if (cnt > 0) {
+ for (let j = 0; j < p; j++) {
+ center[j] = (center[j] ?? 0) / cnt;
+ }
+ } else {
+ // Re-initialize empty cluster to a random point
+ const randIdx = Math.floor(Math.random() * n);
+ newCenters[c] = new Float64Array(X[randIdx] ?? new Float64Array(p));
+ }
+ const shift = euclideanSq(centers[c] ?? new Float64Array(p), newCenters[c] ?? new Float64Array(p));
+ if (shift > maxShift) maxShift = shift;
+ }
+ centers = newCenters;
+ if (maxShift < this.tol ** 2) break;
+ }
+
+ // Compute inertia
+ let inertia = 0;
+ for (let i = 0; i < n; i++) {
+ inertia += euclideanSq(X[i] ?? new Float64Array(p), centers[labels[i] ?? 0] ?? new Float64Array(p));
+ }
+
+ return { centers, labels, inertia };
+ }
+
+ fit(X: Float64Array[]): this {
+ const k = Math.min(this.nClusters, X.length);
+ let best: ReturnType | null = null;
+
+ for (let init = 0; init < this.nInit; init++) {
+ const result = this._run(X, k);
+ if (best === null || result.inertia < best.inertia) {
+ best = result;
+ }
+ }
+
+ this.clusterCenters_ = best?.centers ?? [];
+ this.labels_ = best?.labels ?? new Int32Array(X.length);
+ this.inertia_ = best?.inertia ?? 0;
+ return this;
+ }
+
+ predict(X: Float64Array[]): Int32Array {
+ if (this.clusterCenters_ === null) throw new NotFittedError("KMeans");
+ const centers = this.clusterCenters_;
+ const p = (centers[0] ?? new Float64Array(0)).length;
+ return new Int32Array(
+ X.map((xi) => {
+ let minDist = Number.POSITIVE_INFINITY;
+ let minIdx = 0;
+ for (let c = 0; c < centers.length; c++) {
+ const d = euclideanSq(xi, centers[c] ?? new Float64Array(p));
+ if (d < minDist) {
+ minDist = d;
+ minIdx = c;
+ }
+ }
+ return minIdx;
+ }),
+ );
+ }
+
+ fitPredict(X: Float64Array[]): Int32Array {
+ this.fit(X);
+ return this.labels_ as Int32Array;
+ }
+
+ score(X: Float64Array[]): number {
+ return -this._computeInertia(X, this.clusterCenters_ ?? []);
+ }
+
+ private _computeInertia(X: Float64Array[], centers: Float64Array[]): number {
+ const p = (centers[0] ?? new Float64Array(0)).length;
+ let inertia = 0;
+ for (const xi of X) {
+ let minDist = Number.POSITIVE_INFINITY;
+ for (const c of centers) {
+ const d = euclideanSq(xi, c.length ? c : new Float64Array(p));
+ if (d < minDist) minDist = d;
+ }
+ inertia += minDist;
+ }
+ return inertia;
+ }
+}
+
+export class DBSCAN {
+ eps: number;
+ minSamples: number;
+ metric: string;
+
+ labels_: Int32Array | null = null;
+ coreIndices_: Int32Array | null = null;
+
+ constructor(
+ options: {
+ eps?: number;
+ minSamples?: number;
+ metric?: string;
+ } = {},
+ ) {
+ this.eps = options.eps ?? 0.5;
+ this.minSamples = options.minSamples ?? 5;
+ this.metric = options.metric ?? "euclidean";
+ }
+
+ fitPredict(X: Float64Array[]): Int32Array {
+ const n = X.length;
+ const labels = new Int32Array(n).fill(-2); // -2 = unvisited, -1 = noise
+ let clusterId = 0;
+ const coreIndices: number[] = [];
+
+ function getNeighbors(idx: number): number[] {
+ const neighbors: number[] = [];
+ const xi = X[idx] ?? new Float64Array(0);
+ for (let j = 0; j < n; j++) {
+ if (euclidean(xi, X[j] ?? new Float64Array(0)) <= 0.5) {
+ // placeholder - use eps below
+ }
+ }
+ return neighbors;
+ }
+ void getNeighbors; // suppress unused warning
+
+ const eps = this.eps;
+ const minSamples = this.minSamples;
+
+ function neighbors(idx: number): number[] {
+ const xi = X[idx] ?? new Float64Array(0);
+ const result: number[] = [];
+ for (let j = 0; j < n; j++) {
+ if (euclidean(xi, X[j] ?? new Float64Array(0)) <= eps) {
+ result.push(j);
+ }
+ }
+ return result;
+ }
+
+ for (let i = 0; i < n; i++) {
+ if (labels[i] !== -2) continue;
+ const nb = neighbors(i);
+ if (nb.length < minSamples) {
+ labels[i] = -1;
+ continue;
+ }
+
+ coreIndices.push(i);
+ labels[i] = clusterId;
+ const queue = [...nb.filter((j) => j !== i)];
+
+ while (queue.length > 0) {
+ const j = queue.shift() as number;
+ if (labels[j] === -1) {
+ labels[j] = clusterId;
+ }
+ if (labels[j] !== -2) continue;
+ labels[j] = clusterId;
+ const jNb = neighbors(j);
+ if (jNb.length >= minSamples) {
+ coreIndices.push(j);
+ for (const k of jNb) {
+ if (labels[k] === -2 || labels[k] === -1) {
+ queue.push(k);
+ }
+ }
+ }
+ }
+ clusterId++;
+ }
+
+ // Fix any remaining unvisited (noise)
+ for (let i = 0; i < n; i++) {
+ if (labels[i] === -2) labels[i] = -1;
+ }
+
+ this.labels_ = labels;
+ this.coreIndices_ = new Int32Array(coreIndices);
+ return labels;
+ }
+
+ fit(X: Float64Array[]): this {
+ this.fitPredict(X);
+ return this;
+ }
+}
diff --git a/src/cluster/spectral.ts b/src/cluster/spectral.ts
new file mode 100644
index 0000000..4875131
--- /dev/null
+++ b/src/cluster/spectral.ts
@@ -0,0 +1,549 @@
+/**
+ * SpectralClustering, MeanShift, Birch, and OPTICS clustering.
+ * Mirrors sklearn.cluster SpectralClustering, MeanShift, Birch, OPTICS.
+ */
+
+import { NotFittedError } from "../exceptions.js";
+
+// ─── SpectralClustering ───────────────────────────────────────────────────────
+
+export interface SpectralClusteringOptions {
+ nClusters?: number;
+ nInit?: number;
+ gamma?: number;
+ affinityType?: "rbf" | "nearest_neighbors";
+ nNeighbors?: number;
+ randomState?: number;
+}
+
+function rbfKernel(a: Float64Array, b: Float64Array, gamma: number): number {
+ let d = 0;
+ for (let i = 0; i < a.length; i++) {
+ d += ((a[i] ?? 0) - (b[i] ?? 0)) ** 2;
+ }
+ return Math.exp(-gamma * d);
+}
+
+function computeAffinityMatrix(
+ X: Float64Array[],
+ gamma: number,
+): Float64Array[] {
+ const n = X.length;
+ return X.map((xi, i) =>
+ Float64Array.from(X, (xj, j) => {
+ if (i === j) return 0;
+ return rbfKernel(xi as Float64Array, xj as Float64Array, gamma);
+ }),
+ );
+}
+
+function symmetricNormalizedLaplacian(W: Float64Array[]): Float64Array[] {
+ const n = W.length;
+ const D = W.map((row) => row.reduce((s, v) => s + v, 0));
+ const Dinvhalf = D.map((d) => (d > 0 ? 1 / Math.sqrt(d) : 0));
+ return W.map((row, i) =>
+ Float64Array.from(row, (w, j) => (Dinvhalf[i] ?? 0) * w * (Dinvhalf[j] ?? 0)),
+ );
+}
+
+function powerIterationEigenvectors(
+ L: Float64Array[],
+ k: number,
+ maxIter = 300,
+): Float64Array[] {
+ const n = L.length;
+ const rng = { seed: 42 };
+ const rand = () => {
+ rng.seed = (rng.seed * 1664525 + 1013904223) & 0xffffffff;
+ return (rng.seed >>> 0) / 0xffffffff;
+ };
+ // Initialize random vectors
+ const vecs: Float64Array[] = Array.from({ length: k }, () =>
+ Float64Array.from({ length: n }, () => rand() - 0.5),
+ );
+
+ for (let iter = 0; iter < maxIter; iter++) {
+ // Orthogonalize and normalize via QR (Gram-Schmidt)
+ for (let col = 0; col < k; col++) {
+ const v = vecs[col] as Float64Array;
+ // Multiply: v = L @ v
+ const Lv = new Float64Array(n);
+ for (let i = 0; i < n; i++) {
+ const row = L[i] as Float64Array;
+ let s = 0;
+ for (let j = 0; j < n; j++) s += (row[j] ?? 0) * (v[j] ?? 0);
+ Lv[i] = s;
+ }
+ // Subtract projections of previous vectors
+ for (let prev = 0; prev < col; prev++) {
+ const u = vecs[prev] as Float64Array;
+ let dot = 0;
+ for (let i = 0; i < n; i++) dot += (Lv[i] ?? 0) * (u[i] ?? 0);
+ for (let i = 0; i < n; i++) Lv[i]! -= dot * (u[i] ?? 0);
+ }
+ // Normalize
+ let norm = 0;
+ for (let i = 0; i < n; i++) norm += (Lv[i] ?? 0) ** 2;
+ norm = Math.sqrt(norm) || 1;
+ for (let i = 0; i < n; i++) Lv[i]! /= norm;
+ vecs[col] = Lv;
+ }
+ }
+ return vecs;
+}
+
+function kmeansOnRows(
+ rows: Float64Array[],
+ k: number,
+ maxIter = 100,
+ nInit = 10,
+): Int32Array {
+ const n = rows.length;
+ const d = rows[0]?.length ?? 0;
+ let bestLabels = new Int32Array(n);
+ let bestInertia = Number.POSITIVE_INFINITY;
+
+ const rng = { seed: 0 };
+ const rand = () => {
+ rng.seed = (rng.seed * 1664525 + 1013904223) & 0xffffffff;
+ return (rng.seed >>> 0) / 0xffffffff;
+ };
+
+ for (let init = 0; init < nInit; init++) {
+ rng.seed = init * 1234 + 5678;
+ const centers: Float64Array[] = Array.from({ length: k }, () => {
+ const idx = Math.floor(rand() * n);
+ return Float64Array.from(rows[idx] ?? new Float64Array(d));
+ });
+ const labels = new Int32Array(n);
+
+ for (let iter = 0; iter < maxIter; iter++) {
+ // Assign
+ let changed = false;
+ for (let i = 0; i < n; i++) {
+ const xi = rows[i] as Float64Array;
+ let best = 0;
+ let bestDist = Number.POSITIVE_INFINITY;
+ for (let c = 0; c < k; c++) {
+ const cc = centers[c] as Float64Array;
+ let dd = 0;
+ for (let j = 0; j < d; j++) dd += ((xi[j] ?? 0) - (cc[j] ?? 0)) ** 2;
+ if (dd < bestDist) { bestDist = dd; best = c; }
+ }
+ if (labels[i] !== best) { labels[i]! = best; changed = true; }
+ }
+ if (!changed) break;
+ // Update centers
+ for (const c of centers) c.fill(0);
+ const counts = new Int32Array(k);
+ for (let i = 0; i < n; i++) {
+ const c = labels[i] ?? 0;
+ counts[c]! += 1;
+ const cc = centers[c] as Float64Array;
+ const xi = rows[i] as Float64Array;
+ for (let j = 0; j < d; j++) cc[j]! += xi[j] ?? 0;
+ }
+ for (let c = 0; c < k; c++) {
+ const cnt = counts[c] ?? 1;
+ if (cnt > 0) {
+ const cc = centers[c] as Float64Array;
+ for (let j = 0; j < d; j++) cc[j]! /= cnt;
+ }
+ }
+ }
+
+ let inertia = 0;
+ for (let i = 0; i < n; i++) {
+ const xi = rows[i] as Float64Array;
+ const cc = centers[labels[i] ?? 0] as Float64Array;
+ for (let j = 0; j < d; j++) inertia += ((xi[j] ?? 0) - (cc[j] ?? 0)) ** 2;
+ }
+ if (inertia < bestInertia) {
+ bestInertia = inertia;
+ bestLabels = Int32Array.from(labels);
+ }
+ }
+ return bestLabels;
+}
+
+export class SpectralClustering {
+ nClusters: number;
+ nInit: number;
+ gamma: number;
+
+ labels_: Int32Array | null = null;
+ affinityMatrix_: Float64Array[] | null = null;
+
+ constructor(opts: SpectralClusteringOptions = {}) {
+ this.nClusters = opts.nClusters ?? 8;
+ this.nInit = opts.nInit ?? 10;
+ this.gamma = opts.gamma ?? 1.0;
+ }
+
+ fit(X: Float64Array[]): this {
+ const W = computeAffinityMatrix(X, this.gamma);
+ this.affinityMatrix_ = W;
+ const L = symmetricNormalizedLaplacian(W);
+ const vecs = powerIterationEigenvectors(L, this.nClusters);
+ const n = X.length;
+ const k = this.nClusters;
+ // Assemble rows from eigenvectors
+ const rows: Float64Array[] = Array.from({ length: n }, (_, i) => {
+ const row = new Float64Array(k);
+ for (let c = 0; c < k; c++) {
+ row[c]! = (vecs[c] as Float64Array)[i] ?? 0;
+ }
+ return row;
+ });
+ // Normalize rows to unit norm
+ for (const row of rows) {
+ let norm = 0;
+ for (let j = 0; j < k; j++) norm += (row[j] ?? 0) ** 2;
+ norm = Math.sqrt(norm) || 1;
+ for (let j = 0; j < k; j++) row[j]! /= norm;
+ }
+ this.labels_ = kmeansOnRows(rows, this.nClusters, 100, this.nInit);
+ return this;
+ }
+
+ fitPredict(X: Float64Array[]): Int32Array {
+ this.fit(X);
+ return this.labels_ as Int32Array;
+ }
+}
+
+// ─── MeanShift ────────────────────────────────────────────────────────────────
+
+export interface MeanShiftOptions {
+ bandwidth?: number;
+ maxIter?: number;
+ tol?: number;
+}
+
+function gaussianKernelWeight(dist2: number, bandwidth: number): number {
+ return Math.exp(-dist2 / (2 * bandwidth * bandwidth));
+}
+
+export class MeanShift {
+ bandwidth: number;
+ maxIter: number;
+ tol: number;
+
+ clusterCenters_: Float64Array[] | null = null;
+ labels_: Int32Array | null = null;
+
+ constructor(opts: MeanShiftOptions = {}) {
+ this.bandwidth = opts.bandwidth ?? 1.0;
+ this.maxIter = opts.maxIter ?? 300;
+ this.tol = opts.tol ?? 1e-3;
+ }
+
+ fit(X: Float64Array[]): this {
+ const n = X.length;
+ const d = X[0]?.length ?? 0;
+ // Initialize one seed per point
+ const seeds: Float64Array[] = X.map((x) => Float64Array.from(x));
+
+ for (const seed of seeds) {
+ for (let iter = 0; iter < this.maxIter; iter++) {
+ const newSeed = new Float64Array(d);
+ let totalWeight = 0;
+ for (const xi of X) {
+ let dist2 = 0;
+ for (let j = 0; j < d; j++) dist2 += ((seed[j] ?? 0) - (xi[j] ?? 0)) ** 2;
+ const w = gaussianKernelWeight(dist2, this.bandwidth);
+ totalWeight += w;
+ for (let j = 0; j < d; j++) newSeed[j]! += w * (xi[j] ?? 0);
+ }
+ if (totalWeight > 0) {
+ for (let j = 0; j < d; j++) newSeed[j]! /= totalWeight;
+ }
+ let shift = 0;
+ for (let j = 0; j < d; j++) shift += ((newSeed[j] ?? 0) - (seed[j] ?? 0)) ** 2;
+ for (let j = 0; j < d; j++) seed[j]! = newSeed[j] ?? 0;
+ if (Math.sqrt(shift) < this.tol) break;
+ }
+ }
+
+ // Merge nearby seeds
+ const mergedCenters: Float64Array[] = [];
+ for (const seed of seeds) {
+ let merged = false;
+ for (const center of mergedCenters) {
+ let dist2 = 0;
+ for (let j = 0; j < d; j++) dist2 += ((seed[j] ?? 0) - (center[j] ?? 0)) ** 2;
+ if (Math.sqrt(dist2) < this.bandwidth) { merged = true; break; }
+ }
+ if (!merged) mergedCenters.push(Float64Array.from(seed));
+ }
+
+ this.clusterCenters_ = mergedCenters;
+
+ // Assign labels
+ const labels = new Int32Array(n);
+ for (let i = 0; i < n; i++) {
+ const xi = X[i] as Float64Array;
+ let best = 0;
+ let bestDist = Number.POSITIVE_INFINITY;
+ for (let c = 0; c < mergedCenters.length; c++) {
+ const cc = mergedCenters[c] as Float64Array;
+ let dist2 = 0;
+ for (let j = 0; j < d; j++) dist2 += ((xi[j] ?? 0) - (cc[j] ?? 0)) ** 2;
+ if (dist2 < bestDist) { bestDist = dist2; best = c; }
+ }
+ labels[i]! = best;
+ }
+ this.labels_ = labels;
+ return this;
+ }
+
+ fitPredict(X: Float64Array[]): Int32Array {
+ this.fit(X);
+ return this.labels_ as Int32Array;
+ }
+
+ predict(X: Float64Array[]): Int32Array {
+ if (!this.clusterCenters_) throw new NotFittedError("MeanShift");
+ const n = X.length;
+ const d = X[0]?.length ?? 0;
+ const labels = new Int32Array(n);
+ for (let i = 0; i < n; i++) {
+ const xi = X[i] as Float64Array;
+ let best = 0;
+ let bestDist = Number.POSITIVE_INFINITY;
+ for (let c = 0; c < this.clusterCenters_.length; c++) {
+ const cc = this.clusterCenters_[c] as Float64Array;
+ let dist2 = 0;
+ for (let j = 0; j < d; j++) dist2 += ((xi[j] ?? 0) - (cc[j] ?? 0)) ** 2;
+ if (dist2 < bestDist) { bestDist = dist2; best = c; }
+ }
+ labels[i]! = best;
+ }
+ return labels;
+ }
+}
+
+// ─── Birch ────────────────────────────────────────────────────────────────────
+
+export interface BirchOptions {
+ threshold?: number;
+ branchingFactor?: number;
+ nClusters?: number;
+}
+
+interface CFEntry {
+ n: number;
+ ls: Float64Array;
+ ss: number;
+}
+
+export class Birch {
+ threshold: number;
+ branchingFactor: number;
+ nClusters: number;
+
+ labels_: Int32Array | null = null;
+ subclusterCenters_: Float64Array[] | null = null;
+
+ constructor(opts: BirchOptions = {}) {
+ this.threshold = opts.threshold ?? 0.5;
+ this.branchingFactor = opts.branchingFactor ?? 50;
+ this.nClusters = opts.nClusters ?? 3;
+ }
+
+ fit(X: Float64Array[]): this {
+ const n = X.length;
+ const d = X[0]?.length ?? 0;
+ const entries: CFEntry[] = [];
+
+ for (const xi of X) {
+ let inserted = false;
+ for (const entry of entries) {
+ const centroid = Float64Array.from({ length: d }, (_, j) => (entry.ls[j] ?? 0) / entry.n);
+ let dist2 = 0;
+ for (let j = 0; j < d; j++) dist2 += ((xi[j] ?? 0) - (centroid[j] ?? 0)) ** 2;
+ if (Math.sqrt(dist2) <= this.threshold) {
+ entry.n += 1;
+ for (let j = 0; j < d; j++) entry.ls[j]! += xi[j] ?? 0;
+ entry.ss += xi.reduce((s, v) => s + v * v, 0);
+ inserted = true;
+ break;
+ }
+ }
+ if (!inserted) {
+ entries.push({ n: 1, ls: Float64Array.from(xi), ss: xi.reduce((s, v) => s + v * v, 0) });
+ }
+ }
+
+ const centers: Float64Array[] = entries.map((e) =>
+ Float64Array.from({ length: d }, (_, j) => (e.ls[j] ?? 0) / e.n),
+ );
+ this.subclusterCenters_ = centers;
+
+ // Use k-means on subcluster centers
+ const k = Math.min(this.nClusters, centers.length);
+ const subcluLabels = kmeansOnRows(centers, k, 100, 3);
+
+ // Assign original points to the nearest subcluster then to its k-means label
+ const labels = new Int32Array(n);
+ for (let i = 0; i < n; i++) {
+ const xi = X[i] as Float64Array;
+ let bestIdx = 0;
+ let bestDist = Number.POSITIVE_INFINITY;
+ for (let c = 0; c < centers.length; c++) {
+ const cc = centers[c] as Float64Array;
+ let dist2 = 0;
+ for (let j = 0; j < d; j++) dist2 += ((xi[j] ?? 0) - (cc[j] ?? 0)) ** 2;
+ if (dist2 < bestDist) { bestDist = dist2; bestIdx = c; }
+ }
+ labels[i]! = subcluLabels[bestIdx] ?? 0;
+ }
+ this.labels_ = labels;
+ return this;
+ }
+
+ fitPredict(X: Float64Array[]): Int32Array {
+ this.fit(X);
+ return this.labels_ as Int32Array;
+ }
+
+ predict(X: Float64Array[]): Int32Array {
+ if (!this.subclusterCenters_) throw new NotFittedError("Birch");
+ const n = X.length;
+ const d = X[0]?.length ?? 0;
+ const labels = new Int32Array(n);
+ for (let i = 0; i < n; i++) {
+ const xi = X[i] as Float64Array;
+ let bestIdx = 0;
+ let bestDist = Number.POSITIVE_INFINITY;
+ for (let c = 0; c < this.subclusterCenters_.length; c++) {
+ const cc = this.subclusterCenters_[c] as Float64Array;
+ let dist2 = 0;
+ for (let j = 0; j < d; j++) dist2 += ((xi[j] ?? 0) - (cc[j] ?? 0)) ** 2;
+ if (dist2 < bestDist) { bestDist = dist2; bestIdx = c; }
+ }
+ labels[i]! = bestIdx;
+ }
+ return labels;
+ }
+}
+
+// ─── OPTICS ───────────────────────────────────────────────────────────────────
+
+export interface OPTICSOptions {
+ minSamples?: number;
+ maxEps?: number;
+ xi?: number;
+}
+
+export class OPTICS {
+ minSamples: number;
+ maxEps: number;
+ xi: number;
+
+ labels_: Int32Array | null = null;
+ reachabilityDistances_: Float64Array | null = null;
+ coreDistances_: Float64Array | null = null;
+ ordering_: Int32Array | null = null;
+
+ constructor(opts: OPTICSOptions = {}) {
+ this.minSamples = opts.minSamples ?? 5;
+ this.maxEps = opts.maxEps ?? Number.POSITIVE_INFINITY;
+ this.xi = opts.xi ?? 0.05;
+ }
+
+ fit(X: Float64Array[]): this {
+ const n = X.length;
+ const d = X[0]?.length ?? 0;
+
+ const dist = (a: Float64Array, b: Float64Array): number => {
+ let s = 0;
+ for (let i = 0; i < d; i++) s += ((a[i] ?? 0) - (b[i] ?? 0)) ** 2;
+ return Math.sqrt(s);
+ };
+
+ // Compute all pairwise distances (for small datasets)
+ const dists: Float64Array[] = Array.from({ length: n }, (_, i) =>
+ Float64Array.from({ length: n }, (__, j) =>
+ dist(X[i] as Float64Array, X[j] as Float64Array),
+ ),
+ );
+
+ // Compute core distances
+ const coreDist = new Float64Array(n);
+ for (let i = 0; i < n; i++) {
+ const row = Array.from(dists[i] as Float64Array).sort((a, b) => a - b);
+ coreDist[i]! = row[this.minSamples] ?? Number.POSITIVE_INFINITY;
+ }
+
+ const processed = new Uint8Array(n);
+ const reachDist = new Float64Array(n).fill(Number.POSITIVE_INFINITY);
+ const ordering: number[] = [];
+
+ const seeds: number[] = [];
+ const updateSeeds = (idx: number) => {
+ const cd = coreDist[idx] ?? Number.POSITIVE_INFINITY;
+ for (let j = 0; j < n; j++) {
+ if (processed[j]) continue;
+ const newRD = Math.max(cd, (dists[idx] as Float64Array)[j] ?? Number.POSITIVE_INFINITY);
+ if (newRD < (reachDist[j] ?? Number.POSITIVE_INFINITY)) {
+ reachDist[j]! = newRD;
+ if (!seeds.includes(j)) seeds.push(j);
+ }
+ }
+ };
+
+ for (let start = 0; start < n; start++) {
+ if (processed[start]) continue;
+ processed[start]! = 1;
+ ordering.push(start);
+ if ((coreDist[start] ?? Number.POSITIVE_INFINITY) <= this.maxEps) {
+ updateSeeds(start);
+ while (seeds.length > 0) {
+ // Pick seed with minimum reachability distance
+ let minIdx = 0;
+ let minRD = Number.POSITIVE_INFINITY;
+ for (let s = 0; s < seeds.length; s++) {
+ const sd = seeds[s] ?? 0;
+ const rd = reachDist[sd] ?? Number.POSITIVE_INFINITY;
+ if (rd < minRD) { minRD = rd; minIdx = s; }
+ }
+ const q = seeds[minIdx] ?? 0;
+ seeds.splice(minIdx, 1);
+ if (processed[q]) continue;
+ processed[q]! = 1;
+ ordering.push(q);
+ if ((coreDist[q] ?? Number.POSITIVE_INFINITY) <= this.maxEps) {
+ updateSeeds(q);
+ }
+ }
+ }
+ }
+
+ // Assign labels via xi-cluster extraction (simplified: threshold-based)
+ const labels = new Int32Array(n).fill(-1);
+ let clusterId = 0;
+ const eps = this.xi * (reachDist.reduce((mx, v) => Math.max(mx, isFinite(v) ? v : 0), 0));
+ let currentCluster = -1;
+ for (const idx of ordering) {
+ const rd = reachDist[idx] ?? Number.POSITIVE_INFINITY;
+ if (rd <= eps && (coreDist[idx] ?? Number.POSITIVE_INFINITY) <= this.maxEps) {
+ if (currentCluster === -1) { currentCluster = clusterId++; }
+ labels[idx]! = currentCluster;
+ } else {
+ currentCluster = -1;
+ }
+ }
+
+ this.labels_ = labels;
+ this.reachabilityDistances_ = reachDist;
+ this.coreDistances_ = coreDist;
+ this.ordering_ = Int32Array.from(ordering);
+ return this;
+ }
+
+ fitPredict(X: Float64Array[]): Int32Array {
+ this.fit(X);
+ return this.labels_ as Int32Array;
+ }
+}
diff --git a/src/cluster/ward.ts b/src/cluster/ward.ts
new file mode 100644
index 0000000..de0a6ad
--- /dev/null
+++ b/src/cluster/ward.ts
@@ -0,0 +1,186 @@
+/**
+ * Ward linkage and hierarchical clustering utilities.
+ * Mirrors scipy.cluster.hierarchy (linkage, fcluster, dendrogram helpers)
+ * as used within sklearn.cluster.AgglomerativeClustering.
+ */
+
+export interface LinkageRow {
+ clusterA: number;
+ clusterB: number;
+ distance: number;
+ size: number;
+}
+
+/** Compute the Ward linkage matrix for a dataset (O(n^3) naive implementation). */
+export function wardLinkage(X: Float64Array[]): LinkageRow[] {
+ const n = X.length;
+ if (n < 2) return [];
+
+ // Each point starts as its own cluster
+ const clusterPoints: Map = new Map();
+ for (let i = 0; i < n; i++) clusterPoints.set(i, [i]);
+
+ // Current cluster centroids
+ const centroids: Map = new Map();
+ for (let i = 0; i < n; i++) centroids.set(i, new Float64Array(X[i]!));
+
+ let nextCluster = n;
+ const result: LinkageRow[] = [];
+ const activeClusters = new Set(Array.from({ length: n }, (_, i) => i));
+
+ function centroid(indices: number[]): Float64Array {
+ const d = X[0]!.length;
+ const c = new Float64Array(d);
+ for (const idx of indices) {
+ const pt = X[idx]!;
+ for (let j = 0; j < d; j++) c[j]! += pt[j] ?? 0;
+ }
+ for (let j = 0; j < d; j++) c[j]! /= indices.length;
+ return c;
+ }
+
+ function wardDist(a: number, b: number): number {
+ const pa = clusterPoints.get(a)!;
+ const pb = clusterPoints.get(b)!;
+ const na = pa.length;
+ const nb = pb.length;
+ const ca = centroids.get(a)!;
+ const cb = centroids.get(b)!;
+ let sq = 0;
+ for (let j = 0; j < ca.length; j++) {
+ const diff = (ca[j] ?? 0) - (cb[j] ?? 0);
+ sq += diff * diff;
+ }
+ return Math.sqrt((na * nb) / (na + nb) * sq);
+ }
+
+ while (activeClusters.size > 1) {
+ // Find closest pair
+ const active = [...activeClusters];
+ let minDist = Number.POSITIVE_INFINITY;
+ let bestA = -1;
+ let bestB = -1;
+ for (let i = 0; i < active.length; i++) {
+ for (let j = i + 1; j < active.length; j++) {
+ const d = wardDist(active[i]!, active[j]!);
+ if (d < minDist) { minDist = d; bestA = active[i]!; bestB = active[j]!; }
+ }
+ }
+
+ const pA = clusterPoints.get(bestA)!;
+ const pB = clusterPoints.get(bestB)!;
+ const merged = [...pA, ...pB];
+ clusterPoints.set(nextCluster, merged);
+ centroids.set(nextCluster, centroid(merged));
+
+ result.push({ clusterA: bestA, clusterB: bestB, distance: minDist, size: merged.length });
+ activeClusters.delete(bestA);
+ activeClusters.delete(bestB);
+ activeClusters.add(nextCluster);
+ nextCluster++;
+ }
+
+ return result;
+}
+
+/** Flatten the linkage matrix to cluster labels (fcluster with criterion='maxclust'). */
+export function fcluster(linkage: LinkageRow[], nClusters: number, nPoints: number): Int32Array {
+ const labels = new Int32Array(nPoints);
+ if (nClusters >= nPoints) { for (let i = 0; i < nPoints; i++) labels[i] = i; return labels; }
+
+ // Track which top-level cluster each point belongs to
+ const children: Map = new Map();
+ for (const row of linkage) {
+ children.set(nPoints + children.size, [row.clusterA, row.clusterB]);
+ }
+
+ // The root is the last merged cluster
+ const root = nPoints + linkage.length - 1;
+ // BFS to assign labels — cut the tree to produce nClusters clusters
+ const cutAt = linkage.length - nClusters; // cut after this many merges from the root
+ const mergeCount = linkage.length;
+ const cutThreshold = mergeCount >= nClusters ? linkage[mergeCount - nClusters]?.distance ?? 0 : 0;
+
+ // Assign label by DFS
+ let nextLabel = 0;
+ function assign(node: number, label: number): void {
+ if (node < nPoints) { labels[node] = label; return; }
+ const ch = children.get(node);
+ if (!ch) return;
+ assign(ch[0], label);
+ assign(ch[1], label);
+ }
+
+ // Walk from root, splitting where distance > cutThreshold
+ function split(node: number, rowIdx: number): void {
+ if (node < nPoints) { labels[node] = nextLabel++; return; }
+ const ch = children.get(node);
+ if (!ch) { assign(node, nextLabel++); return; }
+ const row = linkage[rowIdx];
+ if (!row) { assign(node, nextLabel++); return; }
+ if (row.distance > cutThreshold && nextLabel < nClusters) {
+ split(ch[0], rowIdx - 1 - (linkage.length - 1 - rowIdx));
+ split(ch[1], rowIdx - 1);
+ } else {
+ assign(node, nextLabel++);
+ }
+ }
+
+ // Simple BFS approach: top nClusters nodes in the linkage
+ const queue: number[] = [root];
+ const clusters: number[] = [];
+ let label = 0;
+ while (clusters.length < nClusters && queue.length > 0) {
+ const node = queue.shift()!;
+ const ch = children.get(node);
+ if (!ch || clusters.length + queue.length >= nClusters) {
+ clusters.push(node);
+ } else {
+ queue.push(ch[0], ch[1]);
+ }
+ }
+ for (const cl of clusters) assign(cl, label++);
+
+ return labels;
+}
+
+/** Compute cophenetic distances from linkage matrix. */
+export function copheneticDistances(linkage: LinkageRow[], nPoints: number): Float64Array {
+ const n = nPoints;
+ const dist = new Float64Array(n * n);
+ // For each pair of points, find when they first merge
+ function findMerge(a: number, b: number): number {
+ // Walk through linkage in order
+ const clusterOf = new Int32Array(nPoints + linkage.length);
+ for (let i = 0; i < nPoints; i++) clusterOf[i] = i;
+ for (let step = 0; step < linkage.length; step++) {
+ const row = linkage[step]!;
+ const newId = nPoints + step;
+ // Check if a and b are in clusterA and clusterB
+ const inA = isIn(a, row.clusterA, nPoints, linkage, step);
+ const inB = isIn(b, row.clusterB, nPoints, linkage, step);
+ const inBA = isIn(b, row.clusterA, nPoints, linkage, step);
+ const inAB = isIn(a, row.clusterB, nPoints, linkage, step);
+ if ((inA && inB) || (inBA && inAB)) return row.distance;
+ }
+ return 0;
+ }
+ for (let i = 0; i < n; i++) {
+ for (let j = i + 1; j < n; j++) {
+ const d = findMerge(i, j);
+ dist[i * n + j] = d; dist[j * n + i] = d;
+ }
+ }
+ return dist;
+}
+
+function isIn(point: number, cluster: number, nPoints: number, linkage: LinkageRow[], upTo: number): boolean {
+ if (cluster === point) return true;
+ if (cluster < nPoints) return false;
+ const idx = cluster - nPoints;
+ if (idx >= upTo) return false;
+ const row = linkage[idx]!;
+ return isIn(point, row.clusterA, nPoints, linkage, idx) || isIn(point, row.clusterB, nPoints, linkage, idx);
+}
+
+export type { LinkageRow as WardLinkageRow };
diff --git a/src/compose/column_transformer.ts b/src/compose/column_transformer.ts
new file mode 100644
index 0000000..aebbab1
--- /dev/null
+++ b/src/compose/column_transformer.ts
@@ -0,0 +1,102 @@
+/**
+ * ColumnTransformer: applies transformers to columns of an array.
+ * Mirrors sklearn.compose.ColumnTransformer.
+ */
+
+import { NotFittedError } from "../exceptions.js";
+
+export interface Transformer {
+ fit(X: Float64Array[]): this;
+ transform(X: Float64Array[]): Float64Array[];
+ fitTransform?(X: Float64Array[]): Float64Array[];
+}
+
+export type ColumnSpec = number | number[] | "all";
+
+export class ColumnTransformer {
+ transformers: [string, Transformer | "passthrough" | "drop", ColumnSpec][];
+ remainder: "passthrough" | "drop";
+
+ transformers_: [string, Transformer | "passthrough", ColumnSpec][] = [];
+ private _nFeatures = 0;
+ private _allCols = new Set();
+
+ constructor(
+ transformers: [string, Transformer | "passthrough" | "drop", ColumnSpec][],
+ options: { remainder?: "passthrough" | "drop" } = {},
+ ) {
+ this.transformers = transformers;
+ this.remainder = options.remainder ?? "drop";
+ }
+
+ private _getCols(spec: ColumnSpec, nFeatures: number): number[] {
+ if (spec === "all") return Array.from({ length: nFeatures }, (_, i) => i);
+ if (typeof spec === "number") return [spec];
+ return spec;
+ }
+
+ fit(X: Float64Array[]): this {
+ const n = (X[0] ?? new Float64Array(0)).length;
+ this._nFeatures = n;
+ this._allCols.clear();
+
+ this.transformers_ = [];
+ for (const [name, t, spec] of this.transformers) {
+ if (t === "drop") continue;
+ const cols = this._getCols(spec, n);
+ for (const c of cols) this._allCols.add(c);
+
+ if (t === "passthrough") {
+ this.transformers_.push([name, "passthrough", spec]);
+ } else {
+ const Xsub = X.map((row) => new Float64Array(cols.map((c) => row[c] ?? 0)));
+ t.fit(Xsub);
+ this.transformers_.push([name, t, spec]);
+ }
+ }
+ return this;
+ }
+
+ transform(X: Float64Array[]): Float64Array[] {
+ if (this.transformers_.length === 0) throw new NotFittedError("ColumnTransformer");
+ const n = (X[0] ?? new Float64Array(0)).length;
+ const parts: Float64Array[][] = [];
+
+ for (const [, t, spec] of this.transformers_) {
+ const cols = this._getCols(spec, n);
+ const Xsub = X.map((row) => new Float64Array(cols.map((c) => row[c] ?? 0)));
+ if (t === "passthrough") {
+ parts.push(Xsub);
+ } else {
+ parts.push(t.transform(Xsub));
+ }
+ }
+
+ if (this.remainder === "passthrough") {
+ const remainderCols: number[] = [];
+ for (let c = 0; c < n; c++) {
+ if (!this._allCols.has(c)) remainderCols.push(c);
+ }
+ if (remainderCols.length > 0) {
+ parts.push(X.map((row) => new Float64Array(remainderCols.map((c) => row[c] ?? 0))));
+ }
+ }
+
+ // Horizontally concatenate
+ return X.map((_, i) => {
+ const rowParts = parts.map((p) => p[i] ?? new Float64Array(0));
+ const total = rowParts.reduce((s, r) => s + r.length, 0);
+ const result = new Float64Array(total);
+ let offset = 0;
+ for (const part of rowParts) {
+ result.set(part, offset);
+ offset += part.length;
+ }
+ return result;
+ });
+ }
+
+ fitTransform(X: Float64Array[]): Float64Array[] {
+ return this.fit(X).transform(X);
+ }
+}
diff --git a/src/compose/index.ts b/src/compose/index.ts
new file mode 100644
index 0000000..bc4c4e1
--- /dev/null
+++ b/src/compose/index.ts
@@ -0,0 +1,2 @@
+export * from "./column_transformer.js";
+export * from "./transformed_target.js";
diff --git a/src/compose/transformed_target.ts b/src/compose/transformed_target.ts
new file mode 100644
index 0000000..e7b60a5
--- /dev/null
+++ b/src/compose/transformed_target.ts
@@ -0,0 +1,117 @@
+/**
+ * TransformedTargetRegressor.
+ * Mirrors sklearn.compose.TransformedTargetRegressor.
+ */
+
+import { NotFittedError } from "../exceptions.js";
+
+export interface TransformableTarget {
+ fit(y: Float64Array): this;
+ transform(y: Float64Array): Float64Array;
+ inverseTransform(y: Float64Array): Float64Array;
+}
+
+export interface FittableRegressor {
+ fit(X: Float64Array[], y: Float64Array): this;
+ predict(X: Float64Array[]): Float64Array;
+}
+
+export interface TransformedTargetRegressorOptions {
+ regressor?: FittableRegressor;
+ transformer?: TransformableTarget;
+ func?: (y: Float64Array) => Float64Array;
+ inverseFunc?: (y: Float64Array) => Float64Array;
+ checkInverse?: boolean;
+}
+
+export class TransformedTargetRegressor {
+ regressor_: FittableRegressor | null = null;
+ transformer_: TransformableTarget | null = null;
+ func: ((y: Float64Array) => Float64Array) | null;
+ inverseFunc: ((y: Float64Array) => Float64Array) | null;
+
+ private regressorOpt: FittableRegressor | null;
+ private transformerOpt: TransformableTarget | null;
+
+ constructor(opts: TransformedTargetRegressorOptions = {}) {
+ this.regressorOpt = opts.regressor ?? null;
+ this.transformerOpt = opts.transformer ?? null;
+ this.func = opts.func ?? null;
+ this.inverseFunc = opts.inverseFunc ?? null;
+ }
+
+ fit(X: Float64Array[], y: Float64Array): this {
+ let yTrans: Float64Array;
+
+ if (this.func) {
+ yTrans = this.func(y);
+ } else if (this.transformerOpt) {
+ this.transformer_ = this.transformerOpt;
+ this.transformer_.fit(y);
+ yTrans = this.transformer_.transform(y);
+ } else {
+ // Default: identity
+ yTrans = Float64Array.from(y);
+ }
+
+ const reg = this.regressorOpt ?? createDefaultRegressor();
+ this.regressor_ = reg;
+ reg.fit(X, yTrans);
+ return this;
+ }
+
+ predict(X: Float64Array[]): Float64Array {
+ if (!this.regressor_) throw new NotFittedError("TransformedTargetRegressor");
+ const predsTrans = this.regressor_.predict(X);
+
+ if (this.inverseFunc) {
+ return this.inverseFunc(predsTrans);
+ } else if (this.transformer_) {
+ return this.transformer_.inverseTransform(predsTrans);
+ }
+ return predsTrans;
+ }
+
+ score(X: Float64Array[], y: Float64Array): number {
+ const preds = this.predict(X);
+ const mean = y.reduce((s, v) => s + v, 0) / y.length;
+ let ssRes = 0;
+ let ssTot = 0;
+ for (let i = 0; i < y.length; i++) {
+ ssRes += ((y[i] ?? 0) - (preds[i] ?? 0)) ** 2;
+ ssTot += ((y[i] ?? 0) - mean) ** 2;
+ }
+ return ssTot === 0 ? 1 : 1 - ssRes / ssTot;
+ }
+}
+
+function createDefaultRegressor(): FittableRegressor {
+ let coef: Float64Array | null = null;
+ let intercept = 0;
+ return {
+ fit(X: Float64Array[], y: Float64Array) {
+ const n = X.length;
+ const d = X[0]?.length ?? 0;
+ coef = new Float64Array(d);
+ const lr = 0.01;
+ for (let iter = 0; iter < 200; iter++) {
+ for (let i = 0; i < n; i++) {
+ const xi = X[i] as Float64Array;
+ let pred = intercept;
+ for (let j = 0; j < d; j++) pred += (coef![j] ?? 0) * (xi[j] ?? 0);
+ const err = (y[i] ?? 0) - pred;
+ intercept += lr * err;
+ for (let j = 0; j < d; j++) coef![j]! += lr * err * (xi[j] ?? 0);
+ }
+ }
+ return this;
+ },
+ predict(X: Float64Array[]) {
+ return Float64Array.from(X, (xi) => {
+ let pred = intercept;
+ for (let j = 0; j < xi.length; j++) pred += (coef![j] ?? 0) * (xi[j] ?? 0);
+ return pred;
+ });
+ },
+ };
+}
diff --git a/src/covariance/covariance.ts b/src/covariance/covariance.ts
new file mode 100644
index 0000000..534223f
--- /dev/null
+++ b/src/covariance/covariance.ts
@@ -0,0 +1,224 @@
+/**
+ * Covariance estimators: EmpiricalCovariance, ShrunkCovariance, LedoitWolf, OAS.
+ * Mirrors sklearn.covariance.
+ */
+
+import { NotFittedError } from "../exceptions.js";
+
+/** Compute column means of X. */
+function colMeans(X: Float64Array[]): Float64Array {
+ const p = (X[0] ?? new Float64Array(0)).length;
+ const means = new Float64Array(p);
+ const n = X.length;
+ for (const xi of X) {
+ for (let j = 0; j < p; j++) means[j] = (means[j] ?? 0) + (xi[j] ?? 0);
+ }
+ for (let j = 0; j < p; j++) means[j] = (means[j] ?? 0) / n;
+ return means;
+}
+
+/** Compute empirical covariance matrix (biased). */
+function empCov(X: Float64Array[], means: Float64Array): Float64Array[] {
+ const n = X.length;
+ const p = means.length;
+ const C = Array.from({ length: p }, () => new Float64Array(p));
+ for (const xi of X) {
+ for (let i = 0; i < p; i++) {
+ const di = (xi[i] ?? 0) - (means[i] ?? 0);
+ for (let j = i; j < p; j++) {
+ const dj = (xi[j] ?? 0) - (means[j] ?? 0);
+ C[i]![j] = (C[i]![j] ?? 0) + di * dj;
+ }
+ }
+ }
+ for (let i = 0; i < p; i++) {
+ C[i]![i] = (C[i]![i] ?? 0) / n;
+ for (let j = i + 1; j < p; j++) {
+ C[i]![j] = (C[i]![j] ?? 0) / n;
+ C[j]![i] = C[i]![j] ?? 0;
+ }
+ }
+ return C;
+}
+
+/**
+ * Maximum likelihood covariance estimator.
+ * Mirrors sklearn.covariance.EmpiricalCovariance.
+ */
+export class EmpiricalCovariance {
+ assumeCentered: boolean;
+
+ location_: Float64Array | null = null;
+ covariance_: Float64Array[] | null = null;
+
+ constructor(options: { assumeCentered?: boolean } = {}) {
+ this.assumeCentered = options.assumeCentered ?? false;
+ }
+
+ fit(X: Float64Array[]): this {
+ const p = (X[0] ?? new Float64Array(0)).length;
+ if (this.assumeCentered) {
+ this.location_ = new Float64Array(p);
+ } else {
+ this.location_ = colMeans(X);
+ }
+ this.covariance_ = empCov(X, this.location_);
+ return this;
+ }
+
+ score(X: Float64Array[]): number {
+ if (this.covariance_ === null || this.location_ === null) throw new NotFittedError();
+ // Negative log-likelihood
+ const n = X.length;
+ const p = (X[0] ?? new Float64Array(0)).length;
+ let logdet = 0;
+ // Approximate log-det via trace of covariance
+ for (let i = 0; i < p; i++) {
+ logdet += Math.log(Math.abs(this.covariance_[i]![i] ?? 1) + 1e-12);
+ }
+ let trace = 0;
+ for (const xi of X) {
+ const centered = new Float64Array(p);
+ for (let j = 0; j < p; j++) centered[j] = (xi[j] ?? 0) - (this.location_![j] ?? 0);
+ for (let j = 0; j < p; j++) {
+ const cjj = this.covariance_![j]![j] ?? 1e-12;
+ trace += (centered[j] ?? 0) ** 2 / (cjj || 1e-12);
+ }
+ }
+ return -(n * logdet + trace) / 2;
+ }
+
+ mahalanobis(X: Float64Array[]): Float64Array {
+ if (this.covariance_ === null || this.location_ === null) throw new NotFittedError();
+ const p = (X[0] ?? new Float64Array(0)).length;
+ const dists = new Float64Array(X.length);
+ for (let idx = 0; idx < X.length; idx++) {
+ const xi = X[idx] ?? new Float64Array(p);
+ let d = 0;
+ for (let j = 0; j < p; j++) {
+ const diff = (xi[j] ?? 0) - (this.location_![j] ?? 0);
+ const cjj = this.covariance_![j]![j] ?? 1e-12;
+ d += diff ** 2 / (cjj || 1e-12);
+ }
+ dists[idx] = Math.sqrt(d);
+ }
+ return dists;
+ }
+}
+
+/**
+ * Covariance estimator with shrinkage.
+ * Mirrors sklearn.covariance.ShrunkCovariance.
+ */
+export class ShrunkCovariance extends EmpiricalCovariance {
+ shrinkage: number;
+
+ constructor(options: { assumeCentered?: boolean; shrinkage?: number } = {}) {
+ super(options);
+ this.shrinkage = options.shrinkage ?? 0.1;
+ }
+
+ override fit(X: Float64Array[]): this {
+ super.fit(X);
+ if (this.covariance_ !== null) {
+ const p = this.covariance_.length;
+ for (let i = 0; i < p; i++) {
+ for (let j = 0; j < p; j++) {
+ if (i === j) continue;
+ this.covariance_[i]![j] = (this.covariance_![i]![j] ?? 0) * (1 - this.shrinkage);
+ }
+ }
+ }
+ return this;
+ }
+}
+
+/**
+ * Ledoit-Wolf automatic covariance estimator.
+ * Mirrors sklearn.covariance.LedoitWolf.
+ */
+export class LedoitWolf extends EmpiricalCovariance {
+ blockSize: number;
+
+ shrinkage_: number | null = null;
+
+ constructor(options: { assumeCentered?: boolean; blockSize?: number } = {}) {
+ super(options);
+ this.blockSize = options.blockSize ?? 1000;
+ }
+
+ override fit(X: Float64Array[]): this {
+ super.fit(X);
+ const n = X.length;
+ const p = (X[0] ?? new Float64Array(0)).length;
+ if (this.covariance_ !== null) {
+ // Oracle Approximating Shrinkage estimator (simplified Ledoit-Wolf)
+ let mu = 0;
+ for (let i = 0; i < p; i++) mu += this.covariance_![i]![i] ?? 0;
+ mu /= p;
+
+ let delta = 0;
+ for (let i = 0; i < p; i++) {
+ for (let j = 0; j < p; j++) {
+ delta += (this.covariance_![i]![j] ?? 0) ** 2;
+ }
+ }
+
+ const traceS2 = delta;
+ const traceS = p * mu;
+ const beta = (1 / (n * p)) * (traceS2 - traceS ** 2 / p);
+ const alpha = Math.max(0, Math.min(1, beta / delta));
+ this.shrinkage_ = alpha;
+
+ for (let i = 0; i < p; i++) {
+ for (let j = 0; j < p; j++) {
+ this.covariance_![i]![j] =
+ (1 - alpha) * (this.covariance_![i]![j] ?? 0) + (i === j ? alpha * mu : 0);
+ }
+ }
+ }
+ return this;
+ }
+}
+
+/**
+ * Oracle Approximating Shrinkage estimator.
+ * Mirrors sklearn.covariance.OAS.
+ */
+export class OAS extends EmpiricalCovariance {
+ shrinkage_: number | null = null;
+
+ override fit(X: Float64Array[]): this {
+ super.fit(X);
+ const n = X.length;
+ const p = (X[0] ?? new Float64Array(0)).length;
+ if (this.covariance_ !== null) {
+ let trS = 0;
+ let trS2 = 0;
+ for (let i = 0; i < p; i++) {
+ const sii = this.covariance_![i]![i] ?? 0;
+ trS += sii;
+ for (let j = 0; j < p; j++) {
+ trS2 += (this.covariance_![i]![j] ?? 0) ** 2;
+ }
+ }
+ const mu = trS / p;
+ const rho = Math.max(
+ 0,
+ Math.min(
+ 1,
+ ((1 - 2 / p) * trS2 + trS ** 2) /
+ ((n + 1 - 2 / p) * (trS2 - trS ** 2 / p)),
+ ),
+ );
+ this.shrinkage_ = rho;
+ for (let i = 0; i < p; i++) {
+ for (let j = 0; j < p; j++) {
+ this.covariance_![i]![j] =
+ (1 - rho) * (this.covariance_![i]![j] ?? 0) + (i === j ? rho * mu : 0);
+ }
+ }
+ }
+ return this;
+ }
+}
diff --git a/src/covariance/elliptic_envelope.ts b/src/covariance/elliptic_envelope.ts
new file mode 100644
index 0000000..22ad7f2
--- /dev/null
+++ b/src/covariance/elliptic_envelope.ts
@@ -0,0 +1,245 @@
+/**
+ * EllipticEnvelope: outlier detection via robust covariance estimation.
+ * Mirrors sklearn.covariance.EllipticEnvelope.
+ */
+
+import { NotFittedError } from "../exceptions.js";
+
+function colMeans(X: Float64Array[]): Float64Array {
+ const p = (X[0] ?? new Float64Array(0)).length;
+ const means = new Float64Array(p);
+ const n = X.length;
+ for (const xi of X) {
+ for (let j = 0; j < p; j++) means[j] = (means[j] ?? 0) + (xi[j] ?? 0);
+ }
+ for (let j = 0; j < p; j++) means[j] = (means[j] ?? 0) / n;
+ return means;
+}
+
+function empCov(X: Float64Array[], means: Float64Array): Float64Array[] {
+ const n = X.length;
+ const p = means.length;
+ const C = Array.from({ length: p }, () => new Float64Array(p));
+ for (const xi of X) {
+ for (let i = 0; i < p; i++) {
+ const di = (xi[i] ?? 0) - (means[i] ?? 0);
+ for (let j = i; j < p; j++) {
+ const dj = (xi[j] ?? 0) - (means[j] ?? 0);
+ C[i]![j] = (C[i]![j] ?? 0) + di * dj;
+ }
+ }
+ }
+ for (let i = 0; i < p; i++) {
+ C[i]![i] = (C[i]![i] ?? 0) / n;
+ for (let j = i + 1; j < p; j++) {
+ C[i]![j] = (C[i]![j] ?? 0) / n;
+ C[j]![i] = C[i]![j] ?? 0;
+ }
+ }
+ return C;
+}
+
+/** Compute log-determinant of a positive-definite matrix via Cholesky. */
+function logDet(M: Float64Array[]): number {
+ const p = M.length;
+ const L = Array.from({ length: p }, () => new Float64Array(p));
+ for (let i = 0; i < p; i++) {
+ for (let j = 0; j <= i; j++) {
+ let s = M[i]![j] ?? 0;
+ for (let k = 0; k < j; k++) s -= (L[i]![k] ?? 0) * (L[j]![k] ?? 0);
+ if (i === j) {
+ L[i]![j] = Math.sqrt(Math.max(s, 1e-12));
+ } else {
+ L[i]![j] = s / Math.max(L[j]![j] ?? 1e-12, 1e-12);
+ }
+ }
+ }
+ let logd = 0;
+ for (let i = 0; i < p; i++) logd += Math.log(Math.max(L[i]![i] ?? 1e-12, 1e-12));
+ return 2 * logd;
+}
+
+/** Invert a matrix via Gauss-Jordan. Returns null if singular. */
+function invertMatrix(M: Float64Array[]): Float64Array[] | null {
+ const p = M.length;
+ const A = M.map((row) => new Float64Array(row));
+ const I = Array.from({ length: p }, (_, i) => {
+ const r = new Float64Array(p);
+ r[i] = 1;
+ return r;
+ });
+ for (let col = 0; col < p; col++) {
+ let pivotRow = -1;
+ let pivotVal = 0;
+ for (let row = col; row < p; row++) {
+ if (Math.abs(A[row]![col] ?? 0) > Math.abs(pivotVal)) {
+ pivotVal = A[row]![col] ?? 0;
+ pivotRow = row;
+ }
+ }
+ if (pivotRow === -1 || Math.abs(pivotVal) < 1e-12) return null;
+ const tmpA = A[col]!;
+ A[col] = A[pivotRow]!;
+ A[pivotRow] = tmpA;
+ const tmpI = I[col]!;
+ I[col] = I[pivotRow]!;
+ I[pivotRow] = tmpI;
+ const scale = A[col]![col] ?? 1;
+ for (let j = 0; j < p; j++) {
+ A[col]![j] = (A[col]![j] ?? 0) / scale;
+ I[col]![j] = (I[col]![j] ?? 0) / scale;
+ }
+ for (let row = 0; row < p; row++) {
+ if (row === col) continue;
+ const factor = A[row]![col] ?? 0;
+ for (let j = 0; j < p; j++) {
+ A[row]![j] = (A[row]![j] ?? 0) - factor * (A[col]![j] ?? 0);
+ I[row]![j] = (I[row]![j] ?? 0) - factor * (I[col]![j] ?? 0);
+ }
+ }
+ }
+ return I;
+}
+
+/** Mahalanobis distance squared for each row. */
+function mahalanobisDistSq(
+ X: Float64Array[],
+ mean: Float64Array,
+ precisionMat: Float64Array[],
+): Float64Array {
+ const n = X.length;
+ const p = mean.length;
+ const dists = new Float64Array(n);
+ for (let i = 0; i < n; i++) {
+ const xi = X[i] ?? new Float64Array(p);
+ let d = 0;
+ for (let j = 0; j < p; j++) {
+ let row = 0;
+ for (let k = 0; k < p; k++) {
+ row += (precisionMat[j]![k] ?? 0) * ((xi[k] ?? 0) - (mean[k] ?? 0));
+ }
+ d += ((xi[j] ?? 0) - (mean[j] ?? 0)) * row;
+ }
+ dists[i] = d;
+ }
+ return dists;
+}
+
+/**
+ * EllipticEnvelope: fits a robust covariance estimate to detect outliers.
+ * Uses minimum covariance determinant (fast approximation).
+ * Mirrors sklearn.covariance.EllipticEnvelope.
+ */
+export class EllipticEnvelope {
+ contamination: number;
+ supportFraction: number | null;
+ randomState: number;
+
+ location_: Float64Array | null = null;
+ covariance_: Float64Array[] | null = null;
+ precision_: Float64Array[] | null = null;
+ threshold_: number = 0;
+ offset_: number = 0;
+
+ constructor(
+ options: {
+ contamination?: number;
+ supportFraction?: number | null;
+ randomState?: number;
+ } = {},
+ ) {
+ this.contamination = options.contamination ?? 0.1;
+ this.supportFraction = options.supportFraction ?? null;
+ this.randomState = options.randomState ?? 42;
+ }
+
+ fit(X: Float64Array[]): this {
+ const n = X.length;
+ const p = (X[0] ?? new Float64Array(0)).length;
+ const h = this.supportFraction !== null
+ ? Math.floor(this.supportFraction * n)
+ : Math.floor((n + p + 1) / 2);
+
+ // Fast MCD approximation: random subsample + C-step iterations
+ let bestDet = Number.POSITIVE_INFINITY;
+ let bestMean = new Float64Array(p);
+ let bestCov: Float64Array[] = Array.from({ length: p }, () => new Float64Array(p));
+
+ const rng = this.randomState;
+ const nTrials = 10;
+ for (let trial = 0; trial < nTrials; trial++) {
+ // Random subset of h points
+ const indices = Array.from({ length: n }, (_, i) => i);
+ // Pseudo-random shuffle using simple LCG
+ for (let i = n - 1; i > 0; i--) {
+ const j = Math.abs((rng * 1664525 + 1013904223 + i * trial * 31337) % (i + 1));
+ const tmp = indices[i]!;
+ indices[i] = indices[j]!;
+ indices[j] = tmp;
+ }
+ const subset = indices.slice(0, h).map((i) => X[i] ?? new Float64Array(p));
+
+ // C-step iterations
+ let curSubset = subset;
+ for (let cstep = 0; cstep < 30; cstep++) {
+ const mean = colMeans(curSubset);
+ const cov = empCov(curSubset, mean);
+ const inv = invertMatrix(cov);
+ if (!inv) break;
+ const dists = mahalanobisDistSq(X, mean, inv);
+ const sortedIdx = Array.from({ length: n }, (_, i) => i).sort(
+ (a, b) => (dists[a] ?? 0) - (dists[b] ?? 0),
+ );
+ curSubset = sortedIdx.slice(0, h).map((i) => X[i] ?? new Float64Array(p));
+ }
+
+ const mean = colMeans(curSubset);
+ const cov = empCov(curSubset, mean);
+ const det = logDet(cov);
+ if (det < bestDet) {
+ bestDet = det;
+ bestMean = mean;
+ bestCov = cov;
+ }
+ }
+
+ const inv = invertMatrix(bestCov) ?? bestCov;
+ this.location_ = bestMean;
+ this.covariance_ = bestCov;
+ this.precision_ = inv;
+
+ // Compute threshold based on contamination
+ const dists = mahalanobisDistSq(X, bestMean, inv);
+ const sorted = Array.from(dists).sort((a, b) => a - b);
+ const threshIdx = Math.floor((1 - this.contamination) * n);
+ this.threshold_ = sorted[Math.min(threshIdx, n - 1)] ?? 0;
+ this.offset_ = -this.threshold_;
+ return this;
+ }
+
+ mahalanobis(X: Float64Array[]): Float64Array {
+ if (this.location_ === null || this.precision_ === null) {
+ throw new NotFittedError("EllipticEnvelope");
+ }
+ return mahalanobisDistSq(X, this.location_, this.precision_);
+ }
+
+ decisionFunction(X: Float64Array[]): Float64Array {
+ const dists = this.mahalanobis(X);
+ return new Float64Array(dists.map((d) => -d - this.offset_));
+ }
+
+ predict(X: Float64Array[]): Int32Array {
+ const scores = this.decisionFunction(X);
+ return new Int32Array(scores.map((s) => (s >= 0 ? 1 : -1)));
+ }
+
+ score(X: Float64Array[], y: Int32Array): number {
+ const yPred = this.predict(X);
+ let correct = 0;
+ for (let i = 0; i < y.length; i++) {
+ if ((yPred[i] ?? 0) === (y[i] ?? 0)) correct++;
+ }
+ return correct / y.length;
+ }
+}
diff --git a/src/covariance/graphical_lasso.ts b/src/covariance/graphical_lasso.ts
new file mode 100644
index 0000000..00bc9e0
--- /dev/null
+++ b/src/covariance/graphical_lasso.ts
@@ -0,0 +1,252 @@
+/**
+ * GraphicalLasso and MinCovDet (robust covariance).
+ * Mirrors sklearn.covariance.GraphicalLasso and MinCovDet.
+ */
+
+import { NotFittedError } from "../exceptions.js";
+
+function colMeans(X: Float64Array[]): Float64Array {
+ const p = (X[0] ?? new Float64Array(0)).length;
+ const n = X.length;
+ const means = new Float64Array(p);
+ for (const xi of X) for (let j = 0; j < p; j++) means[j] = (means[j] ?? 0) + (xi[j] ?? 0);
+ for (let j = 0; j < p; j++) means[j] = (means[j] ?? 0) / n;
+ return means;
+}
+
+function empiricalCovariance(X: Float64Array[]): Float64Array[] {
+ const n = X.length;
+ const p = (X[0] ?? new Float64Array(0)).length;
+ const means = colMeans(X);
+ const cov: Float64Array[] = Array.from({ length: p }, () => new Float64Array(p));
+ for (const xi of X) {
+ for (let j = 0; j < p; j++) {
+ for (let k = 0; k <= j; k++) {
+ const d = ((xi[j] ?? 0) - (means[j] ?? 0)) * ((xi[k] ?? 0) - (means[k] ?? 0));
+ cov[j]![k] = (cov[j]![k] ?? 0) + d;
+ if (k !== j) cov[k]![j] = (cov[k]![j] ?? 0) + d;
+ }
+ }
+ }
+ for (let j = 0; j < p; j++) for (let k = 0; k < p; k++) cov[j]![k] = (cov[j]![k] ?? 0) / n;
+ return cov;
+}
+
+function matMul(A: Float64Array[], B: Float64Array[]): Float64Array[] {
+ const n = A.length;
+ const m = (B[0] ?? new Float64Array(0)).length;
+ const k = B.length;
+ const C: Float64Array[] = Array.from({ length: n }, () => new Float64Array(m));
+ for (let i = 0; i < n; i++) for (let j = 0; j < m; j++) for (let l = 0; l < k; l++) C[i]![j] = (C[i]![j] ?? 0) + (A[i]![l] ?? 0) * (B[l]![j] ?? 0);
+ return C;
+}
+
+function invertMatrix(A: Float64Array[]): Float64Array[] {
+ const p = A.length;
+ // Augmented matrix [A | I]
+ const M: Float64Array[] = A.map((row, i) => {
+ const r = new Float64Array(2 * p);
+ for (let j = 0; j < p; j++) r[j] = row[j] ?? 0;
+ r[p + i] = 1;
+ return r;
+ });
+
+ for (let col = 0; col < p; col++) {
+ let pivot = col;
+ for (let row = col + 1; row < p; row++) {
+ if (Math.abs(M[row]![col] ?? 0) > Math.abs(M[pivot]![col] ?? 0)) pivot = row;
+ }
+ const tmp = M[col]!; M[col] = M[pivot]!; M[pivot] = tmp;
+ const denom = M[col]![col] ?? 1;
+ for (let j = 0; j < 2 * p; j++) M[col]![j] = (M[col]![j] ?? 0) / denom;
+ for (let row = 0; row < p; row++) {
+ if (row === col) continue;
+ const factor = M[row]![col] ?? 0;
+ for (let j = 0; j < 2 * p; j++) M[row]![j] = (M[row]![j] ?? 0) - factor * (M[col]![j] ?? 0);
+ }
+ }
+
+ return M.map((row) => new Float64Array(Array.from({ length: p }, (_, j) => row[p + j] ?? 0)));
+}
+
+export interface GraphicalLassoOptions {
+ alpha?: number;
+ maxIter?: number;
+ tol?: number;
+}
+
+/**
+ * Sparse inverse covariance estimation with L1 penalty (Graphical Lasso).
+ * Mirrors sklearn.covariance.GraphicalLasso.
+ * Uses the block coordinate descent algorithm (GLASSO).
+ */
+export class GraphicalLasso {
+ alpha: number;
+ maxIter: number;
+ tol: number;
+
+ covariance_: Float64Array[] | null = null;
+ precision_: Float64Array[] | null = null;
+ nIter_: number = 0;
+ location_: Float64Array | null = null;
+
+ constructor(options: GraphicalLassoOptions = {}) {
+ this.alpha = options.alpha ?? 0.01;
+ this.maxIter = options.maxIter ?? 100;
+ this.tol = options.tol ?? 1e-4;
+ }
+
+ fit(X: Float64Array[]): this {
+ const p = (X[0] ?? new Float64Array(0)).length;
+ this.location_ = colMeans(X);
+ const S = empiricalCovariance(X);
+
+ // Initialize with diagonal of S + alpha * I
+ const W: Float64Array[] = Array.from({ length: p }, (_, i) => {
+ const row = new Float64Array(p);
+ for (let j = 0; j < p; j++) row[j] = S[i]![j] ?? 0;
+ row[i] = (row[i] ?? 0) + this.alpha;
+ return row;
+ });
+
+ for (let iter = 0; iter < this.maxIter; iter++) {
+ let maxDelta = 0;
+ for (let j = 0; j < p; j++) {
+ // Partition W into W11 (p-1 x p-1) and w12 (p-1 vector)
+ const idx = Array.from({ length: p }, (_, k) => k).filter((k) => k !== j);
+ const W11: Float64Array[] = idx.map((r) => new Float64Array(idx.map((c) => W[r]![c] ?? 0)));
+ const s12 = new Float64Array(idx.map((r) => S[r]![j] ?? 0));
+
+ // Solve lasso: W11 * beta = s12 with L1 penalty alpha
+ const W11inv = invertMatrix(W11);
+ const q = new Float64Array(p - 1);
+ for (let k = 0; k < p - 1; k++) for (let l = 0; l < p - 1; l++) q[k] = (q[k] ?? 0) + (W11inv[k]![l] ?? 0) * (s12[l] ?? 0);
+
+ // Coordinate descent for lasso subproblem
+ const beta = new Float64Array(p - 1);
+ for (let lasso = 0; lasso < 100; lasso++) {
+ let maxD = 0;
+ for (let k = 0; k < p - 1; k++) {
+ const r = (s12[k] ?? 0) - ((): number => {
+ let s = 0;
+ for (let l = 0; l < p - 1; l++) if (l !== k) s += (W11[k]![l] ?? 0) * (beta[l] ?? 0);
+ return s;
+ })();
+ const wkk = W11[k]![k] ?? 1;
+ const b = r / wkk;
+ const threshold = this.alpha / wkk;
+ const newBeta = b > threshold ? b - threshold : b < -threshold ? b + threshold : 0;
+ maxD = Math.max(maxD, Math.abs(newBeta - (beta[k] ?? 0)));
+ beta[k] = newBeta;
+ }
+ if (maxD < 1e-6) break;
+ }
+
+ // Update W: w12 = W11 * beta
+ for (let k = 0; k < p - 1; k++) {
+ let s = 0;
+ for (let l = 0; l < p - 1; l++) s += (W11[k]![l] ?? 0) * (beta[l] ?? 0);
+ const delta = Math.abs(s - (W[idx[k]!]![j] ?? 0));
+ if (delta > maxDelta) maxDelta = delta;
+ W[idx[k]!]![j] = s;
+ W[j]![idx[k]!] = s;
+ }
+ }
+ this.nIter_ = iter + 1;
+ if (maxDelta < this.tol) break;
+ }
+
+ this.covariance_ = W;
+ this.precision_ = invertMatrix(W);
+ return this;
+ }
+
+ score(X: Float64Array[]): number {
+ if (!this.covariance_) throw new NotFittedError("GraphicalLasso is not fitted yet.");
+ return 0; // Placeholder: log-likelihood requires determinant
+ }
+}
+
+export interface MinCovDetOptions {
+ support?: number | null;
+ randomState?: number;
+}
+
+/**
+ * Minimum Covariance Determinant robust estimator.
+ * Mirrors sklearn.covariance.MinCovDet.
+ * Uses a simplified C-step algorithm.
+ */
+export class MinCovDet {
+ support: number | null;
+ randomState: number;
+
+ location_: Float64Array | null = null;
+ covariance_: Float64Array[] | null = null;
+ precision_: Float64Array[] | null = null;
+ supportFraction_: number = 0;
+ supportIndices_: Int32Array | null = null;
+ rawLocation_: Float64Array | null = null;
+ rawCovariance_: Float64Array[] | null = null;
+
+ private rng_: () => number;
+
+ constructor(options: MinCovDetOptions = {}) {
+ this.support = options.support ?? null;
+ this.randomState = options.randomState ?? 0;
+ let seed = this.randomState + 1;
+ this.rng_ = () => {
+ seed = (seed * 1664525 + 1013904223) & 0xffffffff;
+ return (seed >>> 0) / 0xffffffff;
+ };
+ }
+
+ fit(X: Float64Array[]): this {
+ const n = X.length;
+ const p = (X[0] ?? new Float64Array(0)).length;
+ const h = this.support != null ? Math.floor(this.support * n) : Math.floor((n + p + 1) / 2);
+
+ // Compute Mahalanobis distances from full empirical estimate
+ const fullMeans = colMeans(X);
+ const fullCov = empiricalCovariance(X);
+ let precision: Float64Array[];
+ try { precision = invertMatrix(fullCov); } catch { precision = Array.from({ length: p }, (_, i) => { const r = new Float64Array(p); r[i] = 1; return r; }); }
+
+ // Mahalanobis distance for each point
+ const mDist = X.map((xi) => {
+ const diff = new Float64Array(p);
+ for (let j = 0; j < p; j++) diff[j] = (xi[j] ?? 0) - (fullMeans[j] ?? 0);
+ let d = 0;
+ for (let j = 0; j < p; j++) for (let k = 0; k < p; k++) d += (diff[j] ?? 0) * (precision[j]![k] ?? 0) * (diff[k] ?? 0);
+ return d;
+ });
+
+ // Select h points with smallest Mahalanobis distances
+ const sortedIdx = Array.from({ length: n }, (_, i) => i).sort((a, b) => mDist[a]! - mDist[b]!);
+ const supportIdx = new Int32Array(sortedIdx.slice(0, h));
+
+ const subset = Array.from(supportIdx).map((i) => X[i] ?? new Float64Array(p));
+ this.rawLocation_ = colMeans(subset);
+ this.rawCovariance_ = empiricalCovariance(subset);
+
+ this.location_ = this.rawLocation_;
+ this.covariance_ = this.rawCovariance_;
+ try { this.precision_ = invertMatrix(this.covariance_); } catch { this.precision_ = null; }
+
+ this.supportFraction_ = h / n;
+ this.supportIndices_ = supportIdx;
+ return this;
+ }
+
+ mahalanobis(X: Float64Array[]): Float64Array {
+ if (!this.location_ || !this.precision_) throw new NotFittedError("MinCovDet is not fitted yet.");
+ const p = this.location_.length;
+ return new Float64Array(X.map((xi) => {
+ const diff = new Float64Array(p);
+ for (let j = 0; j < p; j++) diff[j] = (xi[j] ?? 0) - (this.location_![j] ?? 0);
+ let d = 0;
+ for (let j = 0; j < p; j++) for (let k = 0; k < p; k++) d += (diff[j] ?? 0) * (this.precision_![j]![k] ?? 0) * (diff[k] ?? 0);
+ return d;
+ }));
+ }
+}
diff --git a/src/covariance/index.ts b/src/covariance/index.ts
new file mode 100644
index 0000000..30f6f71
--- /dev/null
+++ b/src/covariance/index.ts
@@ -0,0 +1,4 @@
+export * from "./covariance.js";
+export * from "./graphical_lasso.js";
+export * from "./elliptic_envelope.js";
+export * from "./precision.js";
diff --git a/src/covariance/precision.ts b/src/covariance/precision.ts
new file mode 100644
index 0000000..77b6e64
--- /dev/null
+++ b/src/covariance/precision.ts
@@ -0,0 +1,230 @@
+/**
+ * Covariance utilities: precision matrix estimation, covariance selection.
+ * ledoit_wolf() and oas() functional APIs, plus precision/correlation conversion.
+ * Mirrors sklearn.covariance functional API and utility functions.
+ */
+
+import { NotFittedError } from "../exceptions.js";
+
+function colMeans(X: Float64Array[]): Float64Array {
+ const p = (X[0] ?? new Float64Array(0)).length;
+ const m = new Float64Array(p);
+ const n = X.length;
+ for (const xi of X) for (let j = 0; j < p; j++) m[j] = (m[j] ?? 0) + (xi[j] ?? 0);
+ for (let j = 0; j < p; j++) m[j] = (m[j] ?? 0) / n;
+ return m;
+}
+
+function empCovMatrix(X: Float64Array[], means: Float64Array): Float64Array[] {
+ const n = X.length;
+ const p = means.length;
+ const C = Array.from({ length: p }, () => new Float64Array(p));
+ for (const xi of X) {
+ for (let i = 0; i < p; i++) {
+ const di = (xi[i] ?? 0) - (means[i] ?? 0);
+ for (let j = i; j < p; j++) {
+ const dj = (xi[j] ?? 0) - (means[j] ?? 0);
+ C[i]![j] = (C[i]![j] ?? 0) + di * dj;
+ }
+ }
+ }
+ for (let i = 0; i < p; i++) {
+ C[i]![i] = (C[i]![i] ?? 0) / n;
+ for (let j = i + 1; j < p; j++) {
+ C[i]![j] = (C[i]![j] ?? 0) / n;
+ C[j]![i] = C[i]![j] ?? 0;
+ }
+ }
+ return C;
+}
+
+function matTrace(M: Float64Array[]): number {
+ let s = 0;
+ for (let i = 0; i < M.length; i++) s += M[i]![i] ?? 0;
+ return s;
+}
+
+function matFrobSq(M: Float64Array[]): number {
+ let s = 0;
+ for (const row of M) for (let j = 0; j < row.length; j++) s += (row[j] ?? 0) ** 2;
+ return s;
+}
+
+/** Invert diagonal of a matrix (for precision). */
+function invertDiag(M: Float64Array[]): Float64Array[] {
+ return M.map((row, i) => new Float64Array(row.map((v, j) => i === j && v > 0 ? 1 / v : 0)));
+}
+
+/**
+ * Functional API: Ledoit-Wolf analytical shrinkage.
+ * Mirrors sklearn.covariance.ledoit_wolf.
+ */
+export function ledoitWolf(
+ X: Float64Array[],
+ options: { assumeCentered?: boolean } = {},
+): { covariance: Float64Array[]; shrinkage: number } {
+ const n = X.length;
+ const p = (X[0] ?? new Float64Array(0)).length;
+ const location = options.assumeCentered ? new Float64Array(p) : colMeans(X);
+ const S = empCovMatrix(X, location);
+ const trS = matTrace(S);
+ const trS2 = matFrobSq(S);
+ const trSsq = trS ** 2;
+
+ let delta = 0;
+ for (let i = 0; i < p; i++) {
+ for (let k = 0; k < p; k++) {
+ let fourth = 0;
+ for (let t = 0; t < n; t++) {
+ const xt = X[t] ?? new Float64Array(p);
+ fourth += ((xt[i] ?? 0) - (location[i] ?? 0)) ** 2 * ((xt[k] ?? 0) - (location[k] ?? 0)) ** 2;
+ }
+ fourth /= n;
+ delta += fourth - (S[i]![k] ?? 0) ** 2;
+ }
+ }
+ delta /= n;
+
+ const delta2 = trS2 - trSsq / p;
+ const shrinkage = delta2 > 0
+ ? Math.min(1, Math.max(0, (delta + ((n - 2) / n) * delta2) / ((n + 2) * delta2)))
+ : 0;
+
+ const mu = trS / p;
+ const covariance = S.map((row, i) =>
+ new Float64Array(row.map((v, j) => (1 - shrinkage) * v + shrinkage * (i === j ? mu : 0))),
+ );
+ return { covariance, shrinkage };
+}
+
+/**
+ * Functional API: Oracle Approximating Shrinkage (OAS).
+ * Mirrors sklearn.covariance.oas.
+ */
+export function oas(
+ X: Float64Array[],
+ options: { assumeCentered?: boolean } = {},
+): { covariance: Float64Array[]; shrinkage: number } {
+ const n = X.length;
+ const p = (X[0] ?? new Float64Array(0)).length;
+ const location = options.assumeCentered ? new Float64Array(p) : colMeans(X);
+ const S = empCovMatrix(X, location);
+ const trS = matTrace(S);
+ const trS2 = matFrobSq(S);
+ const trSsq = trS ** 2;
+
+ const num = (1 - 2 / p) * trS2 + trSsq;
+ const denom = (n + 1 - 2 / p) * (trS2 - trSsq / p);
+ const shrinkage = denom > 0 ? Math.min(1, Math.max(0, num / denom)) : 0;
+
+ const mu = trS / p;
+ const covariance = S.map((row, i) =>
+ new Float64Array(row.map((v, j) => (1 - shrinkage) * v + shrinkage * (i === j ? mu : 0))),
+ );
+ return { covariance, shrinkage };
+}
+
+/**
+ * Convert a covariance matrix to a correlation matrix.
+ * Mirrors sklearn.covariance.cov_to_corr.
+ */
+export function covToCorr(covariance: Float64Array[]): Float64Array[] {
+ const p = covariance.length;
+ const std = new Float64Array(p).map((_, i) => Math.sqrt(Math.max(covariance[i]![i] ?? 0, 1e-12)));
+ return covariance.map((row, i) =>
+ new Float64Array(row.map((v, j) => v / ((std[i] ?? 1) * (std[j] ?? 1)))),
+ );
+}
+
+/**
+ * Compute the log-likelihood of X under a Gaussian model.
+ * Mirrors sklearn.covariance.empirical_covariance (log_likelihood method).
+ */
+export function gaussianLogLikelihood(
+ X: Float64Array[],
+ mean: Float64Array,
+ covariance: Float64Array[],
+): number {
+ const n = X.length;
+ const p = mean.length;
+
+ // log-det via Cholesky
+ const L = Array.from({ length: p }, () => new Float64Array(p));
+ for (let i = 0; i < p; i++) {
+ for (let j = 0; j <= i; j++) {
+ let s = covariance[i]![j] ?? 0;
+ for (let k = 0; k < j; k++) s -= (L[i]![k] ?? 0) * (L[j]![k] ?? 0);
+ L[i]![j] = i === j ? Math.sqrt(Math.max(s, 1e-12)) : s / Math.max(L[j]![j] ?? 1, 1e-12);
+ }
+ }
+ let logDet = 0;
+ for (let i = 0; i < p; i++) logDet += Math.log(Math.max(L[i]![i] ?? 1e-12, 1e-12));
+ logDet *= 2;
+
+ // trace(S * precision) where S = empirical covariance of X
+ const S = empCovMatrix(X, mean);
+ // Use diagonal approx for precision
+ let trSP = 0;
+ for (let i = 0; i < p; i++) {
+ const cii = covariance[i]![i] ?? 1;
+ trSP += (S[i]![i] ?? 0) / Math.max(cii, 1e-12);
+ }
+
+ return -0.5 * (n * (p * Math.log(2 * Math.PI) + logDet + trSP));
+}
+
+/**
+ * Sparse inverse covariance estimator (precision matrix selector).
+ * Uses a simple soft-threshold approach to zero out small entries.
+ * Mirrors sklearn.covariance sparse precision concepts.
+ */
+export class SparsePrecision {
+ threshold: number;
+ assumeCentered: boolean;
+
+ location_: Float64Array | null = null;
+ covariance_: Float64Array[] | null = null;
+ precision_: Float64Array[] | null = null;
+
+ constructor(options: { threshold?: number; assumeCentered?: boolean } = {}) {
+ this.threshold = options.threshold ?? 0.1;
+ this.assumeCentered = options.assumeCentered ?? false;
+ }
+
+ fit(X: Float64Array[]): this {
+ const p = (X[0] ?? new Float64Array(0)).length;
+ const location = this.assumeCentered ? new Float64Array(p) : colMeans(X);
+ this.location_ = location;
+ const S = empCovMatrix(X, location);
+ this.covariance_ = S;
+
+ // Simple diagonal precision estimate with soft-thresholding
+ const P = invertDiag(S);
+ // Soft-threshold off-diagonal elements
+ this.precision_ = P.map((row, i) =>
+ new Float64Array(row.map((v, j) => {
+ if (i === j) return v;
+ return Math.abs(v) > this.threshold ? v - Math.sign(v) * this.threshold : 0;
+ })),
+ );
+ return this;
+ }
+
+ mahalanobis(X: Float64Array[]): Float64Array {
+ if (this.precision_ === null || this.location_ === null) {
+ throw new NotFittedError("SparsePrecision");
+ }
+ const P = this.precision_;
+ const mu = this.location_;
+ const p = mu.length;
+ return new Float64Array(X.map((xi) => {
+ let d = 0;
+ for (let j = 0; j < p; j++) {
+ let pRow = 0;
+ for (let k = 0; k < p; k++) pRow += (P[j]![k] ?? 0) * ((xi[k] ?? 0) - (mu[k] ?? 0));
+ d += ((xi[j] ?? 0) - (mu[j] ?? 0)) * pRow;
+ }
+ return d;
+ }));
+ }
+}
diff --git a/src/cross_decomposition/cca.ts b/src/cross_decomposition/cca.ts
new file mode 100644
index 0000000..90dbd41
--- /dev/null
+++ b/src/cross_decomposition/cca.ts
@@ -0,0 +1,260 @@
+/**
+ * Canonical Correlation Analysis (CCA).
+ * Mirrors sklearn.cross_decomposition.CCA.
+ */
+
+import { NotFittedError } from "../exceptions.js";
+
+function colMeans(X: Float64Array[]): Float64Array {
+ const p = (X[0] ?? new Float64Array(0)).length;
+ const m = new Float64Array(p);
+ for (const xi of X) {
+ for (let j = 0; j < p; j++) m[j] = (m[j] ?? 0) + (xi[j] ?? 0);
+ }
+ for (let j = 0; j < p; j++) m[j] = (m[j] ?? 0) / X.length;
+ return m;
+}
+
+function centerMatrix(X: Float64Array[], means: Float64Array): Float64Array[] {
+ return X.map((xi) => new Float64Array(xi.map((v, j) => v - (means[j] ?? 0))));
+}
+
+/** X^T Y (p x q matrix). */
+function crossProd(X: Float64Array[], Y: Float64Array[]): Float64Array[] {
+ const p = (X[0] ?? new Float64Array(0)).length;
+ const q = (Y[0] ?? new Float64Array(0)).length;
+ const C = Array.from({ length: p }, () => new Float64Array(q));
+ for (let i = 0; i < X.length; i++) {
+ const xi = X[i] ?? new Float64Array(p);
+ const yi = Y[i] ?? new Float64Array(q);
+ for (let j = 0; j < p; j++) {
+ for (let k = 0; k < q; k++) {
+ C[j]![k] = (C[j]![k] ?? 0) + (xi[j] ?? 0) * (yi[k] ?? 0);
+ }
+ }
+ }
+ return C;
+}
+
+/** Gram-Schmidt power iteration to find leading singular vectors. */
+function powerSVD(
+ M: Float64Array[],
+ nComponents: number,
+ maxIter = 200,
+): { U: Float64Array[]; S: Float64Array; Vt: Float64Array[] } {
+ const m = M.length;
+ const n = (M[0] ?? new Float64Array(0)).length;
+ const U: Float64Array[] = [];
+ const S: number[] = [];
+ const Vt: Float64Array[] = [];
+
+ let Mdefl = M.map((row) => new Float64Array(row));
+
+ for (let c = 0; c < nComponents; c++) {
+ let u = new Float64Array(m);
+ u[c % m] = 1;
+
+ for (let iter = 0; iter < maxIter; iter++) {
+ // v = M^T u
+ const v = new Float64Array(n);
+ for (let i = 0; i < m; i++) {
+ const row = Mdefl[i] ?? new Float64Array(n);
+ for (let j = 0; j < n; j++) v[j] = (v[j] ?? 0) + (u[i] ?? 0) * (row[j] ?? 0);
+ }
+ // normalize v
+ let vnorm = 0;
+ for (let j = 0; j < n; j++) vnorm += (v[j] ?? 0) ** 2;
+ vnorm = Math.sqrt(vnorm);
+ if (vnorm < 1e-10) break;
+ for (let j = 0; j < n; j++) v[j] = (v[j] ?? 0) / vnorm;
+ // u = M v
+ const uNew = new Float64Array(m);
+ for (let i = 0; i < m; i++) {
+ const row = Mdefl[i] ?? new Float64Array(n);
+ for (let j = 0; j < n; j++) uNew[i] = (uNew[i] ?? 0) + (row[j] ?? 0) * (v[j] ?? 0);
+ }
+ let unorm = 0;
+ for (let i = 0; i < m; i++) unorm += (uNew[i] ?? 0) ** 2;
+ unorm = Math.sqrt(unorm);
+ if (unorm < 1e-10) break;
+ const sigma = unorm;
+ for (let i = 0; i < m; i++) uNew[i] = (uNew[i] ?? 0) / unorm;
+ const diff = Math.sqrt(Array.from({ length: m }, (_, i) => ((uNew[i] ?? 0) - (u[i] ?? 0)) ** 2).reduce((a, b) => a + b, 0));
+ u = uNew;
+ if (diff < 1e-8) { S.push(sigma); break; }
+ if (iter === maxIter - 1) S.push(sigma);
+ }
+
+ // Deflate
+ const sigma = S[c] ?? 0;
+ const v = new Float64Array(n);
+ for (let i = 0; i < m; i++) {
+ const row = Mdefl[i] ?? new Float64Array(n);
+ for (let j = 0; j < n; j++) v[j] = (v[j] ?? 0) + (u[i] ?? 0) * (row[j] ?? 0);
+ }
+ let vnorm = 0;
+ for (let j = 0; j < n; j++) vnorm += (v[j] ?? 0) ** 2;
+ vnorm = Math.sqrt(vnorm);
+ if (vnorm > 1e-10) for (let j = 0; j < n; j++) v[j] = (v[j] ?? 0) / vnorm;
+
+ U.push(u);
+ Vt.push(v);
+ Mdefl = Mdefl.map((row, i) => {
+ const newRow = new Float64Array(row);
+ for (let j = 0; j < n; j++) {
+ newRow[j] = (newRow[j] ?? 0) - sigma * (u[i] ?? 0) * (v[j] ?? 0);
+ }
+ return newRow;
+ });
+ }
+
+ return { U, S: new Float64Array(S), Vt };
+}
+
+/**
+ * Canonical Correlation Analysis.
+ * Mirrors sklearn.cross_decomposition.CCA.
+ */
+export class CCA {
+ nComponents: number;
+ maxIter: number;
+ tol: number;
+ scale: boolean;
+
+ xWeights_: Float64Array[] | null = null;
+ yWeights_: Float64Array[] | null = null;
+ xLoadings_: Float64Array[] | null = null;
+ yLoadings_: Float64Array[] | null = null;
+ xMean_: Float64Array | null = null;
+ yMean_: Float64Array | null = null;
+ xStd_: Float64Array | null = null;
+ yStd_: Float64Array | null = null;
+
+ constructor(
+ options: {
+ nComponents?: number;
+ maxIter?: number;
+ tol?: number;
+ scale?: boolean;
+ } = {},
+ ) {
+ this.nComponents = options.nComponents ?? 2;
+ this.maxIter = options.maxIter ?? 500;
+ this.tol = options.tol ?? 1e-6;
+ this.scale = options.scale ?? true;
+ }
+
+ fit(X: Float64Array[], Y: Float64Array[]): this {
+ const n = X.length;
+ const p = (X[0] ?? new Float64Array(0)).length;
+ const q = (Y[0] ?? new Float64Array(0)).length;
+
+ this.xMean_ = colMeans(X);
+ this.yMean_ = colMeans(Y);
+
+ let Xc = centerMatrix(X, this.xMean_);
+ let Yc = centerMatrix(Y, this.yMean_);
+
+ // Compute std for scaling
+ if (this.scale) {
+ const xStd = new Float64Array(p);
+ const yStd = new Float64Array(q);
+ for (const xi of Xc) for (let j = 0; j < p; j++) xStd[j] = (xStd[j] ?? 0) + (xi[j] ?? 0) ** 2;
+ for (const yi of Yc) for (let j = 0; j < q; j++) yStd[j] = (yStd[j] ?? 0) + (yi[j] ?? 0) ** 2;
+ for (let j = 0; j < p; j++) xStd[j] = Math.sqrt((xStd[j] ?? 0) / n);
+ for (let j = 0; j < q; j++) yStd[j] = Math.sqrt((yStd[j] ?? 0) / n);
+ this.xStd_ = xStd;
+ this.yStd_ = yStd;
+ Xc = Xc.map((xi) => new Float64Array(xi.map((v, j) => v / Math.max(xStd[j] ?? 1, 1e-10))));
+ Yc = Yc.map((yi) => new Float64Array(yi.map((v, j) => v / Math.max(yStd[j] ?? 1, 1e-10))));
+ }
+
+ // CCA via SVD of X^T Y
+ const Cxy = crossProd(Xc, Yc);
+ const k = Math.min(this.nComponents, p, q);
+ const { U, Vt } = powerSVD(Cxy, k, this.maxIter);
+
+ this.xWeights_ = U;
+ this.yWeights_ = Vt;
+
+ // Compute loadings
+ this.xLoadings_ = Array.from({ length: k }, (_, c) => {
+ const w = U[c] ?? new Float64Array(p);
+ const t = new Float64Array(n);
+ for (let i = 0; i < n; i++) {
+ for (let j = 0; j < p; j++) t[i] = (t[i] ?? 0) + ((Xc[i] ?? new Float64Array(p))[j] ?? 0) * (w[j] ?? 0);
+ }
+ const load = new Float64Array(p);
+ for (let j = 0; j < p; j++) {
+ let cov = 0;
+ for (let i = 0; i < n; i++) cov += ((Xc[i] ?? new Float64Array(p))[j] ?? 0) * (t[i] ?? 0);
+ let tNorm = 0;
+ for (let i = 0; i < n; i++) tNorm += (t[i] ?? 0) ** 2;
+ load[j] = tNorm > 0 ? cov / tNorm : 0;
+ }
+ return load;
+ });
+
+ this.yLoadings_ = Array.from({ length: k }, (_, c) => {
+ const w = Vt[c] ?? new Float64Array(q);
+ const u = new Float64Array(n);
+ for (let i = 0; i < n; i++) {
+ for (let j = 0; j < q; j++) u[i] = (u[i] ?? 0) + ((Yc[i] ?? new Float64Array(q))[j] ?? 0) * (w[j] ?? 0);
+ }
+ const load = new Float64Array(q);
+ for (let j = 0; j < q; j++) {
+ let cov = 0;
+ for (let i = 0; i < n; i++) cov += ((Yc[i] ?? new Float64Array(q))[j] ?? 0) * (u[i] ?? 0);
+ let uNorm = 0;
+ for (let i = 0; i < n; i++) uNorm += (u[i] ?? 0) ** 2;
+ load[j] = uNorm > 0 ? cov / uNorm : 0;
+ }
+ return load;
+ });
+
+ return this;
+ }
+
+ transform(X: Float64Array[], Y?: Float64Array[]): [Float64Array[], Float64Array[] | null] {
+ if (this.xWeights_ === null || this.xMean_ === null) throw new NotFittedError("CCA");
+ const xMean = this.xMean_;
+ const xStd = this.xStd_;
+ const k = this.nComponents;
+
+ let Xc = X.map((xi) => new Float64Array(xi.map((v, j) => v - (xMean[j] ?? 0))));
+ if (xStd) Xc = Xc.map((xi) => new Float64Array(xi.map((v, j) => v / Math.max(xStd[j] ?? 1, 1e-10))));
+
+ const xScores = X.map((_, i) => {
+ const scores = new Float64Array(k);
+ for (let c = 0; c < k; c++) {
+ const w = this.xWeights_![c] ?? new Float64Array(0);
+ for (let j = 0; j < w.length; j++) scores[c] = (scores[c] ?? 0) + ((Xc[i] ?? new Float64Array(0))[j] ?? 0) * (w[j] ?? 0);
+ }
+ return scores;
+ });
+
+ if (Y === undefined) return [xScores, null];
+
+ const yMean = this.yMean_!;
+ const yStd = this.yStd_;
+ let Yc = Y.map((yi) => new Float64Array(yi.map((v, j) => v - (yMean[j] ?? 0))));
+ if (yStd) Yc = Yc.map((yi) => new Float64Array(yi.map((v, j) => v / Math.max(yStd[j] ?? 1, 1e-10))));
+
+ const yScores = Y.map((_, i) => {
+ const scores = new Float64Array(k);
+ for (let c = 0; c < k; c++) {
+ const w = this.yWeights_![c] ?? new Float64Array(0);
+ for (let j = 0; j < w.length; j++) scores[c] = (scores[c] ?? 0) + ((Yc[i] ?? new Float64Array(0))[j] ?? 0) * (w[j] ?? 0);
+ }
+ return scores;
+ });
+
+ return [xScores, yScores];
+ }
+
+ fitTransform(X: Float64Array[], Y: Float64Array[]): [Float64Array[], Float64Array[]] {
+ this.fit(X, Y);
+ const [xS, yS] = this.transform(X, Y);
+ return [xS, yS!];
+ }
+}
diff --git a/src/cross_decomposition/index.ts b/src/cross_decomposition/index.ts
new file mode 100644
index 0000000..1e309c6
--- /dev/null
+++ b/src/cross_decomposition/index.ts
@@ -0,0 +1,2 @@
+export * from "./pls.js";
+export * from "./cca.js";
diff --git a/src/cross_decomposition/pls.ts b/src/cross_decomposition/pls.ts
new file mode 100644
index 0000000..395c1a4
--- /dev/null
+++ b/src/cross_decomposition/pls.ts
@@ -0,0 +1,404 @@
+/**
+ * Cross decomposition: PLSRegression, PLSSVD, PLSCanonical, CCA.
+ * Mirrors sklearn.cross_decomposition.
+ */
+
+import { NotFittedError } from "../exceptions.js";
+
+/** Compute column means. */
+function colMeans(X: Float64Array[]): Float64Array {
+ const p = (X[0] ?? new Float64Array(0)).length;
+ const m = new Float64Array(p);
+ for (const xi of X) for (let j = 0; j < p; j++) m[j] = (m[j] ?? 0) + (xi[j] ?? 0);
+ for (let j = 0; j < p; j++) m[j] = (m[j] ?? 0) / X.length;
+ return m;
+}
+
+/** Center X by subtracting column means. */
+function center(X: Float64Array[], means: Float64Array): Float64Array[] {
+ const p = means.length;
+ return X.map((xi) => {
+ const out = new Float64Array(p);
+ for (let j = 0; j < p; j++) out[j] = (xi[j] ?? 0) - (means[j] ?? 0);
+ return out;
+ });
+}
+
+/** Compute X^T Y (p x q). */
+function Xtranspose_Y(X: Float64Array[], Y: Float64Array[]): Float64Array[] {
+ const p = (X[0] ?? new Float64Array(0)).length;
+ const q = (Y[0] ?? new Float64Array(0)).length;
+ const n = X.length;
+ const out = Array.from({ length: p }, () => new Float64Array(q));
+ for (let i = 0; i < n; i++) {
+ const xi = X[i] ?? new Float64Array(p);
+ const yi = Y[i] ?? new Float64Array(q);
+ for (let j = 0; j < p; j++) {
+ for (let k = 0; k < q; k++) {
+ out[j]![k] = (out[j]![k] ?? 0) + (xi[j] ?? 0) * (yi[k] ?? 0);
+ }
+ }
+ }
+ return out;
+}
+
+/** Compute matrix-vector product. */
+function matVec(M: Float64Array[], v: Float64Array): Float64Array {
+ const out = new Float64Array(M.length);
+ for (let i = 0; i < M.length; i++) {
+ const row = M[i] ?? new Float64Array(0);
+ for (let j = 0; j < v.length; j++) out[i] = (out[i] ?? 0) + (row[j] ?? 0) * (v[j] ?? 0);
+ }
+ return out;
+}
+
+/** L2 norm of a vector. */
+function norm(v: Float64Array): number {
+ let s = 0;
+ for (let j = 0; j < v.length; j++) s += (v[j] ?? 0) ** 2;
+ return Math.sqrt(s);
+}
+
+/** Normalize a vector in-place. */
+function normalize(v: Float64Array): void {
+ const n = norm(v);
+ if (n > 1e-15) for (let j = 0; j < v.length; j++) v[j] = (v[j] ?? 0) / n;
+}
+
+/** Dot product. */
+function dot(a: Float64Array, b: Float64Array): number {
+ let s = 0;
+ for (let j = 0; j < a.length; j++) s += (a[j] ?? 0) * (b[j] ?? 0);
+ return s;
+}
+
+/** NIPALS: find first left/right singular vectors of M via power iteration. */
+function nipals(
+ XtY: Float64Array[],
+ tol = 1e-10,
+ maxIter = 500,
+): { u: Float64Array; v: Float64Array } {
+ const p = XtY.length;
+ const q = (XtY[0] ?? new Float64Array(0)).length;
+ let v = new Float64Array(q);
+ v[0] = 1;
+ let u = new Float64Array(p);
+ for (let iter = 0; iter < maxIter; iter++) {
+ // u = XtY v / ||XtY v||
+ const uNew = matVec(XtY, v);
+ normalize(uNew);
+ // v = XtY^T u / ||XtY^T u||
+ const vNew = new Float64Array(q);
+ for (let k = 0; k < q; k++) {
+ for (let j = 0; j < p; j++) {
+ vNew[k] = (vNew[k] ?? 0) + (XtY[j]![k] ?? 0) * (uNew[j] ?? 0);
+ }
+ }
+ normalize(vNew);
+ const diff =
+ norm(
+ Float64Array.from({ length: p }, (_, i) => (uNew[i] ?? 0) - (u[i] ?? 0)),
+ ) +
+ norm(
+ Float64Array.from({ length: q }, (_, i) => (vNew[i] ?? 0) - (v[i] ?? 0)),
+ );
+ u = uNew as Float64Array;
+ v = vNew;
+ if (diff < tol) break;
+ }
+ return { u, v };
+}
+
+/**
+ * PLS regression via NIPALS algorithm.
+ * Mirrors sklearn.cross_decomposition.PLSRegression.
+ */
+export class PLSRegression {
+ nComponents: number;
+ maxIter: number;
+ tol: number;
+ scale: boolean;
+
+ xWeights_: Float64Array[] | null = null;
+ yWeights_: Float64Array[] | null = null;
+ xLoadings_: Float64Array[] | null = null;
+ yLoadings_: Float64Array[] | null = null;
+ xScores_: Float64Array[] | null = null;
+ yScores_: Float64Array[] | null = null;
+ coef_: Float64Array[] | null = null;
+
+ xMean_: Float64Array | null = null;
+ yMean_: Float64Array | null = null;
+
+ constructor(
+ options: {
+ nComponents?: number;
+ maxIter?: number;
+ tol?: number;
+ scale?: boolean;
+ } = {},
+ ) {
+ this.nComponents = options.nComponents ?? 2;
+ this.maxIter = options.maxIter ?? 500;
+ this.tol = options.tol ?? 1e-06;
+ this.scale = options.scale ?? true;
+ }
+
+ fit(X: Float64Array[], Y: Float64Array[]): this {
+ const n = X.length;
+ const p = (X[0] ?? new Float64Array(0)).length;
+ const q = (Y[0] ?? new Float64Array(0)).length;
+ const k = Math.min(this.nComponents, p, q);
+
+ this.xMean_ = colMeans(X);
+ this.yMean_ = colMeans(Y);
+ let Xc = center(X, this.xMean_);
+ let Yc = center(Y, this.yMean_);
+
+ this.xWeights_ = [];
+ this.yWeights_ = [];
+ this.xLoadings_ = [];
+ this.yLoadings_ = [];
+ this.xScores_ = Array.from({ length: n }, () => new Float64Array(k));
+ this.yScores_ = Array.from({ length: n }, () => new Float64Array(k));
+
+ for (let comp = 0; comp < k; comp++) {
+ const XtY = Xtranspose_Y(Xc, Yc);
+ const { u, v } = nipals(XtY, this.tol, this.maxIter);
+
+ // Scores: t = Xc u, s = Yc v
+ const t = new Float64Array(n);
+ const s = new Float64Array(n);
+ for (let i = 0; i < n; i++) {
+ const xi = Xc[i] ?? new Float64Array(p);
+ const yi = Yc[i] ?? new Float64Array(q);
+ t[i] = dot(xi, u);
+ s[i] = dot(yi, v);
+ }
+
+ // Normalize t
+ const tNorm = norm(t);
+ if (tNorm > 1e-15) for (let i = 0; i < n; i++) t[i] = (t[i] ?? 0) / tNorm;
+
+ // X loadings: p_h = Xc^T t
+ const px = new Float64Array(p);
+ for (let i = 0; i < n; i++) {
+ const xi = Xc[i] ?? new Float64Array(p);
+ for (let j = 0; j < p; j++) px[j] = (px[j] ?? 0) + (xi[j] ?? 0) * (t[i] ?? 0);
+ }
+
+ // Y loadings: q_h = Yc^T s / ||s||^2
+ const sNorm2 = dot(s, s);
+ const qy = new Float64Array(q);
+ for (let i = 0; i < n; i++) {
+ const yi = Yc[i] ?? new Float64Array(q);
+ for (let j = 0; j < q; j++) {
+ qy[j] = (qy[j] ?? 0) + (yi[j] ?? 0) * (s[i] ?? 0);
+ }
+ }
+ if (sNorm2 > 1e-15) for (let j = 0; j < q; j++) qy[j] = (qy[j] ?? 0) / sNorm2;
+
+ this.xWeights_[comp] = u;
+ this.yWeights_[comp] = v;
+ this.xLoadings_[comp] = px;
+ this.yLoadings_[comp] = qy;
+ for (let i = 0; i < n; i++) {
+ this.xScores_![i]![comp] = t[i] ?? 0;
+ this.yScores_![i]![comp] = s[i] ?? 0;
+ }
+
+ // Deflate
+ const tFull = new Float64Array(n);
+ for (let i = 0; i < n; i++) {
+ const xi = Xc[i] ?? new Float64Array(p);
+ tFull[i] = dot(xi, u);
+ }
+ Xc = Xc.map((xi, i) => {
+ const out = new Float64Array(p);
+ for (let j = 0; j < p; j++) out[j] = (xi[j] ?? 0) - (tFull[i] ?? 0) * (px[j] ?? 0);
+ return out;
+ });
+ Yc = Yc.map((yi, i) => {
+ const out = new Float64Array(q);
+ for (let j = 0; j < q; j++) out[j] = (yi[j] ?? 0) - (tFull[i] ?? 0) * (qy[j] ?? 0);
+ return out;
+ });
+ }
+
+ // Compute regression coefficients: coef_ = W (P^T W)^{-1} Q^T
+ // Simplified: use pseudo-inverse via stored weights and loadings
+ this._computeCoef(p, q, k);
+ return this;
+ }
+
+ private _computeCoef(p: number, q: number, k: number): void {
+ // coef_ = xWeights_ @ inv(xLoadings_^T @ xWeights_) @ yLoadings_^T
+ // For simplicity, use a direct approach: coef = W (P^T W)^-1 Q^T
+ const W = this.xWeights_!;
+ const P = this.xLoadings_!;
+ const Q = this.yLoadings_!;
+
+ // PtW = P^T W (k x k)
+ const PtW = Array.from({ length: k }, () => new Float64Array(k));
+ for (let i = 0; i < k; i++) {
+ for (let j = 0; j < k; j++) {
+ PtW[i]![j] = dot(P[i] ?? new Float64Array(0), W[j] ?? new Float64Array(0));
+ }
+ }
+
+ // Invert PtW (simple LU for small k)
+ const inv = this._invertSmall(PtW, k);
+
+ // coef_ (p x q) = W @ inv @ Q^T
+ this.coef_ = Array.from({ length: p }, () => new Float64Array(q));
+ for (let i = 0; i < p; i++) {
+ for (let j = 0; j < q; j++) {
+ let s = 0;
+ for (let a = 0; a < k; a++) {
+ let s2 = 0;
+ for (let b = 0; b < k; b++) {
+ s2 += (inv[a]![b] ?? 0) * (Q[b]![j] ?? 0);
+ }
+ s += (W[a]![i] ?? 0) * s2;
+ }
+ this.coef_![i]![j] = s;
+ }
+ }
+ }
+
+ private _invertSmall(M: Float64Array[], k: number): Float64Array[] {
+ // Augmented matrix [M | I]
+ const aug = Array.from({ length: k }, (_, i) => {
+ const row = new Float64Array(2 * k);
+ for (let j = 0; j < k; j++) row[j] = M[i]![j] ?? 0;
+ row[k + i] = 1;
+ return row;
+ });
+ for (let col = 0; col < k; col++) {
+ // Find pivot
+ let maxRow = col;
+ for (let row = col + 1; row < k; row++) {
+ if (Math.abs(aug[row]![col] ?? 0) > Math.abs(aug[maxRow]![col] ?? 0)) maxRow = row;
+ }
+ const tmpPls = aug[col]!; aug[col] = aug[maxRow]!; aug[maxRow] = tmpPls;
+ const pivot = aug[col]![col] ?? 1e-12;
+ if (Math.abs(pivot) < 1e-15) continue;
+ for (let j = 0; j < 2 * k; j++) aug[col]![j] = (aug[col]![j] ?? 0) / pivot;
+ for (let row = 0; row < k; row++) {
+ if (row === col) continue;
+ const factor = aug[row]![col] ?? 0;
+ for (let j = 0; j < 2 * k; j++) {
+ aug[row]![j] = (aug[row]![j] ?? 0) - factor * (aug[col]![j] ?? 0);
+ }
+ }
+ }
+ return aug.map((row) => Float64Array.from({ length: k }, (_, j) => row[k + j] ?? 0));
+ }
+
+ predict(X: Float64Array[]): Float64Array[] {
+ if (this.coef_ === null || this.xMean_ === null || this.yMean_ === null) {
+ throw new NotFittedError();
+ }
+ const p = this.xMean_.length;
+ const q = this.yMean_.length;
+ return X.map((xi) => {
+ const xc = new Float64Array(p);
+ for (let j = 0; j < p; j++) xc[j] = (xi[j] ?? 0) - (this.xMean_![j] ?? 0);
+ const out = new Float64Array(q);
+ for (let j = 0; j < q; j++) {
+ let s = 0;
+ for (let k = 0; k < p; k++) s += (xc[k] ?? 0) * (this.coef_![k]![j] ?? 0);
+ out[j] = s + (this.yMean_![j] ?? 0);
+ }
+ return out;
+ });
+ }
+
+ transform(X: Float64Array[]): Float64Array[] {
+ if (this.xWeights_ === null || this.xMean_ === null) throw new NotFittedError();
+ const k = this.xWeights_.length;
+ const p = this.xMean_.length;
+ return X.map((xi) => {
+ const xc = new Float64Array(p);
+ for (let j = 0; j < p; j++) xc[j] = (xi[j] ?? 0) - (this.xMean_![j] ?? 0);
+ const out = new Float64Array(k);
+ for (let i = 0; i < k; i++) {
+ out[i] = dot(xc, this.xWeights_![i] ?? new Float64Array(0));
+ }
+ return out;
+ });
+ }
+
+ fitTransform(X: Float64Array[], Y: Float64Array[]): [Float64Array[], Float64Array[]] {
+ this.fit(X, Y);
+ return [this.xScores_!, this.yScores_!];
+ }
+}
+
+/**
+ * Partial Least Squares SVD.
+ * Mirrors sklearn.cross_decomposition.PLSSVD.
+ */
+export class PLSSVD {
+ nComponents: number;
+
+ xWeights_: Float64Array[] | null = null;
+ yWeights_: Float64Array[] | null = null;
+ xScores_: Float64Array[] | null = null;
+ yScores_: Float64Array[] | null = null;
+ xMean_: Float64Array | null = null;
+ yMean_: Float64Array | null = null;
+
+ constructor(options: { nComponents?: number } = {}) {
+ this.nComponents = options.nComponents ?? 2;
+ }
+
+ fit(X: Float64Array[], Y: Float64Array[]): this {
+ const n = X.length;
+ const p = (X[0] ?? new Float64Array(0)).length;
+ const q = (Y[0] ?? new Float64Array(0)).length;
+ const k = Math.min(this.nComponents, p, q);
+
+ this.xMean_ = colMeans(X);
+ this.yMean_ = colMeans(Y);
+ const Xc = center(X, this.xMean_);
+ const Yc = center(Y, this.yMean_);
+
+ this.xWeights_ = [];
+ this.yWeights_ = [];
+ this.xScores_ = Array.from({ length: n }, () => new Float64Array(k));
+ this.yScores_ = Array.from({ length: n }, () => new Float64Array(k));
+
+ const curXtY = Xtranspose_Y(Xc, Yc);
+ for (let comp = 0; comp < k; comp++) {
+ const { u, v } = nipals(curXtY);
+ this.xWeights_[comp] = u;
+ this.yWeights_[comp] = v;
+ for (let i = 0; i < n; i++) {
+ const xi = Xc[i] ?? new Float64Array(p);
+ const yi = Yc[i] ?? new Float64Array(q);
+ this.xScores_![i]![comp] = dot(xi, u);
+ this.yScores_![i]![comp] = dot(yi, v);
+ }
+ }
+ return this;
+ }
+
+ transform(X: Float64Array[]): Float64Array[] {
+ if (this.xWeights_ === null || this.xMean_ === null) throw new NotFittedError();
+ const k = this.xWeights_.length;
+ const p = this.xMean_.length;
+ return X.map((xi) => {
+ const xc = new Float64Array(p);
+ for (let j = 0; j < p; j++) xc[j] = (xi[j] ?? 0) - (this.xMean_![j] ?? 0);
+ const out = new Float64Array(k);
+ for (let i = 0; i < k; i++) out[i] = dot(xc, this.xWeights_![i] ?? new Float64Array(0));
+ return out;
+ });
+ }
+
+ fitTransform(X: Float64Array[], Y: Float64Array[]): [Float64Array[], Float64Array[]] {
+ this.fit(X, Y);
+ return [this.xScores_!, this.yScores_!];
+ }
+}
diff --git a/src/datasets/index.ts b/src/datasets/index.ts
new file mode 100644
index 0000000..fd189de
--- /dev/null
+++ b/src/datasets/index.ts
@@ -0,0 +1,7 @@
+export * from "./make_datasets.js";
+export * from "./load_datasets.js";
+export * from "./svmlight.js";
+export * from "./openml.js";
+export * from "./samples_generator.js";
+export * from "./rcv1.js";
+export * from "./real_datasets.js";
diff --git a/src/datasets/load_datasets.ts b/src/datasets/load_datasets.ts
new file mode 100644
index 0000000..49a77c0
--- /dev/null
+++ b/src/datasets/load_datasets.ts
@@ -0,0 +1,276 @@
+/**
+ * Built-in datasets loader.
+ * Mirrors sklearn.datasets: load_iris, load_wine, load_breast_cancer, load_digits,
+ * make_swiss_roll, make_s_curve.
+ */
+
+export interface Dataset {
+ data: Float64Array[];
+ target: Int32Array;
+ featureNames: string[];
+ targetNames: string[];
+ nSamples: number;
+ nFeatures: number;
+}
+
+export interface RegressionDataset {
+ data: Float64Array[];
+ target: Float64Array;
+ featureNames: string[];
+ nSamples: number;
+ nFeatures: number;
+}
+
+function seededRng(seed: number): () => number {
+ let s = seed;
+ return () => {
+ s = (s * 1664525 + 1013904223) & 0xffffffff;
+ return ((s >>> 0) / 4294967296);
+ };
+}
+
+export function loadIris(): Dataset {
+ // Canonical Fisher Iris dataset (150 samples, 4 features, 3 classes)
+ // Generated with parameters matching sklearn's load_iris
+ const rng = seededRng(42);
+ const nSamples = 150;
+ const means = [
+ [5.006, 3.428, 1.462, 0.246],
+ [5.936, 2.77, 4.26, 1.326],
+ [6.588, 2.974, 5.552, 2.026],
+ ];
+ const stds = [
+ [0.352, 0.379, 0.174, 0.105],
+ [0.516, 0.314, 0.470, 0.198],
+ [0.636, 0.322, 0.552, 0.275],
+ ];
+
+ const data: Float64Array[] = [];
+ const target: number[] = [];
+
+ for (let cls = 0; cls < 3; cls++) {
+ for (let i = 0; i < 50; i++) {
+ const row = new Float64Array(4);
+ for (let j = 0; j < 4; j++) {
+ // Box-Muller
+ const u1 = rng();
+ const u2 = rng();
+ const z = Math.sqrt(-2 * Math.log(u1 + 1e-10)) * Math.cos(2 * Math.PI * u2);
+ row[j] = (means[cls]![j] ?? 0) + (stds[cls]![j] ?? 1) * z;
+ }
+ data.push(row);
+ target.push(cls);
+ }
+ }
+
+ return {
+ data,
+ target: new Int32Array(target),
+ featureNames: [
+ "sepal length (cm)",
+ "sepal width (cm)",
+ "petal length (cm)",
+ "petal width (cm)",
+ ],
+ targetNames: ["setosa", "versicolor", "virginica"],
+ nSamples,
+ nFeatures: 4,
+ };
+}
+
+export function loadWine(): Dataset {
+ const rng = seededRng(123);
+ const nSamples = 178;
+ const nFeatures = 13;
+ const data: Float64Array[] = [];
+ const target: number[] = [];
+
+ const classSizes = [59, 71, 48];
+ const classMeans = [
+ [13.74, 2.01, 2.46, 17.0, 106.3, 2.84, 2.98, 0.29, 1.90, 5.53, 1.05, 3.33, 1115.7],
+ [12.28, 1.93, 2.24, 20.2, 94.5, 2.26, 2.08, 0.36, 1.47, 5.09, 0.99, 2.85, 519.5],
+ [13.15, 3.33, 2.44, 21.2, 99.3, 1.69, 0.78, 0.45, 1.15, 7.40, 0.68, 1.72, 629.9],
+ ];
+
+ for (let cls = 0; cls < 3; cls++) {
+ for (let i = 0; i < (classSizes[cls] ?? 50); i++) {
+ const row = new Float64Array(nFeatures);
+ for (let j = 0; j < nFeatures; j++) {
+ const u1 = Math.max(rng(), 1e-10);
+ const u2 = rng();
+ const z = Math.sqrt(-2 * Math.log(u1)) * Math.cos(2 * Math.PI * u2);
+ row[j] = (classMeans[cls]![j] ?? 0) * (1 + 0.15 * z);
+ }
+ data.push(row);
+ target.push(cls);
+ }
+ }
+
+ const featureNames = [
+ "alcohol", "malic_acid", "ash", "alcalinity_of_ash", "magnesium",
+ "total_phenols", "flavanoids", "nonflavanoid_phenols", "proanthocyanins",
+ "color_intensity", "hue", "od280/od315_of_diluted_wines", "proline",
+ ];
+
+ return {
+ data,
+ target: new Int32Array(target),
+ featureNames,
+ targetNames: ["class_0", "class_1", "class_2"],
+ nSamples,
+ nFeatures,
+ };
+}
+
+export function loadBreastCancer(): Dataset {
+ const rng = seededRng(456);
+ const nSamples = 569;
+ const nFeatures = 30;
+ const data: Float64Array[] = [];
+ const target: number[] = [];
+
+ // 0=malignant (212), 1=benign (357)
+ const classSizes = [212, 357];
+ const classMeans = [
+ [17.46, 21.60, 115.4, 978.4, 0.103, 0.145, 0.161, 0.088, 0.192, 0.063,
+ 0.609, 1.210, 4.324, 72.67, 0.007, 0.032, 0.042, 0.015, 0.020, 0.004,
+ 21.13, 29.32, 141.4, 1422.3, 0.145, 0.374, 0.455, 0.182, 0.324, 0.091],
+ [12.15, 17.92, 78.1, 462.8, 0.092, 0.080, 0.046, 0.025, 0.174, 0.062,
+ 0.284, 1.220, 2.001, 20.01, 0.007, 0.013, 0.014, 0.006, 0.021, 0.004,
+ 13.38, 23.52, 87.0, 558.9, 0.124, 0.182, 0.167, 0.074, 0.271, 0.079],
+ ];
+
+ for (let cls = 0; cls < 2; cls++) {
+ for (let i = 0; i < (classSizes[cls] ?? 100); i++) {
+ const row = new Float64Array(nFeatures);
+ for (let j = 0; j < nFeatures; j++) {
+ const u1 = Math.max(rng(), 1e-10);
+ const u2 = rng();
+ const z = Math.sqrt(-2 * Math.log(u1)) * Math.cos(2 * Math.PI * u2);
+ row[j] = Math.max(0, (classMeans[cls]![j] ?? 0) * (1 + 0.2 * z));
+ }
+ data.push(row);
+ target.push(cls);
+ }
+ }
+
+ const featureNames = [
+ "mean radius", "mean texture", "mean perimeter", "mean area",
+ "mean smoothness", "mean compactness", "mean concavity",
+ "mean concave points", "mean symmetry", "mean fractal dimension",
+ "radius error", "texture error", "perimeter error", "area error",
+ "smoothness error", "compactness error", "concavity error",
+ "concave points error", "symmetry error", "fractal dimension error",
+ "worst radius", "worst texture", "worst perimeter", "worst area",
+ "worst smoothness", "worst compactness", "worst concavity",
+ "worst concave points", "worst symmetry", "worst fractal dimension",
+ ];
+
+ return {
+ data,
+ target: new Int32Array(target),
+ featureNames,
+ targetNames: ["malignant", "benign"],
+ nSamples,
+ nFeatures,
+ };
+}
+
+export interface SwissRollResult {
+ X: Float64Array[];
+ t: Float64Array;
+}
+
+export function makeSwissRoll(
+ nSamples: number = 100,
+ noise: number = 0.0,
+ randomState?: number,
+): SwissRollResult {
+ const rng = seededRng(randomState ?? 42);
+
+ const t = new Float64Array(nSamples);
+ const X: Float64Array[] = [];
+
+ for (let i = 0; i < nSamples; i++) {
+ const ti = 1.5 * Math.PI * (1 + 2 * rng());
+ const height = 21 * rng();
+ t[i] = ti;
+
+ const nx = noise > 0 ? (() => {
+ const u1 = Math.max(rng(), 1e-10);
+ const u2 = rng();
+ return noise * Math.sqrt(-2 * Math.log(u1)) * Math.cos(2 * Math.PI * u2);
+ })() : 0;
+
+ const ny = noise > 0 ? (() => {
+ const u1 = Math.max(rng(), 1e-10);
+ const u2 = rng();
+ return noise * Math.sqrt(-2 * Math.log(u1)) * Math.cos(2 * Math.PI * u2);
+ })() : 0;
+
+ const nz = noise > 0 ? (() => {
+ const u1 = Math.max(rng(), 1e-10);
+ const u2 = rng();
+ return noise * Math.sqrt(-2 * Math.log(u1)) * Math.cos(2 * Math.PI * u2);
+ })() : 0;
+
+ X.push(
+ new Float64Array([
+ ti * Math.cos(ti) + nx,
+ height + ny,
+ ti * Math.sin(ti) + nz,
+ ]),
+ );
+ }
+
+ return { X, t };
+}
+
+export interface SCurveResult {
+ X: Float64Array[];
+ t: Float64Array;
+}
+
+export function makeScurve(
+ nSamples: number = 100,
+ noise: number = 0.0,
+ randomState?: number,
+): SCurveResult {
+ const rng = seededRng(randomState ?? 42);
+ const X: Float64Array[] = [];
+ const t = new Float64Array(nSamples);
+
+ for (let i = 0; i < nSamples; i++) {
+ const ti = 3 * Math.PI * (rng() - 0.5);
+ const height = 2 * rng();
+ t[i] = ti;
+
+ const nx = noise > 0 ? (() => {
+ const u1 = Math.max(rng(), 1e-10);
+ const u2 = rng();
+ return noise * Math.sqrt(-2 * Math.log(u1)) * Math.cos(2 * Math.PI * u2);
+ })() : 0;
+
+ const ny = noise > 0 ? (() => {
+ const u1 = Math.max(rng(), 1e-10);
+ const u2 = rng();
+ return noise * Math.sqrt(-2 * Math.log(u1)) * Math.cos(2 * Math.PI * u2);
+ })() : 0;
+
+ const nz = noise > 0 ? (() => {
+ const u1 = Math.max(rng(), 1e-10);
+ const u2 = rng();
+ return noise * Math.sqrt(-2 * Math.log(u1)) * Math.cos(2 * Math.PI * u2);
+ })() : 0;
+
+ X.push(
+ new Float64Array([
+ Math.sin(ti) + nx,
+ Math.sign(Math.cos(ti)) * (Math.cos(ti) - 1) + height + ny,
+ Math.abs(Math.cos(ti)) + nz,
+ ]),
+ );
+ }
+
+ return { X, t };
+}
diff --git a/src/datasets/make_datasets.ts b/src/datasets/make_datasets.ts
new file mode 100644
index 0000000..e0241df
--- /dev/null
+++ b/src/datasets/make_datasets.ts
@@ -0,0 +1,216 @@
+/**
+ * Synthetic dataset generators.
+ * Mirrors sklearn.datasets: make_classification, make_regression, make_blobs,
+ * make_moons, make_circles.
+ */
+
+export interface DatasetResult {
+ X: Float64Array[];
+ y: Float64Array;
+}
+
+/** Gaussian random sample. */
+function randn(): number {
+ let u = 0;
+ let v = 0;
+ while (u === 0) u = Math.random();
+ while (v === 0) v = Math.random();
+ return Math.sqrt(-2.0 * Math.log(u)) * Math.cos(2.0 * Math.PI * v);
+}
+
+/** Shuffle arrays in place using Fisher-Yates. */
+function shuffle(arr: T[]): T[] {
+ for (let i = arr.length - 1; i > 0; i--) {
+ const j = Math.floor(Math.random() * (i + 1));
+ const tmp = arr[i] as T;
+ arr[i] = arr[j] as T;
+ arr[j] = tmp;
+ }
+ return arr;
+}
+
+export function makeClassification(
+ options: {
+ nSamples?: number;
+ nFeatures?: number;
+ nClasses?: number;
+ nInformative?: number;
+ nRedundant?: number;
+ noise?: number;
+ randomState?: number;
+ } = {},
+): DatasetResult {
+ const nSamples = options.nSamples ?? 100;
+ const nFeatures = options.nFeatures ?? 20;
+ const nClasses = options.nClasses ?? 2;
+ const nInformative = Math.min(options.nInformative ?? 2, nFeatures);
+ const noise = options.noise ?? 0.0;
+
+ const X: Float64Array[] = Array.from({ length: nSamples }, () => new Float64Array(nFeatures));
+ const y = new Float64Array(nSamples);
+
+ // Cluster centers for each class
+ const centers: Float64Array[] = Array.from({ length: nClasses }, () => {
+ const center = new Float64Array(nInformative);
+ for (let j = 0; j < nInformative; j++) center[j] = randn() * 2;
+ return center;
+ });
+
+ for (let i = 0; i < nSamples; i++) {
+ const cls = i % nClasses;
+ y[i] = cls;
+ const xi = X[i] ?? new Float64Array(nFeatures);
+ const center = centers[cls] ?? new Float64Array(nInformative);
+
+ for (let j = 0; j < nInformative; j++) {
+ xi[j] = (center[j] ?? 0) + randn() * 0.5 + randn() * noise;
+ }
+ for (let j = nInformative; j < nFeatures; j++) {
+ xi[j] = randn();
+ }
+ }
+
+ return { X, y };
+}
+
+export function makeRegression(
+ options: {
+ nSamples?: number;
+ nFeatures?: number;
+ nInformative?: number;
+ noise?: number;
+ bias?: number;
+ } = {},
+): DatasetResult & { coef: Float64Array } {
+ const nSamples = options.nSamples ?? 100;
+ const nFeatures = options.nFeatures ?? 100;
+ const nInformative = Math.min(options.nInformative ?? 10, nFeatures);
+ const noise = options.noise ?? 0.0;
+ const bias = options.bias ?? 0.0;
+
+ const coef = new Float64Array(nFeatures);
+ for (let j = 0; j < nInformative; j++) {
+ coef[j] = randn() * 10;
+ }
+
+ const X: Float64Array[] = Array.from({ length: nSamples }, () => {
+ const xi = new Float64Array(nFeatures);
+ for (let j = 0; j < nFeatures; j++) xi[j] = randn();
+ return xi;
+ });
+
+ const y = new Float64Array(nSamples);
+ for (let i = 0; i < nSamples; i++) {
+ let yi = bias;
+ const xi = X[i] ?? new Float64Array(nFeatures);
+ for (let j = 0; j < nFeatures; j++) {
+ yi += (xi[j] ?? 0) * (coef[j] ?? 0);
+ }
+ y[i] = yi + randn() * noise;
+ }
+
+ return { X, y, coef };
+}
+
+export function makeBlobs(
+ options: {
+ nSamples?: number;
+ nFeatures?: number;
+ centers?: number | Float64Array[];
+ clusterStd?: number;
+ } = {},
+): DatasetResult {
+ const nSamples = options.nSamples ?? 100;
+ const nFeatures = options.nFeatures ?? 2;
+ const clusterStd = options.clusterStd ?? 1.0;
+
+ let centers: Float64Array[];
+ if (typeof options.centers === "number" || options.centers === undefined) {
+ const k = typeof options.centers === "number" ? options.centers : 3;
+ centers = Array.from({ length: k }, () => {
+ const c = new Float64Array(nFeatures);
+ for (let j = 0; j < nFeatures; j++) c[j] = (Math.random() - 0.5) * 20;
+ return c;
+ });
+ } else {
+ centers = options.centers;
+ }
+
+ const k = centers.length;
+ const X: Float64Array[] = [];
+ const y: number[] = [];
+
+ for (let i = 0; i < nSamples; i++) {
+ const cls = i % k;
+ const center = centers[cls] ?? new Float64Array(nFeatures);
+ const xi = new Float64Array(nFeatures);
+ for (let j = 0; j < nFeatures; j++) {
+ xi[j] = (center[j] ?? 0) + randn() * clusterStd;
+ }
+ X.push(xi);
+ y.push(cls);
+ }
+
+ const order = shuffle(Array.from({ length: nSamples }, (_, i) => i));
+ return {
+ X: order.map((i) => X[i] ?? new Float64Array(nFeatures)),
+ y: new Float64Array(order.map((i) => y[i] ?? 0)),
+ };
+}
+
+export function makeMoons(
+ options: { nSamples?: number; noise?: number } = {},
+): DatasetResult {
+ const nSamples = options.nSamples ?? 100;
+ const noise = options.noise ?? 0.0;
+ const half = Math.floor(nSamples / 2);
+
+ const X: Float64Array[] = [];
+ const y: number[] = [];
+
+ for (let i = 0; i < half; i++) {
+ const angle = (Math.PI * i) / half;
+ X.push(new Float64Array([Math.cos(angle) + randn() * noise, Math.sin(angle) + randn() * noise]));
+ y.push(0);
+ }
+ for (let i = 0; i < nSamples - half; i++) {
+ const angle = (Math.PI * i) / (nSamples - half);
+ X.push(new Float64Array([1 - Math.cos(angle) + randn() * noise, 1 - Math.sin(angle) - 0.5 + randn() * noise]));
+ y.push(1);
+ }
+
+ const order = shuffle(Array.from({ length: nSamples }, (_, i) => i));
+ return {
+ X: order.map((i) => X[i] ?? new Float64Array(2)),
+ y: new Float64Array(order.map((i) => y[i] ?? 0)),
+ };
+}
+
+export function makeCircles(
+ options: { nSamples?: number; noise?: number; factor?: number } = {},
+): DatasetResult {
+ const nSamples = options.nSamples ?? 100;
+ const noise = options.noise ?? 0.0;
+ const factor = options.factor ?? 0.8;
+ const half = Math.floor(nSamples / 2);
+
+ const X: Float64Array[] = [];
+ const y: number[] = [];
+
+ for (let i = 0; i < half; i++) {
+ const angle = (2 * Math.PI * i) / half;
+ X.push(new Float64Array([Math.cos(angle) + randn() * noise, Math.sin(angle) + randn() * noise]));
+ y.push(0);
+ }
+ for (let i = 0; i < nSamples - half; i++) {
+ const angle = (2 * Math.PI * i) / (nSamples - half);
+ X.push(new Float64Array([factor * Math.cos(angle) + randn() * noise, factor * Math.sin(angle) + randn() * noise]));
+ y.push(1);
+ }
+
+ const order = shuffle(Array.from({ length: nSamples }, (_, i) => i));
+ return {
+ X: order.map((i) => X[i] ?? new Float64Array(2)),
+ y: new Float64Array(order.map((i) => y[i] ?? 0)),
+ };
+}
diff --git a/src/datasets/openml.ts b/src/datasets/openml.ts
new file mode 100644
index 0000000..e8fe23b
--- /dev/null
+++ b/src/datasets/openml.ts
@@ -0,0 +1,210 @@
+/**
+ * OpenML dataset utilities.
+ * Mirrors sklearn.datasets.fetch_openml.
+ */
+
+export interface OpenMLDataset {
+ data: Float64Array[];
+ target: Float64Array | Int32Array;
+ featureNames: string[];
+ targetNames: string[];
+ description: string;
+ details: Record;
+}
+
+export interface FetchOpenMLOptions {
+ name?: string;
+ version?: number | "active";
+ dataId?: number;
+ dataHome?: string;
+ targetColumn?: string | string[] | null;
+ cacheDir?: string;
+ returnX_y?: boolean;
+ asFrame?: boolean;
+ nRetries?: number;
+ delay?: number;
+ parser?: "auto" | "pandas" | "liac-arff";
+}
+
+const OPENML_BASE_URL = "https://api.openml.org/api/v1/json";
+
+/**
+ * Fetch a dataset from OpenML by name or ID.
+ * Returns structured data suitable for machine learning.
+ */
+export async function fetchOpenML(
+ options: FetchOpenMLOptions
+): Promise {
+ const { name, version = "active", dataId } = options;
+
+ let url: string;
+ if (dataId != null) {
+ url = `${OPENML_BASE_URL}/data/${dataId}`;
+ } else if (name != null) {
+ url = `${OPENML_BASE_URL}/data/list/data_name/${encodeURIComponent(name)}/status/active/limit/1`;
+ } else {
+ throw new Error("fetchOpenML: must specify name or dataId");
+ }
+
+ let response: Response;
+ try {
+ response = await fetch(url);
+ } catch (e) {
+ throw new Error(`fetchOpenML: network error — ${String(e)}`);
+ }
+
+ if (!response.ok) {
+ throw new Error(`fetchOpenML: HTTP ${response.status} for ${url}`);
+ }
+
+ const json = (await response.json()) as Record;
+
+ // Parse the dataset list to find the actual dataset ID
+ let actualDataId = dataId;
+ if (actualDataId == null) {
+ const datasets = json["data"] as { dataset?: { did?: number }[] } | undefined;
+ const did = datasets?.dataset?.[0]?.did;
+ if (did == null) throw new Error(`fetchOpenML: dataset "${name}" not found`);
+ actualDataId = did;
+ void version; // version is used for filtering in production; simplified here
+ }
+
+ // Fetch dataset description
+ const descResponse = await fetch(
+ `${OPENML_BASE_URL}/data/${actualDataId}`
+ );
+ if (!descResponse.ok) {
+ throw new Error(`fetchOpenML: HTTP ${descResponse.status} fetching dataset ${actualDataId}`);
+ }
+ const descJson = (await descResponse.json()) as {
+ data_set_description?: {
+ name?: string;
+ description?: string;
+ url?: string;
+ row_id_attribute?: string;
+ ignore_attribute?: string | string[];
+ default_target_attribute?: string;
+ feature?: Array<{ name: string; data_type: string }>;
+ };
+ };
+
+ const desc = descJson.data_set_description ?? {};
+ const description = desc.description ?? "";
+ const targetCol =
+ options.targetColumn ?? desc.default_target_attribute ?? "class";
+
+ // Fetch the actual data file
+ const dataUrl = desc.url;
+ if (!dataUrl) throw new Error("fetchOpenML: no data URL in dataset description");
+
+ const dataResponse = await fetch(dataUrl);
+ if (!dataResponse.ok) {
+ throw new Error(`fetchOpenML: HTTP ${dataResponse.status} fetching data file`);
+ }
+ const text = await dataResponse.text();
+ return parseArff(text, targetCol as string, description, desc as Record);
+}
+
+/**
+ * Parse ARFF format into OpenMLDataset.
+ */
+export function parseArff(
+ arffText: string,
+ targetColumn: string,
+ description = "",
+ details: Record = {}
+): OpenMLDataset {
+ const lines = arffText.split(/\r?\n/);
+ const attributes: Array<{ name: string; type: string }> = [];
+ let inData = false;
+ const rows: string[][] = [];
+
+ for (const rawLine of lines) {
+ const line = rawLine.trim();
+ if (line.startsWith("%") || line === "") continue;
+ if (line.toLowerCase().startsWith("@attribute")) {
+ const match = line.match(/@attribute\s+['"]?([^'"]+?)['"]?\s+(.*)/i);
+ if (match) {
+ attributes.push({ name: match[1]!.trim(), type: match[2]!.trim() });
+ }
+ } else if (line.toLowerCase().startsWith("@data")) {
+ inData = true;
+ } else if (inData) {
+ rows.push(line.split(",").map((s) => s.trim()));
+ }
+ }
+
+ const targetIdx = attributes.findIndex(
+ (a) => a.name.toLowerCase() === targetColumn.toLowerCase()
+ );
+ const featureIdxs = attributes
+ .map((_, i) => i)
+ .filter((i) => i !== targetIdx);
+
+ const featureNames = featureIdxs.map((i) => attributes[i]?.name ?? `f${i}`);
+ const data: Float64Array[] = rows.map((row) =>
+ new Float64Array(featureIdxs.map((i) => Number.parseFloat(row[i] ?? "0") || 0))
+ );
+
+ const targetAttr = targetIdx >= 0 ? attributes[targetIdx] : null;
+ const targetType = targetAttr?.type ?? "NUMERIC";
+ let target: Float64Array | Int32Array;
+
+ if (
+ targetType.toUpperCase().startsWith("NUMERIC") ||
+ targetType.toUpperCase().startsWith("REAL") ||
+ targetType.toUpperCase().startsWith("INTEGER")
+ ) {
+ target = new Float64Array(
+ rows.map((row) => Number.parseFloat(row[targetIdx] ?? "0") || 0)
+ );
+ } else {
+ // Nominal — encode as integers
+ const vals = new Set(rows.map((row) => row[targetIdx] ?? ""));
+ const valMap = new Map([...vals].map((v, i) => [v, i]));
+ target = new Int32Array(
+ rows.map((row) => valMap.get(row[targetIdx] ?? "") ?? 0)
+ );
+ }
+
+ return {
+ data,
+ target,
+ featureNames,
+ targetNames: targetAttr ? [targetAttr.name] : [],
+ description,
+ details,
+ };
+}
+
+/**
+ * List available OpenML datasets matching the given criteria.
+ */
+export async function listOpenMLDatasets(options: {
+ tag?: string;
+ limit?: number;
+ offset?: number;
+} = {}): Promise> {
+ let url = `${OPENML_BASE_URL}/data/list`;
+ const params: string[] = [];
+ if (options.tag) params.push(`tag/${encodeURIComponent(options.tag)}`);
+ if (params.length > 0) url += "/" + params.join("/");
+
+ const response = await fetch(url);
+ if (!response.ok) throw new Error(`listOpenMLDatasets: HTTP ${response.status}`);
+
+ const json = (await response.json()) as {
+ data?: {
+ dataset?: Array<{ did: number; name: string; version: number; status: string }>;
+ };
+ };
+
+ return (json.data?.dataset ?? [])
+ .slice(0, options.limit ?? 100)
+ .map((d) => ({
+ id: d.did,
+ name: d.name,
+ version: d.version,
+ status: d.status,
+ }));
+}
diff --git a/src/datasets/rcv1.ts b/src/datasets/rcv1.ts
new file mode 100644
index 0000000..f75106d
--- /dev/null
+++ b/src/datasets/rcv1.ts
@@ -0,0 +1,157 @@
+/**
+ * RCV1 dataset utilities and sparse text dataset helpers.
+ * Mirrors sklearn.datasets.rcv1 and related sparse dataset loaders.
+ */
+import type { SparseMatrix } from "../utils/sparsefuncs.js";
+
+export interface RCV1DatasetInfo {
+ nSamples: number;
+ nFeatures: number;
+ nCategories: number;
+ description: string;
+}
+
+/** Metadata about the RCV1 corpus. */
+export const RCV1_INFO: RCV1DatasetInfo = {
+ nSamples: 804414,
+ nFeatures: 47236,
+ nCategories: 103,
+ description:
+ "RCV1 — Reuters Corpus Volume 1. A collection of 804,414 news articles " +
+ "annotated with 103 topic categories. Features are TF-IDF weighted bag-of-words.",
+};
+
+export interface TextDataset {
+ data: SparseMatrix;
+ target: Int32Array;
+ targetNames: string[];
+ featureNames: string[];
+ description: string;
+}
+
+/**
+ * Build a sparse TF-IDF matrix from an array of tokenized documents.
+ * Each document is an array of term strings.
+ */
+export function buildTfIdf(
+ documents: string[][],
+ options: { maxFeatures?: number; sublinearTf?: boolean; smoothIdf?: boolean } = {}
+): { matrix: SparseMatrix; vocabulary: Map; idf: Float64Array } {
+ const { maxFeatures, sublinearTf = false, smoothIdf = true } = options;
+ const nDocs = documents.length;
+
+ // Build vocabulary
+ const df = new Map();
+ for (const doc of documents) {
+ const seen = new Set();
+ for (const term of doc) {
+ if (!seen.has(term)) { df.set(term, (df.get(term) ?? 0) + 1); seen.add(term); }
+ }
+ }
+
+ // Sort by df descending, take top maxFeatures
+ let vocab = [...df.entries()].sort((a, b) => b[1] - a[1]);
+ if (maxFeatures !== undefined) vocab = vocab.slice(0, maxFeatures);
+ const termToIdx = new Map(vocab.map(([t], i) => [t, i]));
+ const nTerms = termToIdx.size;
+
+ // IDF
+ const idf = new Float64Array(nTerms);
+ for (const [term, idx] of termToIdx) {
+ const dfi = df.get(term) ?? 0;
+ idf[idx] = Math.log(((smoothIdf ? 1 : 0) + nDocs) / ((smoothIdf ? 1 : 0) + dfi)) + 1;
+ }
+
+ // Build CSR TF-IDF matrix
+ const dataArr: number[] = [];
+ const indicesArr: number[] = [];
+ const indptrArr: number[] = [0];
+
+ for (const doc of documents) {
+ const tf = new Map();
+ for (const term of doc) {
+ const idx = termToIdx.get(term);
+ if (idx !== undefined) tf.set(idx, (tf.get(idx) ?? 0) + 1);
+ }
+ const docLen = doc.length;
+ const entries = [...tf.entries()].sort((a, b) => a[0] - b[0]);
+ for (const [idx, count] of entries) {
+ const tfVal = sublinearTf ? 1 + Math.log(count) : count / docLen;
+ const val = tfVal * (idf[idx] ?? 0);
+ if (val !== 0) { dataArr.push(val); indicesArr.push(idx); }
+ }
+ indptrArr.push(dataArr.length);
+ }
+
+ const matrix: SparseMatrix = {
+ data: new Float64Array(dataArr),
+ indices: new Int32Array(indicesArr),
+ indptr: new Int32Array(indptrArr),
+ shape: [nDocs, nTerms],
+ };
+
+ return { matrix, vocabulary: termToIdx, idf };
+}
+
+/**
+ * Generate a synthetic sparse text dataset for testing.
+ * Returns documents drawn from `nCategories` topics with `nFeatures` vocabulary.
+ */
+export function makeSparseTextDataset(options: {
+ nSamples?: number;
+ nFeatures?: number;
+ nCategories?: number;
+ avgTermsPerDoc?: number;
+ randomState?: number;
+} = {}): { X: SparseMatrix; y: Int32Array; featureNames: string[]; categoryNames: string[] } {
+ const {
+ nSamples = 200,
+ nFeatures = 500,
+ nCategories = 5,
+ avgTermsPerDoc = 20,
+ randomState = 42,
+ } = options;
+
+ let seed = randomState | 0;
+ const rng = (): number => {
+ seed = (seed ^ (seed << 13)) >>> 0;
+ seed = (seed ^ (seed >>> 17)) >>> 0;
+ seed = (seed ^ (seed << 5)) >>> 0;
+ return (seed >>> 0) / 0xffffffff;
+ };
+
+ const featureNames = Array.from({ length: nFeatures }, (_, i) => `word_${i}`);
+ const categoryNames = Array.from({ length: nCategories }, (_, i) => `category_${i}`);
+
+ const data: number[] = [];
+ const indices: number[] = [];
+ const indptr: number[] = [0];
+ const y = new Int32Array(nSamples);
+
+ for (let i = 0; i < nSamples; i++) {
+ const cat = Math.floor(rng() * nCategories);
+ y[i] = cat;
+ const nTerms = Math.max(1, Math.round(avgTermsPerDoc * (0.5 + rng())));
+ const tfMap = new Map();
+ for (let t = 0; t < nTerms; t++) {
+ // Category-biased term selection
+ const bias = rng() < 0.3 ? cat * Math.floor(nFeatures / nCategories) : 0;
+ const termIdx = (Math.floor(rng() * Math.floor(nFeatures / nCategories)) + bias) % nFeatures;
+ tfMap.set(termIdx, (tfMap.get(termIdx) ?? 0) + 1);
+ }
+ const entries = [...tfMap.entries()].sort((a, b) => a[0] - b[0]);
+ for (const [idx, count] of entries) {
+ data.push(count); indices.push(idx);
+ }
+ indptr.push(data.length);
+ }
+
+ const X: SparseMatrix = {
+ data: new Float64Array(data),
+ indices: new Int32Array(indices),
+ indptr: new Int32Array(indptr),
+ shape: [nSamples, nFeatures],
+ };
+
+ return { X, y, featureNames, categoryNames };
+}
diff --git a/src/datasets/real_datasets.ts b/src/datasets/real_datasets.ts
new file mode 100644
index 0000000..6cf4f44
--- /dev/null
+++ b/src/datasets/real_datasets.ts
@@ -0,0 +1,344 @@
+/**
+ * Real-world dataset generators and synthetic alternatives.
+ * Mirrors sklearn.datasets (california_housing, covtype, kddcup99, etc.)
+ */
+
+export interface RealDataset {
+ data: Float64Array[];
+ target: Float64Array;
+ featureNames: string[];
+ targetNames?: string[];
+ description: string;
+}
+
+export interface RealClassificationDataset extends RealDataset {
+ target: Float64Array; // integer class labels as floats
+ classes: Int32Array;
+}
+
+/**
+ * Generate a synthetic version of the California Housing dataset.
+ * The real dataset has 20,640 instances and 8 features.
+ * This generator produces a statistically similar synthetic dataset.
+ *
+ * Features: MedInc, HouseAge, AveRooms, AveBedrms, Population, AveOccup, Latitude, Longitude
+ * Target: median house value (in $100k)
+ */
+export function makeCaliforniaHousing(options: {
+ nSamples?: number;
+ noise?: number;
+ seed?: number;
+} = {}): RealDataset {
+ const { nSamples = 1000, noise = 0.1, seed = 42 } = options;
+ let rng = seed;
+ const rand = () => {
+ rng = (rng * 1664525 + 1013904223) & 0xffffffff;
+ return ((rng >>> 0) / 0xffffffff);
+ };
+ const randn = () => {
+ const u = rand() || 1e-10;
+ const v = rand() || 1e-10;
+ return Math.sqrt(-2 * Math.log(u)) * Math.cos(2 * Math.PI * v);
+ };
+
+ const featureNames = [
+ "MedInc", "HouseAge", "AveRooms", "AveBedrms",
+ "Population", "AveOccup", "Latitude", "Longitude",
+ ];
+
+ const data: Float64Array[] = [];
+ const target = new Float64Array(nSamples);
+
+ for (let i = 0; i < nSamples; i++) {
+ const medInc = Math.max(0.5, 3.0 + randn() * 2.0);
+ const houseAge = Math.max(1, Math.min(52, 28 + randn() * 12));
+ const aveRooms = Math.max(1, 5.4 + randn() * 2.0);
+ const aveBedrms = Math.max(0.5, 1.1 + randn() * 0.4);
+ const population = Math.max(10, 1400 + randn() * 1100);
+ const aveOccup = Math.max(1, 3.0 + randn() * 1.5);
+ const latitude = 35.6 + randn() * 2.1;
+ const longitude = -119.6 + randn() * 2.0;
+
+ const row = new Float64Array([
+ medInc, houseAge, aveRooms, aveBedrms,
+ population, aveOccup, latitude, longitude,
+ ]);
+ data.push(row);
+
+ // Approximate the California housing formula
+ target[i] = Math.max(0.15, Math.min(5.0,
+ 0.4524 * medInc
+ - 0.0104 * houseAge
+ + 0.0 * aveRooms
+ - 0.0 * aveBedrms
+ - 0.0 * population / 1000
+ - 0.0 * aveOccup
+ - 0.042 * latitude
+ + 0.0 * longitude
+ + 2.1 + randn() * noise,
+ ));
+ }
+
+ return {
+ data,
+ target,
+ featureNames,
+ description: "Synthetic California Housing dataset (sklearn-compatible)",
+ };
+}
+
+/**
+ * Generate a synthetic version of the Forest Covertype dataset.
+ * The real dataset has 581,012 instances and 54 features with 7 cover types.
+ *
+ * Returns integer class labels 1-7 for cover type.
+ */
+export function makeCovtype(options: {
+ nSamples?: number;
+ seed?: number;
+} = {}): RealClassificationDataset {
+ const { nSamples = 500, seed = 42 } = options;
+ let rng = seed;
+ const rand = () => {
+ rng = (rng * 1664525 + 1013904223) & 0xffffffff;
+ return ((rng >>> 0) / 0xffffffff);
+ };
+ const randn = () => {
+ const u = rand() || 1e-10;
+ const v = rand() || 1e-10;
+ return Math.sqrt(-2 * Math.log(u)) * Math.cos(2 * Math.PI * v);
+ };
+
+ // 54 features: 10 continuous, 4 binary wilderness areas, 40 binary soil types
+ const continuousFeatureNames = [
+ "Elevation", "Aspect", "Slope",
+ "Horizontal_Distance_To_Hydrology", "Vertical_Distance_To_Hydrology",
+ "Horizontal_Distance_To_Roadways", "Hillshade_9am", "Hillshade_Noon",
+ "Hillshade_3pm", "Horizontal_Distance_To_Fire_Points",
+ ];
+ const wildernessNames = [
+ "Wilderness_Area1", "Wilderness_Area2", "Wilderness_Area3", "Wilderness_Area4",
+ ];
+ const soilNames = Array.from({ length: 40 }, (_, i) => `Soil_Type${i + 1}`);
+ const featureNames = [...continuousFeatureNames, ...wildernessNames, ...soilNames];
+
+ const data: Float64Array[] = [];
+ const target = new Float64Array(nSamples);
+ const classes = new Int32Array([1, 2, 3, 4, 5, 6, 7]);
+
+ // Cover type priors (approximate): 1=36.5%, 2=48.7%, 3=6.2%, 4=0.5%, 5=1.6%, 6=2.9%, 7=3.5%
+ const priors = [0.365, 0.487, 0.062, 0.005, 0.016, 0.029, 0.035];
+ const cdf = priors.reduce((acc, p, i) => {
+ acc.push((acc[i - 1] ?? 0) + p);
+ return acc;
+ }, []);
+
+ for (let i = 0; i < nSamples; i++) {
+ // Sample class label
+ const u = rand();
+ let cls = 1;
+ for (let c = 0; c < cdf.length; c++) {
+ if (u <= (cdf[c] ?? 1)) { cls = c + 1; break; }
+ }
+ target[i] = cls;
+
+ // Continuous features (mean/std approximate per class)
+ const elevation = 2800 + cls * 50 + randn() * 200;
+ const aspect = 180 + randn() * 90;
+ const slope = 12 + randn() * 8;
+ const horizHydro = 300 + randn() * 250;
+ const vertHydro = 20 + randn() * 50;
+ const horizRoad = 2000 + randn() * 1500;
+ const hillshade9am = Math.max(0, Math.min(255, 200 + randn() * 40));
+ const hillshadeNoon = Math.max(0, Math.min(255, 220 + randn() * 30));
+ const hillshade3pm = Math.max(0, Math.min(255, 135 + randn() * 60));
+ const horizFire = 1500 + randn() * 1200;
+
+ // Binary wilderness area (one-hot)
+ const wArea = Math.floor(rand() * 4);
+ const w = new Float64Array(4);
+ w[wArea] = 1;
+
+ // Binary soil type (one-hot among 40)
+ const sType = Math.floor(rand() * 40);
+ const s = new Float64Array(40);
+ s[sType] = 1;
+
+ const row = new Float64Array([
+ elevation, aspect, slope, horizHydro, vertHydro,
+ horizRoad, hillshade9am, hillshadeNoon, hillshade3pm, horizFire,
+ ...w, ...s,
+ ]);
+ data.push(row);
+ }
+
+ return {
+ data,
+ target,
+ featureNames,
+ targetNames: ["Spruce/Fir", "Lodgepole Pine", "Ponderosa Pine",
+ "Cottonwood/Willow", "Aspen", "Douglas-fir", "Krummholz"],
+ classes,
+ description: "Synthetic Covertype dataset (sklearn-compatible, 7 classes, 54 features)",
+ };
+}
+
+/**
+ * Generate a synthetic version of the KDD Cup 1999 dataset.
+ * Returns a simplified intrusion detection dataset.
+ *
+ * @param subset - 'SA' (small) or 'SF' (larger subset), or '10percent'
+ */
+export function makeKddcup99(options: {
+ nSamples?: number;
+ subset?: "SA" | "SF" | "10percent";
+ percentAnomalies?: number;
+ seed?: number;
+} = {}): RealClassificationDataset {
+ const {
+ nSamples = 500,
+ percentAnomalies = 0.2,
+ seed = 42,
+ } = options;
+
+ let rng = seed;
+ const rand = () => {
+ rng = (rng * 1664525 + 1013904223) & 0xffffffff;
+ return ((rng >>> 0) / 0xffffffff);
+ };
+ const randn = () => {
+ const u = rand() || 1e-10;
+ const v = rand() || 1e-10;
+ return Math.sqrt(-2 * Math.log(u)) * Math.cos(2 * Math.PI * v);
+ };
+
+ const featureNames = [
+ "duration", "protocol_type", "service", "flag",
+ "src_bytes", "dst_bytes", "land", "wrong_fragment",
+ "urgent", "hot", "num_failed_logins", "logged_in",
+ "num_compromised", "root_shell", "su_attempted",
+ "num_root", "num_file_creations", "num_shells",
+ "num_access_files", "num_outbound_cmds", "is_host_login",
+ "is_guest_login", "count", "srv_count",
+ "serror_rate", "srv_serror_rate", "rerror_rate", "srv_rerror_rate",
+ "same_srv_rate", "diff_srv_rate", "srv_diff_host_rate",
+ "dst_host_count", "dst_host_srv_count",
+ "dst_host_same_srv_rate", "dst_host_diff_srv_rate",
+ "dst_host_same_src_port_rate", "dst_host_srv_diff_host_rate",
+ "dst_host_serror_rate", "dst_host_srv_serror_rate",
+ "dst_host_rerror_rate", "dst_host_srv_rerror_rate",
+ ];
+
+ const nAnomalies = Math.floor(nSamples * percentAnomalies);
+ const nNormal = nSamples - nAnomalies;
+
+ const data: Float64Array[] = [];
+ const target = new Float64Array(nSamples);
+ // 0 = normal, 1 = anomaly
+ const classes = new Int32Array([0, 1]);
+
+ for (let i = 0; i < nSamples; i++) {
+ const isAnomaly = i < nAnomalies;
+ target[i] = isAnomaly ? 1 : 0;
+
+ const row = new Float64Array(featureNames.length);
+ if (isAnomaly) {
+ // Anomaly pattern: high src_bytes, high error rates
+ row[0] = Math.max(0, randn() * 2);
+ row[4] = Math.max(0, 100000 + randn() * 50000);
+ row[5] = Math.max(0, randn() * 100);
+ row[24] = Math.max(0, Math.min(1, 0.8 + randn() * 0.2));
+ row[26] = Math.max(0, Math.min(1, 0.7 + randn() * 0.2));
+ } else {
+ // Normal: small transfers, low error
+ row[0] = Math.max(0, randn() * 5);
+ row[4] = Math.max(0, 500 + randn() * 1000);
+ row[5] = Math.max(0, 2000 + randn() * 3000);
+ row[24] = Math.max(0, Math.min(1, 0.02 + randn() * 0.05));
+ row[26] = Math.max(0, Math.min(1, 0.01 + randn() * 0.03));
+ }
+ row[22] = Math.max(0, Math.min(511, Math.abs(randn() * 50 + 10)));
+ row[31] = Math.max(0, Math.min(255, Math.abs(randn() * 50 + 100)));
+ data.push(row);
+ }
+
+ // Shuffle
+ for (let i = nSamples - 1; i > 0; i--) {
+ const j = Math.floor(rand() * (i + 1));
+ const tmp = data[i]!;
+ data[i] = data[j]!;
+ data[j] = tmp;
+ const ttmp = target[i]!;
+ target[i] = target[j]!;
+ target[j] = ttmp;
+ }
+
+ _ = nNormal; // suppress unused var
+
+ return {
+ data,
+ target,
+ featureNames,
+ targetNames: ["normal", "anomaly"],
+ classes,
+ description: "Synthetic KDD Cup 1999 network intrusion detection dataset",
+ };
+}
+
+// Suppress TS unused variable error
+let _: number;
+
+/**
+ * Load a synthetic version of the Olivetti faces dataset.
+ * 400 samples, 64x64 pixel face images (4096 features), 40 subjects.
+ */
+export function makeOlivettiFaces(options: {
+ nSamples?: number;
+ nSubjects?: number;
+ seed?: number;
+} = {}): RealDataset {
+ const { nSamples = 400, nSubjects = 40, seed = 42 } = options;
+ let rng = seed;
+ const rand = () => {
+ rng = (rng * 1664525 + 1013904223) & 0xffffffff;
+ return ((rng >>> 0) / 0xffffffff);
+ };
+ const randn = () => {
+ const u = rand() || 1e-10;
+ const v = rand() || 1e-10;
+ return Math.sqrt(-2 * Math.log(u)) * Math.cos(2 * Math.PI * v);
+ };
+
+ const nFeatures = 4096; // 64x64
+ const data: Float64Array[] = [];
+ const target = new Float64Array(nSamples);
+ const featureNames = Array.from({ length: nFeatures }, (_, i) => `pixel_${i}`);
+
+ // Each subject has a "prototype" face
+ const prototypes: Float64Array[] = Array.from({ length: nSubjects }, () => {
+ const p = new Float64Array(nFeatures);
+ for (let f = 0; f < nFeatures; f++) {
+ p[f] = Math.max(0, Math.min(1, 0.5 + randn() * 0.2));
+ }
+ return p;
+ });
+
+ for (let i = 0; i < nSamples; i++) {
+ const subject = i % nSubjects;
+ target[i] = subject;
+ const proto = prototypes[subject]!;
+ const row = new Float64Array(nFeatures);
+ for (let f = 0; f < nFeatures; f++) {
+ row[f] = Math.max(0, Math.min(1, proto[f]! + randn() * 0.05));
+ }
+ data.push(row);
+ }
+
+ return {
+ data,
+ target,
+ featureNames,
+ targetNames: Array.from({ length: nSubjects }, (_, i) => `subject_${i}`),
+ description: `Synthetic Olivetti faces dataset (${nSubjects} subjects, ${nSamples} samples)`,
+ };
+}
diff --git a/src/datasets/samples_generator.ts b/src/datasets/samples_generator.ts
new file mode 100644
index 0000000..3023de0
--- /dev/null
+++ b/src/datasets/samples_generator.ts
@@ -0,0 +1,228 @@
+/**
+ * Additional synthetic dataset generators.
+ * Mirrors sklearn.datasets: make_hastie_10_2, make_friedman1/2/3,
+ * make_sparse_uncorrelated, make_checkerboard, make_multilabel_classification.
+ */
+
+/** Result type for generated datasets. */
+export interface SamplesDatasetResult {
+ X: Float64Array[];
+ y: Float64Array | Int32Array;
+}
+
+/** Simple seeded Mulberry32 RNG for reproducibility. */
+function makeRng(seed: number): () => number {
+ let s = seed >>> 0;
+ return () => {
+ s = (s + 0x6d2b79f5) >>> 0;
+ let t = Math.imul(s ^ (s >>> 15), s | 1);
+ t ^= t + Math.imul(t ^ (t >>> 7), t | 61);
+ return ((t ^ (t >>> 14)) >>> 0) / 4294967296;
+ };
+}
+
+function randn(rng: () => number): number {
+ const u1 = Math.max(rng(), 1e-14);
+ const u2 = rng();
+ return Math.sqrt(-2 * Math.log(u1)) * Math.cos(2 * Math.PI * u2);
+}
+
+/**
+ * make_hastie_10_2 — 10-feature binary classification problem.
+ * y = sign(sum(X_i^2) - 9.34) where X ~ N(0,1).
+ */
+export function makeHastie10_2(
+ nSamples = 12000,
+ randomState = 0,
+): { X: Float64Array[]; y: Int32Array } {
+ const rng = makeRng(randomState);
+ const X: Float64Array[] = Array.from({ length: nSamples }, () => {
+ const row = new Float64Array(10);
+ for (let j = 0; j < 10; j++) row[j]! = randn(rng);
+ return row;
+ });
+ const y = Int32Array.from(X, (row) => {
+ let s = 0;
+ for (const v of row) s += v * v;
+ return s > 9.34 ? 1 : -1;
+ });
+ return { X, y };
+}
+
+/**
+ * make_friedman1 — regression dataset from Friedman (1991).
+ * y = 10*sin(π*X0*X1) + 20*(X2-0.5)^2 + 10*X3 + 5*X4 + noise
+ */
+export function makeFriedman1(
+ nSamples = 100,
+ nFeatures = 10,
+ noise = 0.0,
+ randomState = 0,
+): SamplesDatasetResult {
+ if (nFeatures < 5) throw new Error("makeFriedman1 requires at least 5 features");
+ const rng = makeRng(randomState);
+ const X: Float64Array[] = Array.from({ length: nSamples }, () => {
+ const row = new Float64Array(nFeatures);
+ for (let j = 0; j < nFeatures; j++) row[j]! = rng();
+ return row;
+ });
+ const y = Float64Array.from(X, (row) => {
+ const x0 = row[0]! ?? 0;
+ const x1 = row[1]! ?? 0;
+ const x2 = row[2]! ?? 0;
+ const x3 = row[3]! ?? 0;
+ const x4 = row[4]! ?? 0;
+ return (
+ 10 * Math.sin(Math.PI * x0 * x1) +
+ 20 * (x2 - 0.5) ** 2 +
+ 10 * x3 +
+ 5 * x4 +
+ (noise > 0 ? noise * randn(rng) : 0)
+ );
+ });
+ return { X, y };
+}
+
+/**
+ * make_friedman2 — regression with nonlinear interactions.
+ * y = sqrt(X0^2 + (X1*X2 - 1/(X1*X3))^2) + noise
+ */
+export function makeFriedman2(
+ nSamples = 100,
+ noise = 0.0,
+ randomState = 0,
+): SamplesDatasetResult {
+ const rng = makeRng(randomState);
+ const bounds: [number, number][] = [[0, 100], [40 * Math.PI, 560 * Math.PI], [0, 1], [1, 11]];
+ const X: Float64Array[] = Array.from({ length: nSamples }, () => {
+ const row = new Float64Array(4);
+ for (let j = 0; j < 4; j++) {
+ const [lo, hi] = bounds[j]!;
+ row[j]! = lo + rng() * (hi - lo);
+ }
+ return row;
+ });
+ const y = Float64Array.from(X, (row) => {
+ const x0 = row[0]! ?? 0;
+ const x1 = row[1]! ?? 0;
+ const x2 = row[2]! ?? 0;
+ const x3 = Math.max(row[3]! ?? 1, 1e-6);
+ const inner = x1 * x2 - 1 / (x1 * x3);
+ return Math.sqrt(x0 ** 2 + inner ** 2) + (noise > 0 ? noise * randn(rng) : 0);
+ });
+ return { X, y };
+}
+
+/**
+ * make_friedman3 — regression with arctan transformation.
+ * y = arctan((X1*X2 - 1/(X1*X3)) / X0) + noise
+ */
+export function makeFriedman3(
+ nSamples = 100,
+ noise = 0.0,
+ randomState = 0,
+): SamplesDatasetResult {
+ const rng = makeRng(randomState);
+ const bounds: [number, number][] = [[0, 100], [40 * Math.PI, 560 * Math.PI], [0, 1], [1, 11]];
+ const X: Float64Array[] = Array.from({ length: nSamples }, () => {
+ const row = new Float64Array(4);
+ for (let j = 0; j < 4; j++) {
+ const [lo, hi] = bounds[j]!;
+ row[j]! = lo + rng() * (hi - lo);
+ }
+ return row;
+ });
+ const y = Float64Array.from(X, (row) => {
+ const x0 = Math.max(Math.abs(row[0]! ?? 0), 1e-6);
+ const x1 = row[1]! ?? 0;
+ const x2 = row[2]! ?? 0;
+ const x3 = Math.max(row[3]! ?? 1, 1e-6);
+ const inner = x1 * x2 - 1 / (x1 * x3);
+ return Math.atan(inner / x0) + (noise > 0 ? noise * randn(rng) : 0);
+ });
+ return { X, y };
+}
+
+/**
+ * make_sparse_uncorrelated — regression dataset with 4 informative features
+ * and `nFeatures - 4` noise features.
+ */
+export function makeSparseUncorrelated(
+ nSamples = 100,
+ nFeatures = 10,
+ randomState = 0,
+): SamplesDatasetResult {
+ const rng = makeRng(randomState);
+ const X: Float64Array[] = Array.from({ length: nSamples }, () =>
+ Float64Array.from({ length: nFeatures }, () => randn(rng)),
+ );
+ const coef = [1, 2, 0.5, -0.5]; // informative coefficients
+ const y = Float64Array.from(X, (row) => {
+ let s = 0;
+ for (let j = 0; j < coef.length; j++) s += (coef[j]! ?? 0) * (row[j]! ?? 0);
+ s += randn(rng);
+ return s;
+ });
+ return { X, y };
+}
+
+/**
+ * make_multilabel_classification — random multilabel dataset.
+ *
+ * @param nSamples - Number of samples.
+ * @param nFeatures - Number of features.
+ * @param nClasses - Number of classes (labels).
+ * @param nLabels - Average number of labels per sample.
+ * @param randomState - Random seed.
+ */
+export function makeMultilabelClassification(
+ nSamples = 100,
+ nFeatures = 20,
+ nClasses = 5,
+ nLabels = 2,
+ randomState = 0,
+): { X: Float64Array[]; y: Int32Array[] } {
+ const rng = makeRng(randomState);
+ const X: Float64Array[] = Array.from({ length: nSamples }, () =>
+ Float64Array.from({ length: nFeatures }, () => rng() > 0.5 ? 1 : 0),
+ );
+ const y: Int32Array[] = Array.from({ length: nSamples }, () => {
+ const row = new Int32Array(nClasses);
+ const nActive = Math.max(1, Math.round(nLabels + (rng() - 0.5) * 2));
+ for (let k = 0; k < nActive && k < nClasses; k++) {
+ row[Math.floor(rng() * nClasses)]! = 1;
+ }
+ return row;
+ });
+ return { X, y };
+}
+
+/**
+ * make_checkerboard — checkerboard pattern for biclustering.
+ *
+ * @param shape - [n_rows, n_cols].
+ * @param nClusters - [n_row_clusters, n_col_clusters].
+ * @param noise - Noise standard deviation.
+ * @param randomState - Random seed.
+ */
+export function makeCheckerboard(
+ shape: [number, number] = [300, 300],
+ nClusters: [number, number] = [4, 3],
+ noise = 0.5,
+ randomState = 0,
+): { data: Float64Array[]; rowLabels: Int32Array; colLabels: Int32Array } {
+ const rng = makeRng(randomState);
+ const [nRows, nCols] = shape;
+ const [nRowC, nColC] = nClusters;
+ const rowLabels = Int32Array.from({ length: nRows }, (_, i) => i % nRowC);
+ const colLabels = Int32Array.from({ length: nCols }, (_, j) => j % nColC);
+ const data: Float64Array[] = Array.from({ length: nRows }, (_, i) => {
+ const row = new Float64Array(nCols);
+ for (let j = 0; j < nCols; j++) {
+ const match = (rowLabels[i]! % 2) === (colLabels[j]! % 2);
+ row[j]! = (match ? 1 : 0) + noise * randn(rng);
+ }
+ return row;
+ });
+ return { data, rowLabels, colLabels };
+}
diff --git a/src/datasets/svmlight.ts b/src/datasets/svmlight.ts
new file mode 100644
index 0000000..3fc6d3d
--- /dev/null
+++ b/src/datasets/svmlight.ts
@@ -0,0 +1,113 @@
+/**
+ * SVMLight format loading and saving utilities.
+ * Ports: load_svmlight_file, dump_svmlight_file
+ */
+
+export interface SVMLightDataset {
+ data: Float64Array[];
+ target: Float64Array;
+ nFeatures: number;
+}
+
+/**
+ * Parse SVMLight / LibSVM format text.
+ * Format: