From 3d3dbc3dea4601c4892989befa7e5037be26dc77 Mon Sep 17 00:00:00 2001 From: godofecht Date: Tue, 29 Sep 2026 12:21:16 +0100 Subject: [PATCH 1/2] Run the algorithm in the four rows that were standing in for one Four implementations said in their own comments that they were simplified, so the registry timed them and showed the numbers without ranking them, and the page carried a paragraph explaining why three of the four would otherwise have been its widest wins. They were wide because they were not doing the work. spectral_coclustering bucketed each row sum between the smallest and the largest. It runs Dhillon's algorithm now: normalize by the square roots of the row and column sums, take the singular vectors after the trivial one, and cluster rows and columns together in the space they span. 0.03 ms against scikit-learn's 8.09. spectral_biclustering bucketed row means and column means the same way. It runs Kluger's: Sinkhorn-Knopp, then rank the singular vectors by how well a one dimensional k-means fits each one, keep the best few, and cluster rows and columns separately. 0.09 ms against 35.17. elliptic_envelope drew random subsets and scored each by the product of the diagonal of its covariance, which is the determinant only when the features do not correlate, and stopped at the draw. The concentration step is the whole content of FastMCD, so it takes the h points closest in Mahalanobis distance, recomputes, and repeats until the determinant stops falling. The determinant comes from a Cholesky factor, and predict uses the full distance rather than the diagonal. 2.03 ms against 19.46. nca computed neither the objective nor its gradient. Its distance summed the squares of the individual products instead of squaring the projected difference, which is not a distance in the projected space, and its gradient kept only the term that pulls same-class points together, with nothing pushing back. Both are written out now. 181 ms against 856. Two test files cover them, by behaviour rather than by shape: the block structure has to come back out of both biclusterings, every planted outlier has to be found with no false alarms, and the projection has to separate the classes better than the identity did. Each was checked by inverting an assertion and confirming a non-zero exit. --- benchmarks/estimator_coverage.json | 15 +- lib/scikit/covariance.flow | 299 ++++++--- lib/scikit/extra_estimators.flow | 704 ++++++++++++++++---- tests/test_outlier_and_metric_learning.flow | 143 ++++ tests/test_spectral_biclustering.flow | 139 ++++ 5 files changed, 1096 insertions(+), 204 deletions(-) create mode 100644 tests/test_outlier_and_metric_learning.flow create mode 100644 tests/test_spectral_biclustering.flow diff --git a/benchmarks/estimator_coverage.json b/benchmarks/estimator_coverage.json index 17e8d94..bd666a8 100644 --- a/benchmarks/estimator_coverage.json +++ b/benchmarks/estimator_coverage.json @@ -2,9 +2,8 @@ "schema_version": 1, "counts": { "estimators": 203, - "runnable": 166, + "runnable": 170, "shaped": 22, - "simplified": 4, "flow_only": 6, "different_shape": 5 }, @@ -2079,8 +2078,7 @@ "contamination", "n_trials" ], - "bucket": "simplified", - "reason": "the implementation's own comments call it simplified, so timing it against scikit-learn's algorithm compares two different things", + "bucket": "runnable", "sklearn_estimator": "EllipticEnvelope" }, { @@ -7395,8 +7393,7 @@ "epochs", "lr" ], - "bucket": "simplified", - "reason": "the implementation's own comments call it simplified, so timing it against scikit-learn's algorithm compares two different things", + "bucket": "runnable", "sklearn_estimator": "NeighborhoodComponentsAnalysis" }, { @@ -11470,8 +11467,7 @@ "n_row_clusters", "n_col_clusters" ], - "bucket": "simplified", - "reason": "the implementation's own comments call it simplified, so timing it against scikit-learn's algorithm compares two different things", + "bucket": "runnable", "sklearn_estimator": "SpectralBiclustering" }, { @@ -11553,8 +11549,7 @@ "X", "n_clusters" ], - "bucket": "simplified", - "reason": "the implementation's own comments call it simplified, so timing it against scikit-learn's algorithm compares two different things", + "bucket": "runnable", "sklearn_estimator": "SpectralCoclustering" }, { diff --git a/lib/scikit/covariance.flow b/lib/scikit/covariance.flow index f863222..b087014 100644 --- a/lib/scikit/covariance.flow +++ b/lib/scikit/covariance.flow @@ -1216,10 +1216,134 @@ export struct EllipticEnvelope { fitted: bool } +# Minimum covariance determinant, by the FastMCD concentration step of +# Rousseeuw and Van Driessen (1999), which is what scikit-learn's MinCovDet +# and EllipticEnvelope run. +# +# What stood here drew random subsets and scored each one by the product of +# the diagonal of its covariance, which is the determinant only when the +# features are uncorrelated. It also stopped at the random draw, where the +# algorithm's whole content is the concentration step: take the h points +# closest in Mahalanobis distance to the current estimate, recompute the +# estimate from those, and repeat until the determinant stops falling. +# +# The helper names carry a prefix because a non-exported function with the +# same name in two modules collides in the generated C. Flow issue #465. + +# Cholesky of a symmetric positive definite p by p matrix, in place into L. +# A ridge is added on failure, which is what a singular subset needs. +function _mcd_cholesky(cov: ptr, L: ptr, p: i32) -> bool { + for i in 0 to p * p { L[i] = 0.0 } + for i in 0 to p { + let mut j: i32 = 0 + while j <= i { + let mut s: f32 = cov[i * p + j] + let mut k: i32 = 0 + while k < j { + s = s - L[i * p + k] * L[j * p + k] + k = k + 1 + } + if i == j { + if s <= 0.0000000001 { return false } + L[i * p + i] = sqrt((s) as f64) as f32 + } else { + L[i * p + j] = s / L[j * p + j] + } + j = j + 1 + } + } + return true +} + +# Twice the sum of the logs of the diagonal, which is the log determinant of +# the matrix the factor came from. Comparing logs keeps a product of p small +# variances from underflowing. +function _mcd_logdet(L: ptr, p: i32) -> f32 { + let mut acc: f32 = 0.0 + for i in 0 to p { + acc = acc + (log((L[i * p + i]) as f64) as f32) + } + return 2.0 * acc +} + +# Squared Mahalanobis distance of every row, from the Cholesky factor: solve +# L z = x - loc by forward substitution and take z^T z. +function _mcd_distances(xdata: ptr, n: i32, p: i32, loc: ptr, L: ptr, out: ptr, z: ptr) -> void { + for i in 0 to n { + let row: ptr = xdata + i * p + let mut a: i32 = 0 + while a < p { + let mut s: f32 = row[a] - loc[a] + let mut b: i32 = 0 + while b < a { + s = s - L[a * p + b] * z[b] + b = b + 1 + } + z[a] = s / L[a * p + a] + a = a + 1 + } + let mut d: f32 = 0.0 + let mut c: i32 = 0 + while c < p { + d = d + z[c] * z[c] + c = c + 1 + } + out[i] = d + } +} + +# Mean and covariance of the rows the index list names. +function _mcd_estimate(xdata: ptr, p: i32, idx: ptr, h: i32, loc: ptr, cov: ptr) -> void { + for j in 0 to p { loc[j] = 0.0 } + for t in 0 to h { + let row: ptr = xdata + idx[t] * p + let mut j2: i32 = 0 + while j2 < p { + loc[j2] = loc[j2] + row[j2] + j2 = j2 + 1 + } + } + let inv_h: f32 = 1.0 / (h as f32) + for j3 in 0 to p { loc[j3] = loc[j3] * inv_h } + + for i in 0 to p * p { cov[i] = 0.0 } + for t2 in 0 to h { + let row2: ptr = xdata + idx[t2] * p + for i2 in 0 to p { + let di: f32 = row2[i2] - loc[i2] + let crow: ptr = cov + i2 * p + let mut j4: i32 = 0 + while j4 < p { + crow[j4] = crow[j4] + di * (row2[j4] - loc[j4]) + j4 = j4 + 1 + } + } + } + for i3 in 0 to p * p { cov[i3] = cov[i3] * inv_h } +} + +# The h smallest distances, by partial selection into idx. +function _mcd_take_closest(dists: ptr, n: i32, h: i32, order: ptr, idx: ptr) -> void { + for i in 0 to n { order[i] = i } + for a in 0 to h { + let mut best: i32 = a + let mut b: i32 = a + 1 + while b < n { + if dists[order[b]] < dists[order[best]] { best = b } + b = b + 1 + } + let t: i32 = order[a] + order[a] = order[best] + order[best] = t + idx[a] = order[a] + } +} + export function elliptic_envelope_fit(X: Matrix, contamination: f32, n_trials: i32) -> EllipticEnvelope { let n: i32 = X.rows let p: i32 = X.cols - let h: i32 = (n * 3) / 4 + let xdata: ptr = X.data + let mut h: i32 = (n * 3) / 4 if h > n { h = n } if h < 2 { h = 2 } @@ -1228,69 +1352,67 @@ export function elliptic_envelope_fit(X: Matrix, contamination: f32, n_trials: i let best_location: ptr = array_new_f32(p) let best_cov: ptr > = malloc((p as i64) * 8) as ptr > for i in 0 to p { best_cov[i] = array_new_f32(p) } - let mut best_det: f32 = 1000000000.0 - - let support_: ptr = array_new_f32(h) as ptr + let support_: ptr = malloc((h as i64) * 4 + 64) as ptr + let mut best_logdet: f32 = 1000000000.0 + let mut have_best: bool = false + + let loc: ptr = array_new_f32(p) + let cov: ptr = array_new_f32(p * p) + let L: ptr = array_new_f32(p * p) + let z: ptr = array_new_f32(p) + let dists: ptr = array_new_f32(n) + let order: ptr = malloc((n as i64) * 4 + 64) as ptr + let idx: ptr = malloc((h as i64) * 4 + 64) as ptr for trial in 0 to n_trials { - # Random subset of size h - let subset: ptr = array_new_f32(h) as ptr - for i in 0 to h { - subset[i] = rand() % n - } - - # Compute mean of subset - let loc: ptr = array_new_f32(p) - for j in 0 to p { - let mut s: f32 = 0.0 - for i in 0 to h { - s = s + matrix_at(X, subset[i], j) + for i2 in 0 to h { idx[i2] = rand() % n } + + let mut trial_logdet: f32 = 1000000000.0 + let mut ok: bool = false + let mut step: i32 = 0 + while step < 30 { + _mcd_estimate(xdata, p, idx, h, loc, cov) + if not _mcd_cholesky(cov, L, p) { + # A degenerate subset gets a ridge rather than being dropped, + # which keeps a trial that started badly from ending the fit. + for r in 0 to p { cov[r * p + r] = cov[r * p + r] + 0.000001 } + if not _mcd_cholesky(cov, L, p) { break } } - loc[j] = s / (h as f32) - } - - # Compute covariance of subset - let cov: ptr > = malloc((p as i64) * 8) as ptr > - for i in 0 to p { cov[i] = array_new_f32(p) } - for i in 0 to p { - for j in 0 to p { - let mut s: f32 = 0.0 - for k in 0 to h { - let di: f32 = matrix_at(X, subset[k], i) - loc[i] - let dj: f32 = matrix_at(X, subset[k], j) - loc[j] - s = s + di * dj - } - cov[i][j] = s / (h as f32) + let ld: f32 = _mcd_logdet(L, p) + ok = true + # The concentration step cannot raise the determinant, so a step + # that fails to lower it means this start has converged. + if ld > trial_logdet - 0.0000001 { + trial_logdet = ld + break } - } - - # Compute determinant (simplified: product of diagonal) - let mut det: f32 = 1.0 - for i in 0 to p { - if cov[i][i] > 0.0000000001 { - det = det * cov[i][i] - } else { - det = det * 0.0000000001 - } - } - - if det < best_det { - best_det = det - for j in 0 to p { best_location[j] = loc[j] } - for i in 0 to p { - for j in 0 to p { - best_cov[i][j] = cov[i][j] + trial_logdet = ld + _mcd_distances(xdata, n, p, loc, L, dists, z) + _mcd_take_closest(dists, n, h, order, idx) + step = step + 1 + } + + if ok { + if trial_logdet < best_logdet || not have_best { + best_logdet = trial_logdet + have_best = true + for j in 0 to p { best_location[j] = loc[j] } + for i3 in 0 to p { + for j2 in 0 to p { best_cov[i3][j2] = cov[i3 * p + j2] } } + for t in 0 to h { support_[t] = idx[t] } } - for i in 0 to h { support_[i] = subset[i] } } - - array_free_f32(loc) - for i in 0 to p { array_free_f32(cov[i]) } - free(cov as ptr) - array_free_f32(subset as ptr) } + array_free_f32(loc) + array_free_f32(cov) + array_free_f32(L) + array_free_f32(z) + array_free_f32(dists) + free(order as ptr) + free(idx as ptr) + return EllipticEnvelope { covariance: best_cov, location: best_location, @@ -1307,34 +1429,55 @@ export function elliptic_envelope_predict(model: EllipticEnvelope, X: Matrix) -> let p: i32 = model.n_features let result: ptr = array_new_f32(n) - # Compute Mahalanobis distance for each sample - # Simplified: use diagonal covariance inverse - let inv_diag: ptr = array_new_f32(p) + # The full Mahalanobis distance, from a Cholesky factor of the fitted + # covariance. It used to invert the diagonal alone, which measures a + # different distance whenever two features are correlated. + let cov: ptr = array_new_f32(p * p) for i in 0 to p { - if model.covariance[i][i] > 0.0000000001 { - inv_diag[i] = 1.0 / model.covariance[i][i] - } else { - inv_diag[i] = 0.0 - } + for j in 0 to p { cov[i * p + j] = model.covariance[i][j] } } - - # Threshold based on contamination - let threshold: f32 = (p as f32) * 2.0 - - for i in 0 to n { - let mut dist: f32 = 0.0 - for j in 0 to p { - let diff: f32 = matrix_at(X, i, j) - model.location[j] - dist = dist + diff * diff * inv_diag[j] - } - if dist > threshold { - result[i] = -1.0 - } else { - result[i] = 1.0 + let L: ptr = array_new_f32(p * p) + let z: ptr = array_new_f32(p) + if not _mcd_cholesky(cov, L, p) { + for r in 0 to p { cov[r * p + r] = cov[r * p + r] + 0.000001 } + if not _mcd_cholesky(cov, L, p) { + for i2 in 0 to n { result[i2] = 1.0 } + array_free_f32(cov) + array_free_f32(L) + array_free_f32(z) + return result } } - array_free_f32(inv_diag) + let dists: ptr = array_new_f32(n) + _mcd_distances(X.data, n, p, model.location, L, dists, z) + + # The cut is the contaminated fraction of the sample, which is how + # scikit-learn sets its own: the points beyond it are the outliers. + let mut n_out: i32 = ((n as f32) * model.contamination) as i32 + if n_out < 0 { n_out = 0 } + if n_out > n { n_out = n } + let order: ptr = malloc((n as i64) * 4 + 64) as ptr + for i3 in 0 to n { order[i3] = i3 } + for a in 0 to n_out { + let mut best: i32 = a + let mut b: i32 = a + 1 + while b < n { + if dists[order[b]] > dists[order[best]] { best = b } + b = b + 1 + } + let t: i32 = order[a] + order[a] = order[best] + order[best] = t + } + for i4 in 0 to n { result[i4] = 1.0 } + for a2 in 0 to n_out { result[order[a2]] = -1.0 } + + free(order as ptr) + array_free_f32(dists) + array_free_f32(cov) + array_free_f32(L) + array_free_f32(z) return result } diff --git a/lib/scikit/extra_estimators.flow b/lib/scikit/extra_estimators.flow index 4c13a99..1b52eab 100644 --- a/lib/scikit/extra_estimators.flow +++ b/lib/scikit/extra_estimators.flow @@ -198,62 +198,373 @@ export struct SpectralBiclustering { fitted: bool } -export function spectral_biclustering_fit(X: Matrix, n_row_clusters: i32, n_col_clusters: i32) -> SpectralBiclustering { - let n: i32 = X.rows - let p: i32 = X.cols +# Modified Gram-Schmidt over the columns of a row-major m by k block. +function _bic_orthonormalize(V: ptr, m: i32, k: i32) -> void { + for c in 0 to k { + for prev in 0 to c { + let mut dot: f32 = 0.0 + let mut i: i32 = 0 + while i < m { + dot = dot + V[i * k + c] * V[i * k + prev] + i = i + 1 + } + let mut i2: i32 = 0 + while i2 < m { + V[i2 * k + c] = V[i2 * k + c] - dot * V[i2 * k + prev] + i2 = i2 + 1 + } + } + let mut norm: f32 = 0.0 + let mut i3: i32 = 0 + while i3 < m { + let v: f32 = V[i3 * k + c] + norm = norm + v * v + i3 = i3 + 1 + } + norm = sqrt((norm) as f64) as f32 + if norm > 0.0000000001 { + let inv: f32 = 1.0 / norm + let mut i4: i32 = 0 + while i4 < m { + V[i4 * k + c] = V[i4 * k + c] * inv + i4 = i4 + 1 + } + } + } +} - # Simplified: assign rows and columns to clusters by k-means on row/col means - let row_means: ptr = array_new_f32(n) - let col_means: ptr = array_new_f32(p) - for i in 0 to n { - let mut m: f32 = 0.0 - for j in 0 to p { m = m + matrix_at(X, i, j) } - row_means[i] = m / (p as f32) +# The k leading right singular vectors of A, as the k leading eigenvectors of +# A^T A by block power iteration. A is row-major n by p, V is p by k. +function _bic_right_singular(A: ptr, n: i32, p: i32, k: i32, V: ptr) -> void { + let M: ptr = array_new_f32(p * p) + cblas_sgemm(101, 112, 111, p, p, n, 1.0, A, p, A, p, 0.0, M, p) + + let Z: ptr = array_new_f32(p * k) + srand(42) + for i in 0 to p * k { V[i] = ((rand() % 1000) as f32) / 500.0 - 1.0 } + _bic_orthonormalize(V, p, k) + + let mut iter: i32 = 0 + while iter < 300 { + blas_matmat(M, V, Z, p, k, p, 1.0, 0.0) + _bic_orthonormalize(Z, p, k) + let mut max_change: f32 = 0.0 + let mut i2: i32 = 0 + let total: i32 = p * k + while i2 < total { + let mut diff: f32 = Z[i2] - V[i2] + if diff < 0.0 { diff = 0.0 - diff } + if diff > max_change { max_change = diff } + V[i2] = Z[i2] + i2 = i2 + 1 + } + if max_change < 0.0000001 { break } + iter = iter + 1 } - for j in 0 to p { - let mut m: f32 = 0.0 - for i in 0 to n { m = m + matrix_at(X, i, j) } - col_means[j] = m / (n as f32) + + array_free_f32(M) + array_free_f32(Z) +} + +# Lloyd's algorithm over m points in d dimensions, seeded by points spread +# evenly through the set so a fit is reproducible. +function _bic_kmeans(Z: ptr, m: i32, d: i32, k: i32, labels: ptr) -> void { + let centers: ptr = array_new_f32(k * d) + let counts: ptr = malloc((k as i64) * 4 + 64) as ptr + for c in 0 to k { + let src: i32 = (c * m) / k + let mut j: i32 = 0 + while j < d { + centers[c * d + j] = Z[src * d + j] + j = j + 1 + } } + for i in 0 to m { labels[i] = 0 } + + let mut iter: i32 = 0 + while iter < 100 { + let mut moved: i32 = 0 + let mut i: i32 = 0 + while i < m { + let zi: ptr = Z + i * d + let mut best: i32 = 0 + let mut best_d: f32 = 999999999999.0 + let mut c: i32 = 0 + while c < k { + let cc: ptr = centers + c * d + let mut acc: f32 = 0.0 + let mut j2: i32 = 0 + while j2 < d { + let diff: f32 = zi[j2] - cc[j2] + acc = acc + diff * diff + j2 = j2 + 1 + } + if acc < best_d { + best_d = acc + best = c + } + c = c + 1 + } + if labels[i] != best { + labels[i] = best + moved = moved + 1 + } + i = i + 1 + } - let row_labels: ptr = array_new_f32(n) as ptr - let col_labels: ptr = array_new_f32(p) as ptr + for c2 in 0 to k { + counts[c2] = 0 + let mut j3: i32 = 0 + while j3 < d { + centers[c2 * d + j3] = 0.0 + j3 = j3 + 1 + } + } + let mut i2: i32 = 0 + while i2 < m { + let lab: i32 = labels[i2] + let zi2: ptr = Z + i2 * d + let cc2: ptr = centers + lab * d + let mut j4: i32 = 0 + while j4 < d { + cc2[j4] = cc2[j4] + zi2[j4] + j4 = j4 + 1 + } + counts[lab] = counts[lab] + 1 + i2 = i2 + 1 + } + for c3 in 0 to k { + if counts[c3] > 0 { + let inv: f32 = 1.0 / (counts[c3] as f32) + let mut j5: i32 = 0 + while j5 < d { + centers[c3 * d + j5] = centers[c3 * d + j5] * inv + j5 = j5 + 1 + } + } + } - # Simple threshold-based clustering - let mut row_min: f32 = row_means[0] - let mut row_max: f32 = row_means[0] - for i in 0 to n { - if row_means[i] < row_min { row_min = row_means[i] } - if row_means[i] > row_max { row_max = row_means[i] } + if moved == 0 { break } + iter = iter + 1 } - let row_range: f32 = row_max - row_min - for i in 0 to n { - if row_range > 0.0000000001 { - row_labels[i] = ((row_means[i] - row_min) / row_range * (n_row_clusters as f32)) as i32 - if row_labels[i] >= n_row_clusters { row_labels[i] = n_row_clusters - 1 } - } else { - row_labels[i] = 0 + + array_free_f32(centers) + free(counts as ptr) +} + +# Spectral biclustering, after Kluger et al. (2003), which is what +# scikit-learn's SpectralBiclustering implements. +# +# The checkerboard structure it looks for shows up in the singular vectors of +# the bistochastically normalized matrix: a vector that fits the structure is +# close to piecewise constant. So the leading vectors are ranked by how well a +# one dimensional k-means fits each of them, the best few are kept, and rows +# and columns are clustered separately in the space they span. +# +# What stood here bucketed row means and column means between their smallest +# and largest values, which finds no checkerboard and came out at thirty +# thousand times scikit-learn for doing almost nothing. + +# Sinkhorn-Knopp, which is the default normalization scikit-learn uses. +function _bic_bistochastic(A: ptr, n: i32, p: i32) -> void { + let r: ptr = array_new_f32(n) + let c: ptr = array_new_f32(p) + let mut iter: i32 = 0 + while iter < 100 { + for i in 0 to n { r[i] = 0.0 } + for j in 0 to p { c[j] = 0.0 } + let mut i2: i32 = 0 + while i2 < n { + let row: ptr = A + i2 * p + let mut s: f32 = 0.0 + let mut j2: i32 = 0 + while j2 < p { + s = s + row[j2] + j2 = j2 + 1 + } + if s < 0.0000000001 { s = 0.0000000001 } + r[i2] = 1.0 / (sqrt((s) as f64) as f32) + i2 = i2 + 1 + } + let mut i3: i32 = 0 + while i3 < n { + let row2: ptr = A + i3 * p + let ri: f32 = r[i3] + let mut j3: i32 = 0 + while j3 < p { + let v: f32 = row2[j3] * ri + row2[j3] = v + c[j3] = c[j3] + v + j3 = j3 + 1 + } + i3 = i3 + 1 + } + let mut max_dev: f32 = 0.0 + let mut j4: i32 = 0 + while j4 < p { + let mut s2: f32 = c[j4] + if s2 < 0.0000000001 { s2 = 0.0000000001 } + let scale: f32 = 1.0 / (sqrt((s2) as f64) as f32) + c[j4] = scale + let dev: f32 = fabs((scale - 1.0) as f64) as f32 + if dev > max_dev { max_dev = dev } + j4 = j4 + 1 + } + let mut i4: i32 = 0 + while i4 < n { + let row3: ptr = A + i4 * p + let mut j5: i32 = 0 + while j5 < p { + row3[j5] = row3[j5] * c[j5] + j5 = j5 + 1 + } + i4 = i4 + 1 } + if max_dev < 0.00001 { break } + iter = iter + 1 } + array_free_f32(r) + array_free_f32(c) +} - let mut col_min: f32 = col_means[0] - let mut col_max: f32 = col_means[0] - for j in 0 to p { - if col_means[j] < col_min { col_min = col_means[j] } - if col_means[j] > col_max { col_max = col_means[j] } +# How well a one dimensional k-means fits a vector, as the sum of squared +# distances to the assigned centroid. A vector that carries the checkerboard +# is close to piecewise constant and scores low. +function _bic_piecewise_error(v: ptr, m: i32, k: i32) -> f32 { + let labels: ptr = malloc((m as i64) * 4 + 64) as ptr + _bic_kmeans(v, m, 1, k, labels) + + let sums: ptr = array_new_f32(k) + let counts: ptr = malloc((k as i64) * 4 + 64) as ptr + for c in 0 to k { counts[c] = 0 } + let mut i: i32 = 0 + while i < m { + sums[labels[i]] = sums[labels[i]] + v[i] + counts[labels[i]] = counts[labels[i]] + 1 + i = i + 1 } - let col_range: f32 = col_max - col_min - for j in 0 to p { - if col_range > 0.0000000001 { - col_labels[j] = ((col_means[j] - col_min) / col_range * (n_col_clusters as f32)) as i32 - if col_labels[j] >= n_col_clusters { col_labels[j] = n_col_clusters - 1 } - } else { - col_labels[j] = 0 + for c2 in 0 to k { + if counts[c2] > 0 { sums[c2] = sums[c2] / (counts[c2] as f32) } + } + let mut err: f32 = 0.0 + let mut i2: i32 = 0 + while i2 < m { + let d: f32 = v[i2] - sums[labels[i2]] + err = err + d * d + i2 = i2 + 1 + } + + free(labels as ptr) + free(counts as ptr) + array_free_f32(sums) + return err +} + +export function spectral_biclustering_fit(X: Matrix, n_row_clusters: i32, n_col_clusters: i32) -> SpectralBiclustering { + let n: i32 = X.rows + let p: i32 = X.cols + let xdata: ptr = X.data + + let An: ptr = array_new_f32(n * p) + for i in 0 to n * p { An[i] = xdata[i] } + _bic_bistochastic(An, n, p) + + # scikit-learn takes six singular vectors and keeps the best three of + # them, which is as many as the matrix has here when it is narrow. + let mut n_vecs: i32 = 6 + if n_vecs > p { n_vecs = p } + if n_vecs < 2 { n_vecs = 2 } + + let V: ptr = array_new_f32(p * n_vecs) + _bic_right_singular(An, n, p, n_vecs, V) + + let U: ptr = array_new_f32(n * n_vecs) + blas_matmat(An, V, U, n, n_vecs, p, 1.0, 0.0) + for c in 0 to n_vecs { + let mut norm: f32 = 0.0 + let mut i2: i32 = 0 + while i2 < n { + let v: f32 = U[i2 * n_vecs + c] + norm = norm + v * v + i2 = i2 + 1 + } + norm = sqrt((norm) as f64) as f32 + if norm > 0.0000000001 { + let inv: f32 = 1.0 / norm + let mut i3: i32 = 0 + while i3 < n { + U[i3 * n_vecs + c] = U[i3 * n_vecs + c] * inv + i3 = i3 + 1 + } } } - array_free_f32(row_means) - array_free_f32(col_means) + # Rank the vectors after the trivial one by how close to piecewise + # constant they are, and keep the best few. + let mut n_best: i32 = 3 + if n_best > n_vecs - 1 { n_best = n_vecs - 1 } + if n_best < 1 { n_best = 1 } + + let scratch: ptr = array_new_f32(n) + let errs: ptr = array_new_f32(n_vecs) + let chosen: ptr = malloc((n_vecs as i64) * 4 + 64) as ptr + for c2 in 0 to n_vecs { chosen[c2] = 0 } + for c3 in 1 to n_vecs { + let mut i4: i32 = 0 + while i4 < n { + scratch[i4] = U[i4 * n_vecs + c3] + i4 = i4 + 1 + } + errs[c3] = _bic_piecewise_error(scratch, n, n_row_clusters) + } + for b in 0 to n_best { + let mut best_c: i32 = -1 + let mut best_err: f32 = 999999999999.0 + for c4 in 1 to n_vecs { + if chosen[c4] == 0 { + if errs[c4] < best_err { + best_err = errs[c4] + best_c = c4 + } + } + } + if best_c >= 0 { chosen[best_c] = 1 } + } + + # Rows and columns are clustered separately, in the space the kept + # vectors span. + let Zr: ptr = array_new_f32(n * n_best) + let Zc: ptr = array_new_f32(p * n_best) + let mut slot: i32 = 0 + for c5 in 1 to n_vecs { + if chosen[c5] == 1 { + let mut i5: i32 = 0 + while i5 < n { + Zr[i5 * n_best + slot] = U[i5 * n_vecs + c5] + i5 = i5 + 1 + } + let mut j2: i32 = 0 + while j2 < p { + Zc[j2 * n_best + slot] = V[j2 * n_vecs + c5] + j2 = j2 + 1 + } + slot = slot + 1 + } + } + + let row_labels: ptr = malloc((n as i64) * 4 + 64) as ptr + let col_labels: ptr = malloc((p as i64) * 4 + 64) as ptr + _bic_kmeans(Zr, n, n_best, n_row_clusters, row_labels) + _bic_kmeans(Zc, p, n_best, n_col_clusters, col_labels) + + array_free_f32(An) + array_free_f32(U) + array_free_f32(V) + array_free_f32(Zr) + array_free_f32(Zc) + array_free_f32(scratch) + array_free_f32(errs) + free(chosen as ptr) return SpectralBiclustering { row_labels: row_labels, @@ -284,49 +595,137 @@ export struct SpectralCoclustering { fitted: bool } +# Spectral co-clustering, after Dhillon (2001), which is what scikit-learn's +# SpectralCoclustering implements. +# +# The data matrix is treated as the adjacency of a bipartite graph. Normalizing +# it by the square roots of its row and column sums turns the partitioning +# problem into a singular value problem: the leading pair is trivial, and the +# next ceil(log2(k)) pairs embed rows and columns in a common space where a +# k-means run recovers the co-clusters. +# +# What stood here scaled each row sum into a bucket between the smallest and +# the largest, which is not that algorithm and is not any other one. It came +# out at three thousand times scikit-learn because it did almost nothing. +# +# The helper names carry a prefix because a non-exported function with the same +# name in two modules collides in the generated C. Flow issue #465. + export function spectral_coclustering_fit(X: Matrix, n_clusters: i32) -> SpectralCoclustering { let n: i32 = X.rows let p: i32 = X.cols + let xdata: ptr = X.data - let row_labels: ptr = array_new_f32(n) as ptr - let col_labels: ptr = array_new_f32(p) as ptr + let row_labels: ptr = malloc((n as i64) * 4 + 64) as ptr + let col_labels: ptr = malloc((p as i64) * 4 + 64) as ptr - # Simplified: cluster rows and columns together by value magnitude - let mut row_sums: ptr = array_new_f32(n) + # Row and column sums, and the inverse square roots the normalization uses. + let r_inv: ptr = array_new_f32(n) + let c_inv: ptr = array_new_f32(p) + let c_sum: ptr = array_new_f32(p) for i in 0 to n { + let row: ptr = xdata + i * p let mut s: f32 = 0.0 - for j in 0 to p { s = s + matrix_at(X, i, j) } - row_sums[i] = s + let mut j: i32 = 0 + while j < p { + let v: f32 = row[j] + s = s + v + c_sum[j] = c_sum[j] + v + j = j + 1 + } + if s < 0.0000000001 { s = 0.0000000001 } + r_inv[i] = 1.0 / (sqrt((s) as f64) as f32) + } + for j2 in 0 to p { + let mut s2: f32 = c_sum[j2] + if s2 < 0.0000000001 { s2 = 0.0000000001 } + c_inv[j2] = 1.0 / (sqrt((s2) as f64) as f32) } - let mut min_val: f32 = row_sums[0] - let mut max_val: f32 = row_sums[0] - for i in 0 to n { - if row_sums[i] < min_val { min_val = row_sums[i] } - if row_sums[i] > max_val { max_val = row_sums[i] } + let An: ptr = array_new_f32(n * p) + for i2 in 0 to n { + let src: ptr = xdata + i2 * p + let dst: ptr = An + i2 * p + let ri: f32 = r_inv[i2] + let mut j3: i32 = 0 + while j3 < p { + dst[j3] = ri * src[j3] * c_inv[j3] + j3 = j3 + 1 + } } - let range: f32 = max_val - min_val - for i in 0 to n { - if range > 0.0000000001 { - row_labels[i] = ((row_sums[i] - min_val) / range * (n_clusters as f32)) as i32 - if row_labels[i] >= n_clusters { row_labels[i] = n_clusters - 1 } - } else { - row_labels[i] = 0 + + # The leading pair is trivial, so ceil(log2(k)) more are needed after it. + let mut n_keep: i32 = 0 + let mut reach: i32 = 1 + while reach < n_clusters { + reach = reach * 2 + n_keep = n_keep + 1 + } + if n_keep < 1 { n_keep = 1 } + let mut n_vecs: i32 = n_keep + 1 + if n_vecs > p { n_vecs = p } + + let V: ptr = array_new_f32(p * n_vecs) + _bic_right_singular(An, n, p, n_vecs, V) + + # The matching left vectors, as A v, which is what the normalized matrix + # maps the right vectors to. + let U: ptr = array_new_f32(n * n_vecs) + blas_matmat(An, V, U, n, n_vecs, p, 1.0, 0.0) + for c in 0 to n_vecs { + let mut norm: f32 = 0.0 + let mut i3: i32 = 0 + while i3 < n { + let v: f32 = U[i3 * n_vecs + c] + norm = norm + v * v + i3 = i3 + 1 + } + norm = sqrt((norm) as f64) as f32 + if norm > 0.0000000001 { + let inv: f32 = 1.0 / norm + let mut i4: i32 = 0 + while i4 < n { + U[i4 * n_vecs + c] = U[i4 * n_vecs + c] * inv + i4 = i4 + 1 + } } } - for j in 0 to p { - let mut s: f32 = 0.0 - for i in 0 to n { s = s + matrix_at(X, i, j) } - if range > 0.0000000001 { - col_labels[j] = ((s - min_val) / range * (n_clusters as f32)) as i32 - if col_labels[j] >= n_clusters { col_labels[j] = n_clusters - 1 } - } else { - col_labels[j] = 0 + # Rows and columns land in one space, scaled back by the same square roots, + # and one k-means run over the stack assigns both. + let d: i32 = n_vecs - 1 + let m: i32 = n + p + let Z: ptr = array_new_f32(m * d) + for i5 in 0 to n { + let ri2: f32 = r_inv[i5] + let mut j4: i32 = 0 + while j4 < d { + Z[i5 * d + j4] = ri2 * U[i5 * n_vecs + j4 + 1] + j4 = j4 + 1 + } + } + for j5 in 0 to p { + let cj: f32 = c_inv[j5] + let mut j6: i32 = 0 + while j6 < d { + Z[(n + j5) * d + j6] = cj * V[j5 * n_vecs + j6 + 1] + j6 = j6 + 1 } } - array_free_f32(row_sums) + let labels: ptr = malloc((m as i64) * 4 + 64) as ptr + _bic_kmeans(Z, m, d, n_clusters, labels) + for i6 in 0 to n { row_labels[i6] = labels[i6] } + for j7 in 0 to p { col_labels[j7] = labels[n + j7] } + + free(labels as ptr) + array_free_f32(Z) + array_free_f32(U) + array_free_f32(V) + array_free_f32(An) + array_free_f32(r_inv) + array_free_f32(c_inv) + array_free_f32(c_sum) return SpectralCoclustering { row_labels: row_labels, @@ -629,76 +1028,149 @@ export struct NeighborhoodComponentsAnalysis { fitted: bool } +# Neighbourhood components analysis, after Goldberger et al. (2004), which is +# what scikit-learn's NeighborhoodComponentsAnalysis maximizes. +# +# The objective is the expected number of correctly classified points under +# stochastic nearest neighbours: with p_ik proportional to exp of minus the +# squared distance between the projected points, the score of a point is the +# probability mass its own class receives, and the gradient of that sum with +# respect to the projection is +# +# 2 * sum_i ( p_i * sum_k p_ik C_ik - sum_{k in class of i} p_ik C_ik ) * T +# +# with C_ik the outer product of the difference of the two rows. +# +# What stood here computed neither. Its distance summed the squares of the +# individual products rather than squaring the projected difference, which is +# not the distance in the projected space unless the projection has one +# column, and its gradient kept only the second of the two terms above, so it +# pulled same-class points together with nothing pushing back. export function nca_fit(X: Matrix, y: ptr, n_components: i32, epochs: i32, lr: f32) -> NeighborhoodComponentsAnalysis { let n: i32 = X.rows let p: i32 = X.cols + let d: i32 = n_components + let xdata: ptr = X.data - let transform: ptr > = malloc((p as i64) * 8) as ptr > + # The projection, held flat while it is being fitted and handed back a row + # at a time at the end. + let T: ptr = array_new_f32(p * d) for i in 0 to p { - transform[i] = array_new_f32(n_components) - for j in 0 to n_components { - if i == j { - transform[i][j] = 1.0 - } else { - transform[i][j] = 0.0 - } + let mut j: i32 = 0 + while j < d { + if i == j { T[i * d + j] = 1.0 } + else { T[i * d + j] = 0.0 } + j = j + 1 } } - # Simplified: gradient descent on leave-one-out NN score + let proj: ptr = array_new_f32(n * d) + let prob: ptr = array_new_f32(n) + let W: ptr = array_new_f32(p * p) + let grad: ptr = array_new_f32(p * d) + let diff: ptr = array_new_f32(p) + for epoch in 0 to epochs { - let grad: ptr > = malloc((p as i64) * 8) as ptr > - for i in 0 to p { - grad[i] = array_new_f32(n_components) - for j in 0 to n_components { grad[i][j] = 0.0 } - } + # One projection of the whole design per epoch. + blas_matmat(xdata, T, proj, n, d, p, 1.0, 0.0) + for i in 0 to p * p { W[i] = 0.0 } for i in 0 to n { - # Compute softmax weights for all other samples + let pi: ptr = proj + i * d + let yi: f32 = y[i] + + # The softmax over the other points, and the mass its own class + # receives. let mut total: f32 = 0.0 - let dists: ptr = array_new_f32(n) - for k in 0 to n { - if k == i { dists[k] = 0.0; continue } - let mut d: f32 = 0.0 - for a in 0 to p { - let diff: f32 = matrix_at(X, i, a) - matrix_at(X, k, a) - for b in 0 to n_components { - d = d + diff * transform[a][b] * diff * transform[a][b] + let mut same: f32 = 0.0 + let mut k: i32 = 0 + while k < n { + if k != i { + let pk: ptr = proj + k * d + let mut dist: f32 = 0.0 + let mut b: i32 = 0 + while b < d { + let dv: f32 = pi[b] - pk[b] + dist = dist + dv * dv + b = b + 1 } + let w: f32 = exp((0.0 - dist) as f64) as f32 + prob[k] = w + total = total + w + if y[k] == yi { same = same + w } + } else { + prob[k] = 0.0 } - let w: f32 = exp((-d) as f64) as f32 - dists[k] = w - total = total + w - } - - if total < 0.0000000001 { array_free_f32(dists); continue } - - # Gradient: encourage same-class neighbors - for k in 0 to n { - if k == i { continue } - let p_ik: f32 = dists[k] / total - if y[i] == y[k] { - for a in 0 to p { - let diff: f32 = matrix_at(X, i, a) - matrix_at(X, k, a) - for b in 0 to n_components { - grad[a][b] = grad[a][b] - p_ik * 2.0 * diff * transform[a][b] * diff + k = k + 1 + } + if total < 0.0000000001 { continue } + + let inv_total: f32 = 1.0 / total + let p_i: f32 = same * inv_total + + # Each pair contributes its outer product, weighted by how much + # the pair helps or hurts the score of this point. + let row_i: ptr = xdata + i * p + let mut k2: i32 = 0 + while k2 < n { + if k2 != i { + let p_ik: f32 = prob[k2] * inv_total + let mut coeff: f32 = p_i * p_ik + if y[k2] == yi { coeff = coeff - p_ik } + if coeff > 0.0000001 || coeff < 0.0 - 0.0000001 { + let row_k: ptr = xdata + k2 * p + let mut a: i32 = 0 + while a < p { + diff[a] = row_i[a] - row_k[a] + a = a + 1 + } + let mut a2: i32 = 0 + while a2 < p { + let da: f32 = coeff * diff[a2] + let wrow: ptr = W + a2 * p + let mut b2: i32 = 0 + while b2 < p { + wrow[b2] = wrow[b2] + da * diff[b2] + b2 = b2 + 1 + } + a2 = a2 + 1 } } } + k2 = k2 + 1 } + } - array_free_f32(dists) + # Ascend the objective. The step is scaled by the sample count so a + # learning rate means the same thing whatever the size of the data. + blas_matmat(W, T, grad, p, d, p, 1.0, 0.0) + let step: f32 = 2.0 * lr / (n as f32) + let mut t: i32 = 0 + let total_t: i32 = p * d + while t < total_t { + T[t] = T[t] + step * grad[t] + t = t + 1 } + } - for a in 0 to p { - for b in 0 to n_components { - transform[a][b] = transform[a][b] - lr * grad[a][b] / (n as f32) - } - array_free_f32(grad[a]) + let transform: ptr > = malloc((p as i64) * 8) as ptr > + for i2 in 0 to p { + let row: ptr = array_new_f32(d) + let mut j2: i32 = 0 + while j2 < d { + row[j2] = T[i2 * d + j2] + j2 = j2 + 1 } - free(grad as ptr) + transform[i2] = row } + array_free_f32(T) + array_free_f32(proj) + array_free_f32(prob) + array_free_f32(W) + array_free_f32(grad) + array_free_f32(diff) + return NeighborhoodComponentsAnalysis { transformation: transform, n_features: p, diff --git a/tests/test_outlier_and_metric_learning.flow b/tests/test_outlier_and_metric_learning.flow new file mode 100644 index 0000000..1a24b29 --- /dev/null +++ b/tests/test_outlier_and_metric_learning.flow @@ -0,0 +1,143 @@ +# EllipticEnvelope has to find the outliers, and NCA has to improve the +# separation it was asked to improve. Both implementations were stand-ins +# before, and nothing in the suite asked either of them for a result. + +import "lib/scikit/scikit.flow" + +extern { + function printf(fmt: string, ...) -> i32 +} + +# A tight cluster with a few points thrown far away from it. +function oml_fixture(n_inliers: i32, n_outliers: i32) -> Matrix { + let n: i32 = n_inliers + n_outliers + let X: Matrix = matrix_new(n, 2) + srand(7) + for i in 0 to n_inliers { + let a: f32 = ((rand() % 200) as f32) / 100.0 - 1.0 + let b: f32 = ((rand() % 200) as f32) / 100.0 - 1.0 + matrix_set(X, i, 0, a) + matrix_set(X, i, 1, b) + } + for k in 0 to n_outliers { + matrix_set(X, n_inliers + k, 0, 40.0 + (k as f32)) + matrix_set(X, n_inliers + k, 1, 40.0 + (k as f32)) + } + return X +} + +function test_envelope_finds_the_outliers() -> i32 { + println("Test: EllipticEnvelope marks the far points") + let n_in: i32 = 60 + let n_out: i32 = 6 + let X: Matrix = oml_fixture(n_in, n_out) + let contamination: f32 = (n_out as f32) / ((n_in + n_out) as f32) + let model: EllipticEnvelope = elliptic_envelope_fit(X, contamination, 10) + let pred: ptr = elliptic_envelope_predict(model, X) + + let mut missed: i32 = 0 + for k in 0 to n_out { + if pred[n_in + k] > 0.0 { missed = missed + 1 } + } + let mut false_alarms: i32 = 0 + for i in 0 to n_in { + if pred[i] < 0.0 { false_alarms = false_alarms + 1 } + } + + let mut failures: i32 = 0 + if missed > 0 { + printf(" FAIL: %d of %d outliers were called inliers\n", missed, n_out) + failures = failures + 1 + } + if false_alarms > 0 { + printf(" FAIL: %d inliers were called outliers\n", false_alarms) + failures = failures + 1 + } + if failures == 0 { + printf(" OK: all %d outliers found, no false alarms\n", n_out) + } + + array_free_f32(pred) + elliptic_envelope_free(model) + matrix_free(X) + return failures +} + +# The mean distance between points of different classes, over the mean +# distance within a class, in whatever space the argument is given in. +function oml_separation(Z: Matrix, y: ptr, n: i32) -> f32 { + let mut between: f32 = 0.0 + let mut n_between: f32 = 0.0 + let mut within: f32 = 0.0 + let mut n_within: f32 = 0.0 + for i in 0 to n { + for j in i + 1 to n { + let mut d: f32 = 0.0 + for c in 0 to Z.cols { + let diff: f32 = matrix_at(Z, i, c) - matrix_at(Z, j, c) + d = d + diff * diff + } + d = sqrt((d) as f64) as f32 + if y[i] == y[j] { + within = within + d + n_within = n_within + 1.0 + } else { + between = between + d + n_between = n_between + 1.0 + } + } + } + if n_within < 1.0 || within < 0.0000001 { return 0.0 } + return (between / n_between) / (within / n_within) +} + +function test_nca_improves_separation() -> i32 { + println("Test: NCA separates the classes better than the identity") + let n: i32 = 60 + let X: Matrix = matrix_new(n, 4) + let y: ptr = array_new_f32(n) + srand(11) + for i in 0 to n { + let cls: i32 = i % 2 + # The first feature carries the class, the rest are noise, so a + # projection that learns anything has to lean on the first. + let centre: f32 = (cls as f32) * 4.0 + matrix_set(X, i, 0, centre + ((rand() % 100) as f32) / 100.0) + matrix_set(X, i, 1, ((rand() % 400) as f32) / 100.0) + matrix_set(X, i, 2, ((rand() % 400) as f32) / 100.0) + matrix_set(X, i, 3, ((rand() % 400) as f32) / 100.0) + y[i] = cls as f32 + } + + let before: f32 = oml_separation(X, y, n) + let model: NeighborhoodComponentsAnalysis = nca_fit(X, y, 2, 30, 0.01) + let Z: Matrix = nca_transform(model, X) + let after: f32 = oml_separation(Z, y, n) + + let mut failures: i32 = 0 + if after <= before { + printf(" FAIL: separation %.4f before, %.4f after\n", before, after) + failures = failures + 1 + } else { + printf(" OK: separation %.4f before, %.4f after\n", before, after) + } + + matrix_free(Z) + nca_free(model) + array_free_f32(y) + matrix_free(X) + return failures +} + +function main() -> i32 { + println("=== Outlier detection and metric learning ===") + let mut failures: i32 = 0 + failures = failures + test_envelope_finds_the_outliers() + failures = failures + test_nca_improves_separation() + if failures > 0 { + printf("%d failure(s)\n", failures) + return 1 + } + println("All outlier and metric learning tests passed!") + return 0 +} diff --git a/tests/test_spectral_biclustering.flow b/tests/test_spectral_biclustering.flow new file mode 100644 index 0000000..4788714 --- /dev/null +++ b/tests/test_spectral_biclustering.flow @@ -0,0 +1,139 @@ +# Spectral co-clustering has to recover a block structure that is obvious by +# eye. The implementation it replaced bucketed row sums between the smallest +# and the largest, which cannot recover a checkerboard at all, and nothing in +# the suite noticed because nothing asked it to. + +import "lib/scikit/scikit.flow" + +extern { + function printf(fmt: string, ...) -> i32 +} + +# Two blocks: the first half of the rows goes with the first half of the +# columns, the second with the second, and the off-blocks are weak. +function coc_fixture(n: i32, p: i32) -> Matrix { + let X: Matrix = matrix_new(n, p) + let half_r: i32 = n / 2 + let half_c: i32 = p / 2 + for i in 0 to n { + for j in 0 to p { + let same: bool = (i < half_r && j < half_c) || (i >= half_r && j >= half_c) + if same { matrix_set(X, i, j, 10.0) } + else { matrix_set(X, i, j, 1.0) } + } + } + return X +} + +# The labels are arbitrary up to renaming, so agreement is checked as "the +# same partition", by requiring every member of a block to share one label and +# the two blocks to differ. +function coc_block_is_uniform(labels: ptr, from_idx: i32, to_idx: i32) -> bool { + let first: i32 = labels[from_idx] + for i in from_idx to to_idx { + if labels[i] != first { return false } + } + return true +} + +function test_recovers_two_blocks() -> i32 { + println("Test: SpectralCoclustering recovers two blocks") + let n: i32 = 12 + let p: i32 = 8 + let X: Matrix = coc_fixture(n, p) + let model: SpectralCoclustering = spectral_coclustering_fit(X, 2) + + let mut failures: i32 = 0 + if not coc_block_is_uniform(model.row_labels, 0, n / 2) { + println(" FAIL: the first row block did not get one label") + failures = failures + 1 + } + if not coc_block_is_uniform(model.row_labels, n / 2, n) { + println(" FAIL: the second row block did not get one label") + failures = failures + 1 + } + if model.row_labels[0] == model.row_labels[n - 1] { + println(" FAIL: the two row blocks landed in the same cluster") + failures = failures + 1 + } + if not coc_block_is_uniform(model.col_labels, 0, p / 2) { + println(" FAIL: the first column block did not get one label") + failures = failures + 1 + } + if not coc_block_is_uniform(model.col_labels, p / 2, p) { + println(" FAIL: the second column block did not get one label") + failures = failures + 1 + } + if model.col_labels[0] == model.col_labels[p - 1] { + println(" FAIL: the two column blocks landed in the same cluster") + failures = failures + 1 + } + + # A row block and the column block it belongs with have to agree, which is + # the whole point of embedding rows and columns in one space. + if model.row_labels[0] != model.col_labels[0] { + println(" FAIL: the first row block and first column block disagree") + failures = failures + 1 + } + + if failures == 0 { + printf(" OK: rows %d %d, columns %d %d\n", model.row_labels[0], model.row_labels[n - 1], model.col_labels[0], model.col_labels[p - 1]) + } + spectral_coclustering_free(model) + matrix_free(X) + return failures +} + +function test_biclustering_recovers_blocks() -> i32 { + println("Test: SpectralBiclustering recovers a checkerboard") + let n: i32 = 12 + let p: i32 = 8 + let X: Matrix = coc_fixture(n, p) + let model: SpectralBiclustering = spectral_biclustering_fit(X, 2, 2) + + let mut failures: i32 = 0 + if not coc_block_is_uniform(model.row_labels, 0, n / 2) { + println(" FAIL: the first row block did not get one label") + failures = failures + 1 + } + if not coc_block_is_uniform(model.row_labels, n / 2, n) { + println(" FAIL: the second row block did not get one label") + failures = failures + 1 + } + if model.row_labels[0] == model.row_labels[n - 1] { + println(" FAIL: the two row blocks landed in the same cluster") + failures = failures + 1 + } + if not coc_block_is_uniform(model.col_labels, 0, p / 2) { + println(" FAIL: the first column block did not get one label") + failures = failures + 1 + } + if not coc_block_is_uniform(model.col_labels, p / 2, p) { + println(" FAIL: the second column block did not get one label") + failures = failures + 1 + } + if model.col_labels[0] == model.col_labels[p - 1] { + println(" FAIL: the two column blocks landed in the same cluster") + failures = failures + 1 + } + + if failures == 0 { + printf(" OK: rows %d %d, columns %d %d\n", model.row_labels[0], model.row_labels[n - 1], model.col_labels[0], model.col_labels[p - 1]) + } + spectral_biclustering_free(model) + matrix_free(X) + return failures +} + +function main() -> i32 { + println("=== Spectral biclustering ===") + let mut failures: i32 = 0 + failures = failures + test_recovers_two_blocks() + failures = failures + test_biclustering_recovers_blocks() + if failures > 0 { + printf("%d failure(s)\n", failures) + return 1 + } + println("All spectral biclustering tests passed!") + return 0 +} From b49b5dd7bb9f27fbf0fc666832cd55fc418121f6 Mon Sep 17 00:00:00 2001 From: godofecht Date: Tue, 29 Sep 2026 12:58:12 +0100 Subject: [PATCH 2/2] Race the five rows that were doing half the job Five Flow functions covered part of what the scikit-learn class they are named after does, so the registry kept them out of the ranking and said why. The missing half is written now, as a second entry point beside each of them, and the five are raced. voting_classifier_full_fit and voting_regressor_full_fit fit the ensemble from the data, where the functions beside them take estimators that are already fitted and their cost is the vote alone. Three trees on both sides, because a vote of one is not what either library is being asked for. stacking_classifier_full_fit cross-validates its base estimators to build the meta-features, fits the meta-learner on those, and refits the base estimators on all of the data. Out of fold predictions are the point: a base estimator's prediction for a row it helped fit would tell the meta-learner more than it will know at prediction time. self_training_classifier_full_fit refits its base estimator on every round, which is the expensive half the function beside it leaves to its caller. incremental_pca_fit walks the design in batches, sized as scikit-learn sizes its own, where the partial fit beside it is one step of that. Racing incremental PCA turned up that it did not work. Its partial fit accumulated the same per-feature quantity for every component, so every component came out as the same vector, and the step it called a power iteration multiplied a component by itself rather than by a covariance. It keeps a running scatter matrix now, which combines across batches exactly once the shift between two means is accounted for, and takes its components from the leading eigenvectors of that. On the test fixture the leading component carried a spread of 640 against the second's 4, where the two were previously identical to the digit. Local timings against scikit-learn on the same data: incremental_pca 0.010 ms against 1.01 voting_classifier_full 0.053 against 2.23 self_training_classifier_full 0.266 against 19.83 stacking_classifier_full 0.529 against 16.07 voting_regressor_full 1.184 against 3.74 A test file covers all five by behaviour: the two ensembles and the stack have to classify or regress the fixture, self training has to label the unlabelled rows and get them right, and the batched decomposition has to see every row and put more spread in the leading component than the second. It was checked as a real gate by inverting an assertion. --- benchmarks/bench_estimators_sklearn.py | 17 +- benchmarks/estimator_coverage.json | 367 +++++++++++++++++- benchmarks/estimator_coverage.py | 69 +++- benchmarks/generate_estimator_bench.py | 4 + benchmarks/generated/bench_estimators_02.flow | 56 +-- benchmarks/generated/bench_estimators_03.flow | 54 +-- benchmarks/generated/bench_estimators_04.flow | 56 ++- benchmarks/generated/bench_estimators_05.flow | 56 +-- benchmarks/generated/bench_estimators_06.flow | 56 ++- benchmarks/generated/bench_estimators_07.flow | 58 +-- benchmarks/generated/bench_estimators_08.flow | 97 ++--- benchmarks/generated/bench_estimators_09.flow | 128 ++++++ lib/scikit/decomposition.flow | 262 +++++++++++-- lib/scikit/ensemble.flow | 228 +++++++++++ lib/scikit/semi_supervised.flow | 113 ++++++ tests/test_full_ensembles.flow | 232 +++++++++++ 16 files changed, 1577 insertions(+), 276 deletions(-) create mode 100644 tests/test_full_ensembles.flow diff --git a/benchmarks/bench_estimators_sklearn.py b/benchmarks/bench_estimators_sklearn.py index a516f87..48b46e4 100644 --- a/benchmarks/bench_estimators_sklearn.py +++ b/benchmarks/bench_estimators_sklearn.py @@ -51,10 +51,16 @@ def _constructors() -> dict: "RFE": lambda c: c(tree_c()), "RFECV": lambda c: c(tree_c()), "SequentialFeatureSelector": lambda c: c(tree_c(), n_features_to_select=2), - "StackingClassifier": lambda c: c(estimators=[("t", tree_c())]), + # Three base estimators on both sides, because a vote or a stack of + # one is not the thing either library is being asked to do. + "StackingClassifier": lambda c: c( + estimators=[("t1", tree_c()), ("t2", tree_c()), ("t3", tree_c())], cv=5), "StackingRegressor": lambda c: c(estimators=[("r", Ridge())]), - "VotingClassifier": lambda c: c(estimators=[("t", tree_c())]), - "VotingRegressor": lambda c: c(estimators=[("r", Ridge())]), + "VotingClassifier": lambda c: c( + estimators=[("t1", tree_c()), ("t2", tree_c()), ("t3", tree_c())]), + "VotingRegressor": lambda c: c( + estimators=[("t1", tree_r()), ("t2", tree_r()), ("t3", tree_r())]), + "SelfTrainingClassifier": lambda c: c(LogisticRegression(max_iter=50)), "SparseCoder": lambda c: c(dictionary=np.eye(4)), # nu=0.5 is infeasible for this class balance. "NuSVC": lambda c: c(nu=0.1), @@ -152,6 +158,11 @@ def main() -> int: elif fit_input == "dicts": first = [{f"f{j}": float(v) for j, v in enumerate(row)} for row in X] second = None + elif fit_input == "semi_labels": + # Every third row keeps its label and the rest are unlabelled, + # which is what the Flow harness hands its own side. + masked = np.array([val if i % 3 == 0 else -1 for i, val in enumerate(y)]) + first, second = X, masked elif fit_input == "labelsets": first = [tuple(int(v) for v in row) for row in y] second = None diff --git a/benchmarks/estimator_coverage.json b/benchmarks/estimator_coverage.json index bd666a8..32ad7c2 100644 --- a/benchmarks/estimator_coverage.json +++ b/benchmarks/estimator_coverage.json @@ -1,9 +1,9 @@ { "schema_version": 1, "counts": { - "estimators": 203, - "runnable": 170, - "shaped": 22, + "estimators": 208, + "runnable": 171, + "shaped": 26, "flow_only": 6, "different_shape": 5 }, @@ -3544,6 +3544,61 @@ "bucket": "runnable", "sklearn_estimator": "HuberRegressor" }, + { + "flow_estimator": "incremental_pca", + "module": "decomposition.flow", + "fit": { + "name": "incremental_pca_fit", + "returns": "IncrementalPCA", + "parameters": [ + { + "name": "X", + "type": "Matrix" + }, + { + "name": "n_components", + "type": "i32" + }, + { + "name": "batch_size", + "type": "i32" + } + ] + }, + "companions": { + "transform": { + "name": "incremental_pca_transform", + "returns": "Matrix", + "parameters": [ + { + "name": "model", + "type": "IncrementalPCA" + }, + { + "name": "X", + "type": "Matrix" + } + ] + }, + "free": { + "name": "incremental_pca_free", + "returns": "void", + "parameters": [ + { + "name": "model", + "type": "IncrementalPCA" + } + ] + } + }, + "parameters": [ + "X", + "n_components", + "batch_size" + ], + "bucket": "runnable", + "sklearn_estimator": "IncrementalPCA" + }, { "flow_estimator": "incremental_pca_partial", "module": "decomposition.flow", @@ -3567,7 +3622,7 @@ "X" ], "bucket": "different_shape", - "reason": "one partial_fit step over one batch, where IncrementalPCA.fit walks the whole design in batches", + "reason": "one partial_fit step over one batch. incremental_pca_fit walks the whole design in batches and is the row that races IncrementalPCA", "sklearn_estimator": "IncrementalPCA" }, { @@ -10773,9 +10828,79 @@ "max_iter" ], "bucket": "different_shape", - "reason": "takes class probabilities as input, so it is the labelling loop alone while SelfTrainingClassifier also fits the base estimator on every round", + "reason": "takes class probabilities as input, so it is the labelling loop alone. self_training_classifier_full_fit refits the base estimator every round and is the row that races the class", "sklearn_estimator": "SelfTrainingClassifier" }, + { + "flow_estimator": "self_training_classifier_full", + "module": "semi_supervised.flow", + "fit": { + "name": "self_training_classifier_full_fit", + "returns": "SelfTrainingClassifier", + "parameters": [ + { + "name": "X", + "type": "Matrix" + }, + { + "name": "y_initial", + "type": "ptr" + }, + { + "name": "n_classes", + "type": "i32" + }, + { + "name": "threshold", + "type": "f32" + }, + { + "name": "max_iter", + "type": "i32" + }, + { + "name": "base_iter", + "type": "i32" + }, + { + "name": "lr", + "type": "f32" + } + ] + }, + "companions": {}, + "parameters": [ + "X", + "y_initial", + "n_classes", + "threshold", + "max_iter", + "base_iter", + "lr" + ], + "bucket": "shaped", + "sklearn_estimator": "SelfTrainingClassifier", + "shape": { + "dataset": "classification", + "flow_preamble": [ + "let st_labels: ptr = malloc((n_c as i64) * 4 + 64) as ptr", + "for i in 0 to n_c {", + " if i % 3 == 0 { st_labels[i] = yi_c[i] }", + " else { st_labels[i] = 0 - 1 }", + "}" + ], + "flow_fit": [ + "X_c", + "st_labels", + "3", + "0.7", + "10", + "50", + "0.1" + ], + "sklearn_input": "semi_labels" + } + }, { "flow_estimator": "sequential_feature_selector", "module": "feature_selection.flow", @@ -11725,9 +11850,109 @@ "n_iter" ], "bucket": "different_shape", - "reason": "takes the base estimators' predictions as input, so it is the meta-learner alone while StackingClassifier also fits the base estimators and cross-validates them", + "reason": "takes the base estimators' predictions as input, so it is the meta-learner alone. stacking_classifier_full_fit does the whole job and is the row that races StackingClassifier", "sklearn_estimator": "StackingClassifier" }, + { + "flow_estimator": "stacking_classifier_full", + "module": "ensemble.flow", + "fit": { + "name": "stacking_classifier_full_fit", + "returns": "StackingClassifierFull", + "parameters": [ + { + "name": "X", + "type": "Matrix" + }, + { + "name": "y", + "type": "ptr" + }, + { + "name": "n_classes", + "type": "i32" + }, + { + "name": "n_base", + "type": "i32" + }, + { + "name": "max_depth", + "type": "i32" + }, + { + "name": "seed", + "type": "i32" + }, + { + "name": "lr", + "type": "f32" + }, + { + "name": "n_iter", + "type": "i32" + } + ] + }, + "companions": { + "predict": { + "name": "stacking_classifier_full_predict", + "returns": "ptr", + "parameters": [ + { + "name": "model", + "type": "StackingClassifierFull" + }, + { + "name": "X", + "type": "Matrix" + } + ] + }, + "free": { + "name": "stacking_classifier_full_free", + "returns": "void", + "parameters": [ + { + "name": "model", + "type": "StackingClassifierFull" + } + ] + } + }, + "parameters": [ + "X", + "y", + "n_classes", + "n_base", + "max_depth", + "seed", + "lr", + "n_iter" + ], + "bucket": "shaped", + "sklearn_estimator": "StackingClassifier", + "shape": { + "dataset": "classification", + "flow_fit": [ + "X_c", + "y_c", + "3", + "3", + "3", + "42", + "0.01", + "50" + ], + "flow_work": [ + "X_c" + ], + "flow_work_fn": "stacking_classifier_full_predict", + "flow_work_returns": "ptr", + "flow_free_fn": "stacking_classifier_full_free", + "sklearn_input": "X" + } + }, { "flow_estimator": "stacking_regressor", "module": "ensemble.flow", @@ -12571,9 +12796,78 @@ "voting" ], "bucket": "different_shape", - "reason": "takes trees that are already fitted, so its fit is the vote alone while VotingClassifier fits every estimator it is given", + "reason": "takes trees that are already fitted, so its fit is the vote alone. voting_classifier_full_fit fits the ensemble and is the row that races VotingClassifier", "sklearn_estimator": "VotingClassifier" }, + { + "flow_estimator": "voting_classifier_full", + "module": "ensemble.flow", + "fit": { + "name": "voting_classifier_full_fit", + "returns": "VotingClassifier", + "parameters": [ + { + "name": "X", + "type": "Matrix" + }, + { + "name": "y", + "type": "ptr" + }, + { + "name": "n_classes", + "type": "i32" + }, + { + "name": "n_estimators", + "type": "i32" + }, + { + "name": "max_depth", + "type": "i32" + }, + { + "name": "seed", + "type": "i32" + }, + { + "name": "voting", + "type": "i32" + } + ] + }, + "companions": {}, + "parameters": [ + "X", + "y", + "n_classes", + "n_estimators", + "max_depth", + "seed", + "voting" + ], + "bucket": "shaped", + "sklearn_estimator": "VotingClassifier", + "shape": { + "dataset": "classification", + "flow_fit": [ + "X_c", + "y_c", + "3", + "3", + "3", + "42", + "0" + ], + "flow_work": [ + "X_c" + ], + "flow_work_fn": "voting_classifier_predict", + "flow_work_returns": "ptr", + "flow_free_fn": "voting_classifier_free", + "sklearn_input": "X" + } + }, { "flow_estimator": "voting_regressor", "module": "ensemble.flow", @@ -12627,8 +12921,65 @@ "weights" ], "bucket": "different_shape", - "reason": "takes regressors that are already fitted, so its fit is the average alone while VotingRegressor fits every estimator it is given", + "reason": "takes regressors that are already fitted, so its fit is the average alone. voting_regressor_full_fit fits the ensemble and is the row that races VotingRegressor", "sklearn_estimator": "VotingRegressor" + }, + { + "flow_estimator": "voting_regressor_full", + "module": "ensemble.flow", + "fit": { + "name": "voting_regressor_full_fit", + "returns": "VotingRegressor", + "parameters": [ + { + "name": "X", + "type": "Matrix" + }, + { + "name": "y", + "type": "ptr" + }, + { + "name": "n_estimators", + "type": "i32" + }, + { + "name": "max_depth", + "type": "i32" + }, + { + "name": "seed", + "type": "i32" + } + ] + }, + "companions": {}, + "parameters": [ + "X", + "y", + "n_estimators", + "max_depth", + "seed" + ], + "bucket": "shaped", + "sklearn_estimator": "VotingRegressor", + "shape": { + "dataset": "regression", + "flow_fit": [ + "X_r", + "y_r", + "3", + "3", + "42" + ], + "flow_work": [ + "X_r" + ], + "flow_work_fn": "voting_regressor_predict", + "flow_work_returns": "ptr", + "flow_free_fn": "voting_regressor_free", + "sklearn_input": "X" + } } ] } diff --git a/benchmarks/estimator_coverage.py b/benchmarks/estimator_coverage.py index bc3499f..8e47ba3 100644 --- a/benchmarks/estimator_coverage.py +++ b/benchmarks/estimator_coverage.py @@ -182,17 +182,20 @@ # same objection the simplified rows carry. The wording says which part. DIFFERENT_JOB: dict[str, str] = { "stacking_classifier": "takes the base estimators' predictions as input, so it is the " - "meta-learner alone while StackingClassifier also fits the base " - "estimators and cross-validates them", + "meta-learner alone. stacking_classifier_full_fit does the whole " + "job and is the row that races StackingClassifier", "self_training_classifier": "takes class probabilities as input, so it is the labelling " - "loop alone while SelfTrainingClassifier also fits the base " - "estimator on every round", - "voting_classifier": "takes trees that are already fitted, so its fit is the vote alone " - "while VotingClassifier fits every estimator it is given", + "loop alone. self_training_classifier_full_fit refits the base " + "estimator every round and is the row that races the class", + "voting_classifier": "takes trees that are already fitted, so its fit is the vote alone. " + "voting_classifier_full_fit fits the ensemble and is the row that " + "races VotingClassifier", "voting_regressor": "takes regressors that are already fitted, so its fit is the average " - "alone while VotingRegressor fits every estimator it is given", - "incremental_pca_partial": "one partial_fit step over one batch, where IncrementalPCA.fit " - "walks the whole design in batches", + "alone. voting_regressor_full_fit fits the ensemble and is the row " + "that races VotingRegressor", + "incremental_pca_partial": "one partial_fit step over one batch. incremental_pca_fit walks " + "the whole design in batches and is the row that races " + "IncrementalPCA", } # One corpus for the two text vectorizers, read by the Flow generator and by @@ -314,6 +317,47 @@ # Two rows whose work function is named for what it returns rather than # predict or transform, so the generic path found nothing to time and the # fit alone fell under the clock's floor. + "voting_classifier_full": { + "dataset": "classification", + "flow_fit": ["X_c", "y_c", "3", "3", "3", "42", "0"], + "flow_work": ["X_c"], + "flow_work_fn": "voting_classifier_predict", + "flow_work_returns": "ptr", + "flow_free_fn": "voting_classifier_free", + "sklearn_input": "X", + }, + "voting_regressor_full": { + "dataset": "regression", + "flow_fit": ["X_r", "y_r", "3", "3", "42"], + "flow_work": ["X_r"], + "flow_work_fn": "voting_regressor_predict", + "flow_work_returns": "ptr", + "flow_free_fn": "voting_regressor_free", + "sklearn_input": "X", + }, + "stacking_classifier_full": { + "dataset": "classification", + "flow_fit": ["X_c", "y_c", "3", "3", "3", "42", "0.01", "50"], + "flow_work": ["X_c"], + "flow_work_fn": "stacking_classifier_full_predict", + "flow_work_returns": "ptr", + "flow_free_fn": "stacking_classifier_full_free", + "sklearn_input": "X", + }, + "self_training_classifier_full": { + # Every third row keeps its label and the rest are unlabelled, on both + # sides, which is the shape the estimator exists for. + "dataset": "classification", + "flow_preamble": [ + "let st_labels: ptr = malloc((n_c as i64) * 4 + 64) as ptr", + "for i in 0 to n_c {", + " if i % 3 == 0 { st_labels[i] = yi_c[i] }", + " else { st_labels[i] = 0 - 1 }", + "}", + ], + "flow_fit": ["X_c", "st_labels", "3", "0.7", "10", "50", "0.1"], + "sklearn_input": "semi_labels", + }, "kernel_density": { "dataset": "classification", "flow_fit": ["X_c", "0.5", "0"], @@ -431,6 +475,13 @@ # Flow names whose scikit-learn counterpart is not a case change away. ALIASES: dict[str, str] = { + # The full fits, which do what the scikit-learn class does rather than the + # part of it the older function covers. Both entries stay in the registry: + # the older one keeps its reason for not being raced. + "voting_classifier_full": "VotingClassifier", + "voting_regressor_full": "VotingRegressor", + "stacking_classifier_full": "StackingClassifier", + "self_training_classifier_full": "SelfTrainingClassifier", # decomposition.flow, and it takes n_topics and eta, so this is Latent # Dirichlet Allocation. discriminant_analysis.flow holds the other one. "lda": "LatentDirichletAllocation", diff --git a/benchmarks/generate_estimator_bench.py b/benchmarks/generate_estimator_bench.py index b4f13b5..0b0eaf8 100644 --- a/benchmarks/generate_estimator_bench.py +++ b/benchmarks/generate_estimator_bench.py @@ -224,6 +224,10 @@ def shaped_block(entry: dict) -> str: ret = entry["fit"]["returns"] call = f"{entry['fit']['name']}({', '.join(shape['flow_fit'])})" free = entry["companions"].get("free") + # A recipe can name the free function, for a row whose model type is not + # named after its own fit. + if shape.get("flow_free_fn"): + free = {"name": shape["flow_free_fn"], "parameters": [{"name": "model"}]} # A fit that takes a composed object built outside the timing mutates that # object and hands it back, so freeing the result on every repeat would free # what the next repeat is about to read. Those rows free once, after the diff --git a/benchmarks/generated/bench_estimators_02.flow b/benchmarks/generated/bench_estimators_02.flow index 03364b2..83b1317 100644 --- a/benchmarks/generated/bench_estimators_02.flow +++ b/benchmarks/generated/bench_estimators_02.flow @@ -440,6 +440,35 @@ function main() -> i32 { printf("ESTIMATOR|huber_regressor|%.9f|%.9f|%d|ok\n", ms_between(t0, t1) / (reps as f32), ms_between(t2, t3) / (reps as f32), reps) fflush(null) + # ---- incremental_pca (unsupervised) ---- + t0 = flow_now_ns() + let probe_incremental_pca: IncrementalPCA = incremental_pca_fit(X_c, 2, 32) + t1 = flow_now_ns() + incremental_pca_free(probe_incremental_pca) + reps = 1 + if (t1 - t0) < 200000 { reps = 200 } + elif (t1 - t0) < 2000000 { reps = 20 } + elif (t1 - t0) < 20000000 { reps = 5 } + t0 = flow_now_ns() + for rep in 0 to reps { + let m_incremental_pca: IncrementalPCA = incremental_pca_fit(X_c, 2, 32) + incremental_pca_free(m_incremental_pca) + } + t1 = flow_now_ns() + let fitted_incremental_pca: IncrementalPCA = incremental_pca_fit(X_c, 2, 32) + t2 = flow_now_ns() + for rep2 in 0 to reps { + let o_incremental_pca: Matrix = incremental_pca_transform(fitted_incremental_pca, X_c) + if o_incremental_pca.rows > 0 { + if o_incremental_pca.cols > 0 { sink = sink + o_incremental_pca.data[0] } + } + matrix_free(o_incremental_pca) + } + t3 = flow_now_ns() + incremental_pca_free(fitted_incremental_pca) + printf("ESTIMATOR|incremental_pca|%.9f|%.9f|%d|ok\n", ms_between(t0, t1) / (reps as f32), ms_between(t2, t3) / (reps as f32), reps) + fflush(null) + # ---- isolation_forest (unsupervised) ---- t0 = flow_now_ns() let probe_isolation_forest: IsolationForest = isolation_forest_fit(X_c, 10, 5, 42) @@ -561,33 +590,6 @@ function main() -> i32 { printf("ESTIMATOR|kbins_discretizer|%.9f|%.9f|%d|ok\n", ms_between(t0, t1) / (reps as f32), ms_between(t2, t3) / (reps as f32), reps) fflush(null) - # ---- kernel_density (classification, written out) ---- - t0 = flow_now_ns() - let probe_kernel_density: KernelDensity = kernel_density_fit(X_c, 0.5, 0) - t1 = flow_now_ns() - kernel_density_free(probe_kernel_density) - reps = 1 - if (t1 - t0) < 200000 { reps = 200 } - elif (t1 - t0) < 2000000 { reps = 20 } - elif (t1 - t0) < 20000000 { reps = 5 } - t0 = flow_now_ns() - for rep in 0 to reps { - let m_kernel_density: KernelDensity = kernel_density_fit(X_c, 0.5, 0) - kernel_density_free(m_kernel_density) - } - t1 = flow_now_ns() - let fitted_kernel_density: KernelDensity = kernel_density_fit(X_c, 0.5, 0) - t2 = flow_now_ns() - for rep2 in 0 to reps { - let o_kernel_density: ptr = kernel_density_score_samples(fitted_kernel_density, X_c) - sink = sink + o_kernel_density[0] - array_free_f32(o_kernel_density) - } - t3 = flow_now_ns() - kernel_density_free(fitted_kernel_density) - printf("ESTIMATOR|kernel_density|%.9f|%.9f|%d|ok\n", ms_between(t0, t1) / (reps as f32), ms_between(t2, t3) / (reps as f32), reps) - fflush(null) - for i in 0 to n_c { array_free_f32(Y_label_rows[i]) } free(Y_label_rows as ptr) matrix_free(Y_labels) diff --git a/benchmarks/generated/bench_estimators_03.flow b/benchmarks/generated/bench_estimators_03.flow index 9086a1b..831d1fd 100644 --- a/benchmarks/generated/bench_estimators_03.flow +++ b/benchmarks/generated/bench_estimators_03.flow @@ -82,6 +82,33 @@ function main() -> i32 { let mut reps: i32 = 1 let mut sink: f32 = 0.0 + # ---- kernel_density (classification, written out) ---- + t0 = flow_now_ns() + let probe_kernel_density: KernelDensity = kernel_density_fit(X_c, 0.5, 0) + t1 = flow_now_ns() + kernel_density_free(probe_kernel_density) + reps = 1 + if (t1 - t0) < 200000 { reps = 200 } + elif (t1 - t0) < 2000000 { reps = 20 } + elif (t1 - t0) < 20000000 { reps = 5 } + t0 = flow_now_ns() + for rep in 0 to reps { + let m_kernel_density: KernelDensity = kernel_density_fit(X_c, 0.5, 0) + kernel_density_free(m_kernel_density) + } + t1 = flow_now_ns() + let fitted_kernel_density: KernelDensity = kernel_density_fit(X_c, 0.5, 0) + t2 = flow_now_ns() + for rep2 in 0 to reps { + let o_kernel_density: ptr = kernel_density_score_samples(fitted_kernel_density, X_c) + sink = sink + o_kernel_density[0] + array_free_f32(o_kernel_density) + } + t3 = flow_now_ns() + kernel_density_free(fitted_kernel_density) + printf("ESTIMATOR|kernel_density|%.9f|%.9f|%d|ok\n", ms_between(t0, t1) / (reps as f32), ms_between(t2, t3) / (reps as f32), reps) + fflush(null) + # ---- kernel_pca (unsupervised) ---- t0 = flow_now_ns() let probe_kernel_pca: KernelPCA = kernel_pca_fit(X_c, 2, 0, 0.1, 2, 0.0) @@ -578,33 +605,6 @@ function main() -> i32 { printf("ESTIMATOR|lasso_lars|%.9f|%.9f|%d|ok\n", ms_between(t0, t1) / (reps as f32), ms_between(t2, t3) / (reps as f32), reps) fflush(null) - # ---- lasso_lars_ic (regression) ---- - t0 = flow_now_ns() - let probe_lasso_lars_ic: LassoLarsIC = lasso_lars_ic_fit(X_r, y_r, 0) - t1 = flow_now_ns() - lasso_lars_ic_free(probe_lasso_lars_ic) - reps = 1 - if (t1 - t0) < 200000 { reps = 200 } - elif (t1 - t0) < 2000000 { reps = 20 } - elif (t1 - t0) < 20000000 { reps = 5 } - t0 = flow_now_ns() - for rep in 0 to reps { - let m_lasso_lars_ic: LassoLarsIC = lasso_lars_ic_fit(X_r, y_r, 0) - lasso_lars_ic_free(m_lasso_lars_ic) - } - t1 = flow_now_ns() - let fitted_lasso_lars_ic: LassoLarsIC = lasso_lars_ic_fit(X_r, y_r, 0) - t2 = flow_now_ns() - for rep2 in 0 to reps { - let o_lasso_lars_ic: ptr = lasso_lars_ic_predict(fitted_lasso_lars_ic, X_r) - sink = sink + o_lasso_lars_ic[0] - array_free_f32(o_lasso_lars_ic) - } - t3 = flow_now_ns() - lasso_lars_ic_free(fitted_lasso_lars_ic) - printf("ESTIMATOR|lasso_lars_ic|%.9f|%.9f|%d|ok\n", ms_between(t0, t1) / (reps as f32), ms_between(t2, t3) / (reps as f32), reps) - fflush(null) - for i in 0 to n_c { array_free_f32(Y_label_rows[i]) } free(Y_label_rows as ptr) matrix_free(Y_labels) diff --git a/benchmarks/generated/bench_estimators_04.flow b/benchmarks/generated/bench_estimators_04.flow index dca7f52..241ee96 100644 --- a/benchmarks/generated/bench_estimators_04.flow +++ b/benchmarks/generated/bench_estimators_04.flow @@ -82,6 +82,33 @@ function main() -> i32 { let mut reps: i32 = 1 let mut sink: f32 = 0.0 + # ---- lasso_lars_ic (regression) ---- + t0 = flow_now_ns() + let probe_lasso_lars_ic: LassoLarsIC = lasso_lars_ic_fit(X_r, y_r, 0) + t1 = flow_now_ns() + lasso_lars_ic_free(probe_lasso_lars_ic) + reps = 1 + if (t1 - t0) < 200000 { reps = 200 } + elif (t1 - t0) < 2000000 { reps = 20 } + elif (t1 - t0) < 20000000 { reps = 5 } + t0 = flow_now_ns() + for rep in 0 to reps { + let m_lasso_lars_ic: LassoLarsIC = lasso_lars_ic_fit(X_r, y_r, 0) + lasso_lars_ic_free(m_lasso_lars_ic) + } + t1 = flow_now_ns() + let fitted_lasso_lars_ic: LassoLarsIC = lasso_lars_ic_fit(X_r, y_r, 0) + t2 = flow_now_ns() + for rep2 in 0 to reps { + let o_lasso_lars_ic: ptr = lasso_lars_ic_predict(fitted_lasso_lars_ic, X_r) + sink = sink + o_lasso_lars_ic[0] + array_free_f32(o_lasso_lars_ic) + } + t3 = flow_now_ns() + lasso_lars_ic_free(fitted_lasso_lars_ic) + printf("ESTIMATOR|lasso_lars_ic|%.9f|%.9f|%d|ok\n", ms_between(t0, t1) / (reps as f32), ms_between(t2, t3) / (reps as f32), reps) + fflush(null) + # ---- lda (unsupervised) ---- t0 = flow_now_ns() let probe_lda: LatentDirichletAllocation = lda_fit(X_c, 3, 50, 1.0, 0.1, 42) @@ -533,35 +560,6 @@ function main() -> i32 { printf("ESTIMATOR|minmax_scaler|%.9f|%.9f|%d|ok\n", ms_between(t0, t1) / (reps as f32), ms_between(t2, t3) / (reps as f32), reps) fflush(null) - # ---- missing_indicator (unsupervised) ---- - t0 = flow_now_ns() - let probe_missing_indicator: MissingIndicator = missing_indicator_fit(X_c, 0.0, 0) - t1 = flow_now_ns() - missing_indicator_free(probe_missing_indicator) - reps = 1 - if (t1 - t0) < 200000 { reps = 200 } - elif (t1 - t0) < 2000000 { reps = 20 } - elif (t1 - t0) < 20000000 { reps = 5 } - t0 = flow_now_ns() - for rep in 0 to reps { - let m_missing_indicator: MissingIndicator = missing_indicator_fit(X_c, 0.0, 0) - missing_indicator_free(m_missing_indicator) - } - t1 = flow_now_ns() - let fitted_missing_indicator: MissingIndicator = missing_indicator_fit(X_c, 0.0, 0) - t2 = flow_now_ns() - for rep2 in 0 to reps { - let o_missing_indicator: Matrix = missing_indicator_transform(fitted_missing_indicator, X_c) - if o_missing_indicator.rows > 0 { - if o_missing_indicator.cols > 0 { sink = sink + o_missing_indicator.data[0] } - } - matrix_free(o_missing_indicator) - } - t3 = flow_now_ns() - missing_indicator_free(fitted_missing_indicator) - printf("ESTIMATOR|missing_indicator|%.9f|%.9f|%d|ok\n", ms_between(t0, t1) / (reps as f32), ms_between(t2, t3) / (reps as f32), reps) - fflush(null) - for i in 0 to n_c { array_free_f32(Y_label_rows[i]) } free(Y_label_rows as ptr) matrix_free(Y_labels) diff --git a/benchmarks/generated/bench_estimators_05.flow b/benchmarks/generated/bench_estimators_05.flow index ae3c9d1..7f440d6 100644 --- a/benchmarks/generated/bench_estimators_05.flow +++ b/benchmarks/generated/bench_estimators_05.flow @@ -82,6 +82,35 @@ function main() -> i32 { let mut reps: i32 = 1 let mut sink: f32 = 0.0 + # ---- missing_indicator (unsupervised) ---- + t0 = flow_now_ns() + let probe_missing_indicator: MissingIndicator = missing_indicator_fit(X_c, 0.0, 0) + t1 = flow_now_ns() + missing_indicator_free(probe_missing_indicator) + reps = 1 + if (t1 - t0) < 200000 { reps = 200 } + elif (t1 - t0) < 2000000 { reps = 20 } + elif (t1 - t0) < 20000000 { reps = 5 } + t0 = flow_now_ns() + for rep in 0 to reps { + let m_missing_indicator: MissingIndicator = missing_indicator_fit(X_c, 0.0, 0) + missing_indicator_free(m_missing_indicator) + } + t1 = flow_now_ns() + let fitted_missing_indicator: MissingIndicator = missing_indicator_fit(X_c, 0.0, 0) + t2 = flow_now_ns() + for rep2 in 0 to reps { + let o_missing_indicator: Matrix = missing_indicator_transform(fitted_missing_indicator, X_c) + if o_missing_indicator.rows > 0 { + if o_missing_indicator.cols > 0 { sink = sink + o_missing_indicator.data[0] } + } + matrix_free(o_missing_indicator) + } + t3 = flow_now_ns() + missing_indicator_free(fitted_missing_indicator) + printf("ESTIMATOR|missing_indicator|%.9f|%.9f|%d|ok\n", ms_between(t0, t1) / (reps as f32), ms_between(t2, t3) / (reps as f32), reps) + fflush(null) + # ---- mlp_classifier (classification) ---- let hidden_sizes_mlp_classifier: array = [8] t0 = flow_now_ns() @@ -611,33 +640,6 @@ function main() -> i32 { printf("ESTIMATOR|oas_estimator|%.9f|%.9f|%d|ok\n", ms_between(t0, t1) / (reps as f32), 0.0, reps) fflush(null) - # ---- omp_cv (regression) ---- - t0 = flow_now_ns() - let probe_omp_cv: OrthogonalMatchingPursuitCV = omp_cv_fit(X_r, y_r, 3) - t1 = flow_now_ns() - omp_cv_free(probe_omp_cv) - reps = 1 - if (t1 - t0) < 200000 { reps = 200 } - elif (t1 - t0) < 2000000 { reps = 20 } - elif (t1 - t0) < 20000000 { reps = 5 } - t0 = flow_now_ns() - for rep in 0 to reps { - let m_omp_cv: OrthogonalMatchingPursuitCV = omp_cv_fit(X_r, y_r, 3) - omp_cv_free(m_omp_cv) - } - t1 = flow_now_ns() - let fitted_omp_cv: OrthogonalMatchingPursuitCV = omp_cv_fit(X_r, y_r, 3) - t2 = flow_now_ns() - for rep2 in 0 to reps { - let o_omp_cv: ptr = omp_cv_predict(fitted_omp_cv, X_r) - sink = sink + o_omp_cv[0] - array_free_f32(o_omp_cv) - } - t3 = flow_now_ns() - omp_cv_free(fitted_omp_cv) - printf("ESTIMATOR|omp_cv|%.9f|%.9f|%d|ok\n", ms_between(t0, t1) / (reps as f32), ms_between(t2, t3) / (reps as f32), reps) - fflush(null) - for i in 0 to n_c { array_free_f32(Y_label_rows[i]) } free(Y_label_rows as ptr) matrix_free(Y_labels) diff --git a/benchmarks/generated/bench_estimators_06.flow b/benchmarks/generated/bench_estimators_06.flow index cf91890..2fada50 100644 --- a/benchmarks/generated/bench_estimators_06.flow +++ b/benchmarks/generated/bench_estimators_06.flow @@ -82,6 +82,33 @@ function main() -> i32 { let mut reps: i32 = 1 let mut sink: f32 = 0.0 + # ---- omp_cv (regression) ---- + t0 = flow_now_ns() + let probe_omp_cv: OrthogonalMatchingPursuitCV = omp_cv_fit(X_r, y_r, 3) + t1 = flow_now_ns() + omp_cv_free(probe_omp_cv) + reps = 1 + if (t1 - t0) < 200000 { reps = 200 } + elif (t1 - t0) < 2000000 { reps = 20 } + elif (t1 - t0) < 20000000 { reps = 5 } + t0 = flow_now_ns() + for rep in 0 to reps { + let m_omp_cv: OrthogonalMatchingPursuitCV = omp_cv_fit(X_r, y_r, 3) + omp_cv_free(m_omp_cv) + } + t1 = flow_now_ns() + let fitted_omp_cv: OrthogonalMatchingPursuitCV = omp_cv_fit(X_r, y_r, 3) + t2 = flow_now_ns() + for rep2 in 0 to reps { + let o_omp_cv: ptr = omp_cv_predict(fitted_omp_cv, X_r) + sink = sink + o_omp_cv[0] + array_free_f32(o_omp_cv) + } + t3 = flow_now_ns() + omp_cv_free(fitted_omp_cv) + printf("ESTIMATOR|omp_cv|%.9f|%.9f|%d|ok\n", ms_between(t0, t1) / (reps as f32), ms_between(t2, t3) / (reps as f32), reps) + fflush(null) + # ---- one_class_svm (unsupervised) ---- t0 = flow_now_ns() let probe_one_class_svm: OneClassSVM = one_class_svm_fit(X_c, 0.5, 0.1, 100) @@ -596,35 +623,6 @@ function main() -> i32 { printf("ESTIMATOR|polynomial_features|%.9f|%.9f|%d|ok\n", ms_between(t0, t1) / (reps as f32), ms_between(t2, t3) / (reps as f32), reps) fflush(null) - # ---- power_transformer (unsupervised) ---- - t0 = flow_now_ns() - let probe_power_transformer: PowerTransformer = power_transformer_fit(X_c, 0) - t1 = flow_now_ns() - power_transformer_free(probe_power_transformer) - reps = 1 - if (t1 - t0) < 200000 { reps = 200 } - elif (t1 - t0) < 2000000 { reps = 20 } - elif (t1 - t0) < 20000000 { reps = 5 } - t0 = flow_now_ns() - for rep in 0 to reps { - let m_power_transformer: PowerTransformer = power_transformer_fit(X_c, 0) - power_transformer_free(m_power_transformer) - } - t1 = flow_now_ns() - let fitted_power_transformer: PowerTransformer = power_transformer_fit(X_c, 0) - t2 = flow_now_ns() - for rep2 in 0 to reps { - let o_power_transformer: Matrix = power_transformer_transform(fitted_power_transformer, X_c) - if o_power_transformer.rows > 0 { - if o_power_transformer.cols > 0 { sink = sink + o_power_transformer.data[0] } - } - matrix_free(o_power_transformer) - } - t3 = flow_now_ns() - power_transformer_free(fitted_power_transformer) - printf("ESTIMATOR|power_transformer|%.9f|%.9f|%d|ok\n", ms_between(t0, t1) / (reps as f32), ms_between(t2, t3) / (reps as f32), reps) - fflush(null) - for i in 0 to n_c { array_free_f32(Y_label_rows[i]) } free(Y_label_rows as ptr) matrix_free(Y_labels) diff --git a/benchmarks/generated/bench_estimators_07.flow b/benchmarks/generated/bench_estimators_07.flow index d96f2eb..d27a4f6 100644 --- a/benchmarks/generated/bench_estimators_07.flow +++ b/benchmarks/generated/bench_estimators_07.flow @@ -82,6 +82,35 @@ function main() -> i32 { let mut reps: i32 = 1 let mut sink: f32 = 0.0 + # ---- power_transformer (unsupervised) ---- + t0 = flow_now_ns() + let probe_power_transformer: PowerTransformer = power_transformer_fit(X_c, 0) + t1 = flow_now_ns() + power_transformer_free(probe_power_transformer) + reps = 1 + if (t1 - t0) < 200000 { reps = 200 } + elif (t1 - t0) < 2000000 { reps = 20 } + elif (t1 - t0) < 20000000 { reps = 5 } + t0 = flow_now_ns() + for rep in 0 to reps { + let m_power_transformer: PowerTransformer = power_transformer_fit(X_c, 0) + power_transformer_free(m_power_transformer) + } + t1 = flow_now_ns() + let fitted_power_transformer: PowerTransformer = power_transformer_fit(X_c, 0) + t2 = flow_now_ns() + for rep2 in 0 to reps { + let o_power_transformer: Matrix = power_transformer_transform(fitted_power_transformer, X_c) + if o_power_transformer.rows > 0 { + if o_power_transformer.cols > 0 { sink = sink + o_power_transformer.data[0] } + } + matrix_free(o_power_transformer) + } + t3 = flow_now_ns() + power_transformer_free(fitted_power_transformer) + printf("ESTIMATOR|power_transformer|%.9f|%.9f|%d|ok\n", ms_between(t0, t1) / (reps as f32), ms_between(t2, t3) / (reps as f32), reps) + fflush(null) + # ---- qda (classification) ---- t0 = flow_now_ns() let probe_qda: QuadraticDiscriminantAnalysis = qda_fit(X_c, y_c, 3) @@ -589,35 +618,6 @@ function main() -> i32 { printf("ESTIMATOR|robust_scaler|%.9f|%.9f|%d|ok\n", ms_between(t0, t1) / (reps as f32), ms_between(t2, t3) / (reps as f32), reps) fflush(null) - # ---- select_fdr (regression) ---- - t0 = flow_now_ns() - let probe_select_fdr: SelectFdr = select_fdr_fit(X_r, y_r, 1.0) - t1 = flow_now_ns() - select_fdr_free(probe_select_fdr) - reps = 1 - if (t1 - t0) < 200000 { reps = 200 } - elif (t1 - t0) < 2000000 { reps = 20 } - elif (t1 - t0) < 20000000 { reps = 5 } - t0 = flow_now_ns() - for rep in 0 to reps { - let m_select_fdr: SelectFdr = select_fdr_fit(X_r, y_r, 1.0) - select_fdr_free(m_select_fdr) - } - t1 = flow_now_ns() - let fitted_select_fdr: SelectFdr = select_fdr_fit(X_r, y_r, 1.0) - t2 = flow_now_ns() - for rep2 in 0 to reps { - let o_select_fdr: Matrix = select_fdr_transform(fitted_select_fdr, X_r) - if o_select_fdr.rows > 0 { - if o_select_fdr.cols > 0 { sink = sink + o_select_fdr.data[0] } - } - matrix_free(o_select_fdr) - } - t3 = flow_now_ns() - select_fdr_free(fitted_select_fdr) - printf("ESTIMATOR|select_fdr|%.9f|%.9f|%d|ok\n", ms_between(t0, t1) / (reps as f32), ms_between(t2, t3) / (reps as f32), reps) - fflush(null) - for i in 0 to n_c { array_free_f32(Y_label_rows[i]) } free(Y_label_rows as ptr) matrix_free(Y_labels) diff --git a/benchmarks/generated/bench_estimators_08.flow b/benchmarks/generated/bench_estimators_08.flow index 7c4579c..7218be0 100644 --- a/benchmarks/generated/bench_estimators_08.flow +++ b/benchmarks/generated/bench_estimators_08.flow @@ -82,6 +82,35 @@ function main() -> i32 { let mut reps: i32 = 1 let mut sink: f32 = 0.0 + # ---- select_fdr (regression) ---- + t0 = flow_now_ns() + let probe_select_fdr: SelectFdr = select_fdr_fit(X_r, y_r, 1.0) + t1 = flow_now_ns() + select_fdr_free(probe_select_fdr) + reps = 1 + if (t1 - t0) < 200000 { reps = 200 } + elif (t1 - t0) < 2000000 { reps = 20 } + elif (t1 - t0) < 20000000 { reps = 5 } + t0 = flow_now_ns() + for rep in 0 to reps { + let m_select_fdr: SelectFdr = select_fdr_fit(X_r, y_r, 1.0) + select_fdr_free(m_select_fdr) + } + t1 = flow_now_ns() + let fitted_select_fdr: SelectFdr = select_fdr_fit(X_r, y_r, 1.0) + t2 = flow_now_ns() + for rep2 in 0 to reps { + let o_select_fdr: Matrix = select_fdr_transform(fitted_select_fdr, X_r) + if o_select_fdr.rows > 0 { + if o_select_fdr.cols > 0 { sink = sink + o_select_fdr.data[0] } + } + matrix_free(o_select_fdr) + } + t3 = flow_now_ns() + select_fdr_free(fitted_select_fdr) + printf("ESTIMATOR|select_fdr|%.9f|%.9f|%d|ok\n", ms_between(t0, t1) / (reps as f32), ms_between(t2, t3) / (reps as f32), reps) + fflush(null) + # ---- select_fpr (regression) ---- t0 = flow_now_ns() let probe_select_fpr: SelectFpr = select_fpr_fit(X_r, y_r, 1.0) @@ -227,6 +256,27 @@ function main() -> i32 { printf("ESTIMATOR|select_percentile|%.9f|%.9f|%d|ok\n", ms_between(t0, t1) / (reps as f32), ms_between(t2, t3) / (reps as f32), reps) fflush(null) + # ---- self_training_classifier_full (classification, written out) ---- + let st_labels: ptr = malloc((n_c as i64) * 4 + 64) as ptr + for i in 0 to n_c { + if i % 3 == 0 { st_labels[i] = yi_c[i] } + else { st_labels[i] = 0 - 1 } + } + t0 = flow_now_ns() + let probe_self_training_classifier_full: SelfTrainingClassifier = self_training_classifier_full_fit(X_c, st_labels, 3, 0.7, 10, 50, 0.1) + t1 = flow_now_ns() + reps = 1 + if (t1 - t0) < 200000 { reps = 200 } + elif (t1 - t0) < 2000000 { reps = 20 } + elif (t1 - t0) < 20000000 { reps = 5 } + t0 = flow_now_ns() + for rep in 0 to reps { + let m_self_training_classifier_full: SelfTrainingClassifier = self_training_classifier_full_fit(X_c, st_labels, 3, 0.7, 10, 50, 0.1) + } + t1 = flow_now_ns() + printf("ESTIMATOR|self_training_classifier_full|%.9f|%.9f|%d|ok\n", ms_between(t0, t1) / (reps as f32), 0.0, reps) + fflush(null) + # ---- sequential_feature_selector (regression) ---- t0 = flow_now_ns() let probe_sequential_feature_selector: SequentialFeatureSelector = sequential_feature_selector_fit(X_r, y_r, 2, 0, n_r, f_r) @@ -554,53 +604,6 @@ function main() -> i32 { printf("ESTIMATOR|spectral_coclustering|%.9f|%.9f|%d|ok\n", ms_between(t0, t1) / (reps as f32), 0.0, reps) fflush(null) - # ---- spectral_embedding (unsupervised) ---- - t0 = flow_now_ns() - let probe_spectral_embedding: SpectralEmbedding = spectral_embedding_fit(X_c, 2, 0.1) - t1 = flow_now_ns() - spectral_embedding_free(probe_spectral_embedding) - reps = 1 - if (t1 - t0) < 200000 { reps = 200 } - elif (t1 - t0) < 2000000 { reps = 20 } - elif (t1 - t0) < 20000000 { reps = 5 } - t0 = flow_now_ns() - for rep in 0 to reps { - let m_spectral_embedding: SpectralEmbedding = spectral_embedding_fit(X_c, 2, 0.1) - spectral_embedding_free(m_spectral_embedding) - } - t1 = flow_now_ns() - printf("ESTIMATOR|spectral_embedding|%.9f|%.9f|%d|ok\n", ms_between(t0, t1) / (reps as f32), 0.0, reps) - fflush(null) - - # ---- spline_transformer (unsupervised) ---- - t0 = flow_now_ns() - let probe_spline_transformer: SplineTransformer = spline_transformer_fit(X_c, 4, 2) - t1 = flow_now_ns() - spline_transformer_free(probe_spline_transformer) - reps = 1 - if (t1 - t0) < 200000 { reps = 200 } - elif (t1 - t0) < 2000000 { reps = 20 } - elif (t1 - t0) < 20000000 { reps = 5 } - t0 = flow_now_ns() - for rep in 0 to reps { - let m_spline_transformer: SplineTransformer = spline_transformer_fit(X_c, 4, 2) - spline_transformer_free(m_spline_transformer) - } - t1 = flow_now_ns() - let fitted_spline_transformer: SplineTransformer = spline_transformer_fit(X_c, 4, 2) - t2 = flow_now_ns() - for rep2 in 0 to reps { - let o_spline_transformer: Matrix = spline_transformer_transform(fitted_spline_transformer, X_c) - if o_spline_transformer.rows > 0 { - if o_spline_transformer.cols > 0 { sink = sink + o_spline_transformer.data[0] } - } - matrix_free(o_spline_transformer) - } - t3 = flow_now_ns() - spline_transformer_free(fitted_spline_transformer) - printf("ESTIMATOR|spline_transformer|%.9f|%.9f|%d|ok\n", ms_between(t0, t1) / (reps as f32), ms_between(t2, t3) / (reps as f32), reps) - fflush(null) - for i in 0 to n_c { array_free_f32(Y_label_rows[i]) } free(Y_label_rows as ptr) matrix_free(Y_labels) diff --git a/benchmarks/generated/bench_estimators_09.flow b/benchmarks/generated/bench_estimators_09.flow index fe79af1..3402a20 100644 --- a/benchmarks/generated/bench_estimators_09.flow +++ b/benchmarks/generated/bench_estimators_09.flow @@ -82,6 +82,80 @@ function main() -> i32 { let mut reps: i32 = 1 let mut sink: f32 = 0.0 + # ---- spectral_embedding (unsupervised) ---- + t0 = flow_now_ns() + let probe_spectral_embedding: SpectralEmbedding = spectral_embedding_fit(X_c, 2, 0.1) + t1 = flow_now_ns() + spectral_embedding_free(probe_spectral_embedding) + reps = 1 + if (t1 - t0) < 200000 { reps = 200 } + elif (t1 - t0) < 2000000 { reps = 20 } + elif (t1 - t0) < 20000000 { reps = 5 } + t0 = flow_now_ns() + for rep in 0 to reps { + let m_spectral_embedding: SpectralEmbedding = spectral_embedding_fit(X_c, 2, 0.1) + spectral_embedding_free(m_spectral_embedding) + } + t1 = flow_now_ns() + printf("ESTIMATOR|spectral_embedding|%.9f|%.9f|%d|ok\n", ms_between(t0, t1) / (reps as f32), 0.0, reps) + fflush(null) + + # ---- spline_transformer (unsupervised) ---- + t0 = flow_now_ns() + let probe_spline_transformer: SplineTransformer = spline_transformer_fit(X_c, 4, 2) + t1 = flow_now_ns() + spline_transformer_free(probe_spline_transformer) + reps = 1 + if (t1 - t0) < 200000 { reps = 200 } + elif (t1 - t0) < 2000000 { reps = 20 } + elif (t1 - t0) < 20000000 { reps = 5 } + t0 = flow_now_ns() + for rep in 0 to reps { + let m_spline_transformer: SplineTransformer = spline_transformer_fit(X_c, 4, 2) + spline_transformer_free(m_spline_transformer) + } + t1 = flow_now_ns() + let fitted_spline_transformer: SplineTransformer = spline_transformer_fit(X_c, 4, 2) + t2 = flow_now_ns() + for rep2 in 0 to reps { + let o_spline_transformer: Matrix = spline_transformer_transform(fitted_spline_transformer, X_c) + if o_spline_transformer.rows > 0 { + if o_spline_transformer.cols > 0 { sink = sink + o_spline_transformer.data[0] } + } + matrix_free(o_spline_transformer) + } + t3 = flow_now_ns() + spline_transformer_free(fitted_spline_transformer) + printf("ESTIMATOR|spline_transformer|%.9f|%.9f|%d|ok\n", ms_between(t0, t1) / (reps as f32), ms_between(t2, t3) / (reps as f32), reps) + fflush(null) + + # ---- stacking_classifier_full (classification, written out) ---- + t0 = flow_now_ns() + let probe_stacking_classifier_full: StackingClassifierFull = stacking_classifier_full_fit(X_c, y_c, 3, 3, 3, 42, 0.01, 50) + t1 = flow_now_ns() + stacking_classifier_full_free(probe_stacking_classifier_full) + reps = 1 + if (t1 - t0) < 200000 { reps = 200 } + elif (t1 - t0) < 2000000 { reps = 20 } + elif (t1 - t0) < 20000000 { reps = 5 } + t0 = flow_now_ns() + for rep in 0 to reps { + let m_stacking_classifier_full: StackingClassifierFull = stacking_classifier_full_fit(X_c, y_c, 3, 3, 3, 42, 0.01, 50) + stacking_classifier_full_free(m_stacking_classifier_full) + } + t1 = flow_now_ns() + let fitted_stacking_classifier_full: StackingClassifierFull = stacking_classifier_full_fit(X_c, y_c, 3, 3, 3, 42, 0.01, 50) + t2 = flow_now_ns() + for rep2 in 0 to reps { + let o_stacking_classifier_full: ptr = stacking_classifier_full_predict(fitted_stacking_classifier_full, X_c) + sink = sink + o_stacking_classifier_full[0] + array_free_f32(o_stacking_classifier_full) + } + t3 = flow_now_ns() + stacking_classifier_full_free(fitted_stacking_classifier_full) + printf("ESTIMATOR|stacking_classifier_full|%.9f|%.9f|%d|ok\n", ms_between(t0, t1) / (reps as f32), ms_between(t2, t3) / (reps as f32), reps) + fflush(null) + # ---- stacking_regressor (regression) ---- t0 = flow_now_ns() let probe_stacking_regressor: StackingRegressor = stacking_regressor_fit(X_r, y_r, 3, 5, 42) @@ -425,6 +499,60 @@ function main() -> i32 { printf("ESTIMATOR|variance_threshold|%.9f|%.9f|%d|ok\n", ms_between(t0, t1) / (reps as f32), ms_between(t2, t3) / (reps as f32), reps) fflush(null) + # ---- voting_classifier_full (classification, written out) ---- + t0 = flow_now_ns() + let probe_voting_classifier_full: VotingClassifier = voting_classifier_full_fit(X_c, y_c, 3, 3, 3, 42, 0) + t1 = flow_now_ns() + voting_classifier_free(probe_voting_classifier_full) + reps = 1 + if (t1 - t0) < 200000 { reps = 200 } + elif (t1 - t0) < 2000000 { reps = 20 } + elif (t1 - t0) < 20000000 { reps = 5 } + t0 = flow_now_ns() + for rep in 0 to reps { + let m_voting_classifier_full: VotingClassifier = voting_classifier_full_fit(X_c, y_c, 3, 3, 3, 42, 0) + voting_classifier_free(m_voting_classifier_full) + } + t1 = flow_now_ns() + let fitted_voting_classifier_full: VotingClassifier = voting_classifier_full_fit(X_c, y_c, 3, 3, 3, 42, 0) + t2 = flow_now_ns() + for rep2 in 0 to reps { + let o_voting_classifier_full: ptr = voting_classifier_predict(fitted_voting_classifier_full, X_c) + sink = sink + o_voting_classifier_full[0] + array_free_f32(o_voting_classifier_full) + } + t3 = flow_now_ns() + voting_classifier_free(fitted_voting_classifier_full) + printf("ESTIMATOR|voting_classifier_full|%.9f|%.9f|%d|ok\n", ms_between(t0, t1) / (reps as f32), ms_between(t2, t3) / (reps as f32), reps) + fflush(null) + + # ---- voting_regressor_full (regression, written out) ---- + t0 = flow_now_ns() + let probe_voting_regressor_full: VotingRegressor = voting_regressor_full_fit(X_r, y_r, 3, 3, 42) + t1 = flow_now_ns() + voting_regressor_free(probe_voting_regressor_full) + reps = 1 + if (t1 - t0) < 200000 { reps = 200 } + elif (t1 - t0) < 2000000 { reps = 20 } + elif (t1 - t0) < 20000000 { reps = 5 } + t0 = flow_now_ns() + for rep in 0 to reps { + let m_voting_regressor_full: VotingRegressor = voting_regressor_full_fit(X_r, y_r, 3, 3, 42) + voting_regressor_free(m_voting_regressor_full) + } + t1 = flow_now_ns() + let fitted_voting_regressor_full: VotingRegressor = voting_regressor_full_fit(X_r, y_r, 3, 3, 42) + t2 = flow_now_ns() + for rep2 in 0 to reps { + let o_voting_regressor_full: ptr = voting_regressor_predict(fitted_voting_regressor_full, X_r) + sink = sink + o_voting_regressor_full[0] + array_free_f32(o_voting_regressor_full) + } + t3 = flow_now_ns() + voting_regressor_free(fitted_voting_regressor_full) + printf("ESTIMATOR|voting_regressor_full|%.9f|%.9f|%d|ok\n", ms_between(t0, t1) / (reps as f32), ms_between(t2, t3) / (reps as f32), reps) + fflush(null) + for i in 0 to n_c { array_free_f32(Y_label_rows[i]) } free(Y_label_rows as ptr) matrix_free(Y_labels) diff --git a/lib/scikit/decomposition.flow b/lib/scikit/decomposition.flow index 45a0bc9..51e17c9 100644 --- a/lib/scikit/decomposition.flow +++ b/lib/scikit/decomposition.flow @@ -1518,72 +1518,207 @@ export struct IncrementalPCA { components: ptr >, mean: ptr, variances: ptr, + # The scatter matrix of everything seen so far, which is what makes a + # partial fit exact: two batches combine into the scatter of their union + # once the shift between their means is accounted for. + scatter: ptr, n_components: i32, n_features: i32, n_samples_seen: i32, fitted: bool } +# The top eigenvectors of a symmetric p by p matrix, by block power iteration +# with modified Gram-Schmidt. The name carries a prefix because a non-exported +# function with the same name in two modules collides in the generated C. +# Flow issue #465. +function _ipca_top_eigenvectors(M: ptr, p: i32, k: i32, V: ptr, eigenvalues: ptr) -> void { + let Z: ptr = array_new_f32(p * k) + srand(42) + for i in 0 to p * k { V[i] = ((rand() % 1000) as f32) / 500.0 - 1.0 } + + let mut iter: i32 = 0 + while iter < 200 { + blas_matmat(M, V, Z, p, k, p, 1.0, 0.0) + # Orthonormalize the block, which keeps the columns from all sliding + # onto the leading eigenvector. + for c in 0 to k { + for prev in 0 to c { + let mut dot: f32 = 0.0 + let mut i2: i32 = 0 + while i2 < p { + dot = dot + Z[i2 * k + c] * Z[i2 * k + prev] + i2 = i2 + 1 + } + let mut i3: i32 = 0 + while i3 < p { + Z[i3 * k + c] = Z[i3 * k + c] - dot * Z[i3 * k + prev] + i3 = i3 + 1 + } + } + let mut norm: f32 = 0.0 + let mut i4: i32 = 0 + while i4 < p { + let v: f32 = Z[i4 * k + c] + norm = norm + v * v + i4 = i4 + 1 + } + norm = sqrt((norm) as f64) as f32 + if norm > 0.0000000001 { + let inv: f32 = 1.0 / norm + let mut i5: i32 = 0 + while i5 < p { + Z[i5 * k + c] = Z[i5 * k + c] * inv + i5 = i5 + 1 + } + } + } + + let mut max_change: f32 = 0.0 + let mut i6: i32 = 0 + let total: i32 = p * k + while i6 < total { + let mut d: f32 = Z[i6] - V[i6] + if d < 0.0 { d = 0.0 - d } + if d > max_change { max_change = d } + V[i6] = Z[i6] + i6 = i6 + 1 + } + if max_change < 0.0000001 { break } + iter = iter + 1 + } + + # The Rayleigh quotient of each column, which is its eigenvalue. + let w: ptr = array_new_f32(p) + for c2 in 0 to k { + let mut i7: i32 = 0 + while i7 < p { + let mrow: ptr = M + i7 * p + let mut s: f32 = 0.0 + let mut j: i32 = 0 + while j < p { + s = s + mrow[j] * V[j * k + c2] + j = j + 1 + } + w[i7] = s + i7 = i7 + 1 + } + let mut lam: f32 = 0.0 + let mut r: i32 = 0 + while r < p { + lam = lam + V[r * k + c2] * w[r] + r = r + 1 + } + eigenvalues[c2] = lam + } + array_free_f32(w) + array_free_f32(Z) +} + +# One batch, folded into what the model has already seen. +# +# What stood here accumulated the same per-feature quantity for every +# component, so every component came out as the same vector, and its power +# iteration multiplied a component by itself rather than by a covariance. The +# estimator was never raced and nothing asserted on its output, so neither +# showed up. export function incremental_pca_partial_fit(model: IncrementalPCA, X: Matrix) -> IncrementalPCA { let n: i32 = X.rows let p: i32 = X.cols + let xdata: ptr = X.data let batch_mean: ptr = array_new_f32(p) - for j in 0 to p { - let mut s: f32 = 0.0 - for i in 0 to n { s = s + matrix_at(X, i, j) } - batch_mean[j] = s / (n as f32) - } - - let total_n: i32 = model.n_samples_seen + n - let mut new_mean: ptr = array_new_f32(p) - for j in 0 to p { - if model.n_samples_seen == 0 { - new_mean[j] = batch_mean[j] - } else { - new_mean[j] = (model.n_samples_seen as f32 * model.mean[j] + n as f32 * batch_mean[j]) / (total_n as f32) + for i in 0 to n { + let row: ptr = xdata + i * p + let mut j: i32 = 0 + while j < p { + batch_mean[j] = batch_mean[j] + row[j] + j = j + 1 } } + let inv_n: f32 = 1.0 / (n as f32) + for j2 in 0 to p { batch_mean[j2] = batch_mean[j2] * inv_n } - for c in 0 to model.n_components { - let Av: ptr = array_new_f32(p) - for i in 0 to p { - let mut s: f32 = 0.0 - for j in 0 to p { - let cov_ij: f32 = 0.0 - if model.n_samples_seen > 0 { - cov_ij = model.components[c][j] - } - s = s + cov_ij * model.components[c][j] + # The batch's own scatter, about its own mean. + let batch_scatter: ptr = array_new_f32(p * p) + let diff: ptr = array_new_f32(p) + for i2 in 0 to n { + let row2: ptr = xdata + i2 * p + let mut j3: i32 = 0 + while j3 < p { + diff[j3] = row2[j3] - batch_mean[j3] + j3 = j3 + 1 + } + let mut a: i32 = 0 + while a < p { + let da: f32 = diff[a] + let srow: ptr = batch_scatter + a * p + let mut b: i32 = 0 + while b < p { + srow[b] = srow[b] + da * diff[b] + b = b + 1 } - Av[i] = s + a = a + 1 } + } - for i in 0 to n { - for j in 0 to p { - let d: f32 = matrix_at(X, i, j) - new_mean[j] - Av[j] = Av[j] + d * (matrix_at(X, i, j) - batch_mean[j]) + let seen: i32 = model.n_samples_seen + let total_n: i32 = seen + n + if seen == 0 { + for i3 in 0 to p * p { model.scatter[i3] = batch_scatter[i3] } + for j4 in 0 to p { model.mean[j4] = batch_mean[j4] } + } else { + # The correction for the two means being in different places, which + # is what makes the sum exact rather than approximate. + let w: f32 = ((seen as f32) * (n as f32)) / (total_n as f32) + for a2 in 0 to p { + let da2: f32 = batch_mean[a2] - model.mean[a2] + let srow2: ptr = model.scatter + a2 * p + let brow: ptr = batch_scatter + a2 * p + let mut b2: i32 = 0 + while b2 < p { + srow2[b2] = srow2[b2] + brow[b2] + w * da2 * (batch_mean[b2] - model.mean[b2]) + b2 = b2 + 1 } } - - let mut norm: f32 = 0.0 - for j in 0 to p { norm = norm + Av[j] * Av[j] } - norm = sqrt((norm as f64)) as f32 - if norm < 0.0000000001 { norm = 0.0000000001 } - - for j in 0 to p { - model.components[c][j] = Av[j] / norm + let inv_total: f32 = 1.0 / (total_n as f32) + for j5 in 0 to p { + model.mean[j5] = ((seen as f32) * model.mean[j5] + (n as f32) * batch_mean[j5]) * inv_total } - array_free_f32(Av) } + model.n_samples_seen = total_n - for j in 0 to p { - model.mean[j] = new_mean[j] + # The components are the leading eigenvectors of the covariance the + # scatter carries. + let k: i32 = model.n_components + let V: ptr = array_new_f32(p * k) + let eigenvalues: ptr = array_new_f32(k) + let cov: ptr = array_new_f32(p * p) + let denom: f32 = 1.0 / ((total_n - 1) as f32) + if total_n > 1 { + for i4 in 0 to p * p { cov[i4] = model.scatter[i4] * denom } + } else { + for i5 in 0 to p * p { cov[i5] = model.scatter[i5] } + } + _ipca_top_eigenvectors(cov, p, k, V, eigenvalues) + for c in 0 to k { + let comp: ptr = model.components[c] + let mut j6: i32 = 0 + while j6 < p { + comp[j6] = V[j6 * k + c] + j6 = j6 + 1 + } + } + for j7 in 0 to p { + if j7 < k { model.variances[j7] = eigenvalues[j7] } } - model.n_samples_seen = total_n + array_free_f32(V) + array_free_f32(eigenvalues) + array_free_f32(cov) + array_free_f32(batch_scatter) array_free_f32(batch_mean) - array_free_f32(new_mean) + array_free_f32(diff) return model } @@ -1605,6 +1740,7 @@ export function incremental_pca_init(n_components: i32, n_features: i32) -> Incr components: components, mean: mean, variances: variances, + scatter: array_new_f32(n_features * n_features), n_components: n_components, n_features: n_features, n_samples_seen: 0, @@ -1612,6 +1748,49 @@ export function incremental_pca_init(n_components: i32, n_features: i32) -> Incr } } +# The whole design, walked in batches, which is what IncrementalPCA.fit does +# with the partial fit its own users call one batch at a time. A single +# partial fit is one step of this, and racing one step against a full fit +# compares two different amounts of work. +# +# scikit-learn sizes its default batch at five times the feature count, so +# this does the same when it is handed a batch size of zero or less. +export function incremental_pca_fit(X: Matrix, n_components: i32, batch_size: i32) -> IncrementalPCA { + let n: i32 = X.rows + let p: i32 = X.cols + let mut step: i32 = batch_size + if step <= 0 { step = 5 * p } + if step < n_components { step = n_components } + if step < 1 { step = 1 } + + let mut model: IncrementalPCA = incremental_pca_init(n_components, p) + let xdata: ptr = X.data + let mut start: i32 = 0 + while start < n { + let mut count: i32 = step + if start + count > n { count = n - start } + # A batch smaller than the component count cannot add one, so it joins + # the batch before it rather than being dropped. + if count < n_components { + start = n + } else { + let batch: Matrix = matrix_new(count, p) + let dst: ptr = batch.data + let src: ptr = xdata + start * p + let mut i: i32 = 0 + let total: i32 = count * p + while i < total { + dst[i] = src[i] + i = i + 1 + } + model = incremental_pca_partial_fit(model, batch) + matrix_free(batch) + start = start + count + } + } + return model +} + export function incremental_pca_transform(model: IncrementalPCA, X: Matrix) -> Matrix { let result: Matrix = matrix_new(X.rows, model.n_components) for i in 0 to X.rows { @@ -1633,6 +1812,7 @@ export function incremental_pca_free(model: IncrementalPCA) -> void { free(model.components as ptr) array_free_f32(model.mean) array_free_f32(model.variances) + array_free_f32(model.scatter) } # ============================================================================ diff --git a/lib/scikit/ensemble.flow b/lib/scikit/ensemble.flow index fb6aa6e..5aa9e28 100644 --- a/lib/scikit/ensemble.flow +++ b/lib/scikit/ensemble.flow @@ -697,6 +697,37 @@ export struct VotingClassifier { fitted: bool } +# The whole ensemble, fitted from the data, which is what scikit-learn's +# VotingClassifier does when it is handed estimators: it clones and fits each +# one, then votes. The function below takes trees that are already fitted, so +# its cost is the vote alone and racing it against the class compares a part +# against the whole. +# +# The trees differ by the seed they draw their feature subsets from, which is +# what makes a vote of identical learners worth taking. +export function voting_classifier_full_fit(X: Matrix, y: ptr, n_classes: i32, n_estimators: i32, max_depth: i32, seed: i32, voting: i32) -> VotingClassifier { + let trees: ptr = malloc((n_estimators as i64) * 256) as ptr + let state: ptr = malloc(8) as ptr + let n_features: i32 = X.cols + for e in 0 to n_estimators { + state[0] = (seed + e * 7919) as u32 + trees[e] = decision_tree_classifier_fit_rf(X, y, n_classes, max_depth, 0, n_features, state) + } + free(state as ptr) + + let classes: ptr = array_new_f32(n_classes) + let found: DistinctF32 = distinct_f32(y, X.rows, n_classes) + for c in 0 to n_classes { + if c < found.count { classes[c] = found.values[c] } + else { classes[c] = c as f32 } + } + array_free_f32(found.values) + + let model: VotingClassifier = voting_classifier_fit(trees, n_estimators, n_classes, classes, voting) + array_free_f32(classes) + return model +} + export function voting_classifier_fit(estimators: ptr, n_estimators: i32, n_classes: i32, classes: ptr, voting: i32) -> VotingClassifier { return VotingClassifier { # The model takes over the fitted sub-estimators and its free @@ -1809,6 +1840,190 @@ export struct StackingClassifier { fitted: bool } +export struct StackingClassifierFull { + trees: ptr, + n_base: i32, + n_classes: i32, + meta_weights: ptr >, + meta_bias: ptr, + fitted: bool +} + +# Stacking with the base estimators fitted and cross-validated here, which is +# what scikit-learn's StackingClassifier does. +# +# The function below takes the base estimators' predictions as input, so its +# cost is the meta-learner alone. The class it is named after fits every base +# estimator once per fold to build the meta-features, fits the meta-learner on +# those, and then refits the base estimators on all of the data. That is +# almost all of the work, and leaving it out compares a part against the whole. +export function stacking_classifier_full_fit(X: Matrix, y: ptr, n_classes: i32, n_base: i32, max_depth: i32, seed: i32, lr: f32, n_iter: i32) -> StackingClassifierFull { + let n: i32 = X.rows + let p: i32 = X.cols + let xdata: ptr = X.data + + # Out of fold predictions, which is what the meta-learner has to be + # trained on: a base estimator's prediction for a row it helped fit would + # tell the meta-learner more than it will know at prediction time. + let meta: ptr = array_new_f32(n * n_base) + let n_folds: i32 = 5 + let fold_size: i32 = n / n_folds + let state: ptr = malloc(8) as ptr + + for fold in 0 to n_folds { + let test_start: i32 = fold * fold_size + let mut test_end: i32 = test_start + fold_size + if fold == n_folds - 1 { test_end = n } + let n_test: i32 = test_end - test_start + let n_train: i32 = n - n_test + if n_train < 2 { continue } + + let X_tr: Matrix = matrix_new(n_train, p) + let y_tr: ptr = array_new_f32(n_train) + let dst: ptr = X_tr.data + let mut at: i32 = 0 + for i in 0 to n { + if i < test_start || i >= test_end { + let src: ptr = xdata + i * p + let drow: ptr = dst + at * p + let mut j: i32 = 0 + while j < p { + drow[j] = src[j] + j = j + 1 + } + y_tr[at] = y[i] + at = at + 1 + } + } + + let X_te: Matrix = matrix_new(n_test, p) + let tdst: ptr = X_te.data + for t in 0 to n_test { + let src2: ptr = xdata + (test_start + t) * p + let trow: ptr = tdst + t * p + let mut j2: i32 = 0 + while j2 < p { + trow[j2] = src2[j2] + j2 = j2 + 1 + } + } + + for b in 0 to n_base { + state[0] = (seed + b * 7919 + fold * 104729) as u32 + let tree: DecisionTreeClassifier = decision_tree_classifier_fit_rf(X_tr, y_tr, n_classes, max_depth, 0, p, state) + let preds: ptr = decision_tree_classifier_predict(tree, X_te) + for t2 in 0 to n_test { + meta[(test_start + t2) * n_base + b] = preds[t2] + } + array_free_f32(preds) + decision_tree_classifier_free(tree) + } + + array_free_f32(y_tr) + matrix_free(X_tr) + matrix_free(X_te) + } + + # The meta-learner, one sigmoid per class over the base predictions. + let meta_weights: ptr > = malloc((n_classes as i64) * 8) as ptr > + let meta_bias: ptr = array_new_f32(n_classes) + for c in 0 to n_classes { + meta_weights[c] = array_new_f32(n_base) + meta_bias[c] = 0.0 + } + for iter in 0 to n_iter { + for i2 in 0 to n { + let row: ptr = meta + i2 * n_base + let yi: i32 = y[i2] as i32 + for c2 in 0 to n_classes { + let wc: ptr = meta_weights[c2] + let mut z: f32 = meta_bias[c2] + let mut j3: i32 = 0 + while j3 < n_base { + z = z + wc[j3] * row[j3] + j3 = j3 + 1 + } + let mut target: f32 = 0.0 + if yi == c2 { target = 1.0 } + let pred: f32 = 1.0 / (1.0 + (exp((0.0 - z) as f64) as f32)) + let err: f32 = (target - pred) * lr + let mut j4: i32 = 0 + while j4 < n_base { + wc[j4] = wc[j4] + err * row[j4] + j4 = j4 + 1 + } + meta_bias[c2] = meta_bias[c2] + err + } + } + } + + # The base estimators are refitted on all of the data, which is what + # predicts for a row nobody has seen. + let trees: ptr = malloc((n_base as i64) * 256) as ptr + for b2 in 0 to n_base { + state[0] = (seed + b2 * 7919) as u32 + trees[b2] = decision_tree_classifier_fit_rf(X, y, n_classes, max_depth, 0, p, state) + } + + free(state as ptr) + array_free_f32(meta) + + return StackingClassifierFull { + trees: trees, + n_base: n_base, + n_classes: n_classes, + meta_weights: meta_weights, + meta_bias: meta_bias, + fitted: true + } +} + +export function stacking_classifier_full_predict(model: StackingClassifierFull, X: Matrix) -> ptr { + let n: i32 = X.rows + let n_base: i32 = model.n_base + let meta: ptr = array_new_f32(n * n_base) + for b in 0 to n_base { + let preds: ptr = decision_tree_classifier_predict(model.trees[b], X) + for i in 0 to n { meta[i * n_base + b] = preds[i] } + array_free_f32(preds) + } + + let result: ptr = array_new_f32(n) + for i2 in 0 to n { + let row: ptr = meta + i2 * n_base + let mut best: i32 = 0 + let mut best_z: f32 = -1000000.0 + for c in 0 to model.n_classes { + let wc: ptr = model.meta_weights[c] + let mut z: f32 = model.meta_bias[c] + let mut j: i32 = 0 + while j < n_base { + z = z + wc[j] * row[j] + j = j + 1 + } + if z > best_z { + best_z = z + best = c + } + } + result[i2] = best as f32 + } + array_free_f32(meta) + return result +} + +export function stacking_classifier_full_free(model: StackingClassifierFull) -> void { + for b in 0 to model.n_base { + decision_tree_classifier_free(model.trees[b]) + } + free(model.trees as ptr) + for c in 0 to model.n_classes { + array_free_f32(model.meta_weights[c]) + } + free(model.meta_weights as ptr) + array_free_f32(model.meta_bias) +} + export function stacking_classifier_fit( base_preds: ptr >, y: ptr, @@ -2700,6 +2915,19 @@ export struct VotingRegressor { fitted: bool } +# The regression counterpart, for the reason the classifier above gives. +export function voting_regressor_full_fit(X: Matrix, y: ptr, n_estimators: i32, max_depth: i32, seed: i32) -> VotingRegressor { + let trees: ptr = malloc((n_estimators as i64) * 256) as ptr + for e in 0 to n_estimators { + trees[e] = decision_tree_regressor_fit(X, y, max_depth, 0) + } + let weights: ptr = array_new_f32(n_estimators) + for e2 in 0 to n_estimators { weights[e2] = 1.0 } + let model: VotingRegressor = voting_regressor_fit(trees, n_estimators, weights) + array_free_f32(weights) + return model +} + export function voting_regressor_fit(estimators: ptr, n_estimators: i32, weights: ptr) -> VotingRegressor { let result: VotingRegressor = VotingRegressor { # The model takes over the fitted sub-estimators and its free diff --git a/lib/scikit/semi_supervised.flow b/lib/scikit/semi_supervised.flow index 05f618c..150eb3b 100644 --- a/lib/scikit/semi_supervised.flow +++ b/lib/scikit/semi_supervised.flow @@ -3,6 +3,7 @@ import "lib/scikit/matrix.flow" import "lib/scikit/blas.flow" +import "lib/scikit/linear.flow" extern { function malloc(size: i64) -> ptr @@ -389,6 +390,118 @@ export struct SelfTrainingClassifier { fitted: bool } +# Self training with the base estimator refitted every round, which is what +# scikit-learn's SelfTrainingClassifier does. +# +# The function below takes class probabilities as input, so its cost is the +# labelling loop alone while the class it is named after also fits its base +# estimator on every round. That is the expensive half, and leaving it out +# compares a part against the whole. +# +# Unlabelled rows carry a negative label, as they do in scikit-learn. +export function self_training_classifier_full_fit(X: Matrix, y_initial: ptr, n_classes: i32, threshold: f32, max_iter: i32, base_iter: i32, lr: f32) -> SelfTrainingClassifier { + let n: i32 = X.rows + let p: i32 = X.cols + let xdata: ptr = X.data + + let labels: ptr = malloc((n as i64) * 4 + 64) as ptr + let confidence: ptr = array_new_f32(n) + for i in 0 to n { + labels[i] = y_initial[i] + confidence[i] = 0.0 + } + + let scores: ptr = array_new_f32(n * n_classes) + let y_lab: ptr = array_new_f32(n) + let mut n_iter: i32 = 0 + + for iter in 0 to max_iter { + n_iter = iter + 1 + + let mut n_lab: i32 = 0 + for i2 in 0 to n { + if labels[i2] >= 0 { n_lab = n_lab + 1 } + } + if n_lab < 2 { break } + + # The labelled rows, gathered for this round's fit. + let X_lab: Matrix = matrix_new(n_lab, p) + let dst: ptr = X_lab.data + let mut at: i32 = 0 + for i3 in 0 to n { + if labels[i3] >= 0 { + let src: ptr = xdata + i3 * p + let drow: ptr = dst + at * p + let mut j: i32 = 0 + while j < p { + drow[j] = src[j] + j = j + 1 + } + y_lab[at] = labels[i3] as f32 + at = at + 1 + } + } + + let base: LogisticRegression = logistic_regression_fit(X_lab, y_lab, n_classes, base_iter, lr, penalty_none()) + + # Per-class scores for every row, then a softmax, which is the + # confidence the threshold is compared against. + let nc: i32 = base.n_classes + blas_matmat_nt(xdata, base.weights, scores, n, nc, p, p, p, 1.0, 0.0) + let mut newly: i32 = 0 + for i4 in 0 to n { + if labels[i4] < 0 { + let srow: ptr = scores + i4 * nc + let mut best: i32 = 0 + let mut best_s: f32 = srow[0] + base.biases[0] + let mut c: i32 = 1 + while c < nc { + let s: f32 = srow[c] + base.biases[c] + if s > best_s { + best_s = s + best = c + } + c = c + 1 + } + let mut total: f32 = 0.0 + let mut c2: i32 = 0 + while c2 < nc { + total = total + (exp((srow[c2] + base.biases[c2] - best_s) as f64) as f32) + c2 = c2 + 1 + } + let prob: f32 = 1.0 / total + if prob >= threshold { + labels[i4] = best + confidence[i4] = prob + newly = newly + 1 + } + } + } + + logistic_regression_free(base) + matrix_free(X_lab) + if newly == 0 { break } + } + + let mut n_labeled: i32 = 0 + for i5 in 0 to n { + if labels[i5] >= 0 { n_labeled = n_labeled + 1 } + } + + array_free_f32(scores) + array_free_f32(y_lab) + + return SelfTrainingClassifier { + labels: labels, + confidence: confidence, + n_iter: n_iter, + n_labeled: n_labeled, + n_samples: n, + threshold: threshold, + fitted: true + } +} + export function self_training_classifier_fit( y_initial: ptr, probabilities: ptr >, diff --git a/tests/test_full_ensembles.flow b/tests/test_full_ensembles.flow new file mode 100644 index 0000000..e3ee76c --- /dev/null +++ b/tests/test_full_ensembles.flow @@ -0,0 +1,232 @@ +# The five estimators whose Flow functions used to cover part of what their +# scikit-learn namesake does. Each has a fit here that does the whole job, and +# each is checked for a result rather than for a shape. + +import "lib/scikit/scikit.flow" + +extern { + function printf(fmt: string, ...) -> i32 +} + +# Two well separated classes in two dimensions, with a third noise feature. +function fe_classes(n: i32) -> Matrix { + let X: Matrix = matrix_new(n, 3) + srand(5) + for i in 0 to n { + let cls: i32 = i % 2 + let centre: f32 = (cls as f32) * 6.0 + matrix_set(X, i, 0, centre + ((rand() % 100) as f32) / 100.0) + matrix_set(X, i, 1, centre + ((rand() % 100) as f32) / 100.0) + matrix_set(X, i, 2, ((rand() % 500) as f32) / 100.0) + } + return X +} + +function fe_labels(n: i32) -> ptr { + let y: ptr = array_new_f32(n) + for i in 0 to n { y[i] = (i % 2) as f32 } + return y +} + +function fe_accuracy(pred: ptr, y: ptr, n: i32) -> f32 { + let mut right: i32 = 0 + for i in 0 to n { + if (pred[i] as i32) == (y[i] as i32) { right = right + 1 } + } + return (right as f32) / (n as f32) +} + +function test_voting_classifier_full() -> i32 { + println("Test: VotingClassifier fits its own ensemble") + let n: i32 = 60 + let X: Matrix = fe_classes(n) + let y: ptr = fe_labels(n) + let model: VotingClassifier = voting_classifier_full_fit(X, y, 2, 3, 3, 42, 0) + let pred: ptr = voting_classifier_predict(model, X) + let acc: f32 = fe_accuracy(pred, y, n) + + let mut failures: i32 = 0 + if acc < 0.9 { + printf(" FAIL: accuracy %.4f, expected at least 0.9\n", acc) + failures = failures + 1 + } else { + printf(" OK: accuracy %.4f\n", acc) + } + array_free_f32(pred) + voting_classifier_free(model) + array_free_f32(y) + matrix_free(X) + return failures +} + +function test_stacking_classifier_full() -> i32 { + println("Test: StackingClassifier fits and cross-validates its base estimators") + let n: i32 = 60 + let X: Matrix = fe_classes(n) + let y: ptr = fe_labels(n) + let model: StackingClassifierFull = stacking_classifier_full_fit(X, y, 2, 3, 3, 42, 0.05, 60) + let pred: ptr = stacking_classifier_full_predict(model, X) + let acc: f32 = fe_accuracy(pred, y, n) + + let mut failures: i32 = 0 + if acc < 0.9 { + printf(" FAIL: accuracy %.4f, expected at least 0.9\n", acc) + failures = failures + 1 + } else { + printf(" OK: accuracy %.4f\n", acc) + } + array_free_f32(pred) + stacking_classifier_full_free(model) + array_free_f32(y) + matrix_free(X) + return failures +} + +function test_voting_regressor_full() -> i32 { + println("Test: VotingRegressor fits its own ensemble") + let n: i32 = 60 + let X: Matrix = matrix_new(n, 2) + let y: ptr = array_new_f32(n) + srand(9) + for i in 0 to n { + let a: f32 = ((rand() % 1000) as f32) / 100.0 + let b: f32 = ((rand() % 1000) as f32) / 100.0 + matrix_set(X, i, 0, a) + matrix_set(X, i, 1, b) + y[i] = 3.0 * a - 2.0 * b + } + let model: VotingRegressor = voting_regressor_full_fit(X, y, 3, 6, 42) + let pred: ptr = voting_regressor_predict(model, X) + + let mut mean: f32 = 0.0 + for i in 0 to n { mean = mean + y[i] } + mean = mean / (n as f32) + let mut ss_res: f32 = 0.0 + let mut ss_tot: f32 = 0.0 + for i in 0 to n { + let d: f32 = y[i] - pred[i] + ss_res = ss_res + d * d + let dm: f32 = y[i] - mean + ss_tot = ss_tot + dm * dm + } + let r2: f32 = 1.0 - ss_res / ss_tot + + let mut failures: i32 = 0 + if r2 < 0.8 { + printf(" FAIL: R2 %.4f, expected at least 0.8\n", r2) + failures = failures + 1 + } else { + printf(" OK: R2 %.4f\n", r2) + } + array_free_f32(pred) + voting_regressor_free(model) + array_free_f32(y) + matrix_free(X) + return failures +} + +function test_self_training_full() -> i32 { + println("Test: SelfTrainingClassifier labels from its own base estimator") + let n: i32 = 60 + let X: Matrix = fe_classes(n) + let y: ptr = fe_labels(n) + let y0: ptr = malloc((n as i64) * 4 + 64) as ptr + for i in 0 to n { + if i % 3 == 0 { y0[i] = y[i] as i32 } + else { y0[i] = 0 - 1 } + } + + let model: SelfTrainingClassifier = self_training_classifier_full_fit(X, y0, 2, 0.7, 10, 60, 0.1) + + let mut labelled: i32 = 0 + let mut right: i32 = 0 + for i in 0 to n { + if model.labels[i] >= 0 { + labelled = labelled + 1 + if model.labels[i] == (y[i] as i32) { right = right + 1 } + } + } + let coverage: f32 = (labelled as f32) / (n as f32) + let acc: f32 = (right as f32) / (labelled as f32) + + let mut failures: i32 = 0 + if coverage < 0.6 { + printf(" FAIL: only %.2f of the rows ended up labelled\n", coverage) + failures = failures + 1 + } + if acc < 0.9 { + printf(" FAIL: %.4f of the labels it assigned are right\n", acc) + failures = failures + 1 + } + if failures == 0 { + printf(" OK: %.2f labelled, %.4f of them right\n", coverage, acc) + } + + free(y0 as ptr) + self_training_classifier_free(model) + array_free_f32(y) + matrix_free(X) + return failures +} + +function test_incremental_pca_batches() -> i32 { + println("Test: IncrementalPCA walks the whole design") + let n: i32 = 60 + let X: Matrix = matrix_new(n, 3) + srand(3) + for i in 0 to n { + let t: f32 = ((rand() % 1000) as f32) / 100.0 + # The first two features carry most of the spread, the third little. + matrix_set(X, i, 0, t) + matrix_set(X, i, 1, 0.5 * t + ((rand() % 100) as f32) / 100.0) + matrix_set(X, i, 2, ((rand() % 20) as f32) / 100.0) + } + + let model: IncrementalPCA = incremental_pca_fit(X, 2, 16) + let Z: Matrix = incremental_pca_transform(model, X) + + # The leading component has to carry more spread than the second, which is + # the one property a principal component decomposition has to have. + let mut var0: f32 = 0.0 + let mut var1: f32 = 0.0 + for i in 0 to n { + let a: f32 = matrix_at(Z, i, 0) + let b: f32 = matrix_at(Z, i, 1) + var0 = var0 + a * a + var1 = var1 + b * b + } + + let mut failures: i32 = 0 + if model.n_samples_seen != n { + printf(" FAIL: saw %d rows of %d\n", model.n_samples_seen, n) + failures = failures + 1 + } + if var0 <= var1 { + printf(" FAIL: leading component spread %.4f, second %.4f\n", var0, var1) + failures = failures + 1 + } + if failures == 0 { + printf(" OK: %d rows seen, spread %.4f then %.4f\n", model.n_samples_seen, var0, var1) + } + + matrix_free(Z) + incremental_pca_free(model) + matrix_free(X) + return failures +} + +function main() -> i32 { + println("=== Ensembles that fit their own estimators ===") + let mut failures: i32 = 0 + failures = failures + test_voting_classifier_full() + failures = failures + test_stacking_classifier_full() + failures = failures + test_voting_regressor_full() + failures = failures + test_self_training_full() + failures = failures + test_incremental_pca_batches() + if failures > 0 { + printf("%d failure(s)\n", failures) + return 1 + } + println("All full ensemble tests passed!") + return 0 +}