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 17e8d94..32ad7c2 100644
--- a/benchmarks/estimator_coverage.json
+++ b/benchmarks/estimator_coverage.json
@@ -1,10 +1,9 @@
 {
   "schema_version": 1,
   "counts": {
-    "estimators": 203,
-    "runnable": 166,
-    "shaped": 22,
-    "simplified": 4,
+    "estimators": 208,
+    "runnable": 171,
+    "shaped": 26,
     "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"
     },
     {
@@ -3546,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",
@@ -3569,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"
     },
     {
@@ -7395,8 +7448,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"
     },
     {
@@ -10776,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<i32>"
+          },
+          {
+            "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<i32> = malloc((n_c as i64) * 4 + 64) as ptr<i32>",
+          "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",
@@ -11470,8 +11592,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 +11674,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"
     },
     {
@@ -11730,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<f32>"
+          },
+          {
+            "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<f32>",
+          "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<f32>",
+        "flow_free_fn": "stacking_classifier_full_free",
+        "sklearn_input": "X"
+      }
+    },
     {
       "flow_estimator": "stacking_regressor",
       "module": "ensemble.flow",
@@ -12576,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<f32>"
+          },
+          {
+            "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<f32>",
+        "flow_free_fn": "voting_classifier_free",
+        "sklearn_input": "X"
+      }
+    },
     {
       "flow_estimator": "voting_regressor",
       "module": "ensemble.flow",
@@ -12632,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<f32>"
+          },
+          {
+            "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<f32>",
+        "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<f32>",
+        "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<f32>",
+        "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<f32>",
+        "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<i32> = malloc((n_c as i64) * 4 + 64) as ptr<i32>",
+            "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<f32> = 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<void>)
     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<f32> = 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<f32> = 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<void>)
     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<f32> = 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<void>)
     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<i32, 1> = [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<f32> = 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<void>)
     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<f32> = 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<void>)
     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<void>)
     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<i32> = malloc((n_c as i64) * 4 + 64) as ptr<i32>
+    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<void>)
     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<f32> = 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<f32> = 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<f32> = 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<void>)
     matrix_free(Y_labels)
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<f32>, L: ptr<f32>, 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<f32>, 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<f32>, n: i32, p: i32, loc: ptr<f32>, L: ptr<f32>, out: ptr<f32>, z: ptr<f32>) -> void {
+    for i in 0 to n {
+        let row: ptr<f32> = 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<f32>, p: i32, idx: ptr<i32>, h: i32, loc: ptr<f32>, cov: ptr<f32>) -> void {
+    for j in 0 to p { loc[j] = 0.0 }
+    for t in 0 to h {
+        let row: ptr<f32> = 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<f32> = xdata + idx[t2] * p
+        for i2 in 0 to p {
+            let di: f32 = row2[i2] - loc[i2]
+            let crow: ptr<f32> = 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<f32>, n: i32, h: i32, order: ptr<i32>, idx: ptr<i32>) -> 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<f32> = 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<f32> = array_new_f32(p)
     let best_cov: ptr<ptr<f32> > = malloc((p as i64) * 8) as ptr<ptr<f32> >
     for i in 0 to p { best_cov[i] = array_new_f32(p) }
-    let mut best_det: f32 = 1000000000.0
-
-    let support_: ptr<i32> = array_new_f32(h) as ptr<i32>
+    let support_: ptr<i32> = malloc((h as i64) * 4 + 64) as ptr<i32>
+    let mut best_logdet: f32 = 1000000000.0
+    let mut have_best: bool = false
+
+    let loc: ptr<f32> = array_new_f32(p)
+    let cov: ptr<f32> = array_new_f32(p * p)
+    let L: ptr<f32> = array_new_f32(p * p)
+    let z: ptr<f32> = array_new_f32(p)
+    let dists: ptr<f32> = array_new_f32(n)
+    let order: ptr<i32> = malloc((n as i64) * 4 + 64) as ptr<i32>
+    let idx: ptr<i32> = malloc((h as i64) * 4 + 64) as ptr<i32>
 
     for trial in 0 to n_trials {
-        # Random subset of size h
-        let subset: ptr<i32> = array_new_f32(h) as ptr<i32>
-        for i in 0 to h {
-            subset[i] = rand() % n
-        }
-
-        # Compute mean of subset
-        let loc: ptr<f32> = 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<ptr<f32> > = malloc((p as i64) * 8) as ptr<ptr<f32> >
-        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<void>)
-        array_free_f32(subset as ptr<f32>)
     }
 
+    array_free_f32(loc)
+    array_free_f32(cov)
+    array_free_f32(L)
+    array_free_f32(z)
+    array_free_f32(dists)
+    free(order as ptr<void>)
+    free(idx as ptr<void>)
+
     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<f32> = array_new_f32(n)
 
-    # Compute Mahalanobis distance for each sample
-    # Simplified: use diagonal covariance inverse
-    let inv_diag: ptr<f32> = 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<f32> = 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<f32> = array_new_f32(p * p)
+    let z: ptr<f32> = 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<f32> = 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<i32> = malloc((n as i64) * 4 + 64) as ptr<i32>
+    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<void>)
+    array_free_f32(dists)
+    array_free_f32(cov)
+    array_free_f32(L)
+    array_free_f32(z)
     return result
 }
 
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<ptr<f32> >,
     mean: ptr<f32>,
     variances: ptr<f32>,
+    # 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<f32>,
     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<f32>, p: i32, k: i32, V: ptr<f32>, eigenvalues: ptr<f32>) -> void {
+    let Z: ptr<f32> = 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<f32> = array_new_f32(p)
+    for c2 in 0 to k {
+        let mut i7: i32 = 0
+        while i7 < p {
+            let mrow: ptr<f32> = 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<f32> = X.data
 
     let batch_mean: ptr<f32> = 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<f32> = 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<f32> = 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<f32> = 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<f32> = array_new_f32(p * p)
+    let diff: ptr<f32> = array_new_f32(p)
+    for i2 in 0 to n {
+        let row2: ptr<f32> = 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<f32> = 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<f32> = model.scatter + a2 * p
+            let brow: ptr<f32> = 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<f32> = array_new_f32(p * k)
+    let eigenvalues: ptr<f32> = array_new_f32(k)
+    let cov: ptr<f32> = 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<f32> = 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<f32> = 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<f32> = batch.data
+            let src: ptr<f32> = 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<void>)
     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<f32>, n_classes: i32, n_estimators: i32, max_depth: i32, seed: i32, voting: i32) -> VotingClassifier {
+    let trees: ptr<DecisionTreeClassifier> = malloc((n_estimators as i64) * 256) as ptr<DecisionTreeClassifier>
+    let state: ptr<u32> = malloc(8) as ptr<u32>
+    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<void>)
+
+    let classes: ptr<f32> = 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<DecisionTreeClassifier>, n_estimators: i32, n_classes: i32, classes: ptr<f32>, 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<DecisionTreeClassifier>,
+    n_base: i32,
+    n_classes: i32,
+    meta_weights: ptr<ptr<f32> >,
+    meta_bias: ptr<f32>,
+    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<f32>, 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<f32> = 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<f32> = array_new_f32(n * n_base)
+    let n_folds: i32 = 5
+    let fold_size: i32 = n / n_folds
+    let state: ptr<u32> = malloc(8) as ptr<u32>
+
+    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<f32> = array_new_f32(n_train)
+        let dst: ptr<f32> = X_tr.data
+        let mut at: i32 = 0
+        for i in 0 to n {
+            if i < test_start || i >= test_end {
+                let src: ptr<f32> = xdata + i * p
+                let drow: ptr<f32> = 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<f32> = X_te.data
+        for t in 0 to n_test {
+            let src2: ptr<f32> = xdata + (test_start + t) * p
+            let trow: ptr<f32> = 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<f32> = 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<ptr<f32> > = malloc((n_classes as i64) * 8) as ptr<ptr<f32> >
+    let meta_bias: ptr<f32> = 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<f32> = meta + i2 * n_base
+            let yi: i32 = y[i2] as i32
+            for c2 in 0 to n_classes {
+                let wc: ptr<f32> = 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<DecisionTreeClassifier> = malloc((n_base as i64) * 256) as ptr<DecisionTreeClassifier>
+    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<void>)
+    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<f32> {
+    let n: i32 = X.rows
+    let n_base: i32 = model.n_base
+    let meta: ptr<f32> = array_new_f32(n * n_base)
+    for b in 0 to n_base {
+        let preds: ptr<f32> = 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<f32> = array_new_f32(n)
+    for i2 in 0 to n {
+        let row: ptr<f32> = 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<f32> = 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<void>)
+    for c in 0 to model.n_classes {
+        array_free_f32(model.meta_weights[c])
+    }
+    free(model.meta_weights as ptr<void>)
+    array_free_f32(model.meta_bias)
+}
+
 export function stacking_classifier_fit(
     base_preds: ptr<ptr<f32> >,
     y: ptr<f32>,
@@ -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<f32>, n_estimators: i32, max_depth: i32, seed: i32) -> VotingRegressor {
+    let trees: ptr<DecisionTreeRegressor> = malloc((n_estimators as i64) * 256) as ptr<DecisionTreeRegressor>
+    for e in 0 to n_estimators {
+        trees[e] = decision_tree_regressor_fit(X, y, max_depth, 0)
+    }
+    let weights: ptr<f32> = 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<DecisionTreeRegressor>, n_estimators: i32, weights: ptr<f32>) -> VotingRegressor {
     let result: VotingRegressor = VotingRegressor {
         # The model takes over the fitted sub-estimators and its free
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<f32>, 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<f32> = array_new_f32(n)
-    let col_means: ptr<f32> = 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<f32>, n: i32, p: i32, k: i32, V: ptr<f32>) -> void {
+    let M: ptr<f32> = 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<f32> = 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<f32>, m: i32, d: i32, k: i32, labels: ptr<i32>) -> void {
+    let centers: ptr<f32> = array_new_f32(k * d)
+    let counts: ptr<i32> = malloc((k as i64) * 4 + 64) as ptr<i32>
+    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<f32> = 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<f32> = 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<i32> = array_new_f32(n) as ptr<i32>
-    let col_labels: ptr<i32> = array_new_f32(p) as ptr<i32>
+        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<f32> = Z + i2 * d
+            let cc2: ptr<f32> = 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<void>)
+}
+
+# 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<f32>, n: i32, p: i32) -> void {
+    let r: ptr<f32> = array_new_f32(n)
+    let c: ptr<f32> = 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<f32> = 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<f32> = 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<f32> = 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<f32>, m: i32, k: i32) -> f32 {
+    let labels: ptr<i32> = malloc((m as i64) * 4 + 64) as ptr<i32>
+    _bic_kmeans(v, m, 1, k, labels)
+
+    let sums: ptr<f32> = array_new_f32(k)
+    let counts: ptr<i32> = malloc((k as i64) * 4 + 64) as ptr<i32>
+    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<void>)
+    free(counts as ptr<void>)
+    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<f32> = X.data
+
+    let An: ptr<f32> = 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<f32> = array_new_f32(p * n_vecs)
+    _bic_right_singular(An, n, p, n_vecs, V)
+
+    let U: ptr<f32> = 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<f32> = array_new_f32(n)
+    let errs: ptr<f32> = array_new_f32(n_vecs)
+    let chosen: ptr<i32> = malloc((n_vecs as i64) * 4 + 64) as ptr<i32>
+    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<f32> = array_new_f32(n * n_best)
+    let Zc: ptr<f32> = 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<i32> = malloc((n as i64) * 4 + 64) as ptr<i32>
+    let col_labels: ptr<i32> = malloc((p as i64) * 4 + 64) as ptr<i32>
+    _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<void>)
 
     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<f32> = X.data
 
-    let row_labels: ptr<i32> = array_new_f32(n) as ptr<i32>
-    let col_labels: ptr<i32> = array_new_f32(p) as ptr<i32>
+    let row_labels: ptr<i32> = malloc((n as i64) * 4 + 64) as ptr<i32>
+    let col_labels: ptr<i32> = malloc((p as i64) * 4 + 64) as ptr<i32>
 
-    # Simplified: cluster rows and columns together by value magnitude
-    let mut row_sums: ptr<f32> = array_new_f32(n)
+    # Row and column sums, and the inverse square roots the normalization uses.
+    let r_inv: ptr<f32> = array_new_f32(n)
+    let c_inv: ptr<f32> = array_new_f32(p)
+    let c_sum: ptr<f32> = array_new_f32(p)
     for i in 0 to n {
+        let row: ptr<f32> = 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<f32> = array_new_f32(n * p)
+    for i2 in 0 to n {
+        let src: ptr<f32> = xdata + i2 * p
+        let dst: ptr<f32> = 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<f32> = 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<f32> = 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<f32> = 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<i32> = malloc((m as i64) * 4 + 64) as ptr<i32>
+    _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<void>)
+    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<f32>, 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<f32> = X.data
 
-    let transform: ptr<ptr<f32> > = malloc((p as i64) * 8) as ptr<ptr<f32> >
+    # The projection, held flat while it is being fitted and handed back a row
+    # at a time at the end.
+    let T: ptr<f32> = 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<f32> = array_new_f32(n * d)
+    let prob: ptr<f32> = array_new_f32(n)
+    let W: ptr<f32> = array_new_f32(p * p)
+    let grad: ptr<f32> = array_new_f32(p * d)
+    let diff: ptr<f32> = array_new_f32(p)
+
     for epoch in 0 to epochs {
-        let grad: ptr<ptr<f32> > = malloc((p as i64) * 8) as ptr<ptr<f32> >
-        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<f32> = 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<f32> = 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<f32> = 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<f32> = 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<f32> = 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<f32> = 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<ptr<f32> > = malloc((p as i64) * 8) as ptr<ptr<f32> >
+    for i2 in 0 to p {
+        let row: ptr<f32> = 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<void>)
+        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/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<void>
@@ -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<i32>, 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<f32> = X.data
+
+    let labels: ptr<i32> = malloc((n as i64) * 4 + 64) as ptr<i32>
+    let confidence: ptr<f32> = array_new_f32(n)
+    for i in 0 to n {
+        labels[i] = y_initial[i]
+        confidence[i] = 0.0
+    }
+
+    let scores: ptr<f32> = array_new_f32(n * n_classes)
+    let y_lab: ptr<f32> = 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<f32> = X_lab.data
+        let mut at: i32 = 0
+        for i3 in 0 to n {
+            if labels[i3] >= 0 {
+                let src: ptr<f32> = xdata + i3 * p
+                let drow: ptr<f32> = 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<f32> = 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<i32>,
     probabilities: ptr<ptr<f32> >,
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<f32> {
+    let y: ptr<f32> = array_new_f32(n)
+    for i in 0 to n { y[i] = (i % 2) as f32 }
+    return y
+}
+
+function fe_accuracy(pred: ptr<f32>, y: ptr<f32>, 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<f32> = fe_labels(n)
+    let model: VotingClassifier = voting_classifier_full_fit(X, y, 2, 3, 3, 42, 0)
+    let pred: ptr<f32> = 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<f32> = fe_labels(n)
+    let model: StackingClassifierFull = stacking_classifier_full_fit(X, y, 2, 3, 3, 42, 0.05, 60)
+    let pred: ptr<f32> = 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<f32> = 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<f32> = 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<f32> = fe_labels(n)
+    let y0: ptr<i32> = malloc((n as i64) * 4 + 64) as ptr<i32>
+    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<void>)
+    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
+}
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<f32> = 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<f32>, 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<f32> = 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<i32>, 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
+}