dataframe-learn 1.0.2.0 → 1.1.0.0

raw patch · 46 files changed

+5735/−791 lines, 46 filesdep +randomdep ~dataframe-coredep ~dataframe-operationsdep ~parallelPVP ok

version bump matches the API change (PVP)

Dependencies added: random

Dependency ranges changed: dataframe-core, dataframe-operations, parallel, text

API changes (from Hackage documentation)

Files

• README.md +473/−0
• dataframe-learn.cabal +43/−8
• src/DataFrame/Boosting.hs +10/−0
• src/DataFrame/Boosting/AdaBoost.hs +211/−0
• src/DataFrame/Boosting/GBM.hs +185/−0
• src/DataFrame/DBSCAN.hs +119/−0
• src/DataFrame/DecisionTree/Cart.hs +90/−29
• src/DataFrame/DecisionTree/Categorical.hs +154/−54
• src/DataFrame/DecisionTree/CondVec.hs +54/−29
• src/DataFrame/DecisionTree/Fit.hs +60/−21
• src/DataFrame/DecisionTree/Linear.hs +57/−29
• src/DataFrame/DecisionTree/Model.hs +101/−0
• src/DataFrame/DecisionTree/Numeric.hs +59/−23
• src/DataFrame/DecisionTree/Pool.hs +82/−32
• src/DataFrame/DecisionTree/Predict.hs +85/−30
• src/DataFrame/DecisionTree/Prune.hs +16/−9
• src/DataFrame/DecisionTree/Regression.hs +144/−0
• src/DataFrame/DecisionTree/Tao.hs +136/−38
• src/DataFrame/DecisionTree/Types.hs +15/−10
• src/DataFrame/Featurize/Internal.hs +132/−0
• src/DataFrame/GMM.hs +310/−0
• src/DataFrame/KMeans.hs +156/−0
• src/DataFrame/LinearAlgebra.hs +122/−0
• src/DataFrame/LinearAlgebra/Eigen.hs +121/−0
• src/DataFrame/LinearAlgebra/Solve.hs +213/−0
• src/DataFrame/LinearModel.hs +11/−0
• src/DataFrame/LinearModel/Logistic.hs +105/−0
• src/DataFrame/LinearModel/Regression.hs +128/−0
• src/DataFrame/LinearSolver.hs +73/−31
• src/DataFrame/LinearSolver/Loss.hs +47/−0
• src/DataFrame/Metrics.hs +237/−0
• src/DataFrame/Metrics/Report.hs +170/−0
• src/DataFrame/Model.hs +69/−0
• src/DataFrame/ModelSelection.hs +91/−0
• src/DataFrame/PCA.hs +131/−0
• src/DataFrame/PCA/Kernel.hs +136/−0
• src/DataFrame/Random.hs +121/−0
• src/DataFrame/SVM.hs +101/−0
• src/DataFrame/SVM/RFF.hs +146/−0
• src/DataFrame/SymbolicRegression.hs +120/−0
• src/DataFrame/SymbolicRegression/Expr.hs +122/−0
• src/DataFrame/SymbolicRegression/GP.hs +207/−0
• src/DataFrame/SymbolicRegression/Optimize.hs +68/−0
• src/DataFrame/SymbolicRegression/Simplify.hs +65/−0
• src/DataFrame/Synthesis.hs +336/−448
• src/DataFrame/Transform.hs +103/−0

@@ -0,0 +1,473 @@+<!--+  This README is a runnable scripths (https://github.com/DataHaskell/scripths)+  notebook. Every ```haskell block runs top-to-bottom in one shared session and+  scripths inserts each block's output beneath it as a blockquote. The+  `-- cabal: packages:` directive builds against the local working tree.+  Regenerate (from the repo root) with:++      scripths dataframe-learn/README.md -o dataframe-learn/README.md+-->++# dataframe-learn++Machine learning for [`dataframe`](https://hackage.haskell.org/package/dataframe)+where **a fitted model is a dataframe expression**. You `fit` a model and+`predict` hands you back an `Expr` over your columns — pretty-print it to read+the formula, apply it with `derive` to score a frame, fold preprocessing into it+with `compileThrough`. The model *is* the prediction, not an opaque blob, and the+scikit-learn-style record (coefficients, centroids, components, support) is right+there too for inspection. Because every prediction is the same kind of `Expr`,+preprocessing, prediction, and deployment all compose the same way — the+[design notes](#design-notes-the-categorical-account) at the end explain why.++## A linear model is a formula++`fit` returns a record (with `regCoef`/`regIntercept` for inspection) and+`predict` compiles it to an `Expr Double` you can read. The `D` import (the+public `DataFrame` umbrella, which also gives `D.col` for the expression DSL)+carries through the rest of the notebook; each later section adds the one model+module it needs:++```haskell+-- cabal: packages: .., ., ../dataframe-core, ../dataframe-operations, ../dataframe-parsing+-- cabal: build-depends: dataframe, dataframe-learn, text+-- cabal: default-extensions: OverloadedStrings, TypeApplications+-- cabal: ghc-options: -w+import qualified DataFrame as D+import DataFrame.LinearModel+import DataFrame.Model (fit, predict)++sales = D.fromNamedColumns+    [ ("x", D.fromList ([1, 2, 3, 4, 5, 6] :: [Double]))+    , ("y", D.fromList ([2 * x + 1 | x <- [1, 2, 3, 4, 5, 6]] :: [Double]))+    ]++model = fit defaultLinearConfig (D.col @Double "y") sales+putStrLn (D.prettyPrint (predict model))+```++> <!-- scripths:mime text/plain -->+> 2.0 * x + 0.9999999999999989++## A decision tree is a readable expression++The tree compiles to nested `if/then/else` over your columns — no special+viewer, it is just an expression:++```haskell+import DataFrame.DecisionTree (defaultTreeConfig)+import DataFrame.DecisionTree.Model ()++flowers = D.fromNamedColumns+    [ ("petal_length", D.fromList ([1.4, 1.3, 1.5, 1.4, 4.5, 4.7, 4.6, 4.4, 5.5, 5.8, 5.6, 5.7] :: [Double]))+    , ("petal_width",  D.fromList ([0.2, 0.2, 0.1, 0.3, 1.5, 1.4, 1.6, 1.3, 2.0, 2.1, 1.9, 2.2] :: [Double]))+    , ("species",      D.fromList ([0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2] :: [Double]))+    ]++tree = fit defaultTreeConfig (D.col @Double "species") flowers+putStrLn (D.prettyPrint (predict tree))+```++> <!-- scripths:mime text/plain -->+> if petal_length .<=. 2.95+>   then 0.0+>   else if petal_length .<=. 5.1+>     then 1.0+>     else 2.0++## Symbolic regression discovers a formula++Genetic programming searches for an expression that fits the data, and returns+it as a dataframe `Expr` plus the accuracy/complexity Pareto front:++```haskell+import DataFrame.SymbolicRegression++curve = D.fromNamedColumns+    [ ("x", D.fromList xs)+    , ("y", D.fromList [x * x + x | x <- xs])+    ]+  where xs = [-3, -2, -1, 0, 1, 2, 3, 4, 5, 6] :: [Double]++sr = fit+        defaultSRConfig { srSeed = 3, srGenerations = 50, srPopSize = 300, srUnaryOps = [] }+        (D.col @Double "y") curve+putStrLn (D.prettyPrint (srBest sr) ++ "   (mse " ++ show (srBestMSE sr) ++ ")")+```++> <!-- scripths:mime text/plain -->+> x + x * x   (mse 0.0)++## When the formula is bigger than a glance++Not every model is a one-liner. A linear model, a small tree, or a symbolic+expression you can *read*; a 40-tree gradient booster you cannot — its `predict`+is an exact sum of forty trees. Counting the characters in each printed formula+shows the gap:++```haskell+import DataFrame.Boosting++gbm = fit defaultGBConfig { gbNEstimators = 40, gbMaxDepth = 2 } (D.col @Double "y") sales+putStr (unlines+    [ "linear prediction: " ++ show (length (D.prettyPrint (predict model))) ++ " characters"+    , "GBM(40 trees):     " ++ show (length (D.prettyPrint (predict gbm)))  ++ " characters" ])+```++> <!-- scripths:mime text/plain -->+> linear prediction: 28 characters+> GBM(40 trees):     7151 characters++Even when it is too big to eyeball, the expression is still the whole story:+a self-contained, dependency-free artifact that scores a frame with `derive` —+no pickled blob, no runtime to ship. For the big ensembles the *interpretability*+comes from `gbFeatureImportances` and pretty-printing individual trees, not from+reading the summed formula.++## Deploy: applying an expression to a frame++Because the model is an `Expr`, deploying it is just `derive` — you add the+prediction as a new column with the ordinary dataframe API:++```haskell+D.columnNames (D.derive "prediction" (predict model) sales)+```++> <!-- scripths:mime text/plain -->+> ["x","y","prediction"]++## A model and its preprocessing compose by substitution++Preprocessing is an expression too, so a model trained in a transformed space and+the transform that produced it *compose* — and composition of expressions is+substitution of one into the other. `compileThrough` performs that composition,+folding a fitted transform into a prediction so the result is a single formula+over the raw inputs. Here we standardize `x`, fit in the scaled space, then fold+the scaler back in to recover a raw-column model:++```haskell+import DataFrame.Transform+import DataFrame.Metrics++scaler      = standardScaler ["x"] sales+scaledSales = applyTransform (scalerTransform scaler) sales+scaledModel = fit defaultLinearConfig (D.col @Double "y") scaledSales++deployed = compileThrough (scalerTransform scaler) (predict scaledModel)+putStr (unlines+    [ "trained in scaled space: " ++ D.prettyPrint (predict scaledModel)+    , "folded to raw columns:   " ++ D.prettyPrint deployed ])+```++> <!-- scripths:mime text/plain -->+> trained in scaled space: 3.4156502553198655 * x + 8.0+> folded to raw columns:   3.4156502553198655 * (x - 3.5) / 1.707825127659933 + 8.0++The folded expression is a function of the raw `x` alone, so it scores the+original frame with no preprocessing step at inference time — and by the+substitution lemma it computes the same result (up to floating point) as+transforming the frame and then predicting:++```haskell+evaluate rmse deployed (D.col @Double "y") sales+```++> <!-- scripths:mime text/plain -->+> 3.6259732146947156e-16++## A realistic run: pick features, split, evaluate held-out, tune++Real frames are noisy and carry columns you must not train on. Here is a noisy+linear signal with a spurious `id` column:++```haskell+realistic = D.fromNamedColumns+    [ ("id", D.fromList [fromIntegral ((i * 7919) `mod` 97) | i <- [1 .. 40 :: Int]])+    , ("x",  D.fromList xs)+    , ("y",  D.fromList [2 * x + 1 + noise i | (i, x) <- zip [0 :: Int ..] xs])+    ]+  where+    xs      = map fromIntegral [1 .. 40 :: Int] :: [Double]+    noise i = fromIntegral ((i * 2654435761 + 12345) `mod` 1000) / 100 - 5+```++> <!-- scripths:mime text/plain -->++**Feature selection.** Supervised `fit` uses *every* non-target column as a+feature, so a naive fit drags `id` into the model. `selectFeatures` restricts to+the columns you mean (mirroring the explicit feature list the unsupervised+fitters take), which is the difference between a leaky model and a clean one:++```haskell+import DataFrame.Model (selectFeatures)++naive   = fit defaultLinearConfig (D.col @Double "y") realistic+guarded = fit defaultLinearConfig (D.col @Double "y")+              (selectFeatures ["x"] (D.col @Double "y") realistic)+putStr (unlines+    [ "all columns:           " ++ D.prettyPrint (predict naive)+    , "selectFeatures [\"x\"]:   " ++ D.prettyPrint (predict guarded) ])+```++> <!-- scripths:mime text/plain -->+> all columns:           -7.746701620642152e-3 * id + 1.9914915483217268 * x + 1.6474622984919354+> selectFeatures ["x"]:   1.9918011257035637 * x + 1.2630769230769452++**Hold-out evaluation.** `trainTestSplit` (seeded, deterministic) keeps the score+honest — evaluate on rows the model never saw, and the metrics are realistic, not+the `1e-15` of an in-sample toy:++```haskell+import DataFrame.ModelSelection++clean        = selectFeatures ["x"] (D.col @Double "y") realistic+(train, test) = trainTestSplit 0.75 7 clean+heldModel     = fit defaultLinearConfig (D.col @Double "y") train+putStr (unlines+    [ "held-out R^2:  " ++ show (evaluate r2   (predict heldModel) (D.col @Double "y") test)+    , "held-out RMSE: " ++ show (evaluate rmse (predict heldModel) (D.col @Double "y") test) ])+```++> <!-- scripths:mime text/plain -->+> held-out R^2:  0.9671190074242891+> held-out RMSE: 3.56674709632647++**Cross-validation.** `crossValidate` is scikit-learn's `cross_val_score`: it+fits on each training fold and scores the prediction expression on the held-out+fold. You pass a `train -> Expr` closure, so it works with any model:++```haskell+cv = crossValidate 5 0 rmse (D.col @Double "y")+         (\tr -> predict (fit defaultLinearConfig (D.col @Double "y") tr))+         clean+putStrLn ("5-fold RMSE: " ++ show (sum cv / fromIntegral (length cv)))+```++> <!-- scripths:mime text/plain -->+> 5-fold RMSE: 3.0325616706245713++`gridSearch` tunes hyperparameters the same way, over a list of configs.++## Reports without hand-rolling metrics++Metrics are plain functions (`rmse`, `mse`, `r2`, `accuracy`, multiclass+`precision`/`recall`/`f1`), and `classificationReport` bundles the common numbers+with a scikit-learn-style layout (per-class precision/recall/F1/support plus+macro/weighted averages):++```haskell+import DataFrame.Metrics.Report++clf = fit defaultLogisticConfig (D.col @Double "species") flowers+putStr (show (classificationReportExpr (predict clf) (D.col @Double "species") flowers))+```++> <!-- scripths:mime text/plain -->+> class       precision   recall      f1          support     +> 0.0         1.0         1.0         1.0         4           +> 1.0         1.0         1.0         1.0         4           +> 2.0         1.0         1.0         1.0         4           +> +> accuracy    = 1.0+> macro f1    = 1.0+> weighted f1 = 1.0++## Pipelines compose as a monoid++A fitted preprocessing step is a `Transform`, and transforms compose with `<>`.+`applyTransform` runs the whole pipeline; `compileThrough` folds it into a single+expression over the raw columns for export:++```haskell+import DataFrame.PCA++features = ["petal_length", "petal_width"]+scalerF  = standardScaler features flowers+pca      = fit (PCAConfig (NComp 2) True) (map (D.col @Double) features) flowers+pipeline = scalerTransform scalerF <> pcaTransform pca++D.columnNames (applyTransform pipeline flowers)+```++> <!-- scripths:mime text/plain -->+> ["petal_length","petal_width","species","pc1","pc2"]++## Unsupervised models are inspectable too++k-means returns `cluster_centers_`-style centroids, and per-cluster distance /+assignment expressions:++```haskell+import DataFrame.KMeans++km = fit defaultKMeansConfig { kmK = 3, kmSeed = 1 } (map (D.col @Double) features) flowers+kmCenters km+```++> <!-- scripths:mime text/plain -->+> [[1.4,0.2],[5.65,2.05],[4.55,1.4500000000000002]]++## Synthesize the feature you would have hand-engineered++`DataFrame.Synthesis` is automated feature engineering: a bottom-up enumerative+search (with observational-equivalence pruning) for a small, interpretable+expression over your columns that tracks the target. Here `y` is the interaction+`a * b`, which a linear model on the raw columns cannot capture; synthesis+discovers the term, and feeding it back as a column lifts the fit from mediocre+to exact — still a formula you can read:++```haskell+import DataFrame.Synthesis++interactions = D.fromNamedColumns+    [ ("a", D.fromList as)+    , ("b", D.fromList bs)+    , ("y", D.fromList (zipWith (*) as bs))+    ]+  where+    as = [-1, -1, 1, 1, -2, 2, -2, 2] :: [Double]+    bs = [-1, 1, -1, 1, -2, -2, 2, 2] :: [Double]++rawModel = fit defaultLinearConfig (D.col @Double "y") interactions+feature  = fit defaultSynthesisConfig (D.col @Double "y") interactions+withFeat = D.derive "synth" (predict feature) interactions+fitModel =+    fit defaultLinearConfig (D.col @Double "y")+        (selectFeatures ["synth"] (D.col @Double "y") withFeat)++putStr (unlines+    [ "discovered feature: " ++ D.prettyPrint (predict feature)+    , "raw linear R^2:     " ++ show (evaluate r2 (predict rawModel) (D.col @Double "y") interactions)+    , "with synth feature: " ++ show (evaluate r2 (predict fitModel) (D.col @Double "y") withFeat)+    ])+```++> <!-- scripths:mime text/plain -->+> discovered feature: a * b+> raw linear R^2:     0.0+> with synth feature: 1.0++`predict feature` is the single best expression; `sfFeatures feature` is the whole+ranked, deduplicated bank, ready to `derive` as a batch of candidate columns.++## What's in the box++| Task | Models |+|------|--------|+| Regression | OLS, ridge, lasso, elastic net, regression trees, gradient boosting, symbolic regression |+| Classification | logistic regression, linear SVC, RFF kernel SVM, decision trees, gradient boosting, AdaBoost |+| Dimensionality reduction | PCA, Nyström kernel PCA |+| Clustering | k-means, Gaussian mixtures, DBSCAN |+| Feature engineering | `DataFrame.Synthesis` (enumerative feature synthesis), symbolic regression |+| Evaluation | `DataFrame.Metrics` (metrics + `evaluate`), `DataFrame.Metrics.Report` (reports) |+| Pipelines & tuning | `DataFrame.Transform` (composable transforms), `DataFrame.ModelSelection` (`trainTestSplit`, `crossValidate`, `gridSearch`) |++Every model is a `Fit` instance, so there is one verb to train — `fit cfg input+df` — and every model with an honest out-of-sample prediction is a `Predict`+instance, so one verb to compile it — `predict model`. Auxiliary outputs+(`gbProbaExpr`, `logisticProbExprs`, `kmeansDistanceExprs`, `pcaTransform`, …)+keep descriptive names; transductive models like DBSCAN deliberately have no+`Predict` instance. Fits that use randomness take a `seed` in their config, so+results are reproducible across Linux, macOS, and Windows. Pure Haskell — the+only extra dependency beyond the dataframe packages is `random`.++## Design notes: the categorical account++The two verbs live in `DataFrame.Model`:++```+class Fit cfg input model | cfg input -> model where+    fit :: cfg -> input -> DataFrame -> model++class Predict model r | model -> r where+    predict :: model -> Expr r+```++They are small on purpose, because the structure they hang on lives in the+expression language, not in the classes. The framing borrows from+[*Seven Sketches in Compositionality*](https://arxiv.org/abs/1803.05316)+(Fong & Spivak) and the Para/Lens account of learners+([Fong, Johnson & Spivak, *Lenses and Learners*](https://arxiv.org/abs/1903.03671);+[Cruttwell et al., *Categorical Foundations of Gradient-Based Learning*](https://arxiv.org/abs/2103.01931)).+What follows is deliberately careful about what is load-bearing and what is only+analogy.++**The row-wise fragment is a category.** Restrict to the row-wise expression+constructors — `Col`, `Lit`, `Unary`, `Binary`, `If`. Take typed column contexts+as objects and, as an arrow `Γ → Δ`, a `Δ`-tuple of such expressions over `Γ`.+Composition is simultaneous substitution (`substituteColumns`, added to+`dataframe-core` for exactly this) and the identities are the column projections+(`Col`). This is the category of contexts (the Lawvere theory) of the column+signature. The restriction is load-bearing for *both* laws, not just composition:+`Agg` and `Over` are column-level/relational, not row-wise maps, and the raw-text+column reference inside `CastWith` is opaque to substitution — so identity-by-`Col`+fails on those constructors too. They are excluded by construction (transforms+reject `Agg`/`Over`), which is why composition and identity stay well defined.++**`predict` gives every model a uniform codomain.** `Predict model r` interprets a+fitted model as an arrow in that category: `predict model :: Expr r` runs from the+model's feature context to the one-column context `{r}`, and the dependency+`model -> r` fixes the codomain object. This is *not* a functor or a denotation in+the technical sense — there is no category of models to be functorial over. The+real, useful property is uniformity: every model's prediction lands in the *same*+expression type (`Expr Double`/`Expr a`/`Expr Int`), so `derive`, the `Transform`+monoid, and `compileThrough` all apply with no per-model glue. That the compiled+`Expr` actually agrees with the fitted record's own parameters is a tested property+(`tests/Learn/Denotation.hs`), not a typeclass law — the class only knows the+symbolic half.++**`fit` is the parametrized-morphism (Para) fragment.** `fit cfg input df` chooses a+parameter — the trained record — and `predict` is the forward map applied at it.+In the Para/Lens picture of learning a learner is a *parametrized lens* carrying a+forward map plus backward update/request maps; we inhabit only the forward (Para)+part and expose no backward maps, because this interface is batch training, not+online gradient exchange. That is a complete, self-contained sub-structure, not a+half-built one — but it does mean the Lens vocabulary is motivation here, not+something the code instantiates. (The functional dependency `cfg input -> model`+fixes the parameter *type*; `fit` is the value-level map that picks the point.)++**`Transform` is a monoid of derived-column lists.** `Transform`'s `<>` keeps the+earlier step's outputs and rewrites the later step's column references through them+by substitution; `mempty` is the empty list. These are context-*extending* maps+(`applyTransform` adds columns), so this is an ordinary algebraic monoid — a monoid+*is* a one-object category (Seven Sketches ch. 3) — not the endomorphism monoid of a+fixed object. Associativity and identity hold for the row-wise fragment **provided+output names do not collide**: the implementation merges output maps with+`Data.Map.fromList`, which keeps the last binding on a clash, so reusing a column+name across steps is the one way to break the law.++**Composition is the point.** `compileThrough t (predict m)` realizes the composite+`predict m ∘ t` (read right-to-left: first `t`, then `predict m`) by substituting+`t`'s definitions into `predict m`. By the substitution lemma it denotes the same+function as transforming the frame and then predicting — equal results up to+floating point, not syntactically identical expressions. That is exactly the+"compose by substitution" example above, and it is why a model trained in a+transformed space deploys as one formula over the raw columns.++**What deliberately has no `predict` — two different reasons.** DBSCAN is+transductive: every clustering *fit* depends on the whole training set, but what+distinguishes the models is whether the *fitted* model induces an out-of-sample+rule. k-means (nearest centroid) and GMM (max posterior) do, so they have honest+`predict` arrows; DBSCAN's density-reachability assignment has no per-row rule, so+we give it no `Predict` instance rather than a fake `Maybe` or a throwing stub.+PCA and kernel PCA are the *opposite* case: they *are* arrows, but multi-output+feature maps with no privileged label column, so their canonical interface is a+`Transform` (`pcaTransform`/`pcaExprs`), not a one-column `predict`.++**A note on classifiers.** A multiclass `predict` is a genuine arrow into the label+object, but it compiles arg-max to a nested-`If` cascade (`argMaxExpr`), quadratic+in the number of classes — so for a 5-class model `prettyPrint (predict m)` is an+If-tree, not a tidy formula. The "model is a readable formula" aesthetic is honest+for affine and tree models; for classifiers and clusterers the value is that the+arrow exists and composes, not that it is short.++**An aside.** A linear or affine model's prediction is a signal-flow graph — a+weighted sum of inputs. `affineExpr` builds the arrow in the prop of affine *maps*+(the single-valued sub-prop of Seven Sketches ch. 5's signal-flow calculus of+affine relations), and dropping zero-weight terms is diagram simplification —+deleting a zero-gain wire.++(The "instance is a functor `C → Set`" slogan from Spivak's functorial data model,+Seven Sketches ch. 3, is sometimes invoked for dataframes; a single flat table is+the degenerate case — a schema with no foreign-key morphisms — so it is an analogy+here, not a structure we use.)

@@ -1,12 +1,14 @@ cabal-version:      2.4 name:               dataframe-learn-version:            1.0.2.0--synopsis:           Decision trees and feature synthesis for the dataframe ecosystem.+version:            1.1.0.0+synopsis:           Interpretable, expression-returning machine learning for the dataframe ecosystem. description:-    @DataFrame.DecisionTree@ — decision-tree training on DataFrames.-    @DataFrame.Synthesis@ — feature synthesis. Built on top of-    @dataframe-operations@.+    A small scikit-learn-style ML library where every model returns both an+    inspectable record and dataframe @Expr@ value(s): linear/ridge/lasso/+    elastic-net and logistic regression, linear and RFF-kernel SVMs, decision+    trees, gradient boosting and AdaBoost, PCA and Nyström kernel PCA, k-means,+    Gaussian mixtures, DBSCAN, and symbolic regression — plus cross-validation+    and grid search. Pure Haskell, built on @dataframe-operations@.  bug-reports:        https://github.com/mchav/dataframe/issues license:            MIT@@ -16,6 +18,7 @@ copyright:          (c) 2024-2026 Michael Chavinda category:           Data tested-with:        GHC ==9.4.8 || ==9.6.7 || ==9.8.4 || ==9.10.3 || ==9.12.2+extra-doc-files:    README.md  common warnings     ghc-options:@@ -41,12 +44,44 @@                         DataFrame.DecisionTree.Tao                         DataFrame.DecisionTree.Fit                         DataFrame.LinearSolver+                        DataFrame.LinearSolver.Loss+                        DataFrame.LinearAlgebra+                        DataFrame.LinearAlgebra.Solve+                        DataFrame.LinearAlgebra.Eigen+                        DataFrame.Random+                        DataFrame.Featurize.Internal+                        DataFrame.Model+                        DataFrame.LinearModel+                        DataFrame.LinearModel.Regression+                        DataFrame.LinearModel.Logistic+                        DataFrame.SVM+                        DataFrame.DecisionTree.Regression+                        DataFrame.DecisionTree.Model+                        DataFrame.PCA+                        DataFrame.PCA.Kernel+                        DataFrame.SVM.RFF+                        DataFrame.KMeans+                        DataFrame.Transform+                        DataFrame.Boosting+                        DataFrame.Boosting.GBM+                        DataFrame.Boosting.AdaBoost+                        DataFrame.GMM+                        DataFrame.DBSCAN+                        DataFrame.Metrics+                        DataFrame.Metrics.Report+                        DataFrame.ModelSelection+                        DataFrame.SymbolicRegression+                        DataFrame.SymbolicRegression.Expr+                        DataFrame.SymbolicRegression.Simplify+                        DataFrame.SymbolicRegression.Optimize+                        DataFrame.SymbolicRegression.GP                         DataFrame.Synthesis     build-depends:      base >= 4 && < 5,                         containers >= 0.6.7 && < 0.9,                         parallel ^>= 3.2,-                        dataframe-core ^>= 1.0,-                        dataframe-operations ^>= 1.1,+                        random >= 1.2 && < 2,+                        dataframe-core ^>= 1.1,+                        dataframe-operations ^>= 1.1.1,                         text >= 2.0 && < 3,                         vector ^>= 0.13,                         vector-algorithms ^>= 0.9

@@ -0,0 +1,10 @@+{- | Tree ensembles: gradient boosting (the recommended default) and AdaBoost+(SAMME). Re-exports the focused submodules.+-}+module DataFrame.Boosting (+    module DataFrame.Boosting.GBM,+    module DataFrame.Boosting.AdaBoost,+) where++import DataFrame.Boosting.AdaBoost+import DataFrame.Boosting.GBM

@@ -0,0 +1,211 @@+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE TypeApplications #-}+{-# LANGUAGE UndecidableInstances #-}++{- | AdaBoost (SAMME) over short, sample-weighted classification trees. The+weighted-Gini stump fitter here is self-contained (it reuses the CART feature+encoding but not the unweighted CART recursion), so the existing decision-tree+path is untouched. 'predict' is the arg-max of weighted votes.+-}+module DataFrame.Boosting.AdaBoost (+    AdaBoostConfig (..),+    defaultAdaBoostConfig,+    AdaBoostModel (..),+) where++import Data.List (sort)+import Data.Maybe (fromMaybe, maybeToList)+import qualified Data.Vector as V+import qualified Data.Vector.Unboxed as VU++import DataFrame.DecisionTree.Cart (+    CartFeature (..),+    cartFeatures,+    sortIndicesByValue,+ )+import DataFrame.DecisionTree.Fit (treeToExpr)+import DataFrame.DecisionTree.Types (Tree (..))+import DataFrame.Featurize.Internal (argMaxExpr, targetValues)+import qualified DataFrame.Functions as F+import DataFrame.Internal.Column (Columnable, TypedColumn (..), toVector)+import DataFrame.Internal.DataFrame (DataFrame)+import DataFrame.Internal.Expression (Expr (..))+import DataFrame.Internal.Interpreter (interpret)+import DataFrame.Model (Fit (..), Predict (..))+import DataFrame.Operators ((.*.), (.+.), (.==.))++data AdaBoostConfig = AdaBoostConfig+    { abNEstimators :: !Int+    , abMaxDepth :: !Int+    }+    deriving (Eq, Show)++defaultAdaBoostConfig :: AdaBoostConfig+defaultAdaBoostConfig = AdaBoostConfig{abNEstimators = 50, abMaxDepth = 1}++-- | A fitted SAMME model: per-stage weights and stumps over the class set.+data AdaBoostModel a = AdaBoostModel+    { abAlphas :: !(VU.Vector Double)+    , abStumps :: !(V.Vector (Tree a))+    , abClasses :: !(V.Vector a)+    }+    deriving (Show)++instance (Columnable a, Ord a) => Fit AdaBoostConfig (Expr a) (AdaBoostModel a) where+    fit = fitAdaBoost++instance (Columnable a, Ord a) => Predict (AdaBoostModel a) a where+    predict = adaBoostExpr++-- | Fit an AdaBoost-SAMME classifier.+fitAdaBoost ::+    (Columnable a, Ord a) =>+    AdaBoostConfig -> Expr a -> DataFrame -> AdaBoostModel a+fitAdaBoost cfg target@(Col name) df =+    AdaBoostModel+        (VU.fromList (reverse alphas))+        (V.fromList (reverse stumps))+        classesV+  where+    feats = V.fromList (cartFeatures name df)+    ys = targetValues target df+    n = V.length ys+    classes = sort (foldr dedup [] (V.toList ys))+    dedup x acc = if x `elem` acc then acc else x : acc+    classesV = V.fromList classes+    kClasses = length classes+    codes = VU.generate n (\i -> classIndex (ys V.! i))+    classIndex v = length (takeWhile (< v) classes)+    (alphas, stumps) = boost 0 (VU.replicate n (1 / fromIntegral (max 1 n))) [] []+    boost !m w as ts+        | m >= abNEstimators cfg = (as, ts)+        | otherwise =+            let stump = fitWeightedTree (abMaxDepth cfg) feats classesV codes kClasses w+                pred = predictCodes df classesV stump+                wrong :: VU.Vector Int+                wrong = VU.generate n (\i -> if pred VU.! i /= codes VU.! i then 1 else 0)+                err = clamp (VU.sum (VU.zipWith (*) w (VU.map fromIntegral wrong)) / VU.sum w)+                alpha = log ((1 - err) / err) + log (fromIntegral (max 1 (kClasses - 1)))+                w' = normalize (VU.zipWith (\wi e -> wi * exp (alpha * fromIntegral e)) w wrong)+             in if err <= 0 || err >= 1 - 1 / fromIntegral kClasses+                    then (alpha : as, stump : ts)+                    else boost (m + 1) w' (alpha : as) (stump : ts)+    clamp e = max 1e-10 (min (1 - 1e-10) e)+    normalize v = let s = VU.sum v in if s == 0 then v else VU.map (/ s) v+fitAdaBoost _ expr _ =+    error ("fitAdaBoost: target must be a column, got " ++ show expr)++predictCodes ::+    forall a.+    (Columnable a, Ord a) =>+    DataFrame -> V.Vector a -> Tree a -> VU.Vector Int+predictCodes df classesV stump =+    VU.fromList (map toCode preds)+  where+    preds :: [a]+    preds = case interpret df (treeToExpr stump) of+        Right (TColumn c) -> either (const []) V.toList (toVector @a @V.Vector c)+        Left e -> error (show e)+    toCode v = fromMaybe 0 (V.findIndex (== v) classesV)++-- | A depth-bounded weighted classification tree (weighted Gini splits).+fitWeightedTree ::+    (Columnable a) =>+    Int ->+    V.Vector CartFeature ->+    V.Vector a ->+    VU.Vector Int ->+    Int ->+    VU.Vector Double ->+    Tree a+fitWeightedTree maxDepth feats classesV codes kClasses weights =+    go 0 (VU.enumFromN 0 (VU.length codes))+  where+    go depth idxs+        | depth >= maxDepth || VU.length idxs < 2 || isPure idxs =+            Leaf (classesV V.! majority idxs)+        | otherwise = case bestSplit idxs of+            Nothing -> Leaf (classesV V.! majority idxs)+            Just (fj, thr) ->+                let vals = cfValues (feats V.! fj)+                    (l, r) = VU.partition (\i -> vals VU.! i <= thr) idxs+                 in if VU.null l || VU.null r+                        then Leaf (classesV V.! majority idxs)+                        else+                            Branch+                                (cfPred (feats V.! fj) thr)+                                (go (depth + 1) l)+                                (go (depth + 1) r)+    classWeights idxs =+        VU.accumulate+            (+)+            (VU.replicate kClasses 0)+            (VU.map (\i -> (codes VU.! i, weights VU.! i)) idxs)+    majority idxs = VU.maxIndex (classWeights idxs)+    isPure idxs = VU.length (VU.filter (> 0) (classWeights idxs)) <= 1+    bestSplit idxs =+        let cands =+                [ (score, fj, thr)+                | fj <- [0 .. V.length feats - 1]+                , (thr, score) <- featureSplits idxs fj+                ]+         in case cands of+                [] -> Nothing+                _ -> let (_, fj, thr) = minimum3 cands in Just (fj, thr)+    featureSplits idxs fj =+        let vals = cfValues (feats V.! fj)+            member =+                VU.replicate (VU.length codes) False+                    VU.// [(i, True) | i <- VU.toList idxs]+            sorted = VU.filter (member VU.!) (sortIndicesByValue vals)+         in sweep vals sorted (classWeights idxs)+    sweep vals sorted totW = go0 0 (VU.replicate kClasses 0) Nothing+      where+        m = VU.length sorted+        totWsum = VU.sum totW+        go0 !k leftW best+            | k >= m - 1 = maybeToList best+            | otherwise =+                let i = sorted VU.! k+                    leftW' = leftW VU.// [(codes VU.! i, leftW VU.! (codes VU.! i) + weights VU.! i)]+                    vCur = vals VU.! i+                    vNext = vals VU.! (sorted VU.! (k + 1))+                    wl = VU.sum leftW'+                    wr = totWsum - wl+                    score = wl * gini leftW' + wr * gini (VU.zipWith (-) totW leftW')+                    valid = vCur /= vNext && wl > 0 && wr > 0+                    best' =+                        if valid && maybe True (\(_, s) -> score < s) best+                            then Just ((vCur + vNext) / 2, score)+                            else best+                 in go0 (k + 1) leftW' best'++gini :: VU.Vector Double -> Double+gini cw =+    let total = VU.sum cw+     in if total == 0+            then 0+            else 1 - VU.sum (VU.map (\c -> (c / total) ^ (2 :: Int)) cw)++minimum3 :: (Ord a) => [(a, b, c)] -> (a, b, c)+minimum3 = foldr1 (\x@(a, _, _) y@(b, _, _) -> if a <= b then x else y)++{- | Compile to an arg-max-of-weighted-votes expression over the class set:+@argmax_c Σ_m αₘ·[stumpₘ = c]@.+-}+adaBoostExpr :: (Columnable a, Ord a) => AdaBoostModel a -> Expr a+adaBoostExpr m = argMaxExpr (zip classes scores)+  where+    classes = V.toList (abClasses m)+    stumpExprs = map treeToExpr (V.toList (abStumps m))+    alphas = VU.toList (abAlphas m)+    scores =+        [ foldr (.+.) (F.lit 0) (zipWith (vote c) alphas stumpExprs)+        | c <- classes+        ]+    vote c a se = F.lit a .*. indicator (se .==. F.lit c)+    indicator cond = F.ifThenElse cond (F.lit 1.0) (F.lit 0.0)

@@ -0,0 +1,185 @@+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE OverloadedStrings #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE TypeApplications #-}++{- | Gradient boosting of regression trees (Friedman). Trees are fitted to the+negative gradient of the loss each round and accumulated with a shrinkage+factor; squared error gives regression, logistic deviance gives binary+classification. 'predict' is the additive score; 'gbProbaExpr' /+'gbDecisionExpr' give the classification probability / decision.+-}+module DataFrame.Boosting.GBM (+    GBLoss (..),+    GBConfig (..),+    defaultGBConfig,+    GBModel (..),+    gbExprAtStage,+    gbProbaExpr,+    gbDecisionExpr,+) where++import Data.Either (fromRight)+import qualified Data.Map.Strict as M+import qualified Data.Text as T+import qualified Data.Vector as V+import qualified Data.Vector.Unboxed as VU++import DataFrame.DecisionTree.Cart (cartFeatures)+import DataFrame.DecisionTree.Fit (treeToExpr)+import DataFrame.DecisionTree.Regression (RegTreeConfig (..), fitRegTreeOn)+import DataFrame.DecisionTree.Types (Tree)+import DataFrame.Featurize.Internal (targetDoubles)+import qualified DataFrame.Functions as F+import DataFrame.Internal.Column (TypedColumn (..), toVector)+import DataFrame.Internal.DataFrame (DataFrame)+import DataFrame.Internal.Expression (Expr (..), getColumns)+import DataFrame.Internal.Interpreter (interpret)+import DataFrame.Model (Fit (..), Predict (..))+import DataFrame.Operators ((.*.), (.+.), (.>.))++-- | The boosting loss.+data GBLoss = SquaredError | LogisticDeviance+    deriving (Eq, Show)++data GBConfig = GBConfig+    { gbLoss :: !GBLoss+    , gbNEstimators :: !Int+    , gbLearningRate :: !Double+    , gbMaxDepth :: !Int+    , gbSeed :: !Int+    }+    deriving (Eq, Show)++defaultGBConfig :: GBConfig+defaultGBConfig =+    GBConfig+        { gbLoss = SquaredError+        , gbNEstimators = 100+        , gbLearningRate = 0.1+        , gbMaxDepth = 3+        , gbSeed = 0+        }++{- | A fitted gradient-boosting model. 'gbInit' is the constant initial score+(mean, or log-odds for classification); 'gbTrees' are the staged regression+trees.+-}+data GBModel = GBModel+    { gbInit :: !Double+    , gbTrees :: !(V.Vector (Tree Double))+    , gbRate :: !Double+    , gbModelLoss :: !GBLoss+    , gbTrainScore :: !(VU.Vector Double)+    , gbFeatureUsage :: !(M.Map T.Text Int)+    }+    deriving (Show)++instance Fit GBConfig (Expr Double) GBModel where+    fit = fitGBM++instance Predict GBModel Double where+    predict = gbExpr++-- | Fit a gradient-boosting ensemble predicting @target@ from the other columns.+fitGBM :: GBConfig -> Expr Double -> DataFrame -> GBModel+fitGBM cfg target@(Col name) df =+    GBModel+        f0+        (V.fromList (reverse trees))+        lr+        (gbLoss cfg)+        (VU.fromList (reverse scores))+        usage+  where+    feats = V.fromList (cartFeatures name df)+    y = targetDoubles target df+    n = VU.length y+    lr = gbLearningRate cfg+    rtCfg =+        RegTreeConfig+            { rtMaxDepth = gbMaxDepth cfg+            , rtMinSamplesSplit = 2+            , rtMinLeafSize = 1+            , rtMinImpurityDecrease = 0.0+            }+    f0 = case gbLoss cfg of+        SquaredError -> VU.sum y / fromIntegral (max 1 n)+        LogisticDeviance ->+            let p = clamp01 (VU.sum y / fromIntegral (max 1 n))+             in log (p / (1 - p))+    (trees, scores, usage) = boost 0 (VU.replicate n f0) [] [] M.empty+    boost !m fScores ts ss usageAcc+        | m >= gbNEstimators cfg = (ts, ss, usageAcc)+        | otherwise =+            let grad = negGradient (gbLoss cfg) y fScores+                tree = fitRegTreeOn rtCfg feats grad Nothing+                pred = predictTree df tree+                fScores' = VU.zipWith (\f p -> f + lr * p) fScores pred+                score = lossValue (gbLoss cfg) y fScores'+                usage' = foldr (\c -> M.insertWith (+) c 1) usageAcc (treeColumns tree)+             in boost (m + 1) fScores' (tree : ts) (score : ss) usage'+fitGBM _ expr _ =+    error ("fitGBM: target must be a column, got " ++ show expr)++negGradient ::+    GBLoss -> VU.Vector Double -> VU.Vector Double -> VU.Vector Double+negGradient SquaredError y f = VU.zipWith (-) y f+negGradient LogisticDeviance y f =+    VU.zipWith (\yi fi -> yi - sigmoid fi) y f++lossValue :: GBLoss -> VU.Vector Double -> VU.Vector Double -> Double+lossValue SquaredError y f =+    VU.sum (VU.zipWith (\yi fi -> (yi - fi) ^ (2 :: Int)) y f)+        / fromIntegral (max 1 (VU.length y))+lossValue LogisticDeviance y f =+    VU.sum+        ( VU.zipWith+            ( \yi fi -> let p = clamp01 (sigmoid fi) in negate (yi * log p + (1 - yi) * log (1 - p))+            )+            y+            f+        )+        / fromIntegral (max 1 (VU.length y))++sigmoid :: Double -> Double+sigmoid z+    | z >= 0 = 1 / (1 + exp (-z))+    | otherwise = let e = exp z in e / (1 + e)++clamp01 :: Double -> Double+clamp01 p = max 1e-12 (min (1 - 1e-12) p)++predictTree :: DataFrame -> Tree Double -> VU.Vector Double+predictTree df t = case interpret @Double df (treeToExpr t) of+    Right (TColumn c) -> fromRight VU.empty (toVector @Double @VU.Vector c)+    Left e -> error (show e)++treeColumns :: Tree Double -> [T.Text]+treeColumns = getColumns . treeToExpr++-- | The full additive prediction expression: @f0 + lr · Σ treeᵢ@.+gbExpr :: GBModel -> Expr Double+gbExpr m = stageExpr (V.length (gbTrees m)) m++-- | The prediction expression using only the first @k@ trees (staged predict).+gbExprAtStage :: Int -> GBModel -> Maybe (Expr Double)+gbExprAtStage k m+    | k < 0 || k > V.length (gbTrees m) = Nothing+    | otherwise = Just (stageExpr k m)++stageExpr :: Int -> GBModel -> Expr Double+stageExpr k m =+    foldr ((.+.) . scaled) (F.lit (gbInit m)) (take k (V.toList (gbTrees m)))+  where+    scaled t = F.lit (gbRate m) .*. treeToExpr t++-- | Probability expression for classification: @sigmoid(score)@.+gbProbaExpr :: GBModel -> Expr Double+gbProbaExpr m = F.lit 1 / (F.lit 1 + exp (negate (gbExpr m)))++-- | Decision expression for classification: positive class when score > 0.+gbDecisionExpr :: GBModel -> Expr Bool+gbDecisionExpr m = gbExpr m .>. F.lit 0

@@ -0,0 +1,119 @@+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE OverloadedStrings #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE TypeApplications #-}++{- | Density-based clustering (DBSCAN). Brute-force @O(n²)@ region queries, no+spatial index — suitable for the in-memory scales this library targets. DBSCAN+is transductive: it has a 'Fit' instance but deliberately no 'Predict' instance+(there is no honest single prediction expression). 'dbscanSurrogateExpr' fits an+interpretable decision-tree surrogate on the cluster labels instead.+-}+module DataFrame.DBSCAN (+    DBSCANConfig (..),+    defaultDBSCANConfig,+    DBSCANModel (..),+    dbscanSurrogateExpr,+) where++import Control.Monad.ST (runST)+import qualified Data.Vector as V+import qualified Data.Vector.Unboxed as VU+import qualified Data.Vector.Unboxed.Mutable as VUM++import DataFrame.DecisionTree.Fit (fitDecisionTree)+import DataFrame.DecisionTree.Types (TreeConfig)+import DataFrame.Featurize.Internal (+    Features (..),+    columnExprName,+    extractFeatures,+    materializeColumn,+ )+import qualified DataFrame.Functions as F+import qualified DataFrame.Internal.Column as DI+import DataFrame.Internal.DataFrame (DataFrame, fromNamedColumns)+import DataFrame.Internal.Expression (Expr)+import DataFrame.LinearAlgebra (epsNeighbors)+import DataFrame.Model (Fit (..))++data DBSCANConfig = DBSCANConfig+    { dbEps :: !Double+    , dbMinSamples :: !Int+    }+    deriving (Eq, Show)++defaultDBSCANConfig :: DBSCANConfig+defaultDBSCANConfig = DBSCANConfig{dbEps = 0.5, dbMinSamples = 5}++{- | A fitted DBSCAN labelling. 'dbLabels' uses @-1@ for noise (sklearn's+@labels_@); 'dbCoreSampleIndices' are the core points.+-}+data DBSCANModel = DBSCANModel+    { dbLabels :: !(VU.Vector Int)+    , dbCoreSampleIndices :: !(VU.Vector Int)+    , dbNClusters :: !Int+    }+    deriving (Eq, Show)++instance Fit DBSCANConfig [Expr Double] DBSCANModel where+    fit = fitDBSCAN++-- | Cluster the feature columns with DBSCAN.+fitDBSCAN :: DBSCANConfig -> [Expr Double] -> DataFrame -> DBSCANModel+fitDBSCAN cfg features df =+    DBSCANModel labels coreIdx nClusters+  where+    Features _ _ rows n _ = extractFeatures features df+    nbrs = V.generate n (epsNeighbors (dbEps cfg) rows)+    isCore i = VU.length (nbrs V.! i) + 1 >= dbMinSamples cfg+    coreIdx = VU.fromList [i | i <- [0 .. n - 1], isCore i]+    labels = clusterLabels n nbrs isCore+    nClusters = if VU.null labels then 0 else 1 + maximum (-1 : VU.toList labels)++clusterLabels ::+    Int -> V.Vector (VU.Vector Int) -> (Int -> Bool) -> VU.Vector Int+clusterLabels n nbrs isCore = runST $ do+    lab <- VUM.replicate n (-2)+    let seedLoop c i+            | i >= n = pure ()+            | otherwise = do+                li <- VUM.read lab i+                if li /= -2+                    then seedLoop c (i + 1)+                    else+                        if not (isCore i)+                            then VUM.write lab i (-1) >> seedLoop c (i + 1)+                            else do+                                VUM.write lab i c+                                expand lab c (VU.toList (nbrs V.! i))+                                seedLoop (c + 1) (i + 1)+        expand _ _ [] = pure ()+        expand lab c (q : qs) = do+            lq <- VUM.read lab q+            if lq == -1+                then VUM.write lab q c >> expand lab c qs+                else+                    if lq /= -2+                        then expand lab c qs+                        else do+                            VUM.write lab q c+                            let extra = if isCore q then VU.toList (nbrs V.! q) else []+                            expand lab c (extra ++ qs)+    seedLoop 0 0+    VU.freeze lab++{- | Fit a decision-tree surrogate on the DBSCAN labels so new rows can be+assigned an (approximate) cluster. Noise (@-1@) is its own class.+-}+dbscanSurrogateExpr ::+    TreeConfig -> [Expr Double] -> DBSCANModel -> DataFrame -> Expr Int+dbscanSurrogateExpr cfg features model df =+    fitDecisionTree cfg (F.col @Int clusterCol) augmented+  where+    clusterCol = "__cluster__"+    cols = map (\e -> (columnExprName e, materializeColumn df e)) features+    augmented =+        fromNamedColumns $+            [(n, DI.fromList (VU.toList v)) | (n, v) <- cols]+                ++ [(clusterCol, DI.fromList (VU.toList (dbLabels model)))]

@@ -4,10 +4,11 @@ {-# LANGUAGE ScopedTypeVariables #-} {-# LANGUAGE TypeApplications #-} --- | sklearn-faithful CART initializer used to seed TAO. One-hot encodes--- categoricals, splits on exact (unsmoothed) Gini over midpoint thresholds--- (@<=@ routes left), and emits a @Tree@ predicting identically to--- @DecisionTreeClassifier(criterion='gini')@ on continuous features.+{- | sklearn-faithful CART initializer used to seed TAO. One-hot encodes+categoricals, splits on exact (unsmoothed) Gini over midpoint thresholds+(@<=@ routes left), and emits a @Tree@ predicting identically to+@DecisionTreeClassifier(criterion='gini')@ on continuous features.+-} module DataFrame.DecisionTree.Cart (     CartFeature (..),     CartNode (..),@@ -39,8 +40,9 @@ import qualified Data.Vector.Unboxed as VU import Type.Reflection (typeRep) --- | A one-hot feature column: per-row Double values plus the sklearn LEFT--- predicate (@x <= threshold@) over the ORIGINAL DataFrame.+{- | A one-hot feature column: per-row Double values plus the sklearn LEFT+predicate (@x <= threshold@) over the ORIGINAL DataFrame.+-} data CartFeature = CartFeature     { cfValues :: !(VU.Vector Double)     , cfPred :: !(Double -> Expr Bool)@@ -59,8 +61,9 @@     , ctxMinLeaf :: !Int     } --- | Indices @0..n-1@ stably sorted by their value (ascending), ties keeping--- ascending index. In-place unboxed merge sort — no boxed-list allocation.+{- | Indices @0..n-1@ stably sorted by their value (ascending), ties keeping+ascending index. In-place unboxed merge sort — no boxed-list allocation.+-} sortIndicesByValue :: VU.Vector Double -> VU.Vector Int sortIndicesByValue vs =     VU.create $ do@@ -68,7 +71,8 @@         VA.sortBy (compare `on` (vs VU.!)) mv         pure mv -buildCartTree :: forall a. (Columnable a, Ord a) => TreeConfig -> T.Text -> DataFrame -> Tree a+buildCartTree ::+    forall a. (Columnable a, Ord a) => TreeConfig -> T.Text -> DataFrame -> Tree a buildCartTree cfg target df =     cartToTree feats classes (buildCartNode ctx 0 (VU.enumFromN 0 nAll) featSorted)   where@@ -109,21 +113,35 @@  classCounts :: CartCtx -> VU.Vector Int -> VU.Vector Int classCounts ctx idxs =-    VU.accumulate (+) (VU.replicate (ctxNClasses ctx) 0) (VU.map (\i -> (ctxCodes ctx VU.! i, 1)) idxs)+    VU.accumulate+        (+)+        (VU.replicate (ctxNClasses ctx) 0)+        (VU.map (\i -> (ctxCodes ctx VU.! i, 1)) idxs)  isPure :: VU.Vector Int -> Bool isPure counts = VU.length (VU.filter (> 0) counts) <= 1 -buildCartNode :: CartCtx -> Int -> VU.Vector Int -> V.Vector (VU.Vector Int) -> CartNode+buildCartNode ::+    CartCtx -> Int -> VU.Vector Int -> V.Vector (VU.Vector Int) -> CartNode buildCartNode ctx depth idxs sortedByFeat     | VU.length idxs < 2 || depth >= ctxMaxDepth ctx || isPure counts = leaf-    | otherwise = maybe leaf (splitNode ctx depth idxs sortedByFeat) (bestSplit ctx sortedByFeat counts n)+    | otherwise =+        maybe+            leaf+            (splitNode ctx depth idxs sortedByFeat)+            (bestSplit ctx sortedByFeat counts n)   where     n = VU.length idxs     counts = classCounts ctx idxs     leaf = CLeaf (VU.maxIndex counts) -splitNode :: CartCtx -> Int -> VU.Vector Int -> V.Vector (VU.Vector Int) -> (Int, Double) -> CartNode+splitNode ::+    CartCtx ->+    Int ->+    VU.Vector Int ->+    V.Vector (VU.Vector Int) ->+    (Int, Double) ->+    CartNode splitNode ctx depth idxs sortedByFeat (fj, thr) =     CSplit fj thr (rec leftIdx leftSorted) (rec rightIdx rightSorted)   where@@ -134,9 +152,15 @@     rightSorted = V.map (VU.filter (\i -> vals VU.! i > thr)) sortedByFeat     rec = buildCartNode ctx (depth + 1) --- | Minimum weighted-child-Gini @(feature, threshold)@; the first feature wins--- ties; 'Nothing' when no feature has a leaf-size-respecting threshold.-bestSplit :: CartCtx -> V.Vector (VU.Vector Int) -> VU.Vector Int -> Int -> Maybe (Int, Double)+{- | Minimum weighted-child-Gini @(feature, threshold)@; the first feature wins+ties; 'Nothing' when no feature has a leaf-size-respecting threshold.+-}+bestSplit ::+    CartCtx ->+    V.Vector (VU.Vector Int) ->+    VU.Vector Int ->+    Int ->+    Maybe (Int, Double) bestSplit ctx sortedByFeat counts n =     fmap (\(_, j, t) -> (j, t)) (foldl' consider Nothing [0 .. ctxNFeats ctx - 1])   where@@ -145,8 +169,9 @@         Just (g, thr) | maybe True (\(gB, _, _) -> g < gB) acc -> Just (g, fj, thr)         _ -> acc --- | Accumulator while sweeping a feature's sorted rows: best @(gini, thr)@ so--- far, per-class left counts, rows moved left, and the previous value seen.+{- | Accumulator while sweeping a feature's sorted rows: best @(gini, thr)@ so+far, per-class left counts, rows moved left, and the previous value seen.+-} data Sweep = Sweep     { swBest :: !(Maybe (Double, Double))     , swLeft :: ![Int]@@ -154,9 +179,20 @@     , swPrev :: !Double     } -sweepFeature :: CartCtx -> [Int] -> VU.Vector Int -> CartFeature -> Int -> Maybe (Double, Double)+sweepFeature ::+    CartCtx ->+    [Int] ->+    VU.Vector Int ->+    CartFeature ->+    Int ->+    Maybe (Double, Double) sweepFeature ctx total si feat n =-    swBest (foldl' step (Sweep Nothing (replicate (ctxNClasses ctx) 0) 0 (0 / 0)) [0 .. VU.length si - 1])+    swBest+        ( foldl'+            step+            (Sweep Nothing (replicate (ctxNClasses ctx) 0) 0 (0 / 0))+            [0 .. VU.length si - 1]+        )   where     vals = cfValues feat     step s k = advance ctx total n (vals VU.! i) (ctxCodes ctx VU.! i) s@@ -165,24 +201,35 @@  advance :: CartCtx -> [Int] -> Int -> Double -> Int -> Sweep -> Sweep advance ctx total n v c s =-    Sweep (considerThreshold ctx total n v s) (bumpClass c (swLeft s)) (swMoved s + 1) v+    Sweep+        (considerThreshold ctx total n v s)+        (bumpClass c (swLeft s))+        (swMoved s + 1)+        v -considerThreshold :: CartCtx -> [Int] -> Int -> Double -> Sweep -> Maybe (Double, Double)+considerThreshold ::+    CartCtx -> [Int] -> Int -> Double -> Sweep -> Maybe (Double, Double) considerThreshold ctx total n v s     | swMoved s >= ctxMinLeaf ctx     , n - swMoved s >= ctxMinLeaf ctx     , v > swPrev s + 1e-7 =-        keepBetter (swBest s) (weightedGini total (swLeft s) (swMoved s) n) ((swPrev s + v) / 2)+        keepBetter+            (swBest s)+            (weightedGini total (swLeft s) (swMoved s) n)+            ((swPrev s + v) / 2)     | otherwise = swBest s -keepBetter :: Maybe (Double, Double) -> Double -> Double -> Maybe (Double, Double)+keepBetter ::+    Maybe (Double, Double) -> Double -> Double -> Maybe (Double, Double) keepBetter best g thr = case best of     Just (wb, _) | wb <= g -> best     _ -> Just (g, thr)  weightedGini :: [Int] -> [Int] -> Int -> Int -> Double weightedGini total leftAcc nl n =-    (fromIntegral nl * giniImpurity leftAcc nl + fromIntegral nr * giniImpurity rightAcc nr)+    ( fromIntegral nl * giniImpurity leftAcc nl+        + fromIntegral nr * giniImpurity rightAcc nr+    )         / fromIntegral n   where     nr = n - nl@@ -204,22 +251,31 @@ featuresOfColumn df c = case unsafeGetColumn c df of     UnboxedColumn _ (v :: VU.Vector b) -> numericFeature @b c v     BoxedColumn _ (v :: V.Vector b) -> oneHotFeatures @b (nRows df) c v+    pt@(PackedText _ _) -> case materializePacked pt of+        BoxedColumn _ (v :: V.Vector b) -> oneHotFeatures @b (nRows df) c v+        _ -> [] -numericFeature :: forall b. (Columnable b, VU.Unbox b) => T.Text -> VU.Vector b -> [CartFeature]+numericFeature ::+    forall b. (Columnable b, VU.Unbox b) => T.Text -> VU.Vector b -> [CartFeature] numericFeature c v = case testEquality (typeRep @b) (typeRep @Double) of     Just Refl -> [CartFeature v (\t -> F.col @Double c .<=. F.lit t)]     Nothing -> case sIntegral @b of-        STrue -> [CartFeature (VU.map fromIntegral v) (\t -> F.toDouble (F.col @b c) .<=. F.lit t)]+        STrue ->+            [ CartFeature (VU.map fromIntegral v) (\t -> F.toDouble (F.col @b c) .<=. F.lit t)+            ]         SFalse -> [] -oneHotFeatures :: forall b. (Columnable b) => Int -> T.Text -> V.Vector b -> [CartFeature]+oneHotFeatures ::+    forall b. (Columnable b) => Int -> T.Text -> V.Vector b -> [CartFeature] oneHotFeatures nAll c v = case testEquality (typeRep @b) (typeRep @T.Text) of     Just Refl -> [oneHot nAll c v cat | cat <- Set.toList (Set.fromList (V.toList v))]     Nothing -> []  oneHot :: Int -> T.Text -> V.Vector T.Text -> T.Text -> CartFeature oneHot nAll c v cat =-    CartFeature (VU.generate nAll (\i -> if v V.! i == cat then 1 else 0)) (const (F.col @T.Text c ./=. F.lit cat))+    CartFeature+        (VU.generate nAll (\i -> if v V.! i == cat then 1 else 0))+        (const (F.col @T.Text c ./=. F.lit cat))  -- | Target column as string labels (matches pandas @y.astype(str)@). cartTargetLabels :: T.Text -> DataFrame -> V.Vector T.Text@@ -228,3 +284,8 @@         Just Refl -> v         Nothing -> V.map (T.pack . show) v     UnboxedColumn _ (v :: VU.Vector b) -> V.map (T.pack . show) (V.convert v)+    pt@(PackedText _ _) -> case materializePacked pt of+        BoxedColumn _ (v :: V.Vector b) -> case testEquality (typeRep @b) (typeRep @T.Text) of+            Just Refl -> v+            Nothing -> V.map (T.pack . show) v+        _ -> V.empty

@@ -5,10 +5,11 @@ {-# LANGUAGE ScopedTypeVariables #-} {-# LANGUAGE TypeApplications #-} --- | Categorical split candidates: Breiman prefixes for binary targets,--- subset/singleton enumeration otherwise, and cross-column equality. Each--- value-list yields an OR-of-equalities condition (as an expression or a--- directly-read membership truth vector).+{- | Categorical split candidates: Breiman prefixes for binary targets,+subset/singleton enumeration otherwise, and cross-column equality. Each+value-list yields an OR-of-equalities condition (as an expression or a+directly-read membership truth vector).+-} module DataFrame.DecisionTree.Categorical (     TargetInfo (..),     mkTargetInfo,@@ -29,7 +30,12 @@ ) where  import DataFrame.DecisionTree.CondVec (CondVec (..), materializeCondVec)-import DataFrame.DecisionTree.Types (ColumnOrdering, SynthConfig (..), TreeConfig (..), withOrdFrom)+import DataFrame.DecisionTree.Types (+    ColumnOrdering,+    SynthConfig (..),+    TreeConfig (..),+    withOrdFrom,+ ) import DataFrame.Internal.Column import DataFrame.Internal.DataFrame (DataFrame, columnNames, unsafeGetColumn) import DataFrame.Internal.Expression (Expr (..))@@ -53,19 +59,25 @@ validBoxedValues :: Bitmap -> V.Vector a -> V.Vector a validBoxedValues bm = V.ifilter (\i _ -> bitmapTestBit bm i) --- | Target-column summary driving the categorical generator: binary vs--- multi-class, the deterministic positive class, and the raw label vector.+{- | Target-column summary driving the categorical generator: binary vs+multi-class, the deterministic positive class, and the raw label vector.+-} data TargetInfo target = TargetInfo     { tiIsBinary :: !Bool     , tiPositiveClass :: !(Maybe target)     , tiValues :: !(V.Vector target)     } --- | Compute 'TargetInfo' once per fit. The positive class for binary targets--- is the lexicographically-first distinct value, for deterministic pools.-mkTargetInfo :: forall target. (Columnable target, Ord target) => T.Text -> DataFrame -> Maybe (TargetInfo target)+{- | Compute 'TargetInfo' once per fit. The positive class for binary targets+is the lexicographically-first distinct value, for deterministic pools.+-}+mkTargetInfo ::+    forall target.+    (Columnable target, Ord target) =>+    T.Text -> DataFrame -> Maybe (TargetInfo target) mkTargetInfo target df = case interpret @target df (Col target) of-    Right (TColumn column) -> either (const Nothing) (Just . targetInfoFromValues) (toVector @target column)+    Right (TColumn column) ->+        either (const Nothing) (Just . targetInfoFromValues) (toVector @target column)     _ -> Nothing  targetInfoFromValues :: (Ord target) => V.Vector target -> TargetInfo target@@ -77,8 +89,9 @@         (p : _) | isBinary -> Just p         _ -> Nothing --- | Distinct values, capped: @Right vs@ (sorted) under the cap, else @Left@--- the count-so-far so the caller routes to the high-cardinality path.+{- | Distinct values, capped: @Right vs@ (sorted) under the cap, else @Left@+the count-so-far so the caller routes to the high-cardinality path.+-} distinctValuesUpTo :: (Ord a) => Int -> V.Vector a -> Either Int [a] distinctValuesUpTo cap values = go Set.empty 0   where@@ -88,8 +101,9 @@         | Set.size s > cap = Left (Set.size s)         | otherwise = go (Set.insert (V.unsafeIndex values i) s) (i + 1) --- | OR-of-equalities for a value-list, shared by the expression and--- truth-vector discrete paths so they stay byte-identical.+{- | OR-of-equalities for a value-list, shared by the expression and+truth-vector discrete paths so they stay byte-identical.+-} orEqs :: (a -> Expr Bool) -> [a] -> Expr Bool orEqs eqLit = foldr1 (.||.) . map eqLit @@ -106,17 +120,28 @@ singletonLists :: [a] -> [[a]] singletonLists = map (: []) -breimanPrefixSplits :: (Ord a, Ord target) => target -> V.Vector a -> V.Vector target -> [a] -> (a -> Expr Bool) -> [Expr Bool]+breimanPrefixSplits ::+    (Ord a, Ord target) =>+    target ->+    V.Vector a ->+    V.Vector target ->+    [a] ->+    (a -> Expr Bool) ->+    [Expr Bool] breimanPrefixSplits pc values targetVals distinctVals eqLit =     map (orEqs eqLit) (breimanPrefixLists pc values targetVals distinctVals) --- | Breiman's binary-target split set: sort levels by Laplace-smoothed--- positive rate, then take every contiguous non-trivial prefix.-breimanPrefixLists :: (Ord a, Ord target) => target -> V.Vector a -> V.Vector target -> [a] -> [[a]]+{- | Breiman's binary-target split set: sort levels by Laplace-smoothed+positive rate, then take every contiguous non-trivial prefix.+-}+breimanPrefixLists ::+    (Ord a, Ord target) => target -> V.Vector a -> V.Vector target -> [a] -> [[a]] breimanPrefixLists pc values targetVals distinctVals =     nonTrivialPrefixes (sortByRate (levelCounts pc values targetVals) distinctVals) -levelCounts :: (Ord a, Eq target) => target -> V.Vector a -> V.Vector target -> M.Map a (Int, Int)+levelCounts ::+    (Ord a, Eq target) =>+    target -> V.Vector a -> V.Vector target -> M.Map a (Int, Int) levelCounts pc values targetVals = V.ifoldl' add M.empty values   where     add acc i v = M.insertWith plus v (indicator (V.unsafeIndex targetVals i == pc), 1) acc@@ -132,17 +157,21 @@ sortByRate counts = sortBy (compare `on` (\v -> (laplaceRate counts v, v)))  nonTrivialPrefixes :: [a] -> [[a]]-nonTrivialPrefixes = tail . init . inits+nonTrivialPrefixes = drop 1 . init . inits --- | Value-lists a categorical column contributes; shared by the expression and--- truth-vector paths so both enumerate identical candidates in the same order.-catValueLists :: (Ord a, Ord target) => Bool -> Maybe target -> V.Vector target -> Int -> V.Vector a -> [[a]]+{- | Value-lists a categorical column contributes; shared by the expression and+truth-vector paths so both enumerate identical candidates in the same order.+-}+catValueLists ::+    (Ord a, Ord target) =>+    Bool -> Maybe target -> V.Vector target -> Int -> V.Vector a -> [[a]] catValueLists isBinary posClass targetVals subsetCap values     | V.null values = []     | isBinary, Just pc <- posClass = binaryLists pc targetVals values     | otherwise = multiclassLists subsetCap values -binaryLists :: (Ord a, Ord target) => target -> V.Vector target -> V.Vector a -> [[a]]+binaryLists ::+    (Ord a, Ord target) => target -> V.Vector target -> V.Vector a -> [[a]] binaryLists pc targetVals values     | length distinct < 2 = []     | otherwise = breimanPrefixLists pc values targetVals distinct@@ -158,15 +187,17 @@ ascDistinct :: (Ord a) => V.Vector a -> [a] ascDistinct = Set.toAscList . Set.fromList . V.toList --- | Truth vector of @col ∈ values@ read directly from the column; equal to--- interpreting @orEqs (== v) values@ because the values are distinct.+{- | Truth vector of @col ∈ values@ read directly from the column; equal to+interpreting @orEqs (== v) values@ because the values are distinct.+-} membershipVec :: (Ord a) => V.Vector a -> [a] -> VU.Vector Bool membershipVec colVals vs =     let !s = Set.fromList vs      in VU.generate (V.length colVals) (\i -> Set.member (colVals `V.unsafeIndex` i) s) --- | Per-fit categorical generation context bundling the target summary and--- the column-ordering registry.+{- | Per-fit categorical generation context bundling the target summary and+the column-ordering registry.+-} data CatCtx target = CatCtx     { ccBinary :: !Bool     , ccPos :: !(Maybe target)@@ -177,7 +208,12 @@  catCtx :: TargetInfo target -> TreeConfig -> CatCtx target catCtx ti cfg =-    CatCtx (tiIsBinary ti) (tiPositiveClass ti) (tiValues ti) (maxCategoricalSubsetCardinality (synthConfig cfg)) (columnOrdering cfg)+    CatCtx+        (tiIsBinary ti)+        (tiPositiveClass ti)+        (tiValues ti)+        (maxCategoricalSubsetCardinality (synthConfig cfg))+        (columnOrdering cfg)  catValueListsFor :: (Ord a, Ord target) => CatCtx target -> V.Vector a -> [[a]] catValueListsFor ctx = catValueLists (ccBinary ctx) (ccPos ctx) (ccTargets ctx) (ccSubsetCap ctx)@@ -190,26 +226,54 @@         STrue -> True         SFalse -> False --- | All equality-based candidate splits from non-numeric columns: per-column--- categorical conditions plus cross-column equality/order conditions.-discreteConditions :: forall target. (Columnable target, Ord target) => TargetInfo target -> TreeConfig -> DataFrame -> [Expr Bool]+{- | All equality-based candidate splits from non-numeric columns: per-column+categorical conditions plus cross-column equality/order conditions.+-}+discreteConditions ::+    forall target.+    (Columnable target, Ord target) =>+    TargetInfo target -> TreeConfig -> DataFrame -> [Expr Bool] discreteConditions targetInfo cfg df =-    concatMap (columnConds (catCtx targetInfo cfg) df) (columnNames df) ++ crossColumnConds cfg df+    concatMap (columnConds (catCtx targetInfo cfg) df) (columnNames df)+        ++ crossColumnConds cfg df -columnConds :: (Columnable target, Ord target) => CatCtx target -> DataFrame -> T.Text -> [Expr Bool]+columnConds ::+    (Columnable target, Ord target) =>+    CatCtx target -> DataFrame -> T.Text -> [Expr Bool] columnConds ctx df colName = case unsafeGetColumn colName df of     BoxedColumn Nothing (column :: V.Vector a) -> nonNullColConds ctx colName column     BoxedColumn (Just bm) (column :: V.Vector a) -> nullableColConds ctx colName bm column     UnboxedColumn _ (_ :: VU.Vector a) -> []+    pt@(PackedText _ _) -> case materializePacked pt of+        BoxedColumn Nothing (column :: V.Vector a) -> nonNullColConds ctx colName column+        BoxedColumn (Just bm) (column :: V.Vector a) -> nullableColConds ctx colName bm column+        _ -> [] -nonNullColConds :: forall a target. (Columnable a, Ord target) => CatCtx target -> T.Text -> V.Vector a -> [Expr Bool]+nonNullColConds ::+    forall a target.+    (Columnable a, Ord target) =>+    CatCtx target -> T.Text -> V.Vector a -> [Expr Bool] nonNullColConds ctx colName column =-    fromMaybe [] (withOrdFrom @a (ccOrds ctx) (map (orEqs (eqExprFor @a colName)) (catValueListsFor ctx column)))+    fromMaybe+        []+        ( withOrdFrom @a+            (ccOrds ctx)+            (map (orEqs (eqExprFor @a colName)) (catValueListsFor ctx column))+        ) -nullableColConds :: forall a target. (Columnable a, Ord target) => CatCtx target -> T.Text -> Bitmap -> V.Vector a -> [Expr Bool]+nullableColConds ::+    forall a target.+    (Columnable a, Ord target) =>+    CatCtx target -> T.Text -> Bitmap -> V.Vector a -> [Expr Bool] nullableColConds ctx colName bm column     | isNumericKind @a || V.null valid = []-    | otherwise = fromMaybe [] (withOrdFrom @a (ccOrds ctx) (map (orEqs (eqJustFor @a colName)) (catValueListsFor ctx valid)))+    | otherwise =+        fromMaybe+            []+            ( withOrdFrom @a+                (ccOrds ctx)+                (map (orEqs (eqJustFor @a colName)) (catValueListsFor ctx valid))+            )   where     valid = validBoxedValues bm column @@ -225,32 +289,50 @@  allowedPairs :: TreeConfig -> DataFrame -> [(T.Text, T.Text)] allowedPairs cfg df =-    [(l, r) | l <- columnNames df, r <- columnNames df, l /= r, not (isDisallowedPair cfg l r)]+    [ (l, r)+    | l <- columnNames df+    , r <- columnNames df+    , l /= r+    , not (isDisallowedPair cfg l r)+    ]  isDisallowedPair :: TreeConfig -> T.Text -> T.Text -> Bool isDisallowedPair cfg l r =-    any (\(l', r') -> sort [l', r'] == sort [l, r]) (disallowedCombinations (synthConfig cfg))+    any+        (\(l', r') -> sort [l', r'] == sort [l, r])+        (disallowedCombinations (synthConfig cfg))  pairConds :: ColumnOrdering -> DataFrame -> (T.Text, T.Text) -> [Expr Bool]-pairConds ords df (l, r) = case (unsafeGetColumn l df, unsafeGetColumn r df) of+pairConds ords df (l, r) = case ( materializePacked (unsafeGetColumn l df)+                                , materializePacked (unsafeGetColumn r df)+                                ) of     (BoxedColumn Nothing (_ :: V.Vector a), BoxedColumn Nothing (_ :: V.Vector b)) -> strictPairConds @a @b l r     (BoxedColumn (Just _) (_ :: V.Vector a), BoxedColumn (Just _) (_ :: V.Vector b)) -> nullablePairConds @a @b ords l r     _ -> [] -strictPairConds :: forall a b. (Columnable a, Columnable b) => T.Text -> T.Text -> [Expr Bool]+strictPairConds ::+    forall a b. (Columnable a, Columnable b) => T.Text -> T.Text -> [Expr Bool] strictPairConds l r = case testEquality (typeRep @a) (typeRep @b) of     Just Refl -> [Col @a l .==. Col @a r]     Nothing -> [] -nullablePairConds :: forall a b. (Columnable a, Columnable b) => ColumnOrdering -> T.Text -> T.Text -> [Expr Bool]+nullablePairConds ::+    forall a b.+    (Columnable a, Columnable b) =>+    ColumnOrdering -> T.Text -> T.Text -> [Expr Bool] nullablePairConds ords l r = case testEquality (typeRep @a) (typeRep @b) of     Nothing -> []     Just Refl -> nullableEqOrLe @a ords l r -nullableEqOrLe :: forall a. (Columnable a) => ColumnOrdering -> T.Text -> T.Text -> [Expr Bool]+nullableEqOrLe ::+    forall a. (Columnable a) => ColumnOrdering -> T.Text -> T.Text -> [Expr Bool] nullableEqOrLe ords l r     | isTextType @a = eqOnly-    | otherwise = maybe eqOnly (++ eqOnly) (withOrdFrom @a ords [Col @(Maybe a) l .<=. Col @(Maybe a) r])+    | otherwise =+        maybe+            eqOnly+            (++ eqOnly)+            (withOrdFrom @a ords [Col @(Maybe a) l .<=. Col @(Maybe a) r])   where     eqOnly = [Col @(Maybe a) l .==. Col @(Maybe a) r] @@ -259,23 +341,41 @@     Just Refl -> True     Nothing -> False --- | 'discreteConditions' materialized with shared per-column reads: the--- non-nullable categorical path builds truth vectors directly from one read--- per column; nullable and cross-column fall back to interpret.-discreteCondVecs :: forall target. (Columnable target, Ord target) => TargetInfo target -> TreeConfig -> DataFrame -> [CondVec]+{- | 'discreteConditions' materialized with shared per-column reads: the+non-nullable categorical path builds truth vectors directly from one read+per column; nullable and cross-column fall back to interpret.+-}+discreteCondVecs ::+    forall target.+    (Columnable target, Ord target) =>+    TargetInfo target -> TreeConfig -> DataFrame -> [CondVec] discreteCondVecs targetInfo cfg df =     concatMap (columnCondVecs (catCtx targetInfo cfg) df) (columnNames df)         ++ mapMaybe (materializeCondVec df) (crossColumnConds cfg df) -columnCondVecs :: (Columnable target, Ord target) => CatCtx target -> DataFrame -> T.Text -> [CondVec]+columnCondVecs ::+    (Columnable target, Ord target) =>+    CatCtx target -> DataFrame -> T.Text -> [CondVec] columnCondVecs ctx df colName = case unsafeGetColumn colName df of     BoxedColumn Nothing (column :: V.Vector a) -> nonNullColCondVecs ctx colName column     BoxedColumn (Just bm) (column :: V.Vector a) -> mapMaybe (materializeCondVec df) (nullableColConds ctx colName bm column)     UnboxedColumn _ (_ :: VU.Vector a) -> []+    pt@(PackedText _ _) -> case materializePacked pt of+        BoxedColumn Nothing (column :: V.Vector a) -> nonNullColCondVecs ctx colName column+        BoxedColumn (Just bm) (column :: V.Vector a) -> mapMaybe (materializeCondVec df) (nullableColConds ctx colName bm column)+        _ -> [] -nonNullColCondVecs :: forall a target. (Columnable a, Ord target) => CatCtx target -> T.Text -> V.Vector a -> [CondVec]+nonNullColCondVecs ::+    forall a target.+    (Columnable a, Ord target) => CatCtx target -> T.Text -> V.Vector a -> [CondVec] nonNullColCondVecs ctx colName column =-    fromMaybe [] (withOrdFrom @a (ccOrds ctx) (map (membershipCondVec colName column) (catValueListsFor ctx column)))+    fromMaybe+        []+        ( withOrdFrom @a+            (ccOrds ctx)+            (map (membershipCondVec colName column) (catValueListsFor ctx column))+        ) -membershipCondVec :: forall a. (Columnable a, Ord a) => T.Text -> V.Vector a -> [a] -> CondVec+membershipCondVec ::+    forall a. (Columnable a, Ord a) => T.Text -> V.Vector a -> [a] -> CondVec membershipCondVec colName column vs = CondVec (orEqs (eqExprFor @a colName) vs) (membershipVec column vs)

@@ -3,9 +3,10 @@ {-# LANGUAGE ScopedTypeVariables #-} {-# LANGUAGE TypeApplications #-} --- | Cached condition truth vectors and the per-fit cache keyed by structural--- form. A condition's truth over a fixed DataFrame is invariant for a whole--- fit, so it is materialized once and reused.+{- | Cached condition truth vectors and the per-fit cache keyed by structural+form. A condition's truth over a fixed DataFrame is invariant for a whole+fit, so it is materialized once and reused.+-} module DataFrame.DecisionTree.CondVec (     CondVec (..),     materializeCondVec,@@ -25,7 +26,12 @@ import qualified DataFrame.Functions as F import DataFrame.Internal.Column (TypedColumn (..), toVector) import DataFrame.Internal.DataFrame (DataFrame)-import DataFrame.Internal.Expression (BinaryOp (binaryName), Expr (..), eqExpr, normalize)+import DataFrame.Internal.Expression (+    BinaryOp (binaryName),+    Expr (..),+    eqExpr,+    normalize,+ ) import DataFrame.Internal.Interpreter (interpret)  import qualified Data.Map.Strict as M@@ -42,8 +48,9 @@     , cvVec :: !(VU.Vector Bool)     } --- | Interpret a condition once over the DataFrame; 'Nothing' on a--- type/interpret failure so the candidate is silently dropped.+{- | Interpret a condition once over the DataFrame; 'Nothing' on a+type/interpret failure so the candidate is silently dropped.+-} materializeCondVec :: DataFrame -> Expr Bool -> Maybe CondVec materializeCondVec df cond = case interpret @Bool df cond of     Left _ -> Nothing@@ -52,13 +59,15 @@ eitherToMaybe :: Either e a -> Maybe a eitherToMaybe = either (const Nothing) Just --- | Full-DataFrame truth vectors keyed by structural form, read-only once--- built. Seeded for free from the candidate pool plus the initial tree so the--- predict/loss passes index a vector instead of re-interpreting per node.+{- | Full-DataFrame truth vectors keyed by structural form, read-only once+built. Seeded for free from the candidate pool plus the initial tree so the+predict/loss passes index a vector instead of re-interpreting per node.+-} type CondCache = M.Map T.Text (VU.Vector Bool) --- | Structural key matching the candidate-dedup key, so a tree branch whose--- condition came from the pool hits the cache (equal keys ⟹ equal vector).+{- | Structural key matching the candidate-dedup key, so a tree branch whose+condition came from the pool hits the cache (equal keys ⟹ equal vector).+-} condCacheKey :: Expr Bool -> T.Text condCacheKey = T.pack . show . normalize @@ -66,8 +75,9 @@ condCacheFromVecs :: [CondVec] -> CondCache condCacheFromVecs cvs = M.fromList [(condCacheKey (cvExpr cv), cvVec cv) | cv <- cvs] --- | Add a tree's branch-condition vectors to a cache (one interpret per--- distinct, not-yet-cached condition).+{- | Add a tree's branch-condition vectors to a cache (one interpret per+distinct, not-yet-cached condition).+-} addTreeCondsToCache :: DataFrame -> Tree a -> CondCache -> CondCache addTreeCondsToCache df = go   where@@ -77,12 +87,14 @@ insertCond :: DataFrame -> Expr Bool -> CondCache -> CondCache insertCond df cond c     | M.member k c = c-    | otherwise = maybe c (\cv -> M.insert k (cvVec cv) c) (materializeCondVec df cond)+    | otherwise =+        maybe c (\cv -> M.insert k (cvVec cv) c) (materializeCondVec df cond)   where     k = condCacheKey cond --- | A condition's truth vector: a cache hit, else interpret over the--- DataFrame. 'Nothing' mirrors the interpret-failure fallback (route left).+{- | A condition's truth vector: a cache hit, else interpret over the+DataFrame. 'Nothing' mirrors the interpret-failure fallback (route left).+-} lookupCondVec :: CondCache -> DataFrame -> Expr Bool -> Maybe (VU.Vector Bool) lookupCondVec cache df cond = case M.lookup (condCacheKey cond) cache of     hit@(Just _) -> hit@@ -98,8 +110,9 @@   where     misrouted cp = (boolVals VU.! cpIndex cp) /= (cpCorrectDir cp == GoLeft) --- | A same-column same-direction Double threshold comparison, with a rebuild--- function to swap in a new threshold.+{- | A same-column same-direction Double threshold comparison, with a rebuild+function to swap in a new threshold.+-} data ThreshCmp = ThreshCmp     { tcCol :: !T.Text     , tcName :: !T.Text@@ -109,7 +122,9 @@  asDoubleThreshold :: Expr Bool -> Maybe ThreshCmp asDoubleThreshold (Binary op (Col c :: Expr cc) (Lit (t :: tt))) =-    case (testEquality (typeRep @cc) (typeRep @Double), testEquality (typeRep @tt) (typeRep @Double)) of+    case ( testEquality (typeRep @cc) (typeRep @Double)+         , testEquality (typeRep @tt) (typeRep @Double)+         ) of         (Just Refl, Just Refl) -> Just (ThreshCmp c (binaryName op) t (Binary op (Col c) . Lit))         _ -> Nothing asDoubleThreshold _ = Nothing@@ -117,8 +132,9 @@ directionalNames :: [T.Text] directionalNames = ["lt", "leq", "gt", "geq"] --- | Tighter (AND) or looser (OR) of two same-direction thresholds: @<@/@<=@--- are left-half-spaces (AND = min), @>@/@>=@ are right-half-spaces (AND = max).+{- | Tighter (AND) or looser (OR) of two same-direction thresholds: @<@/@<=@+are left-half-spaces (AND = min), @>@/@>=@ are right-half-spaces (AND = max).+-} chooseThreshold :: Bool -> T.Text -> Double -> Double -> Double chooseThreshold isAnd name t1 t2     | leftDir = if isAnd then min t1 t2 else max t1 t2@@ -126,8 +142,9 @@   where     leftDir = name == "lt" || name == "leq" --- | Collapse two same-column same-direction strict-Double comparisons into one--- comparison (the @True@ argument selects AND, @False@ OR); 'Nothing' otherwise.+{- | Collapse two same-column same-direction strict-Double comparisons into one+comparison (the @True@ argument selects AND, @False@ OR); 'Nothing' otherwise.+-} consolidateThreshold :: Bool -> Expr Bool -> Expr Bool -> Maybe (Expr Bool) consolidateThreshold isAnd ea eb = do     a <- asDoubleThreshold ea@@ -136,20 +153,28 @@         then Just (tcRebuild a (chooseThreshold isAnd (tcName a) (tcThr a) (tcThr b)))         else Nothing --- | AND-combine two cached conditions: idempotence and threshold consolidation--- first, else the generic @F.and@; the vector is always the elementwise AND.+{- | AND-combine two cached conditions: idempotence and threshold consolidation+first, else the generic @F.and@; the vector is always the elementwise AND.+-} combineAndVec :: CondVec -> CondVec -> CondVec combineAndVec a b     | eqExpr (cvExpr a) (cvExpr b) = a     | otherwise = CondVec expr (VU.zipWith (&&) (cvVec a) (cvVec b))   where-    expr = fromMaybe (F.and (cvExpr a) (cvExpr b)) (consolidateThreshold True (cvExpr a) (cvExpr b))+    expr =+        fromMaybe+            (F.and (cvExpr a) (cvExpr b))+            (consolidateThreshold True (cvExpr a) (cvExpr b)) --- | OR-combine two cached conditions (see 'combineAndVec'; AND/OR direction--- differs in 'consolidateThreshold').+{- | OR-combine two cached conditions (see 'combineAndVec'; AND/OR direction+differs in 'consolidateThreshold').+-} combineOrVec :: CondVec -> CondVec -> CondVec combineOrVec a b     | eqExpr (cvExpr a) (cvExpr b) = a     | otherwise = CondVec expr (VU.zipWith (||) (cvVec a) (cvVec b))   where-    expr = fromMaybe (F.or (cvExpr a) (cvExpr b)) (consolidateThreshold False (cvExpr a) (cvExpr b))+    expr =+        fromMaybe+            (F.or (cvExpr a) (cvExpr b))+            (consolidateThreshold False (cvExpr a) (cvExpr b))

@@ -3,9 +3,10 @@ {-# LANGUAGE ScopedTypeVariables #-} {-# LANGUAGE TypeApplications #-} --- | Top-level fitting: assemble the candidate pool, seed from CART, run TAO,--- and convert the result to an expression. Also the probability-tree variant--- ('fitProbTree') that annotates leaves with class distributions.+{- | Top-level fitting: assemble the candidate pool, seed from CART, run TAO,+and convert the result to an expression. Also the probability-tree variant+('fitProbTree') that annotates leaves with class distributions.+-} module DataFrame.DecisionTree.Fit (     treeToExpr,     fitDecisionTree,@@ -24,7 +25,12 @@ ) where  import DataFrame.DecisionTree.Cart (buildCartTree)-import DataFrame.DecisionTree.Categorical (TargetInfo (..), discreteConditions, discreteCondVecs, mkTargetInfo)+import DataFrame.DecisionTree.Categorical (+    TargetInfo (..),+    discreteCondVecs,+    discreteConditions,+    mkTargetInfo,+ ) import DataFrame.DecisionTree.CondVec (CondVec) import DataFrame.DecisionTree.Numeric (numericCondVecs, numericConditions) import DataFrame.DecisionTree.Pool (dedupCVByExpr, nubByExpr)@@ -54,9 +60,11 @@ treeToExpr (Branch cond left right) = F.ifThenElse cond (treeToExpr left) (treeToExpr right)  -- | Fit a TAO decision tree (CART-seeded) and return it as an expression.-fitDecisionTree :: forall a. (Columnable a, Ord a) => TreeConfig -> Expr a -> DataFrame -> Expr a+fitDecisionTree ::+    forall a. (Columnable a, Ord a) => TreeConfig -> Expr a -> DataFrame -> Expr a fitDecisionTree cfg (Col target) df =-    pruneExpr (treeToExpr (taoOptimizeCV @a cfg target condVecs df indices initialTree))+    pruneExpr+        (treeToExpr (taoOptimizeCV @a cfg target condVecs df indices initialTree))   where     condVecs = candidatePool @a cfg target df     initialTree = buildCartTree @a cfg target df@@ -64,18 +72,24 @@ fitDecisionTree _ expr _ = error ("Cannot create tree for compound expression: " ++ show expr)  -- | The deduplicated numeric + discrete candidate pool for a target column.-candidatePool :: forall a. (Columnable a, Ord a) => TreeConfig -> T.Text -> DataFrame -> [CondVec]+candidatePool ::+    forall a.+    (Columnable a, Ord a) => TreeConfig -> T.Text -> DataFrame -> [CondVec] candidatePool cfg target df = dedupCVByExpr (numericCVs ++ discreteCVs)   where     dfNoTarget = exclude [target] df     numericCVs = numericCondVecs cfg dfNoTarget df     discreteCVs = discreteCondVecs (targetInfoOrEmpty @a target df) cfg dfNoTarget -targetInfoOrEmpty :: forall a. (Columnable a, Ord a) => T.Text -> DataFrame -> TargetInfo a+targetInfoOrEmpty ::+    forall a. (Columnable a, Ord a) => T.Text -> DataFrame -> TargetInfo a targetInfoOrEmpty target df = fromMaybe (TargetInfo False Nothing V.empty) (mkTargetInfo @a target df)  -- | Fit a tree at a given depth from a raw condition list (CART + TAO + prune).-buildTree :: forall a. (Columnable a, Ord a) => TreeConfig -> Int -> T.Text -> [Expr Bool] -> DataFrame -> Expr a+buildTree ::+    forall a.+    (Columnable a, Ord a) =>+    TreeConfig -> Int -> T.Text -> [Expr Bool] -> DataFrame -> Expr a buildTree cfg depth target conds df =     pruneExpr (treeToExpr (taoOptimize @a cfg target conds df indices tree))   where@@ -89,7 +103,8 @@ partitionDataFrame cond df = (filterWhere cond df, filterWhere (F.not cond) df)  -- | Laplace-smoothed Gini impurity of the target distribution.-calculateGini :: forall a. (Columnable a, Ord a) => T.Text -> DataFrame -> Double+calculateGini ::+    forall a. (Columnable a, Ord a) => T.Text -> DataFrame -> Double calculateGini target df     | n == 0 = 0     | otherwise = 1 - sum (map (^ (2 :: Int)) probs)@@ -106,7 +121,8 @@   where     counts = getCounts @a target df -getCounts :: forall a. (Columnable a, Ord a) => T.Text -> DataFrame -> M.Map a Int+getCounts ::+    forall a. (Columnable a, Ord a) => T.Text -> DataFrame -> M.Map a Int getCounts target df = case interpret @a df (Col target) of     Left e -> throw e     Right (TColumn column) -> case toVector @a column of@@ -131,7 +147,9 @@ type ProbTree a = Tree (M.Map a Double)  -- | Normalised class probabilities over a subset of training rows.-probsFromIndices :: forall a. (Columnable a, Ord a) => T.Text -> DataFrame -> V.Vector Int -> M.Map a Double+probsFromIndices ::+    forall a.+    (Columnable a, Ord a) => T.Text -> DataFrame -> V.Vector Int -> M.Map a Double probsFromIndices target df indices = case interpret @a df (Col target) of     Right (TColumn column) -> either (const M.empty) (normaliseCounts indices) (toVector @a column)     _ -> M.empty@@ -139,30 +157,51 @@ normaliseCounts :: (Ord a) => V.Vector Int -> V.Vector a -> M.Map a Double normaliseCounts indices vals = M.map (\c -> fromIntegral c / total) counts   where-    counts = V.foldl' (\acc i -> M.insertWith (+) (vals V.! i) (1 :: Int) acc) M.empty indices+    counts =+        V.foldl'+            (\acc i -> M.insertWith (+) (vals V.! i) (1 :: Int) acc)+            M.empty+            indices     total = fromIntegral (V.length indices) :: Double --- | Re-label a fitted tree's leaves with class distributions, routing the--- training data through the (unchanged) split conditions.-buildProbTree :: forall a. (Columnable a, Ord a) => Tree a -> T.Text -> DataFrame -> V.Vector Int -> ProbTree a+{- | Re-label a fitted tree's leaves with class distributions, routing the+training data through the (unchanged) split conditions.+-}+buildProbTree ::+    forall a.+    (Columnable a, Ord a) =>+    Tree a -> T.Text -> DataFrame -> V.Vector Int -> ProbTree a buildProbTree (Leaf _) target df indices = Leaf (probsFromIndices @a target df indices) buildProbTree (Branch cond left right) target df indices =-    Branch cond (buildProbTree @a left target df l) (buildProbTree @a right target df r)+    Branch+        cond+        (buildProbTree @a left target df l)+        (buildProbTree @a right target df r)   where     (l, r) = partitionIndices cond df indices  -- | Fit a TAO tree and return one probability expression per class.-fitProbTree :: forall a. (Columnable a, Ord a) => TreeConfig -> Expr a -> DataFrame -> M.Map a (Expr Double)+fitProbTree ::+    forall a.+    (Columnable a, Ord a) =>+    TreeConfig -> Expr a -> DataFrame -> M.Map a (Expr Double) fitProbTree cfg (Col target) df = probExprs (buildProbTree @a pruned target df indices)   where-    conds = nubByExpr (numericConditions cfg dfNoTarget ++ discreteConditions (targetInfoOrEmpty @a target df) cfg dfNoTarget)+    conds =+        nubByExpr+            ( numericConditions cfg dfNoTarget+                ++ discreteConditions (targetInfoOrEmpty @a target df) cfg dfNoTarget+            )     dfNoTarget = exclude [target] df     indices = V.enumFromN 0 (nRows df)-    pruned = pruneDead (taoOptimize @a cfg target conds df indices (buildCartTree @a cfg target df))+    pruned =+        pruneDead+            (taoOptimize @a cfg target conds df indices (buildCartTree @a cfg target df)) fitProbTree _ expr _ = error ("Cannot create prob tree for compound expression: " ++ show expr)  -- | Convert a 'ProbTree' into one @Expr Double@ per class.-probExprs :: forall a. (Columnable a, Ord a) => ProbTree a -> M.Map a (Expr Double)+probExprs ::+    forall a. (Columnable a, Ord a) => ProbTree a -> M.Map a (Expr Double) probExprs tree = M.fromList [(c, classExpr c tree) | c <- nub (allClasses tree)]  allClasses :: ProbTree a -> [a]

@@ -1,9 +1,10 @@ {-# LANGUAGE OverloadedStrings #-} {-# LANGUAGE TypeApplications #-} --- | Oblique split candidates: fit an L1-regularised logistic hyperplane to the--- care points (class-balanced) and convert it to a boolean condition, rejecting--- all-zero and degenerate (single-side) hyperplanes.+{- | Oblique split candidates: fit an L1-regularised logistic hyperplane to the+care points (class-balanced) and convert it to a boolean condition, rejecting+all-zero and degenerate (single-side) hyperplanes.+-} module DataFrame.DecisionTree.Linear (     bestLinearCandidate,     fitLinearCandidate,@@ -15,7 +16,11 @@ ) where  import DataFrame.DecisionTree.Numeric (NumExpr (..), numericCols)-import DataFrame.DecisionTree.Types (CarePoint (..), Direction (..), TreeConfig (..))+import DataFrame.DecisionTree.Types (+    CarePoint (..),+    Direction (..),+    TreeConfig (..),+ ) import DataFrame.Internal.Column (TypedColumn (..), toVector) import DataFrame.Internal.DataFrame (DataFrame) import DataFrame.Internal.Expression (Expr, getColumns)@@ -27,18 +32,22 @@ import qualified Data.Vector as V import qualified Data.Vector.Unboxed as VU --- | Best oblique candidate, or 'Nothing' when the linear path is disabled or--- there are too few care points to fit on.-bestLinearCandidate :: TreeConfig -> DataFrame -> [CarePoint] -> Maybe (Expr Bool)+{- | Best oblique candidate, or 'Nothing' when the linear path is disabled or+there are too few care points to fit on.+-}+bestLinearCandidate ::+    TreeConfig -> DataFrame -> [CarePoint] -> Maybe (Expr Bool) bestLinearCandidate cfg df carePoints     | not (useLinearSolver cfg) = Nothing     | length carePoints < minCarePointsForLinear cfg = Nothing     | otherwise = fitLinearCandidate cfg df carePoints --- | Fit an L1 logistic regression to the care points and convert the resulting--- hyperplane to a condition, or 'Nothing' when no numeric features exist or the--- fitted model is all-zero or degenerate.-fitLinearCandidate :: TreeConfig -> DataFrame -> [CarePoint] -> Maybe (Expr Bool)+{- | Fit an L1 logistic regression to the care points and convert the resulting+hyperplane to a condition, or 'Nothing' when no numeric features exist or the+fitted model is all-zero or degenerate.+-}+fitLinearCandidate ::+    TreeConfig -> DataFrame -> [CarePoint] -> Maybe (Expr Bool) fitLinearCandidate cfg df carePoints = case materializedFeatures df carePoints of     [] -> Nothing     mats -> linearFromFeatures cfg carePoints mats@@ -46,7 +55,8 @@ materializedFeatures :: DataFrame -> [CarePoint] -> [(T.Text, VU.Vector Double)] materializedFeatures df carePoints = mapMaybe (materializeFeatureForCare df carePoints) (numericCols df) -linearFromFeatures :: TreeConfig -> [CarePoint] -> [(T.Text, VU.Vector Double)] -> Maybe (Expr Bool)+linearFromFeatures ::+    TreeConfig -> [CarePoint] -> [(T.Text, VU.Vector Double)] -> Maybe (Expr Bool) linearFromFeatures cfg carePoints mats     | VU.all (== 0) weights = Nothing     | degenerateHyperplane rows weights (LS.lmIntercept model) = Nothing@@ -54,14 +64,20 @@   where     rows = careRowsFromFeatures (length carePoints) mats     labels = careLabels carePoints-    model = LS.fitL1Logistic (solverConfigFor cfg labels) rows labels (V.fromList (map fst mats))+    model =+        LS.fitL1Logistic+            (solverConfigFor cfg labels)+            rows+            labels+            (V.fromList (map fst mats))     weights = LS.lmWeights model  solverConfigFor :: TreeConfig -> VU.Vector Double -> LS.SolverConfig solverConfigFor cfg labels = (linearSolverConfig cfg){LS.scSampleWeights = classBalancedWeights labels} --- | Class-balanced sklearn-form weights @w_i = N / (2 · N_class)@ (mean 1), or--- 'Nothing' in the degenerate one-class case (uniform weighting).+{- | Class-balanced sklearn-form weights @w_i = N / (2 · N_class)@ (mean 1), or+'Nothing' in the degenerate one-class case (uniform weighting).+-} classBalancedWeights :: VU.Vector Double -> Maybe (VU.Vector Double) classBalancedWeights labels     | nPos > 0 && nNeg > 0 = Just (VU.generate nCare weightAt)@@ -74,18 +90,25 @@         | VU.unsafeIndex labels i > 0 = fromIntegral nCare / (2 * fromIntegral nPos)         | otherwise = fromIntegral nCare / (2 * fromIntegral nNeg) --- | A hyperplane is degenerate when every care row scores on the same side of--- zero (equivalent to an invalid split, caught upstream).-degenerateHyperplane :: V.Vector (VU.Vector Double) -> VU.Vector Double -> Double -> Bool+{- | A hyperplane is degenerate when every care row scores on the same side of+zero (equivalent to an invalid split, caught upstream).+-}+degenerateHyperplane ::+    V.Vector (VU.Vector Double) -> VU.Vector Double -> Double -> Bool degenerateHyperplane rows weights bias =     nCare > 0 && (VU.minimum scores > 0 || VU.maximum scores < 0)   where     nCare = V.length rows-    scores = VU.generate nCare (\i -> VU.sum (VU.zipWith (*) weights (V.unsafeIndex rows i)) + bias)+    scores =+        VU.generate+            nCare+            (\i -> VU.sum (VU.zipWith (*) weights (V.unsafeIndex rows i)) + bias) --- | Per-care-point feature rows from materialized columns (each of length--- @nCare@, so indexing is in range).-careRowsFromFeatures :: Int -> [(T.Text, VU.Vector Double)] -> V.Vector (VU.Vector Double)+{- | Per-care-point feature rows from materialized columns (each of length+@nCare@, so indexing is in range).+-}+careRowsFromFeatures ::+    Int -> [(T.Text, VU.Vector Double)] -> V.Vector (VU.Vector Double) careRowsFromFeatures nCare mats =     V.generate nCare (\i -> VU.generate nFeat (\j -> snd (matsVec V.! j) VU.! i))   where@@ -94,7 +117,8 @@  -- | Solver labels: @+1@ when 'GoLeft' is correct, @-1@ otherwise. careLabels :: [CarePoint] -> VU.Vector Double-careLabels carePoints = VU.fromList [if cpCorrectDir cp == GoLeft then 1.0 else -1.0 | cp <- carePoints]+careLabels carePoints =+    VU.fromList [if cpCorrectDir cp == GoLeft then 1.0 else -1.0 | cp <- carePoints]  -- | First column referenced by an expression, or a placeholder when none. featName :: Expr b -> T.Text@@ -102,8 +126,9 @@     (c : _) -> c     [] -> "<feat>" --- | Replace missing values with the mean of present ones; 'Nothing' when--- nothing is present so the caller can drop the feature.+{- | Replace missing values with the mean of present ones; 'Nothing' when+nothing is present so the caller can drop the feature.+-} imputeMean :: [Maybe Double] -> Maybe (VU.Vector Double) imputeMean careRaw = case catMaybes careRaw of     [] -> Nothing@@ -116,14 +141,17 @@     Right (TColumn column) -> either (const Nothing) Just (toVector @Double column)     _ -> Nothing -interpretMaybeDoubleVals :: DataFrame -> Expr (Maybe Double) -> Maybe (V.Vector (Maybe Double))+interpretMaybeDoubleVals ::+    DataFrame -> Expr (Maybe Double) -> Maybe (V.Vector (Maybe Double)) interpretMaybeDoubleVals df expr = case interpret @(Maybe Double) df expr of     Right (TColumn column) -> either (const Nothing) Just (toVector @(Maybe Double) column)     _ -> Nothing --- | Materialize a 'NumExpr' over the care rows; 'Nothing' on interpret failure--- or (nullable) when no care point has a present value, else mean-imputed.-materializeFeatureForCare :: DataFrame -> [CarePoint] -> NumExpr -> Maybe (T.Text, VU.Vector Double)+{- | Materialize a 'NumExpr' over the care rows; 'Nothing' on interpret failure+or (nullable) when no care point has a present value, else mean-imputed.+-}+materializeFeatureForCare ::+    DataFrame -> [CarePoint] -> NumExpr -> Maybe (T.Text, VU.Vector Double) materializeFeatureForCare df carePoints (NDouble expr) = do     vals <- interpretDoubleVals df expr     Just (featName expr, VU.fromList [vals V.! cpIndex cp | cp <- carePoints])

@@ -0,0 +1,101 @@+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE UndecidableInstances #-}++{- | sklearn-style standalone tree estimators returning inspectable records+(depth, leaf count, per-feature split usage). 'fit' trains a classifier (from a+'TreeConfig') or a regressor (from a 'RegTreeConfig'); 'predict' is the compiled+tree expression, and the record exposes the raw 'Tree' too. The bare+'DataFrame.DecisionTree.Fit.fitDecisionTree' remains for callers that only want+the classifier @Expr@.+-}+module DataFrame.DecisionTree.Model (+    DecisionTreeClassifier (..),+    DecisionTreeRegressor (..),+) where++import qualified Data.Map.Strict as M+import qualified Data.Text as T++import qualified Data.Vector as V++import DataFrame.DecisionTree.Cart (cartFeatures)+import DataFrame.DecisionTree.Fit (fitDecisionTree, treeToExpr)+import DataFrame.DecisionTree.Regression (RegTreeConfig, fitRegTreeOn)+import DataFrame.DecisionTree.Types (Tree (..), TreeConfig)+import DataFrame.Featurize.Internal (targetDoubles)+import DataFrame.Internal.Column (Columnable)+import DataFrame.Internal.Expression (Expr (..), getColumns)+import DataFrame.Model (Fit (..), Predict (..))++-- | A fitted classification tree with structural diagnostics.+data DecisionTreeClassifier a = DecisionTreeClassifier+    { dtcExpr :: !(Expr a)+    , dtcDepth :: !Int+    , dtcNLeaves :: !Int+    , dtcFeatureUsage :: !(M.Map T.Text Int)+    }+    deriving (Show)++-- | A fitted regression tree with structural diagnostics.+data DecisionTreeRegressor = DecisionTreeRegressor+    { dtrTree :: !(Tree Double)+    , dtrExpr :: !(Expr Double)+    , dtrDepth :: !Int+    , dtrNLeaves :: !Int+    , dtrFeatureUsage :: !(M.Map T.Text Int)+    }+    deriving (Show)++instance (Columnable a, Ord a) => Fit TreeConfig (Expr a) (DecisionTreeClassifier a) where+    fit cfg target df =+        DecisionTreeClassifier+            e+            (exprDepth e)+            (exprLeaves e)+            (usageCounts (exprUsage e))+      where+        e = fitDecisionTree cfg target df++instance Predict (DecisionTreeClassifier a) a where+    predict = dtcExpr++instance Fit RegTreeConfig (Expr Double) DecisionTreeRegressor where+    fit cfg target df =+        DecisionTreeRegressor+            t+            e+            (exprDepth e)+            (exprLeaves e)+            (usageCounts (exprUsage e))+      where+        t = case target of+            Col name ->+                fitRegTreeOn+                    cfg+                    (V.fromList (cartFeatures name df))+                    (targetDoubles target df)+                    Nothing+            _ ->+                error+                    ("fit @DecisionTreeRegressor: target must be a column, got " ++ show target)+        e = treeToExpr t++instance Predict DecisionTreeRegressor Double where+    predict = dtrExpr++usageCounts :: [T.Text] -> M.Map T.Text Int+usageCounts = foldr (\c -> M.insertWith (+) c 1) M.empty++exprUsage :: Expr a -> [T.Text]+exprUsage (If c t e) = getColumns c ++ exprUsage t ++ exprUsage e+exprUsage _ = []++exprLeaves :: Expr a -> Int+exprLeaves (If _ t e) = exprLeaves t + exprLeaves e+exprLeaves _ = 1++exprDepth :: Expr a -> Int+exprDepth (If _ t e) = 1 + max (exprDepth t) (exprDepth e)+exprDepth _ = 0

@@ -5,10 +5,11 @@ {-# LANGUAGE ScopedTypeVariables #-} {-# LANGUAGE TypeApplications #-} --- | Numeric split candidates: per-column Double expressions, arithmetic--- expansion, and threshold conditions. 'numericCondVecs' materializes the--- pool with a single interpret per distinct expression, deriving every--- threshold/operator truth vector by direct comparison.+{- | Numeric split candidates: per-column Double expressions, arithmetic+expansion, and threshold conditions. 'numericCondVecs' materializes the+pool with a single interpret per distinct expression, deriving every+threshold/operator truth vector by direct comparison.+-} module DataFrame.DecisionTree.Numeric (     NumExpr (..),     numExprCols,@@ -67,11 +68,29 @@ combineNumExprs (NDouble e1) (NDouble e2) =     map NDouble [e1 .+ e2, e1 .- e2, e1 .* e2, safeDivD e1 e2] combineNumExprs (NDouble e1) (NMaybeDouble e2) =-    map NMaybeDouble [e1 .+ e2, e1 .- e2, e1 .* e2, safeDivMaybe (F.fromMaybe False (e2 ./= F.lit (0 :: Double))) (e1 ./ e2)]+    map+        NMaybeDouble+        [ e1 .+ e2+        , e1 .- e2+        , e1 .* e2+        , safeDivMaybe (F.fromMaybe False (e2 ./= F.lit (0 :: Double))) (e1 ./ e2)+        ] combineNumExprs (NMaybeDouble e1) (NDouble e2) =-    map NMaybeDouble [e1 .+ e2, e1 .- e2, e1 .* e2, safeDivMaybe (e2 ./= F.lit (0 :: Double)) (e1 ./ e2)]+    map+        NMaybeDouble+        [ e1 .+ e2+        , e1 .- e2+        , e1 .* e2+        , safeDivMaybe (e2 ./= F.lit (0 :: Double)) (e1 ./ e2)+        ] combineNumExprs (NMaybeDouble e1) (NMaybeDouble e2) =-    map NMaybeDouble [e1 .+ e2, e1 .- e2, e1 .* e2, safeDivMaybe (F.fromMaybe False (e2 ./= F.lit (0 :: Double))) (e1 ./ e2)]+    map+        NMaybeDouble+        [ e1 .+ e2+        , e1 .- e2+        , e1 .* e2+        , safeDivMaybe (F.fromMaybe False (e2 ./= F.lit (0 :: Double))) (e1 ./ e2)+        ]  numericConditions :: TreeConfig -> DataFrame -> [Expr Bool] numericConditions = generateNumericConds@@ -92,10 +111,14 @@  condsFromExpr :: NumExpr -> Double -> [Expr Bool] condsFromExpr (NDouble e) t = [e .<= F.lit t, e .>= F.lit t, e .< F.lit t, e .> F.lit t]-condsFromExpr (NMaybeDouble e) t = map (F.fromMaybe False) [e .<= F.lit t, e .>= F.lit t, e .< F.lit t, e .> F.lit t]+condsFromExpr (NMaybeDouble e) t =+    map+        (F.fromMaybe False)+        [e .<= F.lit t, e .>= F.lit t, e .< F.lit t, e .> F.lit t] --- | Percentile thresholds for a value list: sort once, index each percentile.--- Shared by 'generateNumericConds' and 'numericCondVecs' for identical results.+{- | Percentile thresholds for a value list: sort once, index each percentile.+Shared by 'generateNumericConds' and 'numericCondVecs' for identical results.+-} percentilesOf :: [Int] -> [Double] -> [Double] percentilesOf ps valsList     | n == 0 = []@@ -109,15 +132,17 @@     Right (TColumn column) -> either (const Nothing) Just (toVector @Double column)     _ -> Nothing -interpretMaybeDoubleCol :: DataFrame -> Expr (Maybe Double) -> Maybe (V.Vector (Maybe Double))+interpretMaybeDoubleCol ::+    DataFrame -> Expr (Maybe Double) -> Maybe (V.Vector (Maybe Double)) interpretMaybeDoubleCol df e = case interpret @(Maybe Double) df e of     Right (TColumn column) -> either (const Nothing) Just (toVector @(Maybe Double) column)     _ -> Nothing --- | Materialize the numeric pool with one interpret per distinct expression,--- deriving each threshold/operator truth vector by direct comparison.--- Byte-identical to materializing 'numericConditions' one at a time, but--- avoids re-interpreting each LHS per threshold and operator.+{- | Materialize the numeric pool with one interpret per distinct expression,+deriving each threshold/operator truth vector by direct comparison.+Byte-identical to materializing 'numericConditions' one at a time, but+avoids re-interpreting each LHS per threshold and operator.+-} numericCondVecs :: TreeConfig -> DataFrame -> DataFrame -> [CondVec] numericCondVecs cfg dfGen df = concatMap forExpr (numericExprsWithTerms (synthConfig cfg) dfGen)   where@@ -139,12 +164,14 @@   where     gen p = VU.generate n (\i -> p (vals V.! i)) -condsForMaybe :: TreeConfig -> Expr (Maybe Double) -> V.Vector (Maybe Double) -> [CondVec]+condsForMaybe ::+    TreeConfig -> Expr (Maybe Double) -> V.Vector (Maybe Double) -> [CondVec] condsForMaybe cfg e mvals = concatMap (maybeCondsAt e mvals (V.length mvals)) ts   where     ts = percentilesOf (percentiles cfg) (map (fromMaybe 0) (V.toList mvals)) -maybeCondsAt :: Expr (Maybe Double) -> V.Vector (Maybe Double) -> Int -> Double -> [CondVec]+maybeCondsAt ::+    Expr (Maybe Double) -> V.Vector (Maybe Double) -> Int -> Double -> [CondVec] maybeCondsAt e mvals n t =     [ CondVec (F.fromMaybe False (e .<= F.lit t)) (gen (<= t))     , CondVec (F.fromMaybe False (e .>= F.lit t)) (gen (>= t))@@ -154,13 +181,15 @@   where     gen p = VU.generate n (\i -> maybe False p (mvals V.! i)) --- | Arithmetic candidate expansion, generated already-deduped: each round--- combines @frontier × base@ and admits only normalized-novel candidates.--- Produces @base@ plus @maxExprDepth-1@ combination rounds.+{- | Arithmetic candidate expansion, generated already-deduped: each round+combines @frontier × base@ and admits only normalized-novel candidates.+Produces @base@ plus @maxExprDepth-1@ combination rounds.+-} numericExprsWithTerms :: SynthConfig -> DataFrame -> [NumExpr] numericExprsWithTerms cfg df     | not (enableArithOps cfg) = base-    | otherwise = base ++ expandRounds cfg base (max 0 (maxExprDepth cfg - 1)) base seen0+    | otherwise =+        base ++ expandRounds cfg base (max 0 (maxExprDepth cfg - 1)) base seen0   where     base = numericCols df     seen0 = Set.fromList (map keyNum base)@@ -177,9 +206,16 @@  roundProducts :: SynthConfig -> [NumExpr] -> [NumExpr] -> [NumExpr] roundProducts cfg frontier base =-    [c | e1 <- frontier, e2 <- base, not (numExprEq e1 e2), not (isDisallowed cfg e1 e2), c <- combineNumExprs e1 e2]+    [ c+    | e1 <- frontier+    , e2 <- base+    , not (numExprEq e1 e2)+    , not (isDisallowed cfg e1 e2)+    , c <- combineNumExprs e1 e2+    ] -expandRounds :: SynthConfig -> [NumExpr] -> Int -> [NumExpr] -> Set.Set String -> [NumExpr]+expandRounds ::+    SynthConfig -> [NumExpr] -> Int -> [NumExpr] -> Set.Set String -> [NumExpr] expandRounds _ _ 0 _ _ = [] expandRounds cfg base d frontier seen     | null admitted = []

@@ -1,9 +1,10 @@ {-# LANGUAGE BangPatterns #-} {-# LANGUAGE OverloadedStrings #-} --- | Candidate-pool scoring and boolean expansion: penalized scoring, diverse--- top-K selection, AND/OR saturation, and structural/truth-vector dedup. The--- per-node scoring scans run in parallel chunks.+{- | Candidate-pool scoring and boolean expansion: penalized scoring, diverse+top-K selection, AND/OR saturation, and structural/truth-vector dedup. The+per-node scoring scans run in parallel chunks.+-} module DataFrame.DecisionTree.Pool (     evalWithPenaltyVec,     primaryColExpr,@@ -21,9 +22,25 @@     nubByExpr, ) where -import DataFrame.DecisionTree.CondVec (CondVec (..), combineAndVec, combineOrVec, countErrorsByVec)-import DataFrame.DecisionTree.Types (CarePoint, SynthConfig (..), TreeConfig (..))-import DataFrame.Internal.Expression (Expr, compareExpr, eSize, eqExpr, getColumns, normalize)+import DataFrame.DecisionTree.CondVec (+    CondVec (..),+    combineAndVec,+    combineOrVec,+    countErrorsByVec,+ )+import DataFrame.DecisionTree.Types (+    CarePoint,+    SynthConfig (..),+    TreeConfig (..),+ )+import DataFrame.Internal.Expression (+    Expr,+    compareExpr,+    eSize,+    eqExpr,+    getColumns,+    normalize,+ )  import Control.Parallel.Strategies (parListChunk, rdeepseq, using) import Data.Function (on)@@ -33,16 +50,18 @@ import qualified Data.Text as T import qualified Data.Vector.Unboxed as VU --- | Penalized score of a candidate: care-point errors plus a complexity--- penalty, tie-broken by expression size.+{- | Penalized score of a candidate: care-point errors plus a complexity+penalty, tie-broken by expression size.+-} evalWithPenaltyVec :: TreeConfig -> [CarePoint] -> CondVec -> (Int, Int) evalWithPenaltyVec cfg carePoints cv = (countErrorsByVec (cvVec cv) carePoints + penalty, sz)   where     sz = eSize (cvExpr cv)     penalty = floor (complexityPenalty (synthConfig cfg) * fromIntegral sz) --- | First referenced column of a condition (a sentinel for literal-only ones),--- used by 'takeDiverse' to enforce per-column diversity.+{- | First referenced column of a condition (a sentinel for literal-only ones),+used by 'takeDiverse' to enforce per-column diversity.+-} primaryColExpr :: Expr Bool -> T.Text primaryColExpr e = case getColumns e of     [] -> "<noncol>"@@ -51,8 +70,9 @@ primaryColCV :: CondVec -> T.Text primaryColCV = primaryColExpr . cvExpr --- | Keep the first @k@ of an already-sorted list, admitting at most @quota@ per--- primary column (@Nothing@ disables the per-column cap).+{- | Keep the first @k@ of an already-sorted list, admitting at most @quota@ per+primary column (@Nothing@ disables the per-column cap).+-} takeDiverse :: Int -> Maybe Int -> (a -> T.Text) -> [a] -> [a] takeDiverse k Nothing _ = take k takeDiverse k (Just quota) primary = go M.empty 0@@ -65,36 +85,51 @@       where         !col = primary x --- | Chunk size for the parallel per-node candidate scans; tuned by an -N--- sweep, not correctness-affecting.+{- | Chunk size for the parallel per-node candidate scans; tuned by an -N+sweep, not correctness-affecting.+-} candidateParChunk :: Int candidateParChunk = 64 --- | Decorate candidates with their penalty in parallel chunks, forcing only--- the @(Int, Int)@ key so the order (hence later sorts/minima) is preserved.+{- | Decorate candidates with their penalty in parallel chunks, forcing only+the @(Int, Int)@ key so the order (hence later sorts/minima) is preserved.+-} decorate :: (CondVec -> (Int, Int)) -> [CondVec] -> [((Int, Int), CondVec)] decorate penaltyCV xs = zip (map penaltyCV xs `using` parListChunk candidateParChunk rdeepseq) xs  -- | The diverse top-@expressionPairs@ valid candidates by penalty. sortedTopK :: TreeConfig -> (CondVec -> (Int, Int)) -> [CondVec] -> [CondVec] sortedTopK cfg penaltyCV validCondVecs =-    map snd (takeDiverse (expressionPairs cfg) (perColumnQuota (synthConfig cfg)) (primaryColCV . snd) sorted)+    map+        snd+        ( takeDiverse+            (expressionPairs cfg)+            (perColumnQuota (synthConfig cfg))+            (primaryColCV . snd)+            sorted+        )   where     sorted = sortBy (compare `on` fst) (decorate penaltyCV validCondVecs)  -- | Lowest-penalty candidate after boolean saturation of the diverse top-K.-bestDiscreteCandidate :: TreeConfig -> (CondVec -> (Int, Int)) -> [CondVec] -> Maybe CondVec+bestDiscreteCandidate ::+    TreeConfig -> (CondVec -> (Int, Int)) -> [CondVec] -> Maybe CondVec bestDiscreteCandidate _ _ [] = Nothing bestDiscreteCandidate cfg penaltyCV validCondVecs =-    case saturateCandidates Structural (boolExpansion (synthConfig cfg)) (sortedTopK cfg penaltyCV validCondVecs) of+    case saturateCandidates+        Structural+        (boolExpansion (synthConfig cfg))+        (sortedTopK cfg penaltyCV validCondVecs) of         [] -> Nothing         xs -> Just (snd (minimumBy (compare `on` fst) (decorate penaltyCV xs))) --- | AND/OR expansion of cached conditions to depth @maxDepth@ (each--- combination is a single vector op, not an interpret).+{- | AND/OR expansion of cached conditions to depth @maxDepth@ (each+combination is a single vector op, not an interpret).+-} boolExprsVec :: [CondVec] -> [CondVec] -> Int -> Int -> [CondVec] boolExprsVec baseExprs prevExprs depth maxDepth-    | depth == 0 = baseExprs ++ boolExprsVec baseExprs prevExprs (depth + 1) maxDepth+    | depth == 0 =+        baseExprs ++ boolExprsVec baseExprs prevExprs (depth + 1) maxDepth     | depth >= maxDepth = []     | otherwise = combined ++ boolExprsVec baseExprs combined (depth + 1) maxDepth   where@@ -103,27 +138,38 @@ data DedupMode = Structural | TruthVector     deriving (Eq, Show) --- | Saturate the pool with AND/OR combinations, deduplicating structurally--- (byte-identical, first occurrence kept) or by truth vector (opt-in).+{- | Saturate the pool with AND/OR combinations, deduplicating structurally+(byte-identical, first occurrence kept) or by truth vector (opt-in).+-} saturateCandidates :: DedupMode -> Int -> [CondVec] -> [CondVec] saturateCandidates Structural maxDepth base = base' ++ go 1 base' seen0   where     (base', seen0) = admitKeys Set.empty base     go !depth frontier seen         | depth >= maxDepth || null frontier = []-        | otherwise = let (admitted, seen') = admitKeys seen (roundProducts frontier base) in admitted ++ go (depth + 1) admitted seen'+        | otherwise =+            let (admitted, seen') = admitKeys seen (roundProducts frontier base)+             in admitted ++ go (depth + 1) admitted seen' saturateCandidates TruthVector maxDepth base = M.elems (go 1 frontier0 reps0)   where     (reps0, frontier0) = admitVecs M.empty base     go !depth frontier reps         | depth >= maxDepth || null frontier = reps-        | otherwise = let (reps', admitted) = admitVecs reps (roundProducts frontier base) in go (depth + 1) admitted reps'+        | otherwise =+            let (reps', admitted) = admitVecs reps (roundProducts frontier base)+             in go (depth + 1) admitted reps' --- | One combination round: @frontier × base@ via AND then OR, skipping--- self-pairs (mirrors 'boolExprsVec' for byte-identical structural output).+{- | One combination round: @frontier × base@ via AND then OR, skipping+self-pairs (mirrors 'boolExprsVec' for byte-identical structural output).+-} roundProducts :: [CondVec] -> [CondVec] -> [CondVec] roundProducts frontier base =-    [c | e1 <- frontier, e2 <- base, not (eqExpr (cvExpr e1) (cvExpr e2)), c <- [combineAndVec e1 e2, combineOrVec e1 e2]]+    [ c+    | e1 <- frontier+    , e2 <- base+    , not (eqExpr (cvExpr e1) (cvExpr e2))+    , c <- [combineAndVec e1 e2, combineOrVec e1 e2]+    ]  -- | Admit candidates with a not-yet-seen normalized form, preserving order. admitKeys :: Set.Set String -> [CondVec] -> ([CondVec], Set.Set String)@@ -137,9 +183,13 @@ structuralKey :: CondVec -> String structuralKey = show . normalize . cvExpr --- | Admit candidates by distinct truth vector, keeping the smallest-expression--- representative per vector.-admitVecs :: M.Map (VU.Vector Bool) CondVec -> [CondVec] -> (M.Map (VU.Vector Bool) CondVec, [CondVec])+{- | Admit candidates by distinct truth vector, keeping the smallest-expression+representative per vector.+-}+admitVecs ::+    M.Map (VU.Vector Bool) CondVec ->+    [CondVec] ->+    (M.Map (VU.Vector Bool) CondVec, [CondVec]) admitVecs = go []   where     go acc reps [] = (reps, reverse acc)

@@ -2,9 +2,10 @@ {-# LANGUAGE ScopedTypeVariables #-} {-# LANGUAGE TypeApplications #-} --- | Prediction, care-point identification, node validity, and tree loss. The--- batched, cache-aware variants resolve each branch condition's truth vector--- once per call instead of once per row.+{- | Prediction, care-point identification, node validity, and tree loss. The+batched, cache-aware variants resolve each branch condition's truth vector+once per call instead of once per row.+-} module DataFrame.DecisionTree.Predict (     predictWithTree,     predictManyWithTree,@@ -20,8 +21,17 @@     isValidAtNode, ) where -import DataFrame.DecisionTree.CondVec (CondCache, countErrorsByVec, lookupCondVec)-import DataFrame.DecisionTree.Types (CarePoint (..), Direction (..), Tree (..), TreeConfig (..))+import DataFrame.DecisionTree.CondVec (+    CondCache,+    countErrorsByVec,+    lookupCondVec,+ )+import DataFrame.DecisionTree.Types (+    CarePoint (..),+    Direction (..),+    Tree (..),+    TreeConfig (..),+ ) import DataFrame.Internal.Column (Columnable, TypedColumn (..), toVector) import DataFrame.Internal.DataFrame (DataFrame) import DataFrame.Internal.Expression (Expr (..))@@ -37,21 +47,24 @@ import qualified Data.Vector.Mutable as VM import qualified Data.Vector.Unboxed as VU --- | A condition's truth vector over the DataFrame, or 'Nothing' on a--- type/interpret failure (callers default such rows to the left child).+{- | A condition's truth vector over the DataFrame, or 'Nothing' on a+type/interpret failure (callers default such rows to the left child).+-} branchBool :: DataFrame -> Expr Bool -> Maybe (VU.Vector Bool) branchBool df cond = case interpret @Bool df cond of     Right (TColumn column) -> either (const Nothing) Just (toVector @Bool @VU.Vector column)     _ -> Nothing  -- | The target column as a label vector, or 'Nothing' on failure.-interpretLabelCol :: forall a. (Columnable a) => DataFrame -> T.Text -> Maybe (V.Vector a)+interpretLabelCol ::+    forall a. (Columnable a) => DataFrame -> T.Text -> Maybe (V.Vector a) interpretLabelCol df target = case interpret @a df (Col target) of     Right (TColumn column) -> either (const Nothing) Just (toVector @a column)     _ -> Nothing  -- | Predict the label for a single row by walking a fixed tree (@True@ → left).-predictWithTree :: forall a. (Columnable a) => T.Text -> DataFrame -> Int -> Tree a -> a+predictWithTree ::+    forall a. (Columnable a) => T.Text -> DataFrame -> Int -> Tree a -> a predictWithTree _ _ _ (Leaf v) = v predictWithTree target df idx (Branch cond left right) =     predictWithTree @a target df idx (childFor cond left right idx df)@@ -61,12 +74,16 @@     Nothing -> left     Just boolVals -> if boolVals VU.! idx then left else right -predictManyWithTree :: forall a. (Columnable a) => Tree a -> DataFrame -> V.Vector Int -> V.Vector a+predictManyWithTree ::+    forall a. (Columnable a) => Tree a -> DataFrame -> V.Vector Int -> V.Vector a predictManyWithTree = predictManyWithTreeCached @a M.empty --- | 'predictManyWithTree' resolving each branch condition through a 'CondCache'.--- Each condition is read at most once per call rather than once per row.-predictManyWithTreeCached :: forall a. (Columnable a) => CondCache -> Tree a -> DataFrame -> V.Vector Int -> V.Vector a+{- | 'predictManyWithTree' resolving each branch condition through a 'CondCache'.+Each condition is read at most once per call rather than once per row.+-}+predictManyWithTreeCached ::+    forall a.+    (Columnable a) => CondCache -> Tree a -> DataFrame -> V.Vector Int -> V.Vector a predictManyWithTreeCached cache tree df indices = V.create $ do     mv <- VM.new (V.length indices)     fill mv (V.zip (V.enumFromN 0 (V.length indices)) indices) tree@@ -78,15 +95,33 @@         Nothing -> fill mv prs left         Just boolVals -> fillSplit mv (V.partition (\(_, i) -> boolVals VU.! i) prs) left right -    fillSplit :: VM.MVector s a -> (V.Vector (Int, Int), V.Vector (Int, Int)) -> Tree a -> Tree a -> ST s ()+    fillSplit ::+        VM.MVector s a ->+        (V.Vector (Int, Int), V.Vector (Int, Int)) ->+        Tree a ->+        Tree a ->+        ST s ()     fillSplit mv (leftPrs, rightPrs) left right = fill mv leftPrs left >> fill mv rightPrs right -identifyCarePoints :: forall a. (Columnable a) => T.Text -> DataFrame -> V.Vector Int -> Tree a -> Tree a -> [CarePoint]+identifyCarePoints ::+    forall a.+    (Columnable a) =>+    T.Text -> DataFrame -> V.Vector Int -> Tree a -> Tree a -> [CarePoint] identifyCarePoints = identifyCarePointsCached @a M.empty --- | Rows the parent must route to a specific child for the (fixed) subtrees to--- classify correctly; a 'CondCache' avoids re-interpreting subtree conditions.-identifyCarePointsCached :: forall a. (Columnable a) => CondCache -> T.Text -> DataFrame -> V.Vector Int -> Tree a -> Tree a -> [CarePoint]+{- | Rows the parent must route to a specific child for the (fixed) subtrees to+classify correctly; a 'CondCache' avoids re-interpreting subtree conditions.+-}+identifyCarePointsCached ::+    forall a.+    (Columnable a) =>+    CondCache ->+    T.Text ->+    DataFrame ->+    V.Vector Int ->+    Tree a ->+    Tree a ->+    [CarePoint] identifyCarePointsCached cache target df indices leftTree rightTree =     maybe [] carePoints (interpretLabelCol @a df target)   where@@ -94,7 +129,9 @@     rightPreds = predictManyWithTreeCached cache rightTree df indices     carePoints targetVals = V.toList (V.imapMaybe (checkPoint targetVals leftPreds rightPreds) indices) -checkPoint :: (Eq a) => V.Vector a -> V.Vector a -> V.Vector a -> Int -> Int -> Maybe CarePoint+checkPoint ::+    (Eq a) =>+    V.Vector a -> V.Vector a -> V.Vector a -> Int -> Int -> Maybe CarePoint checkPoint targetVals leftPreds rightPreds k idx =     case (leftPreds V.! k == trueLabel, rightPreds V.! k == trueLabel) of         (True, False) -> Just (CarePoint idx GoLeft)@@ -108,12 +145,19 @@ countCarePointErrors cond df carePoints =     maybe (length carePoints) (`countErrorsByVec` carePoints) (branchBool df cond) -partitionIndices :: Expr Bool -> DataFrame -> V.Vector Int -> (V.Vector Int, V.Vector Int)+partitionIndices ::+    Expr Bool -> DataFrame -> V.Vector Int -> (V.Vector Int, V.Vector Int) partitionIndices = partitionIndicesCached M.empty --- | 'partitionIndices' resolving the condition through a 'CondCache'; a miss--- routes every index left (matching the uncached fallback).-partitionIndicesCached :: CondCache -> Expr Bool -> DataFrame -> V.Vector Int -> (V.Vector Int, V.Vector Int)+{- | 'partitionIndices' resolving the condition through a 'CondCache'; a miss+routes every index left (matching the uncached fallback).+-}+partitionIndicesCached ::+    CondCache ->+    Expr Bool ->+    DataFrame ->+    V.Vector Int ->+    (V.Vector Int, V.Vector Int) partitionIndicesCached cache cond df indices = case lookupCondVec cache df cond of     Nothing -> (indices, V.empty)     Just boolVals -> V.partition (boolVals VU.!) indices@@ -125,7 +169,8 @@   where     (t, f) = partitionIndices c df indices -majorityValueFromIndices :: forall a. (Columnable a, Ord a) => T.Text -> DataFrame -> V.Vector Int -> a+majorityValueFromIndices ::+    forall a. (Columnable a, Ord a) => T.Text -> DataFrame -> V.Vector Int -> a majorityValueFromIndices target df indices = majorityOf (countLabels (labelColOrThrow @a df target) indices)  labelColOrThrow :: forall a. (Columnable a) => DataFrame -> T.Text -> V.Vector a@@ -141,21 +186,31 @@     | M.null counts = error "Empty indices in majorityValueFromIndices"     | otherwise = fst (maximumBy (compare `on` snd) (M.toList counts)) -computeTreeLoss :: forall a. (Columnable a) => T.Text -> DataFrame -> V.Vector Int -> Tree a -> Double+computeTreeLoss ::+    forall a.+    (Columnable a) => T.Text -> DataFrame -> V.Vector Int -> Tree a -> Double computeTreeLoss = computeTreeLossCached @a M.empty  -- | 0/1 loss of a tree over @indices@, with a 'CondCache' for the predictions.-computeTreeLossCached :: forall a. (Columnable a) => CondCache -> T.Text -> DataFrame -> V.Vector Int -> Tree a -> Double+computeTreeLossCached ::+    forall a.+    (Columnable a) =>+    CondCache -> T.Text -> DataFrame -> V.Vector Int -> Tree a -> Double computeTreeLossCached cache target df indices tree     | V.null indices = 0-    | otherwise = maybe 1.0 (treeLoss cache tree df indices) (interpretLabelCol @a df target)+    | otherwise =+        maybe 1.0 (treeLoss cache tree df indices) (interpretLabelCol @a df target) -treeLoss :: (Columnable a) => CondCache -> Tree a -> DataFrame -> V.Vector Int -> V.Vector a -> Double+treeLoss ::+    (Columnable a) =>+    CondCache -> Tree a -> DataFrame -> V.Vector Int -> V.Vector a -> Double treeLoss cache tree df indices targetVals =-    fromIntegral (countMismatches targetVals indices preds) / fromIntegral (V.length indices)+    fromIntegral (countMismatches targetVals indices preds)+        / fromIntegral (V.length indices)   where     preds = predictManyWithTreeCached cache tree df indices  countMismatches :: (Eq a) => V.Vector a -> V.Vector Int -> V.Vector a -> Int countMismatches targetVals indices preds =-    V.length (V.ifilter (\k _ -> targetVals V.! (indices V.! k) /= preds V.! k) preds)+    V.length+        (V.ifilter (\k _ -> targetVals V.! (indices V.! k) /= preds V.! k) preds)

@@ -2,9 +2,10 @@ {-# LANGUAGE GADTs #-} {-# LANGUAGE ScopedTypeVariables #-} --- | Post-convergence simplification of a fitted tree and its expression form:--- drop branches forced by path-condition entailment, collapse identical--- siblings, and fold redundant nested conditionals.+{- | Post-convergence simplification of a fitted tree and its expression form:+drop branches forced by path-condition entailment, collapse identical+siblings, and fold redundant nested conditionals.+-} module DataFrame.DecisionTree.Prune (     pruneDead,     treeEq,@@ -16,9 +17,10 @@ import DataFrame.Internal.Expression (Expr (..), eqExpr) import DataFrame.Internal.Simplify (PredFact, entails, factFalse, factTrue) --- | Drop branches whose test is forced by the path conditions reaching them,--- and collapse @Branch c t t@ to @t@. Sound for the decidable threshold subset;--- other tests are left untouched.+{- | Drop branches whose test is forced by the path conditions reaching them,+and collapse @Branch c t t@ to @t@. Sound for the decidable threshold subset;+other tests are left untouched.+-} pruneDead :: forall a. (Columnable a) => Tree a -> Tree a pruneDead = go []   where@@ -27,7 +29,11 @@     go facts (Branch cond left right) = case entails facts cond of         Just True -> go facts left         Just False -> go facts right-        Nothing -> reconcile cond (go (addFact (factTrue cond) facts) left) (go (addFact (factFalse cond) facts) right)+        Nothing ->+            reconcile+                cond+                (go (addFact (factTrue cond) facts) left)+                (go (addFact (factFalse cond) facts) right)  reconcile :: (Columnable a) => Expr Bool -> Tree a -> Tree a -> Tree a reconcile cond left right@@ -43,8 +49,9 @@ treeEq (Branch c1 l1 r1) (Branch c2 l2 r2) = eqExpr c1 c2 && treeEq l1 l2 && treeEq r1 r2 treeEq _ _ = False --- | Recursively fold @If@ expressions whose branches coincide or nest the same--- condition; leave other expressions structurally unchanged.+{- | Recursively fold @If@ expressions whose branches coincide or nest the same+condition; leave other expressions structurally unchanged.+-} pruneExpr :: forall a. (Columnable a) => Expr a -> Expr a pruneExpr (If cond t0 f0) = collapseIf cond (pruneExpr t0) (pruneExpr f0) pruneExpr (Unary op e) = Unary op (pruneExpr e)

@@ -0,0 +1,144 @@+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE ScopedTypeVariables #-}++{- | Variance-reduction (weighted-SSE) regression trees, reusing the CART+feature machinery. Leaves predict the (weighted) mean of their rows. The+matrix-level 'fitRegTreeOn' lets gradient boosting refit on residuals without+re-extracting features each round.+-}+module DataFrame.DecisionTree.Regression (+    RegTreeConfig (..),+    defaultRegTreeConfig,+    fitRegTreeOn,+) where++import Data.Function (on)+import Data.Maybe (maybeToList)+import qualified Data.Vector as V+import qualified Data.Vector.Algorithms.Merge as VA+import qualified Data.Vector.Unboxed as VU++import DataFrame.DecisionTree.Cart (CartFeature (..))+import DataFrame.DecisionTree.Types (Tree (..))++-- | Stopping criteria for the regression tree.+data RegTreeConfig = RegTreeConfig+    { rtMaxDepth :: !Int+    , rtMinSamplesSplit :: !Int+    , rtMinLeafSize :: !Int+    , rtMinImpurityDecrease :: !Double+    }+    deriving (Eq, Show)++defaultRegTreeConfig :: RegTreeConfig+defaultRegTreeConfig =+    RegTreeConfig+        { rtMaxDepth = 3+        , rtMinSamplesSplit = 2+        , rtMinLeafSize = 1+        , rtMinImpurityDecrease = 0.0+        }++{- | Fit on pre-extracted features, a target vector, and optional per-row+weights (length @n@). Used by gradient boosting on residual targets.+-}+fitRegTreeOn ::+    RegTreeConfig ->+    V.Vector CartFeature ->+    VU.Vector Double ->+    Maybe (VU.Vector Double) ->+    Tree Double+fitRegTreeOn cfg feats y mw = go 0 (VU.enumFromN 0 n)+  where+    n = VU.length y+    wt i = maybe 1 (VU.! i) mw+    go depth idxs+        | VU.length idxs < rtMinSamplesSplit cfg+            || depth >= rtMaxDepth cfg =+            Leaf (nodeMean idxs)+        | otherwise = case bestSplit idxs of+            Nothing -> Leaf (nodeMean idxs)+            Just (fj, thr) ->+                let vals = cfValues (feats V.! fj)+                    (lefts, rights) = VU.partition (\i -> vals VU.! i <= thr) idxs+                 in if VU.null lefts || VU.null rights+                        then Leaf (nodeMean idxs)+                        else+                            Branch+                                (cfPred (feats V.! fj) thr)+                                (go (depth + 1) lefts)+                                (go (depth + 1) rights)+    nodeMean idxs =+        let (sw, sy) = VU.foldl' (\(!a, !b) i -> (a + wt i, b + wt i * (y VU.! i))) (0, 0) idxs+         in if sw == 0 then 0 else sy / sw+    bestSplit idxs =+        let (totW, totSY, totSY2) = moments idxs+            nodeSSE = totSY2 - safeDiv (totSY * totSY) totW+            candidates =+                [ (red, fj, thr)+                | fj <- [0 .. V.length feats - 1]+                , (thr, red) <- featureSplits idxs fj totW totSY totSY2 nodeSSE+                ]+         in case candidates of+                [] -> Nothing+                _ ->+                    let (red, fj, thr) = maximumByFst candidates+                     in if red >= rtMinImpurityDecrease cfg && red > 0+                            then Just (fj, thr)+                            else Nothing+    featureSplits idxs fj totW totSY totSY2 nodeSSE =+        let vals = cfValues (feats V.! fj)+            sorted = sortByVal vals idxs+         in sweep sorted vals totW totSY totSY2 nodeSSE+    sweep sorted vals totW totSY totSY2 nodeSSE = go0 0 0 0 0 Nothing+      where+        m = VU.length sorted+        go0 !k !wl !syl !syl2 best+            | k >= m - 1 = maybeToList best+            | otherwise =+                let i = sorted VU.! k+                    wi = wt i+                    yi = y VU.! i+                    wl' = wl + wi+                    syl' = syl + wi * yi+                    syl2' = syl2 + wi * yi * yi+                    vCur = vals VU.! i+                    vNext = vals VU.! (sorted VU.! (k + 1))+                    nl = k + 1+                    nr = m - nl+                    wr = totW - wl'+                    valid =+                        vCur /= vNext+                            && nl >= rtMinLeafSize cfg+                            && nr >= rtMinLeafSize cfg+                            && wl' > 0+                            && wr > 0+                    red =+                        nodeSSE+                            - ( (syl2' - safeDiv (syl' * syl') wl')+                                    + ( (totSY2 - syl2')+                                            - safeDiv ((totSY - syl') * (totSY - syl')) wr+                                      )+                              )+                    thr = (vCur + vNext) / 2+                    best' =+                        if valid && maybe True (\(_, b) -> red > b) best+                            then Just (thr, red)+                            else best+                 in go0 (k + 1) wl' syl' syl2' best'+    moments =+        VU.foldl'+            ( \(!w, !sy, !sy2) i ->+                let wi = wt i; yi = y VU.! i+                 in (w + wi, sy + wi * yi, sy2 + wi * yi * yi)+            )+            (0, 0, 0)++safeDiv :: Double -> Double -> Double+safeDiv a b = if b == 0 then 0 else a / b++sortByVal :: VU.Vector Double -> VU.Vector Int -> VU.Vector Int+sortByVal vals = VU.modify (VA.sortBy (compare `on` (vals VU.!)))++maximumByFst :: (Ord a) => [(a, b, c)] -> (a, b, c)+maximumByFst = foldr1 (\x@(a, _, _) y@(b, _, _) -> if a >= b then x else y)

@@ -3,9 +3,10 @@ {-# LANGUAGE ScopedTypeVariables #-} {-# LANGUAGE TypeApplications #-} --- | Tree Alternating Optimization: hold the tree fixed and re-optimize one node--- at a time, bottom-up, minimizing care-point misroutes. Sibling subtrees at a--- depth level are independent and optimized in parallel.+{- | Tree Alternating Optimization: hold the tree fixed and re-optimize one node+at a time, bottom-up, minimizing care-point misroutes. Sibling subtrees at a+depth level are independent and optimized in parallel.+-} module DataFrame.DecisionTree.Tao (     taoOptimize,     taoOptimizeCV,@@ -17,7 +18,11 @@  import DataFrame.DecisionTree.CondVec import DataFrame.DecisionTree.Linear (bestLinearCandidate)-import DataFrame.DecisionTree.Pool (bestDiscreteCandidate, candidateParChunk, evalWithPenaltyVec)+import DataFrame.DecisionTree.Pool (+    bestDiscreteCandidate,+    candidateParChunk,+    evalWithPenaltyVec,+ ) import DataFrame.DecisionTree.Predict import DataFrame.DecisionTree.Prune (pruneDead) import DataFrame.DecisionTree.Types@@ -34,8 +39,9 @@ import qualified Data.Vector as V import qualified Data.Vector.Unboxed as VU --- | The constant per-fit context threaded through the node-optimization--- recursion (the cache is rebuilt each iteration).+{- | The constant per-fit context threaded through the node-optimization+recursion (the cache is rebuilt each iteration).+-} data TaoEnv = TaoEnv     { teCache :: !CondCache     , teCfg :: !TreeConfig@@ -45,13 +51,32 @@     }  -- | Public TAO entry point over raw conditions; materializes each once.-taoOptimize :: forall a. (Columnable a, Ord a) => TreeConfig -> T.Text -> [Expr Bool] -> DataFrame -> V.Vector Int -> Tree a -> Tree a+taoOptimize ::+    forall a.+    (Columnable a, Ord a) =>+    TreeConfig ->+    T.Text ->+    [Expr Bool] ->+    DataFrame ->+    V.Vector Int ->+    Tree a ->+    Tree a taoOptimize cfg target conds df =     taoOptimizeCV @a cfg target (mapMaybe (materializeCondVec df) conds) df --- | TAO outer loop over pre-evaluated candidates: iterate until the iteration--- budget or convergence tolerance is reached, then prune dead branches.-taoOptimizeCV :: forall a. (Columnable a, Ord a) => TreeConfig -> T.Text -> [CondVec] -> DataFrame -> V.Vector Int -> Tree a -> Tree a+{- | TAO outer loop over pre-evaluated candidates: iterate until the iteration+budget or convergence tolerance is reached, then prune dead branches.+-}+taoOptimizeCV ::+    forall a.+    (Columnable a, Ord a) =>+    TreeConfig ->+    T.Text ->+    [CondVec] ->+    DataFrame ->+    V.Vector Int ->+    Tree a ->+    Tree a taoOptimizeCV cfg target condVecs df rootIndices initialTree =     go 0 initialTree (lossWith baseCache initialTree)   where@@ -67,32 +92,63 @@         newLoss = lossWith cache tree'  -- | Public single-iteration entry point.-taoIteration :: forall a. (Columnable a, Ord a) => TreeConfig -> T.Text -> [Expr Bool] -> DataFrame -> V.Vector Int -> Tree a -> Tree a+taoIteration ::+    forall a.+    (Columnable a, Ord a) =>+    TreeConfig ->+    T.Text ->+    [Expr Bool] ->+    DataFrame ->+    V.Vector Int ->+    Tree a ->+    Tree a taoIteration cfg target conds df rootIndices tree =     let condVecs = mapMaybe (materializeCondVec df) conds         cache = addTreeCondsToCache df tree (condCacheFromVecs condVecs)      in taoIterationCV @a cache cfg target condVecs df rootIndices tree  -- | One bottom-to-top sweep: re-optimize every node level by level.-taoIterationCV :: forall a. (Columnable a, Ord a) => CondCache -> TreeConfig -> T.Text -> [CondVec] -> DataFrame -> V.Vector Int -> Tree a -> Tree a+taoIterationCV ::+    forall a.+    (Columnable a, Ord a) =>+    CondCache ->+    TreeConfig ->+    T.Text ->+    [CondVec] ->+    DataFrame ->+    V.Vector Int ->+    Tree a ->+    Tree a taoIterationCV cache cfg target condVecs df rootIndices tree =-    foldl' (optimizeDepthLevel env rootIndices) tree [treeDepth tree, treeDepth tree - 1 .. 0]+    foldl'+        (optimizeDepthLevel env rootIndices)+        tree+        [treeDepth tree, treeDepth tree - 1 .. 0]   where     env = TaoEnv cache cfg target condVecs df -optimizeDepthLevel :: forall a. (Columnable a, Ord a) => TaoEnv -> V.Vector Int -> Tree a -> Int -> Tree a+optimizeDepthLevel ::+    forall a.+    (Columnable a, Ord a) => TaoEnv -> V.Vector Int -> Tree a -> Int -> Tree a optimizeDepthLevel env rootIndices tree = optimizeAtDepth @a env rootIndices tree 0 -optimizeAtDepth :: forall a. (Columnable a, Ord a) => TaoEnv -> V.Vector Int -> Tree a -> Int -> Int -> Tree a+optimizeAtDepth ::+    forall a.+    (Columnable a, Ord a) =>+    TaoEnv -> V.Vector Int -> Tree a -> Int -> Int -> Tree a optimizeAtDepth env indices tree currentDepth targetDepth     | currentDepth == targetDepth = optimizeNode @a env indices tree     | otherwise = case tree of         Leaf v -> Leaf v         Branch cond left right -> optimizeChildren @a env indices cond left right currentDepth targetDepth --- | Optimize the two subtrees over their disjoint index sets, scoring the left--- in parallel with the right (the cache is read-only, so this is pure).-optimizeChildren :: forall a. (Columnable a, Ord a) => TaoEnv -> V.Vector Int -> Expr Bool -> Tree a -> Tree a -> Int -> Int -> Tree a+{- | Optimize the two subtrees over their disjoint index sets, scoring the left+in parallel with the right (the cache is read-only, so this is pure).+-}+optimizeChildren ::+    forall a.+    (Columnable a, Ord a) =>+    TaoEnv -> V.Vector Int -> Expr Bool -> Tree a -> Tree a -> Int -> Int -> Tree a optimizeChildren env indices cond left right currentDepth targetDepth =     forceTreeWork left' `par` (forceTreeWork right' `pseq` Branch cond left' right')   where@@ -100,15 +156,18 @@     left' = optimizeAtDepth @a env indicesL left (currentDepth + 1) targetDepth     right' = optimizeAtDepth @a env indicesR right (currentDepth + 1) targetDepth --- | Force a subtree's optimization work to WHNF so the parallel scheduler has--- something substantial to evaluate; pure and value-preserving.+{- | Force a subtree's optimization work to WHNF so the parallel scheduler has+something substantial to evaluate; pure and value-preserving.+-} forceTreeWork :: Tree a -> () forceTreeWork (Leaf v) = v `seq` () forceTreeWork (Branch c l r) = c `seq` forceTreeWork l `seq` forceTreeWork r --- | Re-optimize one node: pick its best split, or collapse to a leaf when the--- node is empty or the chosen split underflows 'minLeafSize'.-optimizeNode :: forall a. (Columnable a, Ord a) => TaoEnv -> V.Vector Int -> Tree a -> Tree a+{- | Re-optimize one node: pick its best split, or collapse to a leaf when the+node is empty or the chosen split underflows 'minLeafSize'.+-}+optimizeNode ::+    forall a. (Columnable a, Ord a) => TaoEnv -> V.Vector Int -> Tree a -> Tree a optimizeNode env indices tree     | V.null indices = tree     | otherwise = case tree of@@ -117,7 +176,10 @@   where     leaf = Leaf (majorityValueFromIndices @a (teTarget env) (teDf env) indices) -rebuiltBranch :: forall a. (Columnable a, Ord a) => TaoEnv -> V.Vector Int -> Expr Bool -> Tree a -> Tree a -> Tree a -> Tree a+rebuiltBranch ::+    forall a.+    (Columnable a, Ord a) =>+    TaoEnv -> V.Vector Int -> Expr Bool -> Tree a -> Tree a -> Tree a -> Tree a rebuiltBranch env indices oldCond left right leaf     | underflows = leaf     | otherwise = Branch newCond left right@@ -126,43 +188,79 @@     (l, r) = partitionIndicesCached (teCache env) newCond (teDf env) indices     underflows = V.length l < minLeafSize (teCfg env) || V.length r < minLeafSize (teCfg env) --- | The lowest-penalty replacement condition for a node, falling back to the--- current condition when no valid candidate beats it.-findBestSplitTAO :: forall a. (Columnable a) => TaoEnv -> V.Vector Int -> Tree a -> Tree a -> Expr Bool -> Expr Bool+{- | The lowest-penalty replacement condition for a node, falling back to the+current condition when no valid candidate beats it.+-}+findBestSplitTAO ::+    forall a.+    (Columnable a) =>+    TaoEnv -> V.Vector Int -> Tree a -> Tree a -> Expr Bool -> Expr Bool findBestSplitTAO env indices leftTree rightTree currentCond     | V.null indices || null carePoints = currentCond-    | pureReplacementLinear cfg, Just c <- linearCandidate, isValidAtNode cfg (teDf env) indices c = c+    | pureReplacementLinear cfg+    , Just c <- linearCandidate+    , isValidAtNode cfg (teDf env) indices c =+        c     | otherwise = bestOfPool penaltyCV currentCond pool   where     cfg = teCfg env-    carePoints = identifyCarePointsCached @a (teCache env) (teTarget env) (teDf env) indices leftTree rightTree+    carePoints =+        identifyCarePointsCached @a+            (teCache env)+            (teTarget env)+            (teDf env)+            indices+            leftTree+            rightTree     penaltyCV = evalWithPenaltyVec cfg carePoints     linearCandidate = bestLinearCandidate cfg (teDf env) carePoints     valid = filterValidCandidates cfg indices (teConds env)-    pool = candidatePool env indices currentCond (bestDiscreteCandidate cfg penaltyCV valid) linearCandidate+    pool =+        candidatePool+            env+            indices+            currentCond+            (bestDiscreteCandidate cfg penaltyCV valid)+            linearCandidate  bestOfPool :: (CondVec -> (Int, Int)) -> Expr Bool -> [CondVec] -> Expr Bool bestOfPool _ currentCond [] = currentCond bestOfPool penaltyCV _ pool = cvExpr (minimumBy (compare `on` penaltyCV) pool) --- | Validity-filtered candidates the node could split on: both children must--- keep at least 'minLeafSize'. Scored in parallel chunks, order preserved.+{- | Validity-filtered candidates the node could split on: both children must+keep at least 'minLeafSize'. Scored in parallel chunks, order preserved.+-} filterValidCandidates :: TreeConfig -> V.Vector Int -> [CondVec] -> [CondVec] filterValidCandidates cfg indices condVecs = map snd (filter fst (zip validity condVecs))   where-    validity = map (validAtNode cfg indices) condVecs `using` parListChunk candidateParChunk rdeepseq+    validity =+        map (validAtNode cfg indices) condVecs+            `using` parListChunk candidateParChunk rdeepseq  validAtNode :: TreeConfig -> V.Vector Int -> CondVec -> Bool validAtNode cfg indices cv = nTrue >= minLeaf && (V.length indices - nTrue) >= minLeaf   where     minLeaf = minLeafSize cfg-    nTrue = V.foldl' (\ !acc i -> if cvVec cv VU.! i then acc + 1 else acc) (0 :: Int) indices+    nTrue =+        V.foldl'+            (\ !acc i -> if cvVec cv VU.! i then acc + 1 else acc)+            (0 :: Int)+            indices --- | The candidate pool to minimize over: the current condition, the best--- discrete candidate, and the linear candidate, each kept only if valid.-candidatePool :: TaoEnv -> V.Vector Int -> Expr Bool -> Maybe CondVec -> Maybe (Expr Bool) -> [CondVec]+{- | The candidate pool to minimize over: the current condition, the best+discrete candidate, and the linear candidate, each kept only if valid.+-}+candidatePool ::+    TaoEnv ->+    V.Vector Int ->+    Expr Bool ->+    Maybe CondVec ->+    Maybe (Expr Bool) ->+    [CondVec] candidatePool env indices currentCond discreteCV linearCandidate =-    filter (isValidAtNode (teCfg env) (teDf env) indices . cvExpr) (catMaybes [currentCV, discreteCV, linearCV])+    filter+        (isValidAtNode (teCfg env) (teDf env) indices . cvExpr)+        (catMaybes [currentCV, discreteCV, linearCV])   where     currentCV = CondVec currentCond <$> lookupCondVec (teCache env) (teDf env) currentCond     linearCV = linearCandidate >>= materializeCondVec (teDf env)

@@ -5,8 +5,9 @@ {-# LANGUAGE ScopedTypeVariables #-} {-# LANGUAGE TypeApplications #-} --- | Shared types, configuration and ordering machinery for the decision-tree--- learner. Imported by every other @DataFrame.DecisionTree.*@ module.+{- | Shared types, configuration and ordering machinery for the decision-tree+learner. Imported by every other @DataFrame.DecisionTree.*@ module.+-} module DataFrame.DecisionTree.Types (     Tree (..),     treeDepth,@@ -38,8 +39,9 @@ import System.IO.Unsafe (unsafePerformIO) import Type.Reflection (SomeTypeRep (..), typeRep) --- | A fitted tree: a leaf value, or an internal node testing a boolean--- expression with @True@ routing left.+{- | A fitted tree: a leaf value, or an internal node testing a boolean+expression with @True@ routing left.+-} data Tree a     = Leaf !a     | Branch !(Expr Bool) !(Tree a) !(Tree a)@@ -49,8 +51,9 @@ treeDepth (Leaf _) = 0 treeDepth (Branch _ l r) = 1 + max (treeDepth l) (treeDepth r) --- | A row the parent node must route to a specific child for the subtrees to--- classify it correctly (the TAO objective is the count of misroutes).+{- | A row the parent node must route to a specific child for the subtrees to+classify it correctly (the TAO objective is the count of misroutes).+-} data CarePoint = CarePoint     { cpIndex :: !Int     , cpCorrectDir :: !Direction@@ -121,8 +124,9 @@         , pureReplacementLinear = False         } --- | Which column types support ordering for splits. Register a type with--- 'orderable' and combine with @<>@.+{- | Which column types support ordering for splits. Register a type with+'orderable' and combine with @<>@.+-} newtype ColumnOrdering = ColumnOrdering (M.Map SomeTypeRep OrdDict)  instance Semigroup ColumnOrdering where@@ -164,8 +168,9 @@ data OrdDict where     OrdDict :: (Columnable a, Ord a) => Proxy a -> OrdDict --- | Run @k@ with the @Ord a@ instance recovered from the ordering registry,--- or 'Nothing' when @a@ is not registered.+{- | Run @k@ with the @Ord a@ instance recovered from the ordering registry,+or 'Nothing' when @a@ is not registered.+-} withOrdFrom ::     forall a r. (Columnable a) => ColumnOrdering -> ((Ord a) => r) -> Maybe r withOrdFrom (ColumnOrdering m) k = case M.lookup (SomeTypeRep (typeRep @a)) m of

@@ -0,0 +1,132 @@+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE TypeApplications #-}++{- | Shared internal helpers used across the model fitters: turning a 'DataFrame'+plus a target/feature 'Expr' into the numeric matrices the algorithms consume,+and the common expression builders (affine score, arg-max/arg-min over named+scores) that every linear model and classifier would otherwise re-implement.+-}+module DataFrame.Featurize.Internal (+    -- * Supervised extraction+    featureNames,+    numericMatrix,+    targetDoubles,+    targetValues,++    -- * Unsupervised extraction+    Features (..),+    extractFeatures,+    columnExprName,+    materializeColumn,++    -- * Expression builders+    affineExpr,+    argMaxExpr,+    argMinExpr,+) where++import Control.Exception (throw)+import qualified Data.Text as T+import qualified Data.Vector as V+import qualified Data.Vector.Unboxed as VU++import qualified DataFrame.Functions as F+import DataFrame.Internal.Column (Columnable)+import DataFrame.Internal.DataFrame (DataFrame, columnNames)+import DataFrame.Internal.Expression (Expr (..))+import DataFrame.LinearAlgebra (Matrix, transposeM)+import DataFrame.Operations.Core (columnAsDoubleVector, columnAsVector)+import DataFrame.Operators ((.&&.), (.*.), (.+.), (.<=.), (.>=.))++-- | Every column name except the supervised target's.+featureNames :: Expr a -> DataFrame -> [T.Text]+featureNames (Col target) df = filter (/= target) (columnNames df)+featureNames _ df = columnNames df++{- | The named columns as a row-major @n×d@ matrix of doubles (non-numeric+columns are coerced through 'columnAsDoubleVector'), paired with the names.+-}+numericMatrix :: [T.Text] -> DataFrame -> (V.Vector T.Text, Matrix)+numericMatrix names df = (V.fromList names, transposeM colMajor)+  where+    colMajor = V.fromList (map column names)+    column name = case columnAsDoubleVector (F.col @Double name) df of+        Right v -> v+        Left e -> throw e++-- | The target column as a vector of doubles.+targetDoubles :: Expr Double -> DataFrame -> VU.Vector Double+targetDoubles expr df = case columnAsDoubleVector expr df of+    Right v -> v+    Left e -> throw e++-- | The target column as a vector of its own type (for classifiers).+targetValues :: (Columnable a) => Expr a -> DataFrame -> V.Vector a+targetValues expr df = case columnAsVector expr df of+    Right v -> v+    Left e -> throw e++{- | The extracted feature columns of an unsupervised fit, in the shapes the+algorithms need: names, column-major vectors, the row-major matrix, and the+@(n, d)@ dimensions.+-}+data Features = Features+    { ftNames :: ![T.Text]+    , ftCols :: ![VU.Vector Double]+    , ftRows :: !Matrix+    , ftN :: !Int+    , ftD :: !Int+    }++-- | Extract the given feature columns once, in every shape the fitters use.+extractFeatures :: [Expr Double] -> DataFrame -> Features+extractFeatures features df = Features names cols rows n d+  where+    names = map columnExprName features+    cols = map (materializeColumn df) features+    n = if null cols then 0 else VU.length (head cols)+    d = length cols+    rows = V.generate n (\i -> VU.generate d (\j -> (cols !! j) VU.! i))++-- | The column name behind a @Col@ feature expression.+columnExprName :: Expr Double -> T.Text+columnExprName (Col n) = n+columnExprName e = error ("expected a column expression, got " ++ show e)++-- | Interpret a @Col@ (or numeric) expression to a @Double@ vector.+materializeColumn :: DataFrame -> Expr Double -> VU.Vector Double+materializeColumn df e = case columnAsDoubleVector e df of+    Right v -> v+    Left err -> throw err++{- | An affine score @b + Σ wⱼ·colⱼ@ over named columns, dropping zero weights+(the shared core of linear/logistic/SVM margins).+-}+affineExpr :: Double -> [(Double, T.Text)] -> Expr Double+affineExpr b terms =+    foldr+        (.+.)+        (F.lit b)+        [F.lit w .*. (Col n :: Expr Double) | (w, n) <- terms, w /= 0]++{- | The class whose score is greatest, as a nested-@If@ expression; ties go to+the earlier class.+-}+argMaxExpr :: (Columnable a) => [(a, Expr Double)] -> Expr a+argMaxExpr = argExtreme (.>=.)++-- | The class whose score is smallest (e.g. nearest cluster by distance).+argMinExpr :: (Columnable a) => [(a, Expr Double)] -> Expr a+argMinExpr = argExtreme (.<=.)++argExtreme ::+    (Columnable a) =>+    (Expr Double -> Expr Double -> Expr Bool) -> [(a, Expr Double)] -> Expr a+argExtreme _ [] = error "argExtreme: no classes"+argExtreme _ [(c, _)] = Lit c+argExtreme cmp ((c, sc) : rest) =+    If+        (foldr ((\o acc -> cmp sc o .&&. acc) . snd) (F.lit True) rest)+        (Lit c)+        (argExtreme cmp rest)

@@ -0,0 +1,310 @@+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE ScopedTypeVariables #-}++{- | Gaussian mixture models fitted by EM. Full covariance by default (with a+diagonal option and an automatic fall-back when a covariance is not positive+definite), log-space responsibilities, and Cholesky-based densities for+stability. 'predict' is the hard (arg-max) component assignment; per-component+log-densities are available via 'gmmLogDensityExprs'.+-}+module DataFrame.GMM (+    CovType (..),+    GMMConfig (..),+    defaultGMMConfig,+    GMMModel (..),+    gmmLogDensityExprs,+    gmmBIC,+    gmmAIC,+) where++import Data.List (sortBy)+import qualified Data.Map.Strict as M+import Data.Ord (comparing)+import qualified Data.Text as T+import qualified Data.Vector as V+import qualified Data.Vector.Unboxed as VU++import DataFrame.Featurize.Internal (Features (..), argMaxExpr, extractFeatures)+import qualified DataFrame.Functions as F+import DataFrame.Internal.DataFrame (DataFrame)+import DataFrame.Internal.Expression (Expr (..))+import DataFrame.LinearAlgebra (Matrix, logSumExp)+import DataFrame.LinearAlgebra.Solve (backSubst, cholesky, forwardSubst)+import DataFrame.Model (Fit (..), Predict (..))+import DataFrame.Operators ((.*.), (.+.), (.-.))+import DataFrame.Random (mkGen, sampleIndices)++data CovType = FullCov | DiagCov+    deriving (Eq, Show)++data GMMConfig = GMMConfig+    { gmmK :: !Int+    , gmmCovType :: !CovType+    , gmmMaxIter :: !Int+    , gmmTol :: !Double+    , gmmRegCovar :: !Double+    , gmmSeed :: !Int+    }+    deriving (Eq, Show)++defaultGMMConfig :: GMMConfig+defaultGMMConfig =+    GMMConfig+        { gmmK = 2+        , gmmCovType = FullCov+        , gmmMaxIter = 100+        , gmmTol = 1.0e-3+        , gmmRegCovar = 1.0e-6+        , gmmSeed = 0+        }++-- | A fitted mixture. 'gmmCovariances' are the per-component covariance matrices.+data GMMModel = GMMModel+    { gmmWeights :: !(VU.Vector Double)+    , gmmMeans :: !(V.Vector (VU.Vector Double))+    , gmmCovariances :: !(V.Vector Matrix)+    , gmmConverged :: !Bool+    , gmmNIter :: !Int+    , gmmLogLikelihood :: !Double+    , gmmNObs :: !Int+    , gmmFeatureNames :: !(V.Vector T.Text)+    }+    deriving (Eq, Show)++instance Fit GMMConfig [Expr Double] GMMModel where+    fit = fitGMM++instance Predict GMMModel Int where+    predict = gmmAssignExpr++-- | Fit a Gaussian mixture over the given feature columns.+fitGMM :: GMMConfig -> [Expr Double] -> DataFrame -> GMMModel+fitGMM cfg features df = canonical finalModel+  where+    Features names _ rows n d = extractFeatures features df+    k = min (gmmK cfg) (max 1 n)+    reg = gmmRegCovar cfg+    (initIdx, _) = sampleIndices k n (mkGen (gmmSeed cfg))+    means0 = V.map (rows V.!) (V.convert initIdx)+    varDiag =+        VU.generate d $ \j ->+            let mu = sum [(rows V.! i) VU.! j | i <- [0 .. n - 1]] / fromIntegral (max 1 n)+             in ( sum [((rows V.! i) VU.! j - mu) ^ (2 :: Int) | i <- [0 .. n - 1]]+                    / fromIntegral (max 1 n)+                )+                    + reg+    cov0 = diagMatrix varDiag+    covs0 = V.replicate k cov0+    weights0 = VU.replicate k (1 / fromIntegral k)+    finalModel = em 0 weights0 means0 covs0 (-(1 / 0)) False+    em !iter weights means covs prevLL converged+        | iter >= gmmMaxIter cfg || converged =+            GMMModel weights means covs converged iter prevLL n (V.fromList names)+        | otherwise =+            let (logResp, ll) = eStep cfg rows weights means covs+                (weights', means', covs') = mStep cfg rows logResp d reg+                done = abs (ll - prevLL) < gmmTol cfg+             in em (iter + 1) weights' means' covs' ll done++-- | Per-component log-density expressions (log weight + Gaussian log pdf).+gmmLogDensityExprs :: GMMModel -> M.Map Int (Expr Double)+gmmLogDensityExprs m =+    M.fromList+        [ ( c+          , logDensityExpr+                (gmmWeights m VU.! c)+                (gmmMeans m V.! c)+                (gmmCovariances m V.! c)+                names+          )+        | c <- [0 .. V.length (gmmMeans m) - 1]+        ]+  where+    names = V.toList (gmmFeatureNames m)++-- | The hard-assignment expression: arg-max of component log-densities.+gmmAssignExpr :: GMMModel -> Expr Int+gmmAssignExpr m =+    argMaxExpr (M.toList (gmmLogDensityExprs m))++-- | Bayesian information criterion (lower is better).+gmmBIC :: GMMModel -> Double+gmmBIC m =+    negate (2 * gmmLogLikelihood m)+        + fromIntegral (nParams m) * log (fromIntegral (max 1 (gmmNObs m)))++-- | Akaike information criterion (lower is better).+gmmAIC :: GMMModel -> Double+gmmAIC m = negate (2 * gmmLogLikelihood m) + 2 * fromIntegral (nParams m)++nParams :: GMMModel -> Int+nParams m =+    let k = VU.length (gmmWeights m)+        d = if V.null (gmmMeans m) then 0 else VU.length (V.head (gmmMeans m))+     in (k - 1) + k * d + k * (d * (d + 1) `div` 2)++eStep ::+    GMMConfig ->+    Matrix ->+    VU.Vector Double ->+    V.Vector (VU.Vector Double) ->+    V.Vector Matrix ->+    (V.Vector (VU.Vector Double), Double)+eStep _ rows weights means covs = (logResp, totalLL)+  where+    k = VU.length weights+    comps =+        V.generate k $ \c ->+            ( log (max 1e-300 (weights VU.! c))+            , means V.! c+            , gaussianLogPdf (covs V.! c) (means V.! c)+            )+    perRow x =+        let lps = VU.generate k (\c -> let (lw, _, f) = comps V.! c in lw + f x)+            lse = logSumExp lps+         in (VU.map (subtract lse) lps, lse)+    results = V.map perRow rows+    logResp = V.map fst results+    totalLL = V.sum (V.map snd results)++mStep ::+    GMMConfig ->+    Matrix ->+    V.Vector (VU.Vector Double) ->+    Int ->+    Double ->+    (VU.Vector Double, V.Vector (VU.Vector Double), V.Vector Matrix)+mStep cfg rows logResp d reg = (weights, means, covs)+  where+    n = V.length rows+    k = if V.null logResp then 0 else VU.length (V.head logResp)+    resp = V.map (VU.map exp) logResp+    nk = VU.generate k (\c -> sum [resp V.! i VU.! c | i <- [0 .. n - 1]])+    weights = VU.map (/ fromIntegral (max 1 n)) nk+    means =+        V.generate k $ \c ->+            let s =+                    foldr+                        (VU.zipWith (+))+                        (VU.replicate d 0)+                        [VU.map (* (resp V.! i VU.! c)) (rows V.! i) | i <- [0 .. n - 1]]+             in VU.map (/ max 1e-12 (nk VU.! c)) s+    covs =+        V.generate k $ \c ->+            let mu = means V.! c+                acc = foldr addOuter (zeroMatrix d) [0 .. n - 1]+                addOuter i m =+                    let diff = VU.zipWith (-) (rows V.! i) mu+                        w = resp V.! i VU.! c+                     in addScaledOuter w diff m+                scaled = scaleMatrix (1 / max 1e-12 (nk VU.! c)) acc+                regd = addDiagScalar reg scaled+             in case gmmCovType cfg of+                    FullCov -> regd+                    DiagCov -> diagOnly regd++gaussianLogPdf :: Matrix -> VU.Vector Double -> VU.Vector Double -> Double+gaussianLogPdf cov mu =+    case cholesky cov of+        Just l ->+            let logdet = 2 * sum [log ((l V.! i) VU.! i) | i <- [0 .. d - 1]]+             in \x ->+                    let diff = VU.zipWith (-) x mu+                        z = forwardSubst l diff+                        quad = VU.sum (VU.map (^ (2 :: Int)) z)+                     in negate 0.5 * (fromIntegral d * log (2 * pi) + logdet + quad)+        Nothing ->+            let var = VU.generate d (\i -> max 1e-12 ((cov V.! i) VU.! i))+             in \x ->+                    negate 0.5+                        * VU.sum+                            ( VU.generate d $ \j ->+                                let diff = x VU.! j - mu VU.! j+                                 in diff * diff / var VU.! j + log (2 * pi * var VU.! j)+                            )+  where+    d = VU.length mu++logDensityExpr ::+    Double -> VU.Vector Double -> Matrix -> [T.Text] -> Expr Double+logDensityExpr weight mu cov names =+    case precisionAndLogdet cov of+        Just (prec, logdet) ->+            let constTerm =+                    log (max 1e-300 weight)+                        - 0.5 * (fromIntegral d * log (2 * pi) + logdet)+                quad =+                    foldr (.+.) (F.lit 0) $+                        [ F.lit (-(0.5 * prec V.! a VU.! b)) .*. (centered a .*. centered b)+                        | a <- [0 .. d - 1]+                        , b <- [0 .. d - 1]+                        ]+             in F.lit constTerm .+. quad+        Nothing -> F.lit (log (max 1e-300 weight))+  where+    d = VU.length mu+    centered j = (Col (names !! j) :: Expr Double) .-. F.lit (mu VU.! j)++precisionAndLogdet :: Matrix -> Maybe (Matrix, Double)+precisionAndLogdet cov = do+    l <- cholesky cov+    let d = V.length cov+        logdet = 2 * sum [log ((l V.! i) VU.! i) | i <- [0 .. d - 1]]+        cols = [forwardThenBack l (unitVec d i) | i <- [0 .. d - 1]]+        prec =+            V.fromList+                [VU.fromList [cols !! j VU.! i | j <- [0 .. d - 1]] | i <- [0 .. d - 1]]+    pure (prec, logdet)++forwardThenBack :: Matrix -> VU.Vector Double -> VU.Vector Double+forwardThenBack l b = backSubst l (forwardSubst l b)++canonical :: GMMModel -> GMMModel+canonical m =+    let order =+            map snd $+                sortBy+                    (comparing fst)+                    [ (firstCoord (gmmMeans m V.! c), c)+                    | c <- [0 .. V.length (gmmMeans m) - 1]+                    ]+        firstCoord v = if VU.null v then 0 else VU.head v+     in m+            { gmmWeights = VU.fromList [gmmWeights m VU.! c | c <- order]+            , gmmMeans = V.fromList [gmmMeans m V.! c | c <- order]+            , gmmCovariances = V.fromList [gmmCovariances m V.! c | c <- order]+            }++diagMatrix :: VU.Vector Double -> Matrix+diagMatrix v =+    let d = VU.length v+     in V.generate d (\i -> VU.generate d (\j -> if i == j then v VU.! i else 0))++diagOnly :: Matrix -> Matrix+diagOnly m =+    let d = V.length m+     in V.generate+            d+            (\i -> VU.generate d (\j -> if i == j then (m V.! i) VU.! j else 0))++zeroMatrix :: Int -> Matrix+zeroMatrix d = V.replicate d (VU.replicate d 0)++scaleMatrix :: Double -> Matrix -> Matrix+scaleMatrix s = V.map (VU.map (* s))++addDiagScalar :: Double -> Matrix -> Matrix+addDiagScalar s = V.imap (\i row -> row VU.// [(i, row VU.! i + s)])++addScaledOuter :: Double -> VU.Vector Double -> Matrix -> Matrix+addScaledOuter w diff m =+    let d = VU.length diff+     in V.generate d $ \i ->+            VU.generate d $ \j ->+                (m V.! i) VU.! j + w * (diff VU.! i) * (diff VU.! j)++unitVec :: Int -> Int -> VU.Vector Double+unitVec d i = VU.generate d (\j -> if i == j then 1 else 0)

@@ -0,0 +1,156 @@+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE OverloadedStrings #-}++{- | k-means clustering (Lloyd's algorithm with k-means++ seeding and multiple+restarts). 'fit' trains a 'KMeansModel' (inspectable centres); 'predict' is the+arg-min cluster assignment. Per-cluster distance features are available via+'kmeansDistanceExprs' / 'kmeansTransform'.+-}+module DataFrame.KMeans (+    KMeansConfig (..),+    defaultKMeansConfig,+    KMeansModel (..),+    kmeansDistanceExprs,+    kmeansTransform,+) where++import Data.List (minimumBy)+import Data.Ord (comparing)+import qualified Data.Text as T+import qualified Data.Vector as V+import qualified Data.Vector.Unboxed as VU++import DataFrame.Featurize.Internal (Features (..), argMinExpr, extractFeatures)+import qualified DataFrame.Functions as F+import DataFrame.Internal.DataFrame (DataFrame)+import DataFrame.Internal.Expression (Expr (..), UExpr (..))+import DataFrame.LinearAlgebra (Matrix, nearestCenter, sqDist)+import DataFrame.Model (Fit (..), Predict (..))+import DataFrame.Operators ((.*.), (.+.), (.-.))+import DataFrame.Random (Gen, mkGen, nextDouble, nextIntR, splitGen)+import DataFrame.Transform (Transform (..))++data KMeansConfig = KMeansConfig+    { kmK :: !Int+    , kmNInit :: !Int+    , kmMaxIter :: !Int+    , kmTol :: !Double+    , kmSeed :: !Int+    }+    deriving (Eq, Show)++defaultKMeansConfig :: KMeansConfig+defaultKMeansConfig =+    KMeansConfig{kmK = 8, kmNInit = 10, kmMaxIter = 300, kmTol = 1.0e-4, kmSeed = 0}++-- | A fitted k-means model. 'kmCenters' are sklearn's @cluster_centers_@.+data KMeansModel = KMeansModel+    { kmCenters :: !(V.Vector (VU.Vector Double))+    , kmLabels :: !(VU.Vector Int)+    , kmInertia :: !Double+    , kmNIter :: !Int+    , kmFeatureNames :: !(V.Vector T.Text)+    }+    deriving (Eq, Show)++instance Fit KMeansConfig [Expr Double] KMeansModel where+    fit = fitKMeans++instance Predict KMeansModel Int where+    predict m = argMinExpr (zip [0 :: Int ..] (map snd (kmeansDistanceExprs m)))++-- | Fit k-means over the given feature columns.+fitKMeans :: KMeansConfig -> [Expr Double] -> DataFrame -> KMeansModel+fitKMeans cfg features df = best+  where+    Features names _ rows n _ = extractFeatures features df+    k = min (kmK cfg) (max 1 n)+    seeds = take (max 1 (kmNInit cfg)) (genSeeds (mkGen (kmSeed cfg)))+    runs = map (lloyd cfg k rows) seeds+    best =+        let (centers, labels, inertia, iters) =+                minimumBy (comparing (\(_, _, i, _) -> i)) runs+         in KMeansModel centers labels inertia iters (V.fromList names)++genSeeds :: Gen -> [Gen]+genSeeds g = let (g1, g2) = splitGen g in g1 : genSeeds g2++{- | One k-means run: k-means++ seeding then Lloyd iterations. Returns+@(centers, labels, inertia, nIter)@.+-}+lloyd ::+    KMeansConfig ->+    Int ->+    Matrix ->+    Gen ->+    (V.Vector (VU.Vector Double), VU.Vector Int, Double, Int)+lloyd cfg k rows g0+    | V.null rows = (V.empty, VU.empty, 0, 0)+    | otherwise = iterate' 0 initCenters+  where+    initCenters = kmeansPP k rows g0+    iterate' !iter centers =+        let labels = VU.generate (V.length rows) (\i -> fst (nearestCenter centers (rows V.! i)))+            newCenters = recompute centers labels+            shift = V.sum (V.zipWith sqDist centers newCenters)+         in if iter + 1 >= kmMaxIter cfg || shift <= kmTol cfg+                then (newCenters, labels, inertiaOf newCenters labels, iter + 1)+                else iterate' (iter + 1) newCenters+    recompute centers labels =+        V.generate k $ \c ->+            let members = [rows V.! i | i <- [0 .. V.length rows - 1], labels VU.! i == c]+             in if null members then centers V.! c else meanOf members+    inertiaOf centers labels =+        sum+            [ sqDist (rows V.! i) (centers V.! (labels VU.! i))+            | i <- [0 .. V.length rows - 1]+            ]++meanOf :: [VU.Vector Double] -> VU.Vector Double+meanOf vs =+    let d = VU.length (head vs)+        s = foldr (VU.zipWith (+)) (VU.replicate d 0) vs+     in VU.map (/ fromIntegral (length vs)) s++-- | k-means++ seeding: spread initial centres by squared-distance sampling.+kmeansPP :: Int -> Matrix -> Gen -> V.Vector (VU.Vector Double)+kmeansPP k rows g0 = V.fromList (reverse (pick [first] g1))+  where+    n = V.length rows+    (i0, g1) = nextIntR (0, n - 1) g0+    first = rows V.! i0+    pick chosen g+        | length chosen >= k = chosen+        | otherwise =+            let dists = VU.generate n (\i -> minimum [sqDist (rows V.! i) c | c <- chosen])+                (u, g') = nextDouble g+                idx = sampleCumulative dists (u * VU.sum dists)+             in pick (rows V.! idx : chosen) g'++sampleCumulative :: VU.Vector Double -> Double -> Int+sampleCumulative dists target = go 0 0+  where+    go i acc+        | i >= VU.length dists - 1 = VU.length dists - 1+        | acc + dists VU.! i >= target = i+        | otherwise = go (i + 1) (acc + dists VU.! i)++-- | Per-cluster squared-distance expressions, named @dist1@, @dist2@, …+kmeansDistanceExprs :: KMeansModel -> [(T.Text, Expr Double)]+kmeansDistanceExprs m =+    [ ("dist" <> T.pack (show (c + 1)), distExpr (kmCenters m V.! c))+    | c <- [0 .. V.length (kmCenters m) - 1]+    ]+  where+    names = V.toList (kmFeatureNames m)+    distExpr center =+        foldr (.+.) (F.lit 0) $+            [ let diff = (Col n :: Expr Double) .-. F.lit ci in diff .*. diff+            | (n, ci) <- zip names (VU.toList center)+            ]++-- | The per-cluster distance features as a composable fitted 'Transform'.+kmeansTransform :: KMeansModel -> Transform+kmeansTransform m = Transform [(n, UExpr e) | (n, e) <- kmeansDistanceExprs m]

@@ -0,0 +1,122 @@+{-# LANGUAGE BangPatterns #-}++{- | Dependency-free dense linear algebra over row-major matrices, shared by the+models in @dataframe-learn@. Basic vector/matrix operations plus stability and+distance helpers; solvers live in "DataFrame.LinearAlgebra.Solve" and+eigenproblems in "DataFrame.LinearAlgebra.Eigen".+-}+module DataFrame.LinearAlgebra (+    Matrix,+    dot,+    axpy,+    scaleV,+    matVec,+    tMatVec,+    gram,+    transposeM,+    identityM,+    logSumExp,+    sqDist,+    nearestCenter,+    epsNeighbors,+) where++import qualified Data.Vector as V+import qualified Data.Vector.Unboxed as VU++{- | Row-major dense matrix: an outer boxed vector of equal-length rows. An+@n×d@ matrix has @n@ rows of length @d@.+-}+type Matrix = V.Vector (VU.Vector Double)++-- | Inner product of two equal-length vectors.+dot :: VU.Vector Double -> VU.Vector Double -> Double+dot a b = VU.sum (VU.zipWith (*) a b)+{-# INLINE dot #-}++-- | @axpy a x y = a*x + y@.+axpy :: Double -> VU.Vector Double -> VU.Vector Double -> VU.Vector Double+axpy a = VU.zipWith (\xi yi -> a * xi + yi)+{-# INLINE axpy #-}++-- | Scalar-vector product.+scaleV :: Double -> VU.Vector Double -> VU.Vector Double+scaleV a = VU.map (* a)+{-# INLINE scaleV #-}++-- | @matVec A v@ for @A@ of shape @n×d@ and @v@ of length @d@; result length @n@.+matVec :: Matrix -> VU.Vector Double -> VU.Vector Double+matVec a v = VU.convert (V.map (`dot` v) a)++-- | @tMatVec A v = Aᵀ v@ for @A@ of shape @n×d@, @v@ of length @n@; result length @d@.+tMatVec :: Matrix -> VU.Vector Double -> VU.Vector Double+tMatVec a v+    | V.null a = VU.empty+    | otherwise = V.foldl' step (VU.replicate d 0) (V.zipWith (,) vBoxed a)+  where+    d = VU.length (V.head a)+    vBoxed = V.generate (V.length a) (v VU.!)+    step !acc (vi, row) = axpy vi row acc++-- | @gram A = Aᵀ A@, the @d×d@ symmetric matrix of column inner products.+gram :: Matrix -> Matrix+gram a+    | V.null a = V.empty+    | otherwise =+        V.generate d $ \i ->+            VU.generate d $ \j ->+                V.foldl' (\ !acc row -> acc + (row VU.! i) * (row VU.! j)) 0 a+  where+    d = VU.length (V.head a)++-- | Transpose an @n×d@ matrix to @d×n@.+transposeM :: Matrix -> Matrix+transposeM a+    | V.null a = V.empty+    | otherwise = V.generate d $ \j -> VU.generate n $ \i -> (a V.! i) VU.! j+  where+    n = V.length a+    d = VU.length (V.head a)++-- | @d×d@ identity matrix.+identityM :: Int -> Matrix+identityM d = V.generate d $ \i -> VU.generate d $ \j -> if i == j then 1 else 0++-- | Numerically stable @log Σ exp xᵢ@.+logSumExp :: VU.Vector Double -> Double+logSumExp xs+    | VU.null xs = negate (1 / 0)+    | otherwise = m + log (VU.sum (VU.map (\x -> exp (x - m)) xs))+  where+    m = VU.maximum xs++-- | Squared Euclidean distance.+sqDist :: VU.Vector Double -> VU.Vector Double -> Double+sqDist a b = VU.sum (VU.zipWith (\x y -> let z = x - y in z * z) a b)+{-# INLINE sqDist #-}++-- | Index of and squared distance to the nearest centre.+nearestCenter ::+    V.Vector (VU.Vector Double) -> VU.Vector Double -> (Int, Double)+nearestCenter centers p =+    V.ifoldl'+        ( \(!bi, !bd) i c ->+            let dd = sqDist c p in if dd < bd then (i, dd) else (bi, bd)+        )+        (-1, 1 / 0)+        centers++{- | Indices @j@ (excluding @i@) within squared radius @eps²@ of row @i@, by+brute force; @O(n d)@ per query.+-}+epsNeighbors :: Double -> Matrix -> Int -> VU.Vector Int+epsNeighbors eps rows i =+    VU.fromList+        [ j+        | j <- [0 .. n - 1]+        , j /= i+        , sqDist (rows V.! i) (rows V.! j) <= eps2+        ]+  where+    n = V.length rows+    eps2 = eps * eps

@@ -0,0 +1,121 @@+{-# LANGUAGE BangPatterns #-}++{- | Symmetric eigenproblems in pure Haskell: cyclic Jacobi for full+decomposition (PCA covariance, @m×m@ kernels) and power iteration for the+dominant eigenpair (FISTA step sizes). Deterministic, sign-canonicalised output.+-}+module DataFrame.LinearAlgebra.Eigen (+    jacobiEigenSym,+    powerIterTop,+) where++import Control.Monad (forM_, when)+import Control.Monad.ST (runST)+import Data.List (sortBy)+import Data.Ord (Down (..), comparing)+import qualified Data.Vector as V+import qualified Data.Vector.Unboxed as VU+import qualified Data.Vector.Unboxed.Mutable as VUM+import DataFrame.LinearAlgebra (Matrix, dot, matVec, scaleV)++{- | Cyclic Jacobi eigendecomposition of a symmetric matrix. Eigenvalues are+returned in descending order paired with eigenvectors as rows; each eigenvector+is sign-canonicalised (largest-magnitude component positive) so the output is+unique.+-}+jacobiEigenSym :: Matrix -> (VU.Vector Double, Matrix)+jacobiEigenSym a0+    | V.null a0 = (VU.empty, V.empty)+    | otherwise = runST $ do+        a <- VUM.new (d * d)+        forM_ [0 .. d - 1] $ \i ->+            forM_ [0 .. d - 1] $ \j ->+                VUM.write a (i * d + j) ((a0 V.! i) VU.! j)+        v <- VUM.replicate (d * d) 0+        forM_ [0 .. d - 1] $ \i -> VUM.write v (i * d + i) 1+        sweep a v 0+        afrozen <- VU.freeze a+        vmat <- VU.freeze v+        let diag = VU.generate d (\i -> afrozen VU.! (i * d + i))+            vecs =+                V.generate d $ \col ->+                    VU.generate d $ \row -> vmat VU.! (row * d + col)+            paired =+                sortBy+                    (comparing (Down . fst))+                    (zip (VU.toList diag) (V.toList vecs))+        pure+            ( VU.fromList (map fst paired)+            , V.fromList (map (canonicalSign . snd) paired)+            )+  where+    d = V.length a0+    maxSweeps = 100+    tol = 1e-12+    sweep a v s+        | s >= maxSweeps = pure ()+        | otherwise = do+            off <- offNorm a+            when (off >= tol) $ do+                forM_ [0 .. d - 2] $ \p ->+                    forM_ [p + 1 .. d - 1] $ \q -> rotate a v p q+                sweep a v (s + 1)+    offNorm a = go 0 0+      where+        go i !acc+            | i >= d = pure acc+            | otherwise = do+                r <- goRow i (i + 1) acc+                go (i + 1) r+        goRow i j !acc+            | j >= d = pure acc+            | otherwise = do+                x <- VUM.read a (i * d + j)+                goRow i (j + 1) (acc + x * x)+    rotate a v p q = do+        apq <- VUM.read a (p * d + q)+        when (abs apq > 1e-300) $ do+            app <- VUM.read a (p * d + p)+            aqq <- VUM.read a (q * d + q)+            let theta = (aqq - app) / (2 * apq)+                s' = if theta == 0 then 1 else signum theta+                t = s' / (abs theta + sqrt (theta * theta + 1))+                c = 1 / sqrt (t * t + 1)+                sn = t * c+            forM_ [0 .. d - 1] $ \i -> do+                aip <- VUM.read a (i * d + p)+                aiq <- VUM.read a (i * d + q)+                VUM.write a (i * d + p) (c * aip - sn * aiq)+                VUM.write a (i * d + q) (sn * aip + c * aiq)+            forM_ [0 .. d - 1] $ \j -> do+                apj <- VUM.read a (p * d + j)+                aqj <- VUM.read a (q * d + j)+                VUM.write a (p * d + j) (c * apj - sn * aqj)+                VUM.write a (q * d + j) (sn * apj + c * aqj)+            forM_ [0 .. d - 1] $ \i -> do+                vip <- VUM.read v (i * d + p)+                viq <- VUM.read v (i * d + q)+                VUM.write v (i * d + p) (c * vip - sn * viq)+                VUM.write v (i * d + q) (sn * vip + c * viq)++canonicalSign :: VU.Vector Double -> VU.Vector Double+canonicalSign vec =+    let idx = VU.maxIndex (VU.map abs vec)+     in if vec VU.! idx < 0 then VU.map negate vec else vec++{- | Dominant eigenvalue and eigenvector of a symmetric PSD matrix via power+iteration with a deterministic all-ones start.+-}+powerIterTop :: Int -> Matrix -> (Double, VU.Vector Double)+powerIterTop iters a+    | V.null a = (0, VU.empty)+    | otherwise = go iters (normalize (VU.replicate d 1))+  where+    d = V.length a+    normalize v =+        let nrm = sqrt (dot v v) in if nrm == 0 then v else scaleV (1 / nrm) v+    go 0 v = (dot v (matVec a v), v)+    go k v =+        let av = matVec a v+            nrm = sqrt (dot av av)+         in if nrm < 1e-300 then (0, v) else go (k - 1) (scaleV (1 / nrm) av)

@@ -0,0 +1,213 @@+{-# LANGUAGE BangPatterns #-}++{- | Householder QR (for ordinary least squares) and Cholesky factorisation (for+ridge normal equations and Gaussian log-densities). Pure, deterministic, no+LAPACK; sound at the @d@ ≤ low-hundreds scales this library targets.+-}+module DataFrame.LinearAlgebra.Solve (+    qrLeastSquares,+    cholesky,+    choleskySolve,+    logDetFromChol,+    forwardSubst,+    backSubst,+) where++import Control.Monad (forM_)+import Control.Monad.ST (ST, runST)+import qualified Data.Vector as V+import qualified Data.Vector.Unboxed as VU+import qualified Data.Vector.Unboxed.Mutable as VUM+import DataFrame.LinearAlgebra (Matrix)++{- | Solve @min ‖A x − b‖₂@ for an @n×d@ matrix @A@ (@n ≥ d@) via Householder QR.+@Left cols@ reports rank deficiency (near-zero @R@ diagonal) with the offending+column indices; @Right x@ is the least-squares solution.+-}+qrLeastSquares :: Matrix -> VU.Vector Double -> Either [Int] (VU.Vector Double)+qrLeastSquares a b+    | V.null a = Right VU.empty+    | n < d = Left [0 .. d - 1]+    | otherwise = runST $ do+        mat <- VUM.new (n * d)+        forM_ [0 .. n - 1] $ \i ->+            forM_ [0 .. d - 1] $ \j ->+                VUM.write mat (j * n + i) ((a V.! i) VU.! j)+        rhs <- VUM.new n+        forM_ [0 .. n - 1] $ \i -> VUM.write rhs i (b VU.! i)+        deficient <- householder mat rhs n d+        if not (null deficient)+            then pure (Left deficient)+            else do+                x <- VUM.new d+                backSubstQR mat rhs n d x+                Right <$> VU.freeze x+  where+    n = V.length a+    d = VU.length (V.head a)++householder ::+    VUM.STVector s Double -> VUM.STVector s Double -> Int -> Int -> ST s [Int]+householder mat rhs n d = go 0 []+  where+    tol = 1e-10+    go k acc+        | k >= d = pure (reverse acc)+        | otherwise = do+            normSq <- sumSq k+            let alphaMag = sqrt normSq+            akk <- VUM.read mat (k * n + k)+            let alpha = if akk > 0 then negate alphaMag else alphaMag+            if alphaMag < tol+                then go (k + 1) (k : acc)+                else do+                    VUM.write mat (k * n + k) (akk - alpha)+                    vNormSq <- sumSq k+                    if vNormSq < tol * tol+                        then go (k + 1) acc+                        else do+                            forM_ [k + 1 .. d - 1] $ \j -> reflectColumn k j vNormSq+                            reflectRhs k vNormSq+                            VUM.write mat (k * n + k) alpha+                            go (k + 1) acc+    sumSq k = foldRows k+      where+        foldRows i+            | i >= n = pure 0+            | otherwise = do+                x <- VUM.read mat (k * n + i)+                rest <- foldRows (i + 1)+                pure (x * x + rest)+    reflectColumn k j vNormSq = do+        dotv <- dotV k j k+        let beta = 2 * dotv / vNormSq+        forM_ [k .. n - 1] $ \i -> do+            vi <- VUM.read mat (k * n + i)+            aij <- VUM.read mat (j * n + i)+            VUM.write mat (j * n + i) (aij - beta * vi)+    reflectRhs k vNormSq = do+        dotv <- dotRhs k k+        let beta = 2 * dotv / vNormSq+        forM_ [k .. n - 1] $ \i -> do+            vi <- VUM.read mat (k * n + i)+            bi <- VUM.read rhs i+            VUM.write rhs i (bi - beta * vi)+    dotV k j i+        | i >= n = pure 0+        | otherwise = do+            vi <- VUM.read mat (k * n + i)+            aij <- VUM.read mat (j * n + i)+            rest <- dotV k j (i + 1)+            pure (vi * aij + rest)+    dotRhs k i+        | i >= n = pure 0+        | otherwise = do+            vi <- VUM.read mat (k * n + i)+            bi <- VUM.read rhs i+            rest <- dotRhs k (i + 1)+            pure (vi * bi + rest)++backSubstQR ::+    VUM.STVector s Double ->+    VUM.STVector s Double ->+    Int ->+    Int ->+    VUM.STVector s Double ->+    ST s ()+backSubstQR mat rhs n d x = forM_ [d - 1, d - 2 .. 0] $ \i -> do+    bi <- VUM.read rhs i+    s <- sumAbove i (i + 1) 0+    rii <- VUM.read mat (i * n + i)+    VUM.write x i ((bi - s) / rii)+  where+    sumAbove i j !acc+        | j >= d = pure acc+        | otherwise = do+            rij <- VUM.read mat (j * n + i)+            xj <- VUM.read x j+            sumAbove i (j + 1) (acc + rij * xj)++{- | Cholesky factor @L@ (lower-triangular, @A = L Lᵀ@) of a symmetric+positive-definite matrix, or 'Nothing' if a non-positive pivot is hit.+-}+cholesky :: Matrix -> Maybe Matrix+cholesky a+    | V.null a = Just V.empty+    | otherwise = runST $ do+        l <- VUM.replicate (d * d) 0+        ok <- buildL l+        if ok then Just <$> freezeLower l else pure Nothing+  where+    d = V.length a+    buildL l = go 0+      where+        go j+            | j >= d = pure True+            | otherwise = do+                s <- sumLk l j j (j - 1) 0+                let ajj = (a V.! j) VU.! j+                    diag = ajj - s+                if diag <= 0+                    then pure False+                    else do+                        let ljj = sqrt diag+                        VUM.write l (j * d + j) ljj+                        forM_ [j + 1 .. d - 1] $ \i -> do+                            sij <- sumLk l i j (j - 1) 0+                            let aij = (a V.! i) VU.! j+                            VUM.write l (i * d + j) ((aij - sij) / ljj)+                        go (j + 1)+    sumLk l i j k !acc+        | k < 0 = pure acc+        | otherwise = do+            lik <- VUM.read l (i * d + k)+            ljk <- VUM.read l (j * d + k)+            sumLk l i j (k - 1) (acc + lik * ljk)+    freezeLower l = do+        frozen <- VU.freeze l+        pure $ V.generate d $ \i -> VU.slice (i * d) d frozen++-- | Solve @L y = b@ for lower-triangular @L@.+forwardSubst :: Matrix -> VU.Vector Double -> VU.Vector Double+forwardSubst l b = runST $ do+    y <- VUM.new d+    forM_ [0 .. d - 1] $ \i -> do+        let row = l V.! i+        s <- sumKnown y row i 0 0+        VUM.write y i ((b VU.! i - s) / (row VU.! i))+    VU.freeze y+  where+    d = V.length l+    sumKnown y row i j !acc+        | j >= i = pure acc+        | otherwise = do+            yj <- VUM.read y j+            sumKnown y row i (j + 1) (acc + (row VU.! j) * yj)++-- | Solve @Lᵀ x = y@ for lower-triangular @L@.+backSubst :: Matrix -> VU.Vector Double -> VU.Vector Double+backSubst l y = runST $ do+    x <- VUM.new d+    forM_ [d - 1, d - 2 .. 0] $ \i -> do+        s <- sumKnown x i (i + 1) 0+        VUM.write x i ((y VU.! i - s) / ((l V.! i) VU.! i))+    VU.freeze x+  where+    d = V.length l+    sumKnown x i j !acc+        | j >= d = pure acc+        | otherwise = do+            xj <- VUM.read x j+            sumKnown x i (j + 1) (acc + ((l V.! j) VU.! i) * xj)++{- | Solve the SPD system @A x = b@ via Cholesky; 'Nothing' when @A@ is not+positive-definite.+-}+choleskySolve :: Matrix -> VU.Vector Double -> Maybe (VU.Vector Double)+choleskySolve a b = do+    l <- cholesky a+    pure (backSubst l (forwardSubst l b))++-- | @log det A = 2 Σ log Lᵢᵢ@ from a Cholesky factor @L@.+logDetFromChol :: Matrix -> Double+logDetFromChol l = 2 * V.sum (V.imap (\i row -> log (row VU.! i)) l)

@@ -0,0 +1,11 @@+{- | Linear models for the dataframe ecosystem: regression (OLS, ridge, lasso,+elastic net) and one-vs-rest logistic classification. Re-exports the focused+submodules.+-}+module DataFrame.LinearModel (+    module DataFrame.LinearModel.Regression,+    module DataFrame.LinearModel.Logistic,+) where++import DataFrame.LinearModel.Logistic+import DataFrame.LinearModel.Regression

@@ -0,0 +1,105 @@+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE UndecidableInstances #-}++{- | Logistic regression: binary and one-vs-rest multiclass over the FISTA+solver. 'fit' trains a 'LogisticModel'; 'predict' is the arg-max class decision.+Per-class margins and (normalized) probabilities stay available as named+auxiliary expressions.+-}+module DataFrame.LinearModel.Logistic (+    LogisticConfig (..),+    defaultLogisticConfig,+    LogisticModel (..),+    logisticMarginExprs,+    logisticProbExprs,+) where++import Data.List (sort)+import qualified Data.Map.Strict as M+import qualified Data.Vector as V+import qualified Data.Vector.Unboxed as VU++import DataFrame.Featurize.Internal (+    affineExpr,+    argMaxExpr,+    featureNames,+    numericMatrix,+    targetValues,+ )+import qualified DataFrame.Functions as F+import DataFrame.Internal.Column (Columnable)+import DataFrame.Internal.DataFrame (DataFrame)+import DataFrame.Internal.Expression (Expr)+import DataFrame.LinearSolver (+    LinearModel (..),+    SolverConfig,+    defaultSolverConfig,+    fitL1Logistic,+ )+import DataFrame.Model (Fit (..), Predict (..))++-- | Hyperparameters for logistic regression (the FISTA solver config).+newtype LogisticConfig = LogisticConfig {lgSolver :: SolverConfig}+    deriving (Eq, Show)++defaultLogisticConfig :: LogisticConfig+defaultLogisticConfig = LogisticConfig defaultSolverConfig++{- | A fitted (one-vs-rest) logistic model: parallel vectors of class labels and+their binary sub-models. 'lgModels' carries sklearn's per-class @coef_@.+-}+data LogisticModel a = LogisticModel+    { lgClasses :: !(V.Vector a)+    , lgModels :: !(V.Vector LinearModel)+    }+    deriving (Eq, Show)++instance (Columnable a, Ord a) => Fit LogisticConfig (Expr a) (LogisticModel a) where+    fit = fitLogistic++instance (Columnable a, Ord a) => Predict (LogisticModel a) a where+    predict m = argMaxExpr (labelledMargins m)++-- | Fit one-vs-rest logistic regression; the target column supplies the classes.+fitLogistic ::+    (Columnable a, Ord a) =>+    LogisticConfig -> Expr a -> DataFrame -> LogisticModel a+fitLogistic (LogisticConfig cfg) target df =+    LogisticModel (V.fromList classes) (V.fromList (map fitOne classes))+  where+    names = featureNames target df+    (nameVec, mat) = numericMatrix names df+    ys = targetValues target df+    classes = sort (foldr dedup [] (V.toList ys))+    dedup x acc = if x `elem` acc then acc else x : acc+    fitOne c =+        let labels =+                VU.generate (V.length ys) (\i -> if ys V.! i == c then 1 else -1)+         in fitL1Logistic cfg mat labels nameVec++-- | The raw margin @Expr@ for each class.+logisticMarginExprs ::+    (Columnable a, Ord a) => LogisticModel a -> M.Map a (Expr Double)+logisticMarginExprs m = M.fromList (labelledMargins m)++-- | Per-class probability expressions: @1 / (1 + exp(-margin))@.+logisticProbExprs ::+    (Columnable a, Ord a) => LogisticModel a -> M.Map a (Expr Double)+logisticProbExprs = M.map sigmoidExpr . logisticMarginExprs++labelledMargins :: LogisticModel a -> [(a, Expr Double)]+labelledMargins m =+    [ (lgClasses m V.! i, marginOf (lgModels m V.! i))+    | i <- [0 .. V.length (lgClasses m) - 1]+    ]++marginOf :: LinearModel -> Expr Double+marginOf m =+    affineExpr+        (lmIntercept m)+        (zip (VU.toList (lmWeights m)) (V.toList (lmFeatureNames m)))++sigmoidExpr :: Expr Double -> Expr Double+sigmoidExpr z = F.lit 1 / (F.lit 1 + exp (negate z))

@@ -0,0 +1,128 @@+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE MultiParamTypeClasses #-}++{- | Linear regression with the standard penalties: ordinary least squares+(Householder QR), ridge (Cholesky on the regularized normal equations), and+lasso / elastic net (FISTA). 'fit' produces a 'LinearRegressor' record;+'predict' compiles it to an @Expr Double@ over the raw feature columns.+-}+module DataFrame.LinearModel.Regression (+    Penalty (..),+    LinearConfig (..),+    defaultLinearConfig,+    LinearRegressor (..),+    predictLinear,+) where++import qualified Data.Text as T+import qualified Data.Vector as V+import qualified Data.Vector.Unboxed as VU++import DataFrame.Featurize.Internal (+    affineExpr,+    featureNames,+    numericMatrix,+    targetDoubles,+ )+import DataFrame.Internal.Expression (Expr)+import DataFrame.LinearAlgebra (Matrix, dot, gram, tMatVec)+import DataFrame.LinearAlgebra.Solve (choleskySolve, qrLeastSquares)+import DataFrame.LinearSolver (+    LinearModel (..),+    SolverConfig (..),+    defaultSolverConfig,+    fitProx,+ )+import DataFrame.LinearSolver.Loss (squaredLoss)+import DataFrame.Model (Fit (..), Predict (..))++-- | Regularization choice. @alpha@ is the penalty strength; @l1Ratio@ mixes L1/L2.+data Penalty+    = OLS+    | Ridge !Double+    | Lasso !Double+    | ElasticNet !Double !Double+    deriving (Eq, Show)++-- | Hyperparameters for linear regression: the penalty and the FISTA solver config.+data LinearConfig = LinearConfig+    { lcPenalty :: !Penalty+    , lcSolver :: !SolverConfig+    }+    deriving (Eq, Show)++defaultLinearConfig :: LinearConfig+defaultLinearConfig = LinearConfig{lcPenalty = OLS, lcSolver = defaultSolverConfig}++{- | A fitted linear regressor. @regCoef@ and @regIntercept@ are sklearn's+@coef_@ / @intercept_@ in raw feature space.+-}+data LinearRegressor = LinearRegressor+    { regCoef :: !(VU.Vector Double)+    , regIntercept :: !Double+    , regFeatureNames :: !(V.Vector T.Text)+    , regPenalty :: !Penalty+    }+    deriving (Eq, Show)++instance Fit LinearConfig (Expr Double) LinearRegressor where+    fit (LinearConfig penalty cfg) target df =+        case penalty of+            OLS -> closedForm (olsSolve mat y)+            Ridge alpha -> closedForm (ridgeSolve alpha mat y)+            Lasso alpha -> proxFit alpha 1.0+            ElasticNet alpha l1r -> proxFit alpha l1r+      where+        names = featureNames target df+        (nameVec, mat) = numericMatrix names df+        y = targetDoubles target df+        closedForm (coef, intercept) =+            LinearRegressor coef intercept nameVec penalty+        proxFit alpha l1r =+            let proxCfg =+                    cfg{scL1Lambda = alpha * l1r, scL2Lambda = alpha * (1 - l1r)}+                m = fitProx squaredLoss proxCfg mat y nameVec+             in LinearRegressor (lmWeights m) (lmIntercept m) nameVec penalty++instance Predict LinearRegressor Double where+    predict m =+        affineExpr+            (regIntercept m)+            (zip (VU.toList (regCoef m)) (V.toList (regFeatureNames m)))++-- | OLS via QR on the intercept-augmented design matrix.+olsSolve :: Matrix -> VU.Vector Double -> (VU.Vector Double, Double)+olsSolve mat y =+    case qrLeastSquares augmented y of+        Right sol -> (VU.drop 1 sol, sol VU.! 0)+        Left _ -> ridgeSolve 1e-8 mat y+  where+    augmented = V.map (VU.cons 1) mat++{- | Ridge via Cholesky on @(XcᵀXc + αI) w = Xcᵀ yc@ over centred data; the+intercept is recovered from the column/target means.+-}+ridgeSolve :: Double -> Matrix -> VU.Vector Double -> (VU.Vector Double, Double)+ridgeSolve alpha mat y =+    case choleskySolve a rhs of+        Just w -> (w, meanY - dot w meansX)+        Nothing -> (VU.replicate d 0, meanY)+  where+    n = V.length mat+    d = if n == 0 then 0 else VU.length (V.head mat)+    meansX =+        VU.generate d $ \j ->+            sum [(mat V.! i) VU.! j | i <- [0 .. n - 1]] / fromIntegral n+    meanY = VU.sum y / fromIntegral n+    centered = V.map (\row -> VU.zipWith (-) row meansX) mat+    yc = VU.map (subtract meanY) y+    a = addDiag alpha (gram centered)+    rhs = tMatVec centered yc++-- | Add @alpha@ to the diagonal of a square matrix.+addDiag :: Double -> Matrix -> Matrix+addDiag alpha = V.imap (\i row -> row VU.// [(i, row VU.! i + alpha)])++-- | Predict the target for each row of a feature matrix.+predictLinear :: LinearRegressor -> Matrix -> VU.Vector Double+predictLinear m = VU.convert . V.map (\x -> regIntercept m + dot (regCoef m) x)

@@ -1,9 +1,12 @@ {-# LANGUAGE BangPatterns #-} {-# LANGUAGE ScopedTypeVariables #-} -{- | L1-regularized logistic regression used as the per-node split solver in-'DataFrame.DecisionTree'. Produces a sparse oblique hyperplane that can be-compiled to an 'Expr Bool' over numeric columns.+{- | Proximal-gradient (FISTA) solver for L1/L2-regularized generalized linear+models. 'fitL1Logistic' is the binary logistic split solver used by+'DataFrame.DecisionTree'; 'fitProx' generalizes it to any 'SmoothLoss'+(squared loss for lasso/elastic-net, squared hinge for LinearSVC). Features are+standardized internally and weights de-standardized, so the model applies to+raw column values. -} module DataFrame.LinearSolver (     -- * Model@@ -13,14 +16,16 @@     SolverConfig (..),     defaultSolverConfig, -    -- * Solver+    -- * Solvers     fitL1Logistic,+    fitProx,      -- * Expr conversion     modelToExpr,      -- * Internals (exposed for testing)     standardize,+    columnStats,     softThreshold,     sigmoid,     dotProduct,@@ -28,6 +33,13 @@  import qualified DataFrame.Functions as F import DataFrame.Internal.Expression (Expr (..))+import DataFrame.LinearAlgebra (Matrix, gram, scaleV)+import DataFrame.LinearAlgebra.Eigen (powerIterTop)+import DataFrame.LinearSolver.Loss (+    SmoothLoss (..),+    logisticLoss,+    sigmoid,+ ) import DataFrame.Operators ((.*.), (.+.), (.>.))  import Control.Monad.ST (ST, runST)@@ -95,7 +107,44 @@     V.Vector T.Text ->     LinearModel {-# INLINEABLE fitL1Logistic #-}-fitL1Logistic cfg rows labels featureNames+fitL1Logistic = runFista logisticLoss logisticLipschitz+  where+    logisticLipschitz _ keepN = fromIntegral (keepN + 1) / 4++{- | Fit any 'SmoothLoss' with the elastic-net proximal-gradient engine. The+Lipschitz constant uses the spectral norm of the standardized Gram matrix+(power iteration), which is tight for squared and squared-hinge losses where+the logistic trace bound would be far too loose.+-}+fitProx ::+    SmoothLoss ->+    SolverConfig ->+    V.Vector (VU.Vector Double) ->+    VU.Vector Double ->+    V.Vector T.Text ->+    LinearModel+fitProx loss = runFista loss specNormLipschitz+  where+    specNormLipschitz xKept _ =+        let n = V.length xKept+            gramN = V.map (scaleV (1 / fromIntegral n)) (gram xKept)+            (specNorm, _) = powerIterTop 50 gramN+         in slCurvBound loss * (specNorm + 1)++{- | Shared FISTA scaffolding: standardize, drop near-constant columns, run the+inner loop, de-standardize. @lipschitzOf@ receives the standardized kept-feature+matrix and the number of kept columns and returns the smooth-part Lipschitz+bound; the L2 penalty contribution @λ₂@ is added here.+-}+runFista ::+    SmoothLoss ->+    (Matrix -> Int -> Double) ->+    SolverConfig ->+    V.Vector (VU.Vector Double) ->+    VU.Vector Double ->+    V.Vector T.Text ->+    LinearModel+runFista loss lipschitzOf cfg rows labels featureNames     | n == 0 || d == 0 = zeroModel     | otherwise =         let (!means, !stds, !variances) = columnStats rows@@ -106,14 +155,11 @@                     let !meansKept = gatherBy keep means                         !stdsKept = gatherBy keep stds                         !xKept = V.map (standardizeRowKept keep means stds) rows-                        -- Elastic-Net Lipschitz: standard logistic bound-                        -- @(d+1)/4@ plus the L2 part's Hessian-norm-                        -- contribution @λ₂·I@ (operator norm @λ₂@).                         !lipschitz =-                            fromIntegral (VU.length keep + 1) / 4-                                + scL2Lambda cfg+                            lipschitzOf xKept (VU.length keep) + scL2Lambda cfg                         (!wStdKept, !bStd) =                             fistaLoop+                                loss                                 (scL1Lambda cfg)                                 (scL2Lambda cfg)                                 lipschitz@@ -270,48 +316,42 @@     | v < -lambda = v + lambda     | otherwise = 0 --- | Numerically stable logistic sigmoid.-sigmoid :: Double -> Double-sigmoid z-    | z >= 0 = 1 / (1 + exp (-z))-    | otherwise = let ez = exp z in ez / (1 + ez)- {- | Dot product of two unboxed vectors. Caller must ensure equal length; lengths are not checked. -} dotProduct :: VU.Vector Double -> VU.Vector Double -> Double dotProduct u v = VU.sum (VU.zipWith (*) u v) -{- | Gradient of the average binary logistic loss at @(w, b)@ for labels in-@{\-1,+1}@. Returns @(gradW, gradB)@.+{- | Gradient of the average loss at @(w, b)@. Returns @(gradW, gradB)@.  Sample-weighted variant: when @sampleWeights@ is @Just ws@ the per-row contribution is multiplied by @ws[i]@. With weights of mean 1 (i.e. @Σ w_i = N@; the class-balanced convention used by 'fitLinearCandidate'), the @1/N@ normalisation is preserved exactly. -}-logisticGradient ::+lossGradient ::+    SmoothLoss ->     Maybe (VU.Vector Double) ->     V.Vector (VU.Vector Double) ->     VU.Vector Double ->     VU.Vector Double ->     Double ->     (VU.Vector Double, Double)-logisticGradient sampleWeights features labels w b = (gradW, gradB)+lossGradient loss sampleWeights features labels w b = (gradW, gradB)   where     !invN = 1 / fromIntegral (V.length features)-    !coeffs = rowCoeffs sampleWeights features labels w b invN+    !coeffs = rowCoeffs loss sampleWeights features labels w b invN     !gradW = accumulateGradW (VU.length w) features coeffs     !gradB = VU.sum coeffs -{- | Per-row loss coefficient. Without sample weights:-@c_i = -y_i * sigmoid(-y_i * margin_i) / N@. With @Just ws@, each row's-contribution is additionally multiplied by @ws[i]@.+{- | Per-row loss coefficient @c_i = ℓ'(y_i, z_i) / N@ at margin+@z_i = w·x_i + b@, optionally scaled by @ws[i]@.  unsafeIndex is safe: @i@ ranges over @[0,n-1]@ and @labels@ / @sampleWeights@ both have length @n@ by construction. -} rowCoeffs ::+    SmoothLoss ->     Maybe (VU.Vector Double) ->     V.Vector (VU.Vector Double) ->     VU.Vector Double ->@@ -319,12 +359,12 @@     Double ->     Double ->     VU.Vector Double-rowCoeffs sampleWeights features labels w b invN =+rowCoeffs loss sampleWeights features labels w b invN =     VU.generate (V.length features) $ \i ->         let !yi = VU.unsafeIndex labels i             !row = V.unsafeIndex features i-            !margin = yi * (dotProduct w row + b)-            !base = -(yi * sigmoid (-margin) * invN)+            !z = dotProduct w row + b+            !base = slGradZ loss yi z * invN          in case sampleWeights of                 Nothing -> base                 Just ws -> base * VU.unsafeIndex ws i@@ -350,6 +390,7 @@ @prox(z) = softThreshold(z, λ₁/lp) / (1 + λ₂/lp)@. -} fistaLoop ::+    SmoothLoss ->     Double ->     Double ->     Double ->@@ -361,13 +402,13 @@     VU.Vector Double ->     Double ->     (VU.Vector Double, Double)-fistaLoop lambda1 lambda2 lp maxIter tol sampleWeights features labels w0 b0 =+fistaLoop loss lambda1 lambda2 lp maxIter tol sampleWeights features labels w0 b0 =     go 0 w0 b0 w0 b0 1.0   where     !shrink = lambda1 / lp     !ridgeDenom = 1 + lambda2 / lp     !stepInv = 1 / lp-    proxStep = fistaProxStep sampleWeights features labels shrink ridgeDenom stepInv+    proxStep = fistaProxStep loss sampleWeights features labels shrink ridgeDenom stepInv     go !iter !xWPrev !xBPrev !yW !yB !t         | iter >= maxIter = (xWPrev, xBPrev)         | iter > 0 && delta < tol = (xW, xB)@@ -386,6 +427,7 @@ Teboulle 2009 §4). The intercept is unregularised (no L1 or L2 applied). -} fistaProxStep ::+    SmoothLoss ->     Maybe (VU.Vector Double) ->     V.Vector (VU.Vector Double) ->     VU.Vector Double ->@@ -395,8 +437,8 @@     VU.Vector Double ->     Double ->     (VU.Vector Double, Double)-fistaProxStep sampleWeights features labels shrink ridgeDenom stepInv yW yB =-    let (gW, gB) = logisticGradient sampleWeights features labels yW yB+fistaProxStep loss sampleWeights features labels shrink ridgeDenom stepInv yW yB =+    let (gW, gB) = lossGradient loss sampleWeights features labels yW yB         !wNew =             VU.zipWith                 (\yi gi -> softThreshold shrink (yi - gi * stepInv) / ridgeDenom)

@@ -0,0 +1,47 @@+{-# LANGUAGE OverloadedStrings #-}++{- | Smooth losses for the proximal-gradient engine. Each carries its+derivative @∂ℓ/∂z@ at @z = w·x + b@ and a global bound on the curvature+@∂²ℓ/∂z²@ (used for the FISTA step size).+-}+module DataFrame.LinearSolver.Loss (+    SmoothLoss (..),+    sigmoid,+    logisticLoss,+    squaredLoss,+    sqHingeLoss,+) where++import qualified Data.Text as T++{- | A convex, @C¹@ per-sample loss @ℓ(y, z)@. 'slGradZ' is @∂ℓ/∂z@;+'slCurvBound' bounds @∂²ℓ/∂z²@ over all @(y, z)@.+-}+data SmoothLoss = SmoothLoss+    { slName :: !T.Text+    , slGradZ :: Double -> Double -> Double+    , slCurvBound :: !Double+    }++-- | Numerically stable logistic sigmoid.+sigmoid :: Double -> Double+sigmoid z+    | z >= 0 = 1 / (1 + exp (-z))+    | otherwise = let ez = exp z in ez / (1 + ez)++-- | Binary logistic loss for labels in @{\-1,+1}@: @ℓ = log(1 + exp(-y z))@.+logisticLoss :: SmoothLoss+logisticLoss =+    SmoothLoss "logistic" (\y z -> negate (y * sigmoid (negate (y * z)))) 0.25++-- | Squared error for regression: @ℓ = ½ (z - y)²@.+squaredLoss :: SmoothLoss+squaredLoss = SmoothLoss "squared" (flip (-)) 1.0++-- | Squared hinge for classification (LinearSVC default), labels @{\-1,+1}@.+sqHingeLoss :: SmoothLoss+sqHingeLoss =+    SmoothLoss+        "squared_hinge"+        (\y z -> let m = 1 - y * z in if m > 0 then negate (2 * y * m) else 0)+        2.0

@@ -0,0 +1,237 @@+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE TypeApplications #-}++{- | Evaluation metrics for fitted models. The everyday entry point is+'evaluate', which applies a model's prediction expression and a truth column to+a frame and folds a metric — no manual @interpret@/extract plumbing. Metrics are+plain functions (@type Metric = Vector -> Vector -> Double@), so you pass @mse@+or @accuracy@ directly. Classification metrics handle multiclass via 'Average';+'classificationReport' / 'regressionReport' bundle the common numbers with a+scikit-learn-style 'Show'.+-}+module DataFrame.Metrics (+    -- * Metric type + evaluation+    Metric,+    evaluate,+    predictColumn,+    columnOf,++    -- * Regression metrics+    mse,+    rmse,+    mae,+    r2,++    -- * Classification metrics+    accuracy,+    logLoss,+    Average (..),+    precision,+    recall,+    f1,+    rocAuc,++    -- * Per-class helpers (for reports)+    classCounts,+    precOf,+    recOf,+    f1Of,+) where++import Data.Either (fromRight)+import Data.List (nub, sort, sortBy)+import Data.Ord (comparing)+import qualified Data.Text as T+import qualified Data.Vector.Unboxed as VU++import DataFrame.Internal.Column (TypedColumn (..), toVector)+import DataFrame.Internal.DataFrame (DataFrame)+import DataFrame.Internal.Expression (Expr)+import DataFrame.Internal.Interpreter (interpret)+import DataFrame.Operations.Transformations (derive)++-- | A metric maps predictions and ground truth to a scalar score.+type Metric = VU.Vector Double -> VU.Vector Double -> Double++{- | Evaluate a model's prediction expression against a truth column on a frame.++> evaluate rmse (linearExpr model) (F.col @Double "target") df+> evaluate accuracy (logisticDecisionExpr model) (F.col @Double "label") df+-}+evaluate :: Metric -> Expr Double -> Expr Double -> DataFrame -> Double+evaluate metric predExpr truthExpr df =+    metric (columnOf df predExpr) (columnOf df truthExpr)++-- | Add a model's prediction expression to a frame as a named column.+predictColumn :: T.Text -> Expr Double -> DataFrame -> DataFrame+predictColumn = derive++-- | Interpret an expression to a @Double@ vector over a frame.+columnOf :: DataFrame -> Expr Double -> VU.Vector Double+columnOf df e = case interpret @Double df e of+    Right (TColumn c) -> fromRight VU.empty (toVector @Double @VU.Vector c)+    Left err -> error (show err)++n2 :: VU.Vector Double -> Double+n2 = fromIntegral . VU.length++-- | Mean squared error.+mse :: Metric+mse preds truth+    | VU.null truth = 0+    | otherwise =+        VU.sum (VU.zipWith (\p t -> (p - t) ^ (2 :: Int)) preds truth) / n2 truth++-- | Root mean squared error.+rmse :: Metric+rmse preds truth = sqrt (mse preds truth)++-- | Mean absolute error.+mae :: Metric+mae preds truth+    | VU.null truth = 0+    | otherwise = VU.sum (VU.zipWith (\p t -> abs (p - t)) preds truth) / n2 truth++-- | Coefficient of determination @R²@.+r2 :: Metric+r2 preds truth+    | VU.null truth || ssTot == 0 = 0+    | otherwise = 1 - ssRes / ssTot+  where+    mean = VU.sum truth / n2 truth+    ssRes = VU.sum (VU.zipWith (\p t -> (t - p) ^ (2 :: Int)) preds truth)+    ssTot = VU.sum (VU.map (\t -> (t - mean) ^ (2 :: Int)) truth)++-- | Fraction of exact matches.+accuracy :: Metric+accuracy preds truth+    | VU.null truth = 0+    | otherwise =+        fromIntegral (VU.length (VU.filter id (VU.zipWith (==) preds truth))) / n2 truth++-- | Binary log loss; probabilities clamped away from @0@/@1@.+logLoss :: Metric+logLoss probs truth+    | VU.null truth = 0+    | otherwise =+        negate+            ( VU.sum+                ( VU.zipWith+                    (\p y -> let q = clampP p in y * log q + (1 - y) * log (1 - q))+                    probs+                    truth+                )+            )+            / n2 truth+  where+    clampP p = max 1e-15 (min (1 - 1e-15) p)++-- | Averaging strategy for multiclass precision/recall/F1.+data Average+    = -- | one class is positive; the rest negative+      Binary Double+    | -- | unweighted mean over classes+      Macro+    | -- | pool per-class counts (equals accuracy for single-label)+      Micro+    | -- | support-weighted mean over classes+      Weighted+    deriving (Eq, Show)++-- | Per-class @(tp, fp, fn, support)@ over the class set of @truth ∪ preds@.+classCounts ::+    VU.Vector Double -> VU.Vector Double -> [(Double, (Int, Int, Int, Int))]+classCounts preds truth =+    [(c, countsFor c) | c <- classes]+  where+    classes = sort (nub (VU.toList truth ++ VU.toList preds))+    countsFor c =+        VU.foldl'+            ( \(tp, fp, fn, sup) (p, y) ->+                ( if p == c && y == c then tp + 1 else tp+                , if p == c && y /= c then fp + 1 else fp+                , if p /= c && y == c then fn + 1 else fn+                , if y == c then sup + 1 else sup+                )+            )+            (0, 0, 0, 0)+            (VU.zip preds truth)++safeDiv :: Int -> Int -> Double+safeDiv a b = if b == 0 then 0 else fromIntegral a / fromIntegral b++precOf, recOf :: (Int, Int, Int, Int) -> Double+precOf (tp, fp, _, _) = safeDiv tp (tp + fp)+recOf (tp, _, fn, _) = safeDiv tp (tp + fn)++f1Of :: (Int, Int, Int, Int) -> Double+f1Of cs =+    let p = precOf cs; r = recOf cs in if p + r == 0 then 0 else 2 * p * r / (p + r)++averaged ::+    ((Int, Int, Int, Int) -> Double) ->+    Average ->+    VU.Vector Double ->+    VU.Vector Double ->+    Double+averaged stat avg preds truth =+    case avg of+        Binary pos -> maybe 0 stat (lookup pos cc)+        Macro -> meanOf [stat c | (_, c) <- cc]+        Weighted ->+            let total = sum [sup | (_, (_, _, _, sup)) <- cc]+             in if total == 0+                    then 0+                    else+                        sum [fromIntegral sup * stat c | (_, c@(_, _, _, sup)) <- cc]+                            / fromIntegral total+        Micro ->+            let tp = sum [t | (_, (t, _, _, _)) <- cc]+                fp = sum [x | (_, (_, x, _, _)) <- cc]+                fn = sum [x | (_, (_, _, x, _)) <- cc]+             in stat (tp, fp, fn, 0)+  where+    cc = classCounts preds truth+    meanOf xs = if null xs then 0 else sum xs / fromIntegral (length xs)++-- | Precision with the given averaging.+precision :: Average -> VU.Vector Double -> VU.Vector Double -> Double+precision = averaged precOf++-- | Recall with the given averaging.+recall :: Average -> VU.Vector Double -> VU.Vector Double -> Double+recall = averaged recOf++-- | F1 with the given averaging.+f1 :: Average -> VU.Vector Double -> VU.Vector Double -> Double+f1 = averaged f1Of++{- | Binary ROC-AUC (Mann–Whitney). @scores@ are predicted probabilities, @truth@+is @0@/@1@.+-}+rocAuc :: VU.Vector Double -> VU.Vector Double -> Double+rocAuc scores truth+    | nPos == 0 || nNeg == 0 = 0.5+    | otherwise = (rankSum - nPos * (nPos + 1) / 2) / (nPos * nNeg)+  where+    ranked = rankAverages (VU.toList scores)+    pairs = zip (VU.toList truth) ranked+    rankSum = sum [r | (y, r) <- pairs, y == 1]+    nPos = fromIntegral (length (filter (== 1) (VU.toList truth)))+    nNeg = fromIntegral (VU.length truth) - nPos++-- | Average ranks (ties share the mean rank), returned in input order.+rankAverages :: [Double] -> [Double]+rankAverages xs =+    let indexed = zip [0 :: Int ..] xs+        sorted = sortBy (comparing snd) indexed+        ranked = assignRanks 1 sorted+     in map snd (sortBy (comparing fst) ranked)+  where+    assignRanks _ [] = []+    assignRanks start grp =+        let v = snd (head grp)+            (tied, rest) = span ((== v) . snd) grp+            k = length tied+            avgRank = fromIntegral (sum [start .. start + k - 1]) / fromIntegral k+         in [(i, avgRank) | (i, _) <- tied] ++ assignRanks (start + k) rest

@@ -0,0 +1,170 @@+{- | Bundled, pretty-printing evaluation summaries: a labelled confusion matrix+and scikit-learn-style regression / classification reports. The @*Expr@ variants+take a model's prediction expression and a truth column directly, so a full+report is a one-liner after fitting.+-}+module DataFrame.Metrics.Report (+    ConfusionMatrix (..),+    confusionMatrix,+    confusionMatrixExpr,+    RegressionReport (..),+    regressionReport,+    regressionReportExpr,+    ClassStats (..),+    ClassificationReport (..),+    classificationReport,+    classificationReportExpr,+) where++import Data.List (nub, sort, sortBy)+import Data.Ord (Down (..), comparing)+import qualified Data.Vector.Unboxed as VU++import DataFrame.Internal.DataFrame (DataFrame)+import DataFrame.Internal.Expression (Expr)+import DataFrame.Metrics (+    Average (..),+    accuracy,+    classCounts,+    columnOf,+    f1,+    f1Of,+    mae,+    mse,+    precOf,+    r2,+    recOf,+    rmse,+ )++-- | A labelled confusion matrix: class order plus row-major @actual×predicted@.+data ConfusionMatrix = ConfusionMatrix+    { cmClasses :: ![Double]+    , cmCounts :: ![[Int]]+    }+    deriving (Eq)++-- | Confusion matrix over the class set of @truth ∪ preds@.+confusionMatrix :: VU.Vector Double -> VU.Vector Double -> ConfusionMatrix+confusionMatrix preds truth = ConfusionMatrix classes counts+  where+    classes = sort (nub (VU.toList truth ++ VU.toList preds))+    counts =+        [ [ VU.length (VU.filter id (VU.zipWith (\p t -> t == a && p == c) preds truth))+          | c <- classes+          ]+        | a <- classes+        ]++-- | Confusion matrix from a prediction expression and a truth column.+confusionMatrixExpr ::+    Expr Double -> Expr Double -> DataFrame -> ConfusionMatrix+confusionMatrixExpr predExpr truthExpr df =+    confusionMatrix (columnOf df predExpr) (columnOf df truthExpr)++instance Show ConfusionMatrix where+    show (ConfusionMatrix classes counts) =+        unlines (header : rows)+      where+        lbls = map show classes+        w = maximum (8 : map length lbls) + 2+        cell s = replicate (max 1 (w - length s)) ' ' ++ s+        header = cell "a\\p" ++ concatMap cell lbls+        rows = [cell a ++ concatMap (cell . show) row | (a, row) <- zip lbls counts]++-- | Regression metrics bundle.+data RegressionReport = RegressionReport+    { rrMSE :: !Double+    , rrRMSE :: !Double+    , rrMAE :: !Double+    , rrR2 :: !Double+    }+    deriving (Eq)++instance Show RegressionReport where+    show r =+        unlines+            [ "Regression report"+            , "  mse  = " ++ show (rrMSE r)+            , "  rmse = " ++ show (rrRMSE r)+            , "  mae  = " ++ show (rrMAE r)+            , "  r2   = " ++ show (rrR2 r)+            ]++-- | Regression report from prediction/truth vectors.+regressionReport :: VU.Vector Double -> VU.Vector Double -> RegressionReport+regressionReport preds truth =+    RegressionReport+        (mse preds truth)+        (rmse preds truth)+        (mae preds truth)+        (r2 preds truth)++-- | Regression report from a prediction expression and a truth column.+regressionReportExpr ::+    Expr Double -> Expr Double -> DataFrame -> RegressionReport+regressionReportExpr predExpr truthExpr df =+    regressionReport (columnOf df predExpr) (columnOf df truthExpr)++-- | Per-class precision/recall/F1/support.+data ClassStats = ClassStats+    { csPrecision :: !Double+    , csRecall :: !Double+    , csF1 :: !Double+    , csSupport :: !Int+    }+    deriving (Eq, Show)++{- | A scikit-learn-style classification report: per-class stats plus accuracy+and macro/weighted F1.+-}+data ClassificationReport = ClassificationReport+    { crPerClass :: ![(Double, ClassStats)]+    , crAccuracy :: !Double+    , crMacroF1 :: !Double+    , crWeightedF1 :: !Double+    }+    deriving (Eq)++instance Show ClassificationReport where+    show r =+        unlines $+            (pad "class" ++ pad "precision" ++ pad "recall" ++ pad "f1" ++ pad "support")+                : [ pad (show c)+                        ++ pad (num (csPrecision s))+                        ++ pad (num (csRecall s))+                        ++ pad (num (csF1 s))+                        ++ pad (show (csSupport s))+                  | (c, s) <- crPerClass r+                  ]+                ++ [ ""+                   , "accuracy    = " ++ num (crAccuracy r)+                   , "macro f1    = " ++ num (crMacroF1 r)+                   , "weighted f1 = " ++ num (crWeightedF1 r)+                   ]+      where+        pad s = let w = 12 in s ++ replicate (max 1 (w - length s)) ' '+        num x = show (fromIntegral (round (x * 1000) :: Int) / 1000 :: Double)++-- | Classification report from prediction/truth vectors.+classificationReport ::+    VU.Vector Double -> VU.Vector Double -> ClassificationReport+classificationReport preds truth =+    ClassificationReport+        perClass+        (accuracy preds truth)+        (f1 Macro preds truth)+        (f1 Weighted preds truth)+  where+    perClass =+        sortBy (comparing (Down . csSupport . snd)) $+            [ (c, ClassStats (precOf s) (recOf s) (f1Of s) (supOf s))+            | (c, s) <- classCounts preds truth+            ]+    supOf (_, _, _, sup) = sup++-- | Classification report from a prediction expression and a truth column.+classificationReportExpr ::+    Expr Double -> Expr Double -> DataFrame -> ClassificationReport+classificationReportExpr predExpr truthExpr df =+    classificationReport (columnOf df predExpr) (columnOf df truthExpr)

@@ -0,0 +1,69 @@+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE FunctionalDependencies #-}++{- | The two verbs every model speaks. Instead of a per-model @fitX@ / @xExpr@+zoo, every estimator is an instance of these classes:++  * 'fit' trains a model from a hyperparameter config, an @input@ (the supervised+    target @Expr a@ or the unsupervised feature list @[Expr Double]@), and a frame.+  * 'predict' compiles the model's canonical prediction to an @Expr@ over the raw+    columns (regressors give @Expr Double@, classifiers @Expr a@, clusterers+    @Expr Int@). Models with no honest out-of-sample prediction (e.g. DBSCAN) simply+    have no instance — @predict@ on them is a compile error, not a fake.++Every prediction lands in the /same/ expression type, so a fitted model composes+with 'DataFrame.Operations.Transformations.derive', the 'DataFrame.Transform'+monoid, and 'DataFrame.Transform.compileThrough' with no per-model glue.++Auxiliary outputs (class probabilities, per-cluster distances, component+loadings, the @*Transform@ pipeline pieces) keep their own descriptive+functions — they are not the one canonical prediction, so they are not forced+through 'predict'.++Supervised 'fit' treats every non-target column as a feature; use+'selectFeatures' to restrict to an explicit set when the frame also carries ids,+timestamps, or a second candidate target.+-}+module DataFrame.Model (+    Fit (..),+    Predict (..),+    selectFeatures,+) where++import qualified Data.Text as T++import DataFrame.Internal.DataFrame (DataFrame)+import DataFrame.Internal.Expression (Expr (..))+import DataFrame.Operations.Subset (select)++{- | Train a model. @cfg@ is the hyperparameter config; @input@ is the supervised+target @Expr a@ or the unsupervised feature list @[Expr Double]@. The config and+input together determine the model, so @fit cfg target df@ needs no annotation+(a classifier's label type comes from its @Expr a@ target).+-}+class Fit cfg input model | cfg input -> model where+    fit :: cfg -> input -> DataFrame -> model++{- | Compile a fitted model's canonical prediction to an expression over the raw+columns. The result type @r@ is determined by the model.+-}+class Predict model r | model -> r where+    predict :: model -> Expr r++{- | Restrict @df@ to exactly the named feature columns plus the supervised+target (when the target is a column), so a following 'fit' trains on those+features only.++Supervised 'fit' otherwise uses /every/ non-target column as a feature —+convenient on a clean frame, a leakage hazard when the frame carries ids,+timestamps, or a second candidate target. This mirrors the explicit+@[Expr Double]@ feature list the unsupervised fitters already take:++> model = fit defaultLinearConfig target (selectFeatures ["age", "income"] target df)+-}+selectFeatures :: [T.Text] -> Expr a -> DataFrame -> DataFrame+selectFeatures cols target = select (cols ++ targetCols target)+  where+    targetCols :: Expr b -> [T.Text]+    targetCols (Col n) = [n]+    targetCols _ = []

@@ -0,0 +1,91 @@+{- | Cross-validation and grid search for hyperparameter tuning. The model+fitters have heterogeneous types, so these helpers are parameterized by a+user-supplied @train -> test -> score@ closure; the search maximizes the mean+cross-validated score (use a negated error metric to minimize). Splitting reuses+the deterministic 'kFolds' / 'randomSplit' from @dataframe-operations@.+-}+module DataFrame.ModelSelection (+    trainTestSplit,+    crossValScore,+    crossValidate,+    GridSearchResult (..),+    gridSearch,+) where++import Data.List (maximumBy)+import Data.Ord (comparing)+import System.Random (mkStdGen)++import DataFrame.Internal.DataFrame (DataFrame)+import DataFrame.Internal.Expression (Expr)+import DataFrame.Metrics (Metric, evaluate)+import DataFrame.Operations.Merge ()+import DataFrame.Operations.Subset (kFolds, randomSplit)++{- | Split into @(train, test)@ with the given training fraction and seed+(deterministic).+-}+trainTestSplit :: Double -> Int -> DataFrame -> (DataFrame, DataFrame)+trainTestSplit trainFrac seed = randomSplit (mkStdGen seed) trainFrac++{- | Per-fold scores from k-fold cross-validation. @scoreFn train test@ fits on+the training rows and returns a score on the held-out fold.+-}+crossValScore ::+    Int -> Int -> (DataFrame -> DataFrame -> Double) -> DataFrame -> [Double]+crossValScore folds seed scoreFn df =+    [ scoreFn (combine (others i)) (fs !! i)+    | i <- [0 .. length fs - 1]+    , not (null (others i))+    ]+  where+    fs = kFolds (mkStdGen seed) folds df+    others i = [f | (j, f) <- zip [0 ..] fs, j /= i]+    combine = foldr1 (<>)++{- | scikit-learn @cross_val_score@: fit a model on each training fold and score+its prediction expression against a truth column on the held-out fold.++@fitPredict train@ fits on the training frame and returns the prediction+expression; @truth@ is the target column. Returns the per-fold metric values.++> crossValidate 5 0 rmse (F.col @Double "target")+>   (\tr -> predict (fit defaultLinearConfig (F.col @Double "target") tr)) df+-}+crossValidate ::+    Int ->+    Int ->+    Metric ->+    Expr Double ->+    (DataFrame -> Expr Double) ->+    DataFrame ->+    [Double]+crossValidate folds seed metric truth fitPredict =+    crossValScore folds seed score+  where+    score train = evaluate metric (fitPredict train) truth++-- | The outcome of a grid search: the best config, its score, and all results.+data GridSearchResult c = GridSearchResult+    { gsBest :: !c+    , gsBestScore :: !Double+    , gsAll :: ![(c, Double)]+    }+    deriving (Show)++{- | Search configurations by mean cross-validated score, returning the+maximizer. @scoreFn cfg train test@ fits @cfg@ on @train@ and scores on @test@.+-}+gridSearch ::+    Int ->+    Int ->+    [c] ->+    (c -> DataFrame -> DataFrame -> Double) ->+    DataFrame ->+    GridSearchResult c+gridSearch folds seed configs scoreFn df =+    GridSearchResult bestC bestS scored+  where+    scored = [(c, mean (crossValScore folds seed (scoreFn c) df)) | c <- configs]+    (bestC, bestS) = maximumBy (comparing snd) scored+    mean xs = if null xs then -(1 / 0) else sum xs / fromIntegral (length xs)

@@ -0,0 +1,131 @@+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE OverloadedStrings #-}++{- | Principal component analysis via the symmetric Jacobi eigensolver on the+covariance of the (optionally standardized) feature columns. 'fit' trains a+'PCAModel' (components + explained variance); the projection is exposed as+'pcaExprs' / 'pcaTransform' (PCA is a transformer, so it has no 'Predict').+-}+module DataFrame.PCA (+    NComponents (..),+    PCAConfig (..),+    defaultPCAConfig,+    PCAModel (..),+    pcaExprs,+    pcaTransform,+) where++import qualified Data.Text as T+import qualified Data.Vector as V+import qualified Data.Vector.Unboxed as VU++import DataFrame.Featurize.Internal (Features (..), extractFeatures)+import qualified DataFrame.Functions as F+import DataFrame.Internal.DataFrame (DataFrame)+import DataFrame.Internal.Expression (Expr (..), UExpr (..))+import DataFrame.LinearAlgebra (gram)+import DataFrame.LinearAlgebra.Eigen (jacobiEigenSym)+import DataFrame.Model (Fit (..))+import DataFrame.Operators ((.*.), (.+.), (.-.))+import DataFrame.Transform (Transform (..))++-- | How many components to keep.+data NComponents = NComp !Int | VarianceCovered !Double+    deriving (Eq, Show)++data PCAConfig = PCAConfig+    { pcaNComponents :: !NComponents+    , pcaStandardize :: !Bool+    }+    deriving (Eq, Show)++defaultPCAConfig :: PCAConfig+defaultPCAConfig = PCAConfig{pcaNComponents = NComp 2, pcaStandardize = False}++{- | A fitted PCA. 'pcaComponents' are sklearn's @components_@ (row @i@ is the+@i@-th loading vector); 'pcaScale' is @Just@ the per-column std when+standardizing.+-}+data PCAModel = PCAModel+    { pcaComponents :: !(V.Vector (VU.Vector Double))+    , pcaExplainedVariance :: !(VU.Vector Double)+    , pcaExplainedVarianceRatio :: !(VU.Vector Double)+    , pcaMean :: !(VU.Vector Double)+    , pcaScale :: !(Maybe (VU.Vector Double))+    , pcaFeatureNames :: !(V.Vector T.Text)+    }+    deriving (Eq, Show)++instance Fit PCAConfig [Expr Double] PCAModel where+    fit = fitPCA++-- | Fit PCA on the given feature columns (each must be a @Col@).+fitPCA :: PCAConfig -> [Expr Double] -> DataFrame -> PCAModel+fitPCA cfg features df =+    PCAModel+        { pcaComponents = V.take k vecs+        , pcaExplainedVariance = VU.take k evar+        , pcaExplainedVarianceRatio = VU.take k ratio+        , pcaMean = means+        , pcaScale = if pcaStandardize cfg then Just scales else Nothing+        , pcaFeatureNames = V.fromList names+        }+  where+    Features names cols _ n d = extractFeatures features df+    means = VU.fromList [VU.sum c / fromIntegral (max 1 n) | c <- cols]+    scales =+        VU.fromList+            [ let mu = means VU.! j+                  v = VU.sum (VU.map (\x -> (x - mu) ^ (2 :: Int)) c) / fromIntegral (max 1 n)+                  s = sqrt v+               in if s < 1e-12 then 1 else s+            | (j, c) <- zip [0 ..] cols+            ]+    scaled =+        V.generate n $ \i ->+            VU.generate d $ \j ->+                let mu = means VU.! j+                    s = if pcaStandardize cfg then scales VU.! j else 1+                 in ((cols !! j) VU.! i - mu) / s+    denom = fromIntegral (max 1 (n - 1))+    cov = V.map (VU.map (/ denom)) (gram scaled)+    (evals, vecs) = jacobiEigenSym cov+    evar = VU.map (max 0) evals+    total = VU.sum evar+    ratio = if total == 0 then evar else VU.map (/ total) evar+    k = resolveK (pcaNComponents cfg) d ratio++-- | Per-component projection expressions, named @pc1@, @pc2@, …+pcaExprs :: PCAModel -> [(T.Text, Expr Double)]+pcaExprs m =+    [ ("pc" <> T.pack (show i), componentExpr (pcaComponents m V.! (i - 1)))+    | i <- [1 .. V.length (pcaComponents m)]+    ]+  where+    names = V.toList (pcaFeatureNames m)+    means = VU.toList (pcaMean m)+    scales = maybe (repeat 1) VU.toList (pcaScale m)+    componentExpr vec =+        foldr (.+.) (F.lit 0) $+            [ F.lit (w / s) .*. ((Col n :: Expr Double) .-. F.lit mu)+            | (w, n, mu, s) <- zip4 (VU.toList vec) names means scales+            ]++-- | The PCA projection as a composable fitted 'Transform'.+pcaTransform :: PCAModel -> Transform+pcaTransform m = Transform [(n, UExpr e) | (n, e) <- pcaExprs m]++resolveK :: NComponents -> Int -> VU.Vector Double -> Int+resolveK (NComp k) d _ = max 1 (min k d)+resolveK (VarianceCovered frac) d ratio = max 1 (min d (go 0 0 1))+  where+    go !acc !cum !i+        | i > VU.length ratio = VU.length ratio+        | cum >= frac = acc+        | otherwise = go (acc + 1) (cum + ratio VU.! (i - 1)) (i + 1)++zip4 :: [a] -> [b] -> [c] -> [d] -> [(a, b, c, d)]+zip4 (a : as) (b : bs) (c : cs) (d : ds) = (a, b, c, d) : zip4 as bs cs ds+zip4 _ _ _ _ = []

@@ -0,0 +1,136 @@+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE OverloadedStrings #-}++{- | Kernel PCA with an RBF kernel, solved on a set of landmark points (Nyström).+Exact kernel PCA when the landmark count covers every row, a principled+approximation otherwise. 'fit' trains the model; the projection is exposed as+'kernelPCAExprs' / 'kernelPcaTransform' (a transformer, so no 'Predict').+-}+module DataFrame.PCA.Kernel (+    KernelPCAConfig (..),+    defaultKernelPCAConfig,+    KernelPCAModel (..),+    kernelPCAExprs,+    kernelPcaTransform,+) where++import qualified Data.Text as T+import qualified Data.Vector as V+import qualified Data.Vector.Unboxed as VU++import DataFrame.Featurize.Internal (Features (..), extractFeatures)+import qualified DataFrame.Functions as F+import DataFrame.Internal.DataFrame (DataFrame)+import DataFrame.Internal.Expression (Expr (..), UExpr (..))+import DataFrame.LinearAlgebra (sqDist)+import DataFrame.LinearAlgebra.Eigen (jacobiEigenSym)+import DataFrame.Model (Fit (..))+import DataFrame.Operators ((.*.), (.+.), (.-.))+import DataFrame.Random (mkGen, sampleIndices)+import DataFrame.Transform (Transform (..))++data KernelPCAConfig = KernelPCAConfig+    { kpcaNComponents :: !Int+    , kpcaGamma :: !(Maybe Double)+    , kpcaNLandmarks :: !Int+    , kpcaSeed :: !Int+    }+    deriving (Eq, Show)++defaultKernelPCAConfig :: KernelPCAConfig+defaultKernelPCAConfig =+    KernelPCAConfig+        { kpcaNComponents = 2+        , kpcaGamma = Nothing+        , kpcaNLandmarks = 128+        , kpcaSeed = 0+        }++{- | A fitted kernel PCA. Each component is @Σ_l βₗ·K(x, landmarkₗ) + cᵢ@ with an+RBF kernel of bandwidth 'kpcaGammaUsed'.+-}+data KernelPCAModel = KernelPCAModel+    { kpcaLandmarks :: !(V.Vector (VU.Vector Double))+    , kpcaBetas :: !(V.Vector (VU.Vector Double))+    , kpcaConsts :: !(VU.Vector Double)+    , kpcaEigenvalues :: !(VU.Vector Double)+    , kpcaGammaUsed :: !Double+    , kpcaFeatureNames :: !(V.Vector T.Text)+    }+    deriving (Eq, Show)++instance Fit KernelPCAConfig [Expr Double] KernelPCAModel where+    fit = fitKernelPCA++-- | Fit kernel PCA over the given feature columns.+fitKernelPCA :: KernelPCAConfig -> [Expr Double] -> DataFrame -> KernelPCAModel+fitKernelPCA cfg features df =+    KernelPCAModel+        { kpcaLandmarks = landmarks+        , kpcaBetas = betas+        , kpcaConsts = consts+        , kpcaEigenvalues = VU.take k evals+        , kpcaGammaUsed = gamma+        , kpcaFeatureNames = V.fromList names+        }+  where+    Features names _ rows n d = extractFeatures features df+    m = min (max 1 (kpcaNLandmarks cfg)) n+    (idx, _) = sampleIndices m n (mkGen (kpcaSeed cfg))+    landmarks = V.map (rows V.!) (V.convert idx)+    gamma = case kpcaGamma cfg of+        Just g -> g+        Nothing -> 1 / fromIntegral (max 1 d)+    kmat =+        V.generate m $ \i ->+            VU.generate m $ \j ->+                exp (negate gamma * sqDist (landmarks V.! i) (landmarks V.! j))+    rowMean i = VU.sum (kmat V.! i) / fromIntegral m+    totalMean = sum [rowMean i | i <- [0 .. m - 1]] / fromIntegral m+    centered =+        V.generate m $ \i ->+            VU.generate m $ \j ->+                (kmat V.! i) VU.! j - rowMean i - rowMean j + totalMean+    (evals, vecs) = jacobiEigenSym centered+    k = min (kpcaNComponents cfg) m+    alphas =+        V.generate k $ \i ->+            let lam = max 1e-12 (evals VU.! i)+             in VU.map (/ sqrt lam) (vecs V.! i)+    betas =+        V.map+            (\a -> let s = VU.sum a / fromIntegral m in VU.map (subtract s) a)+            alphas+    consts =+        VU.generate k $ \i ->+            let a = alphas V.! i+                sA = VU.sum a+             in negate (sum [a VU.! l * rowMean l | l <- [0 .. m - 1]])+                    + totalMean * sA++-- | Per-component projection expressions, named @kpc1@, @kpc2@, …+kernelPCAExprs :: KernelPCAModel -> [(T.Text, Expr Double)]+kernelPCAExprs m =+    [ ("kpc" <> T.pack (show (i + 1)), componentExpr i)+    | i <- [0 .. V.length (kpcaBetas m) - 1]+    ]+  where+    names = V.toList (kpcaFeatureNames m)+    gamma = kpcaGammaUsed m+    componentExpr i =+        foldr (.+.) (F.lit (kpcaConsts m VU.! i)) $+            [ F.lit (kpcaBetas m V.! i VU.! l) .*. kernelExpr (kpcaLandmarks m V.! l)+            | l <- [0 .. V.length (kpcaLandmarks m) - 1]+            ]+    kernelExpr landmark =+        exp (F.lit (negate gamma) .*. sqDistExpr landmark)+    sqDistExpr landmark =+        foldr (.+.) (F.lit 0) $+            [ let diff = (Col n :: Expr Double) .-. F.lit lj in diff .*. diff+            | (n, lj) <- zip names (VU.toList landmark)+            ]++-- | The kernel-PCA projection as a composable fitted 'Transform'.+kernelPcaTransform :: KernelPCAModel -> Transform+kernelPcaTransform m = Transform [(n, UExpr e) | (n, e) <- kernelPCAExprs m]

@@ -0,0 +1,121 @@+{-# LANGUAGE CPP #-}++{- | Deterministic, platform-independent random sampling for the stochastic+fitters. Built on @random@'s SplitMix 'StdGen' (only 'genWord64' and a+version-bridged split are used, the stable surface across @random@ versions);+the distributions here are our own so a seeded fit is bit-reproducible on+Linux, macOS, and Windows.+-}+module DataFrame.Random (+    Gen,+    mkGen,+    splitGen,+    nextWord64,+    nextDouble,+    nextIntR,+    gaussianPair,+    gaussianVector,+    shuffleInts,+    sampleIndices,+) where++import Control.Monad (forM_)+import Control.Monad.ST (runST)+import Data.Bits (shiftR)+import qualified Data.Vector.Unboxed as VU+import qualified Data.Vector.Unboxed.Mutable as VUM+import Data.Word (Word64)+import System.Random (StdGen, genWord64, mkStdGen)+import qualified System.Random as R++-- | The pure splittable generator. A fit is a function of @(seed, data)@.+type Gen = StdGen++-- | Seed a generator from an 'Int'.+mkGen :: Int -> Gen+mkGen = mkStdGen++-- | Split into two independent generators.+splitGen :: Gen -> (Gen, Gen)+#if MIN_VERSION_random(1,3,0)+splitGen = R.splitGen+#else+splitGen = R.split+#endif++-- | Raw 64-bit draw.+nextWord64 :: Gen -> (Word64, Gen)+nextWord64 = genWord64++-- | Uniform 'Double' in @[0, 1)@ from the top 53 bits (exact mantissa).+nextDouble :: Gen -> (Double, Gen)+nextDouble g =+    let (w, g') = genWord64 g+        d = fromIntegral (w `shiftR` 11) * (1 / 9007199254740992)+     in (d, g')++{- | Uniform 'Int' in the inclusive range @[lo, hi]@ by rejection sampling+(unbiased). Returns @lo@ when @hi <= lo@.+-}+nextIntR :: (Int, Int) -> Gen -> (Int, Gen)+nextIntR (lo, hi) g+    | hi <= lo = (lo, g)+    | otherwise = loop g+  where+    range = fromIntegral (hi - lo + 1) :: Word64+    threshold = negate range `mod` range+    loop gg =+        let (w, gg') = genWord64 gg+         in if w >= threshold+                then (lo + fromIntegral (w `mod` range), gg')+                else loop gg'++{- | A pair of independent standard normals via Box-Muller, consuming exactly two+uniforms so stream offsets stay data-independent.+-}+gaussianPair :: Gen -> ((Double, Double), Gen)+gaussianPair g =+    let (u1, g1) = nextDouble g+        (u2, g2) = nextDouble g1+        u1' = if u1 <= 0 then 2.220446049250313e-16 else u1+        r = sqrt (-(2 * log u1'))+        a = 2 * pi * u2+     in ((r * cos a, r * sin a), g2)++-- | A length-@n@ vector of standard normals.+gaussianVector :: Int -> Gen -> (VU.Vector Double, Gen)+gaussianVector n g0 = go n g0 []+  where+    go k g acc+        | k <= 0 = (VU.fromList (take n (reverse acc)), g)+        | otherwise =+            let ((z0, z1), g') = gaussianPair g+             in go (k - 2) g' (z1 : z0 : acc)++{- | A uniformly random permutation of @[0 .. n-1]@ (Fisher-Yates), threading the+generator purely.+-}+shuffleInts :: Int -> Gen -> (VU.Vector Int, Gen)+shuffleInts n g0+    | n <= 1 = (VU.enumFromN 0 (max 0 n), g0)+    | otherwise =+        let (swaps, g1) = genSwaps (n - 1) g0 []+            v = runST $ do+                m <- VU.thaw (VU.enumFromN 0 n)+                forM_ swaps $ uncurry (VUM.swap m)+                VU.freeze m+         in (v, g1)+  where+    genSwaps i g acc+        | i < 1 = (reverse acc, g)+        | otherwise =+            let (j, g') = nextIntR (0, i) g+             in genSwaps (i - 1) g' ((i, j) : acc)++{- | @sampleIndices k n@ draws @k@ distinct indices from @[0 .. n-1]@ (the first+@k@ of a full shuffle); returns all @n@ when @k >= n@.+-}+sampleIndices :: Int -> Int -> Gen -> (VU.Vector Int, Gen)+sampleIndices k n g =+    let (perm, g') = shuffleInts n g+     in (VU.take (min k n) perm, g')

@@ -0,0 +1,101 @@+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE UndecidableInstances #-}++{- | Linear support vector classification: L2-regularized squared hinge fitted+with the FISTA engine (sklearn's LinearSVC default loss). 'fit' trains a+one-vs-rest 'LinearSVCModel'; 'predict' is the arg-max class margin. There is no+@predict_proba@, matching sklearn's LinearSVC.+-}+module DataFrame.SVM (+    LinearSVCModel (..),+    SVCConfig (..),+    defaultSVCConfig,+    svcMarginExprs,+) where++import Data.List (sort)+import qualified Data.Map.Strict as M+import qualified Data.Vector as V+import qualified Data.Vector.Unboxed as VU++import DataFrame.Featurize.Internal (+    affineExpr,+    argMaxExpr,+    featureNames,+    numericMatrix,+    targetValues,+ )+import DataFrame.Internal.Column (Columnable)+import DataFrame.Internal.DataFrame (DataFrame)+import DataFrame.Internal.Expression (Expr)+import DataFrame.LinearSolver (LinearModel (..), SolverConfig (..), fitProx)+import DataFrame.LinearSolver.Loss (sqHingeLoss)+import DataFrame.Model (Fit (..), Predict (..))++-- | Hyper-parameters. @svcC@ is the inverse regularization strength (sklearn @C@).+data SVCConfig = SVCConfig+    { svcC :: !Double+    , svcMaxIter :: !Int+    , svcTol :: !Double+    }+    deriving (Eq, Show)++defaultSVCConfig :: SVCConfig+defaultSVCConfig = SVCConfig{svcC = 1.0, svcMaxIter = 1000, svcTol = 1.0e-4}++-- | A fitted one-vs-rest linear SVC: class labels and their margin sub-models.+data LinearSVCModel a = LinearSVCModel+    { svcClasses :: !(V.Vector a)+    , svcModels :: !(V.Vector LinearModel)+    }+    deriving (Eq, Show)++instance (Columnable a, Ord a) => Fit SVCConfig (Expr a) (LinearSVCModel a) where+    fit = fitLinearSVC++instance (Columnable a, Ord a) => Predict (LinearSVCModel a) a where+    predict m = argMaxExpr (labelledMargins m)++-- | Fit a one-vs-rest linear SVC.+fitLinearSVC ::+    (Columnable a, Ord a) =>+    SVCConfig -> Expr a -> DataFrame -> LinearSVCModel a+fitLinearSVC cfg target df =+    LinearSVCModel (V.fromList classes) (V.fromList (map fitOne classes))+  where+    names = featureNames target df+    (nameVec, mat) = numericMatrix names df+    ys = targetValues target df+    classes = sort (foldr dedup [] (V.toList ys))+    dedup x acc = if x `elem` acc then acc else x : acc+    solverCfg =+        SolverConfig+            { scL1Lambda = 0+            , scL2Lambda = 1 / svcC cfg+            , scMaxIter = svcMaxIter cfg+            , scTol = svcTol cfg+            , scSampleWeights = Nothing+            }+    fitOne c =+        let labels =+                VU.generate (V.length ys) (\i -> if ys V.! i == c then 1 else -1)+         in fitProx sqHingeLoss solverCfg mat labels nameVec++-- | The raw margin expression for each class.+svcMarginExprs ::+    (Columnable a, Ord a) => LinearSVCModel a -> M.Map a (Expr Double)+svcMarginExprs m = M.fromList (labelledMargins m)++labelledMargins :: LinearSVCModel a -> [(a, Expr Double)]+labelledMargins m =+    [ (svcClasses m V.! i, marginOf (svcModels m V.! i))+    | i <- [0 .. V.length (svcClasses m) - 1]+    ]++marginOf :: LinearModel -> Expr Double+marginOf m =+    affineExpr+        (lmIntercept m)+        (zip (VU.toList (lmWeights m)) (V.toList (lmFeatureNames m)))

@@ -0,0 +1,146 @@+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE OverloadedStrings #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE UndecidableInstances #-}++{- | Approximate RBF-kernel SVM via Random Fourier Features (Rahimi & Recht): map+each row through @z(x) = √(2/D)·cos(W x + b)@ with @W ~ N(0, 2γI)@ (seeded), then+fit a linear SVC in the random-feature space. 'predict' compiles to a closed+@Σ_r β_r·cos(…)@ expression of size @O(D·d)@, independent of the row count.+-}+module DataFrame.SVM.RFF (+    RFFConfig (..),+    defaultRFFConfig,+    RFFSVMModel (..),+) where++import Data.List (sort)+import qualified Data.Text as T+import qualified Data.Vector as V+import qualified Data.Vector.Unboxed as VU++import DataFrame.Featurize.Internal (featureNames, numericMatrix, targetValues)+import qualified DataFrame.Functions as F+import DataFrame.Internal.Column (Columnable)+import DataFrame.Internal.DataFrame (DataFrame)+import DataFrame.Internal.Expression (Expr (..))+import DataFrame.LinearAlgebra (dot)+import DataFrame.LinearSolver (LinearModel (..), SolverConfig (..), fitProx)+import DataFrame.LinearSolver.Loss (sqHingeLoss)+import DataFrame.Model (Fit (..), Predict (..))+import DataFrame.Operators ((.*.), (.+.), (.>.))+import DataFrame.Random (Gen, gaussianVector, mkGen, nextDouble)++data RFFConfig = RFFConfig+    { rffD :: !Int+    , rffGamma :: !Double+    , rffC :: !Double+    , rffMaxIter :: !Int+    , rffTol :: !Double+    , rffSeed :: !Int+    }+    deriving (Eq, Show)++defaultRFFConfig :: RFFConfig+defaultRFFConfig =+    RFFConfig+        { rffD = 100+        , rffGamma = 0.1+        , rffC = 1.0+        , rffMaxIter = 1000+        , rffTol = 1.0e-4+        , rffSeed = 0+        }++{- | A fitted RFF SVM (binary). 'rffW' / 'rffB' are the random projection;+'rffCoef' / 'rffIntercept' the linear SVC in feature space.+-}+data RFFSVMModel a = RFFSVMModel+    { rffW :: !(V.Vector (VU.Vector Double))+    , rffB :: !(VU.Vector Double)+    , rffCoef :: !(VU.Vector Double)+    , rffIntercept :: !Double+    , rffScale :: !Double+    , rffNegClass :: !a+    , rffPosClass :: !a+    , rffFeatureNames :: !(V.Vector T.Text)+    }+    deriving (Show)++instance (Columnable a, Ord a) => Fit RFFConfig (Expr a) (RFFSVMModel a) where+    fit = fitRFFSVM++instance (Columnable a) => Predict (RFFSVMModel a) a where+    predict m =+        If (margin .>. F.lit 0) (Lit (rffPosClass m)) (Lit (rffNegClass m))+      where+        names = V.toList (rffFeatureNames m)+        margin =+            foldr (.+.) (F.lit (rffIntercept m)) $+                [ F.lit (rffCoef m VU.! r * rffScale m) .*. cosTerm r+                | r <- [0 .. V.length (rffW m) - 1]+                , rffCoef m VU.! r /= 0+                ]+        cosTerm r = cos (linComb (rffW m V.! r) (rffB m VU.! r))+        linComb w b =+            foldr (.+.) (F.lit b) $+                [ F.lit (w VU.! j) .*. (Col n :: Expr Double)+                | (j, n) <- zip [0 ..] names+                ]++-- | Fit a binary RFF SVM. Targets with more than two classes are rejected.+fitRFFSVM ::+    (Columnable a, Ord a) => RFFConfig -> Expr a -> DataFrame -> RFFSVMModel a+fitRFFSVM cfg target df =+    case classes of+        [neg, pos] -> build neg pos+        _ -> error "fitRFFSVM: binary classification only (got /= 2 classes)"+  where+    names = featureNames target df+    (nameVec, mat) = numericMatrix names df+    ys = targetValues target df+    classes = sort (foldr dedup [] (V.toList ys))+    dedup x acc = if x `elem` acc then acc else x : acc+    d = if V.null mat then 0 else VU.length (V.head mat)+    bigD = max 1 (rffD cfg)+    (ws, bs) = sampleRFF bigD d (rffGamma cfg) (mkGen (rffSeed cfg))+    scale = sqrt (2 / fromIntegral bigD)+    z = V.map (featureRow ws bs scale) mat+    featNames = V.fromList ["rff" <> T.pack (show r) | r <- [0 .. bigD - 1]]+    build neg pos =+        let labels = VU.generate (V.length ys) (\i -> if ys V.! i == pos then 1 else -1)+            solverCfg =+                SolverConfig+                    { scL1Lambda = 0+                    , scL2Lambda = 1 / rffC cfg+                    , scMaxIter = rffMaxIter cfg+                    , scTol = rffTol cfg+                    , scSampleWeights = Nothing+                    }+            model = fitProx sqHingeLoss solverCfg z labels featNames+         in RFFSVMModel ws bs (lmWeights model) (lmIntercept model) scale neg pos nameVec++sampleRFF ::+    Int -> Int -> Double -> Gen -> (V.Vector (VU.Vector Double), VU.Vector Double)+sampleRFF bigD d gamma g0 = (V.fromList ws, VU.fromList bs)+  where+    sigma = sqrt (2 * gamma)+    (ws, g1) = goW bigD g0 []+    goW 0 g acc = (reverse acc, g)+    goW k g acc =+        let (vec, g') = gaussianVector d g+         in goW (k - 1) g' (VU.map (* sigma) vec : acc)+    bs = take bigD (goB g1)+    goB g = let (u, g') = nextDouble g in (u * 2 * pi) : goB g'++featureRow ::+    V.Vector (VU.Vector Double) ->+    VU.Vector Double ->+    Double ->+    VU.Vector Double ->+    VU.Vector Double+featureRow ws bs scale x =+    VU.generate (V.length ws) $ \r ->+        scale * cos (dot (ws V.! r) x + bs VU.! r)

@@ -0,0 +1,120 @@+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE ScopedTypeVariables #-}++{- | Symbolic regression by genetic programming (modelled on the+@symbolic-regression@ library, ported dependency-light: no e-graphs, no NLOPT).+'predict' is the best discovered @Expr Double@; the search also returns the+accuracy-vs-complexity Pareto front. Deterministic given the seed.+-}+module DataFrame.SymbolicRegression (+    UnOp (..),+    SRConfig (..),+    defaultSRConfig,+    SRModel (..),+) where++import Control.Exception (throw)+import qualified Data.Vector as V+import qualified Data.Vector.Unboxed as VU++import DataFrame.Featurize.Internal (featureNames, targetDoubles)+import DataFrame.Internal.DataFrame (DataFrame)+import DataFrame.Internal.Expression (Expr (..))+import DataFrame.Model (Fit (..), Predict (..))+import DataFrame.Operations.Core (columnAsDoubleVector)+import DataFrame.Random (mkGen)+import DataFrame.SymbolicRegression.Expr (+    UnOp (..),+    allUnOps,+    toDataFrameExpr,+ )+import DataFrame.SymbolicRegression.GP (GPParams (..), runGP)+import DataFrame.SymbolicRegression.Simplify (simplify)++data SRConfig = SRConfig+    { srSeed :: !Int+    , srPopSize :: !Int+    , srGenerations :: !Int+    , srMaxSize :: !Int+    , srTournament :: !Int+    , srCrossoverP :: !Double+    , srMutationP :: !Double+    , srOptimizeP :: !Double+    , srParsimony :: !Double+    , srUnaryOps :: ![UnOp]+    }+    deriving (Eq, Show)++defaultSRConfig :: SRConfig+defaultSRConfig =+    SRConfig+        { srSeed = 42+        , srPopSize = 200+        , srGenerations = 40+        , srMaxSize = 25+        , srTournament = 5+        , srCrossoverP = 0.9+        , srMutationP = 0.3+        , srOptimizeP = 0.15+        , srParsimony = 1.0e-3+        , srUnaryOps = allUnOps+        }++{- | A fitted symbolic regressor. 'srBest' is the lowest-error expression;+'srPareto' is the @(complexity, mse, expr)@ frontier.+-}+data SRModel = SRModel+    { srBest :: !(Expr Double)+    , srBestMSE :: !Double+    , srPareto :: ![(Int, Double, Expr Double)]+    , srGenerationsRun :: !Int+    }++instance Fit SRConfig (Expr Double) SRModel where+    fit = fitSymbolicRegression++instance Predict SRModel Double where+    predict = srBest++-- | Search for an expression predicting @target@ from the other columns.+fitSymbolicRegression :: SRConfig -> Expr Double -> DataFrame -> SRModel+fitSymbolicRegression cfg target df =+    SRModel+        { srBest = translate best+        , srBestMSE = bestMse+        , srPareto = [(sz, mse, translate e) | (sz, mse, e) <- front]+        , srGenerationsRun = gens+        }+  where+    names = featureNames target df+    nameVec = V.fromList names+    cols = V.fromList (map (materialize df . Col) names)+    target' = targetDoubles target df+    n = VU.length target'+    params =+        GPParams+            { gpFeats = cols+            , gpN = n+            , gpTarget = target'+            , gpNVars = length names+            , gpUnOps = srUnaryOps cfg+            , gpPopSize = srPopSize cfg+            , gpGenerations = srGenerations cfg+            , gpMaxSize = srMaxSize cfg+            , gpTournament = srTournament cfg+            , gpCrossoverP = srCrossoverP cfg+            , gpMutationP = srMutationP cfg+            , gpOptimizeP = srOptimizeP cfg+            , gpParsimony = srParsimony cfg+            }+    (best, front, gens) = runGP params (mkGen (srSeed cfg))+    bestMse = case [m | (_, m, e) <- front, e == best] of+        (m : _) -> m+        [] -> 1 / 0+    translate = toDataFrameExpr nameVec . simplify++materialize :: DataFrame -> Expr Double -> VU.Vector Double+materialize df e = case columnAsDoubleVector e df of+    Right v -> v+    Left err -> throw err

@@ -0,0 +1,122 @@+{-# LANGUAGE FlexibleContexts #-}++{- | The symbolic-regression expression tree: a small first-order ADT (no hegg,+no 'Fix'). Vectorized evaluation over a feature matrix and a total translation+to a dataframe 'Expr Double' (the SR result IS a dataframe expression). Division,+log, and sqrt are protected so evaluation never produces @NaN@.+-}+module DataFrame.SymbolicRegression.Expr (+    SRExpr (..),+    BinOp (..),+    UnOp (..),+    evalSR,+    toDataFrameExpr,+    srSize,+    constants,+    setConstants,+    allBinOps,+    allUnOps,+) where++import qualified Data.Text as T+import qualified Data.Vector as V+import qualified Data.Vector.Unboxed as VU++import qualified DataFrame.Functions as F+import DataFrame.Internal.Expression (Expr (..))+import DataFrame.Operators ((.*.), (.+.), (.-.), (./.))++data BinOp = SAdd | SSub | SMul | SDiv+    deriving (Eq, Ord, Show, Enum, Bounded)++data UnOp = SNeg | SSin | SCos | SExp | SLog | SSqrt+    deriving (Eq, Ord, Show, Enum, Bounded)++-- | A symbolic-regression expression over feature variables and constants.+data SRExpr+    = SVar !Int+    | SConst !Double+    | SUn !UnOp SRExpr+    | SBin !BinOp SRExpr SRExpr+    deriving (Eq, Ord, Show)++allBinOps :: [BinOp]+allBinOps = [minBound .. maxBound]++allUnOps :: [UnOp]+allUnOps = [minBound .. maxBound]++{- | Evaluate over a feature matrix given column-major (@feats ! j@ is feature+@j@ across all rows). Protected operators keep results finite.+-}+evalSR :: V.Vector (VU.Vector Double) -> Int -> SRExpr -> VU.Vector Double+evalSR feats n = go+  where+    go (SVar j)+        | j < V.length feats = feats V.! j+        | otherwise = VU.replicate n 0+    go (SConst c) = VU.replicate n c+    go (SUn op e) = VU.map (unFn op) (go e)+    go (SBin op a b) = VU.zipWith (binFn op) (go a) (go b)++binFn :: BinOp -> Double -> Double -> Double+binFn SAdd a b = a + b+binFn SSub a b = a - b+binFn SMul a b = a * b+binFn SDiv a b = if abs b < 1e-9 then 1 else a / b++unFn :: UnOp -> Double -> Double+unFn SNeg = negate+unFn SSin = sin+unFn SCos = cos+unFn SExp = exp . min 50+unFn SLog = \x -> log (abs x + 1e-9)+unFn SSqrt = sqrt . abs++-- | Translate to a dataframe expression over the named feature columns.+toDataFrameExpr :: V.Vector T.Text -> SRExpr -> Expr Double+toDataFrameExpr names = go+  where+    go (SVar j)+        | j < V.length names = Col (names V.! j)+        | otherwise = F.lit 0+    go (SConst c) = F.lit c+    go (SUn op e) = unExpr op (go e)+    go (SBin op a b) = binExpr op (go a) (go b)+    unExpr SNeg = negate+    unExpr SSin = sin+    unExpr SCos = cos+    unExpr SExp = exp+    unExpr SLog = log+    unExpr SSqrt = sqrt+    binExpr SAdd = (.+.)+    binExpr SSub = (.-.)+    binExpr SMul = (.*.)+    binExpr SDiv = (./.)++srSize :: SRExpr -> Int+srSize (SVar _) = 1+srSize (SConst _) = 1+srSize (SUn _ e) = 1 + srSize e+srSize (SBin _ a b) = 1 + srSize a + srSize b++-- | The constant values in left-to-right traversal order.+constants :: SRExpr -> [Double]+constants (SConst c) = [c]+constants (SVar _) = []+constants (SUn _ e) = constants e+constants (SBin _ a b) = constants a ++ constants b++-- | Replace the constants in traversal order; extra values are ignored.+setConstants :: [Double] -> SRExpr -> SRExpr+setConstants vals e = fst (go vals e)+  where+    go vs (SConst _) = case vs of+        (v : rest) -> (SConst v, rest)+        [] -> (SConst 0, [])+    go vs (SVar j) = (SVar j, vs)+    go vs (SUn op a) = let (a', vs') = go vs a in (SUn op a', vs')+    go vs (SBin op a b) =+        let (a', vs') = go vs a+            (b', vs'') = go vs' b+         in (SBin op a' b', vs'')

@@ -0,0 +1,207 @@+{- | A compact generational genetic-programming search over 'SRExpr': ramped+random initialization, tournament selection, subtree crossover and mutation,+elitism, and a complexity-keyed Pareto archive. Deterministic given the seed+(the splitmix generator from "DataFrame.Random").+-}+module DataFrame.SymbolicRegression.GP (+    GPParams (..),+    runGP,+) where++import Data.List (foldl', minimumBy, sortBy)+import qualified Data.Map.Strict as M+import Data.Ord (comparing)+import qualified Data.Vector as V+import qualified Data.Vector.Unboxed as VU++import DataFrame.Random (Gen, nextDouble, nextIntR)+import DataFrame.SymbolicRegression.Expr+import DataFrame.SymbolicRegression.Optimize (+    meanSquaredError,+    optimizeConstants,+ )+import DataFrame.SymbolicRegression.Simplify (simplify)++-- | GP hyper-parameters resolved from the public config.+data GPParams = GPParams+    { gpFeats :: !(V.Vector (VU.Vector Double))+    , gpN :: !Int+    , gpTarget :: !(VU.Vector Double)+    , gpNVars :: !Int+    , gpUnOps :: ![UnOp]+    , gpPopSize :: !Int+    , gpGenerations :: !Int+    , gpMaxSize :: !Int+    , gpTournament :: !Int+    , gpCrossoverP :: !Double+    , gpMutationP :: !Double+    , gpOptimizeP :: !Double+    , gpParsimony :: !Double+    }++type Scored = (SRExpr, Double)++-- | Run the search; returns @(best, pareto front, generations run)@.+runGP :: GPParams -> Gen -> (SRExpr, [(Int, Double, SRExpr)], Int)+runGP p g0 =+    let (pop0, g1) = initPop p g0+        scored0 = map (scoreOf p) pop0+        arch0 = foldl' (archiveInsert p) M.empty scored0+        (_, finalArch, gN, _) =+            iterate' 0 scored0 arch0 g1+        best = bestOfArchive finalArch+        front =+            [ (sz, mse, e)+            | (sz, (mse, e)) <- M.toList finalArch+            ]+     in (snd3 best, sortBy (comparing fst3) front, gN)+  where+    iterate' gen pop arch g+        | gen >= gpGenerations p = (pop, arch, gen, g)+        | otherwise =+            let (pop', g') = nextGen p pop g+                arch' = foldl' (archiveInsert p) arch pop'+             in iterate' (gen + 1) pop' arch' g'+    fst3 (a, _, _) = a+    snd3 (_, b, _) = b+    bestOfArchive arch =+        case M.toList arch of+            [] -> (0 :: Int, SConst 0, 1 / 0)+            xs ->+                let (sz, (mse, e)) = minimumBy (comparing (fst . snd)) xs+                 in (sz, e, mse)++scoreOf :: GPParams -> SRExpr -> Scored+scoreOf p e = (e, meanSquaredError (gpFeats p) (gpN p) (gpTarget p) e)++fitness :: GPParams -> Scored -> Double+fitness p (e, mse) = mse + gpParsimony p * fromIntegral (srSize e)++archiveInsert ::+    GPParams -> M.Map Int (Double, SRExpr) -> Scored -> M.Map Int (Double, SRExpr)+archiveInsert _ arch (e, mse)+    | isNaN mse || isInfinite mse = arch+    | otherwise =+        let key = srSize (simplify e)+         in M.insertWith better key (mse, e) arch+  where+    better newv@(m1, _) oldv@(m2, _) = if m1 < m2 then newv else oldv++initPop :: GPParams -> Gen -> ([SRExpr], Gen)+initPop p = go (gpPopSize p) []+  where+    go 0 acc g = (acc, g)+    go k acc g =+        let (depth, g1) = nextIntR (1, 4) g+            (e, g2) = randomExpr p depth g1+         in go (k - 1) (e : acc) g2++randomExpr :: GPParams -> Int -> Gen -> (SRExpr, Gen)+randomExpr p depth g+    | depth <= 1 = randomLeaf p g+    | otherwise =+        let (r, g1) = nextDouble g+         in if r < 0.3+                then randomLeaf p g1+                else+                    let (isUn, g2) = nextDouble g1+                     in if isUn < 0.3 && not (null (gpUnOps p))+                            then+                                let (oi, g3) = nextIntR (0, length (gpUnOps p) - 1) g2+                                    (e, g4) = randomExpr p (depth - 1) g3+                                 in (SUn (gpUnOps p !! oi) e, g4)+                            else+                                let (oi, g3) = nextIntR (0, length allBinOps - 1) g2+                                    (a, g4) = randomExpr p (depth - 1) g3+                                    (b, g5) = randomExpr p (depth - 1) g4+                                 in (SBin (allBinOps !! oi) a b, g5)++randomLeaf :: GPParams -> Gen -> (SRExpr, Gen)+randomLeaf p g =+    let (r, g1) = nextDouble g+     in if r < 0.6 && gpNVars p > 0+            then let (j, g2) = nextIntR (0, gpNVars p - 1) g1 in (SVar j, g2)+            else let (c, g2) = nextDouble g1 in (SConst (c * 4 - 2), g2)++nextGen :: GPParams -> [Scored] -> Gen -> ([Scored], Gen)+nextGen p pop g0 =+    let elite = minimumBy (comparing (fitness p)) pop+        (rest, g1) = go (gpPopSize p - 1) [] g0+     in (elite : rest, g1)+  where+    go 0 acc g = (acc, g)+    go k acc g =+        let (child, g') = breed p pop g+            scored = optimizeMaybe p child g'+         in go (k - 1) (fst scored : acc) (snd scored)++optimizeMaybe :: GPParams -> SRExpr -> Gen -> (Scored, Gen)+optimizeMaybe p e g =+    let (r, g1) = nextDouble g+        e' =+            if r < gpOptimizeP p+                then optimizeConstants (gpFeats p) (gpN p) (gpTarget p) 15 e+                else e+     in (scoreOf p e', g1)++breed :: GPParams -> [Scored] -> Gen -> (SRExpr, Gen)+breed p pop g0 =+    let (pa, g1) = tournament p pop g0+        (doX, g2) = nextDouble g1+        (child, g3) =+            if doX < gpCrossoverP p+                then+                    let (pb, g2') = tournament p pop g2+                        (c, g3') = crossover pa pb g2'+                     in (c, g3')+                else (pa, g2)+        (doM, g4) = nextDouble g3+        (child', g5) =+            if doM < gpMutationP p then mutate p child g4 else (child, g4)+        capped = if srSize child' > gpMaxSize p then pa else child'+     in (simplify capped, g5)++tournament :: GPParams -> [Scored] -> Gen -> (SRExpr, Gen)+tournament p pop g0 =+    let (picks, g1) = pickN (gpTournament p) g0+        chosen = map (pop !!) picks+     in (fst (minimumBy (comparing (fitness p)) chosen), g1)+  where+    n = length pop+    pickN 0 g = ([], g)+    pickN k g =+        let (i, g') = nextIntR (0, n - 1) g+            (is, g'') = pickN (k - 1) g'+         in (i : is, g'')++crossover :: SRExpr -> SRExpr -> Gen -> (SRExpr, Gen)+crossover a b g0 =+    let (ia, g1) = nextIntR (0, srSize a - 1) g0+        (ib, g2) = nextIntR (0, srSize b - 1) g1+        sub = subtreeAt ib b+     in (replaceAt ia a sub, g2)++mutate :: GPParams -> SRExpr -> Gen -> (SRExpr, Gen)+mutate p e g0 =+    let (i, g1) = nextIntR (0, srSize e - 1) g0+        (depth, g2) = nextIntR (1, 3) g1+        (newSub, g3) = randomExpr p depth g2+     in (replaceAt i e newSub, g3)++subtreeAt :: Int -> SRExpr -> SRExpr+subtreeAt 0 e = e+subtreeAt i (SUn _ e) = subtreeAt (i - 1) e+subtreeAt i (SBin _ a b) =+    let sa = srSize a+     in if i <= sa then subtreeAt (i - 1) a else subtreeAt (i - 1 - sa) b+subtreeAt _ e = e++replaceAt :: Int -> SRExpr -> SRExpr -> SRExpr+replaceAt 0 _ new = new+replaceAt i (SUn op e) new = SUn op (replaceAt (i - 1) e new)+replaceAt i (SBin op a b) new =+    let sa = srSize a+     in if i <= sa+            then SBin op (replaceAt (i - 1) a new) b+            else SBin op a (replaceAt (i - 1 - sa) b new)+replaceAt _ e _ = e

@@ -0,0 +1,68 @@+{-# LANGUAGE BangPatterns #-}++{- | Constant optimization for symbolic-regression candidates: finite-difference+gradient descent with backtracking line search on the embedded constants. Pure+and dependency-free; effective at the one-to-few constants a tree carries, which+is where random restarts plus a quasi-Newton step would otherwise be used.+-}+module DataFrame.SymbolicRegression.Optimize (+    optimizeConstants,+    meanSquaredError,+) where++import qualified Data.Vector as V+import qualified Data.Vector.Unboxed as VU++import DataFrame.SymbolicRegression.Expr (+    SRExpr,+    constants,+    evalSR,+    setConstants,+ )++-- | Mean squared error of an expression's predictions against the target.+meanSquaredError ::+    V.Vector (VU.Vector Double) -> Int -> VU.Vector Double -> SRExpr -> Double+meanSquaredError feats n target e =+    let pred = evalSR feats n e+        diff = VU.zipWith (-) pred target+     in VU.sum (VU.map (\x -> x * x) diff) / fromIntegral (max 1 n)++-- | Refine an expression's constants to reduce MSE (no-op when constant-free).+optimizeConstants ::+    V.Vector (VU.Vector Double) ->+    Int ->+    VU.Vector Double ->+    Int ->+    SRExpr ->+    SRExpr+optimizeConstants feats n target iters expr+    | null theta0 = expr+    | otherwise = setConstants (descend iters theta0) expr+  where+    theta0 = constants expr+    eps = 1e-6+    mseAt theta = meanSquaredError feats n target (setConstants theta expr)+    descend 0 theta = theta+    descend k theta =+        let f0 = mseAt theta+            g = numGrad theta+            gn = sqrt (sum (map (\x -> x * x) g))+         in if gn < 1e-10+                then theta+                else+                    let theta' = lineSearch theta g f0+                     in if theta' == theta then theta else descend (k - 1) theta'+    numGrad theta =+        [ (mseAt (bump i eps theta) - mseAt (bump i (negate eps) theta)) / (2 * eps)+        | i <- [0 .. length theta - 1]+        ]+    bump i delta theta =+        [if j == i then t + delta else t | (j, t) <- zip [0 ..] theta]+    lineSearch theta g f0 = go (1.0 :: Double)+      where+        go !step+            | step < 1e-8 = theta+            | otherwise =+                let theta' = zipWith (\t gi -> t - step * gi) theta g+                 in if mseAt theta' < f0 then theta' else go (step / 2)

@@ -0,0 +1,65 @@+{- | A fuel-bounded, deterministic algebraic simplifier — the dependency-light+stand-in for equality saturation. It is used as a canonical dedup key for the+Pareto archive and to tidy reported expressions; no confluence is claimed, only+that it is a total, size-non-increasing, idempotent function whose rewrites+preserve evaluation.+-}+module DataFrame.SymbolicRegression.Simplify (+    simplify,+) where++import DataFrame.SymbolicRegression.Expr (BinOp (..), SRExpr (..), UnOp (..))++-- | Simplify to a fixed point (bounded by a fuel counter).+simplify :: SRExpr -> SRExpr+simplify = go (10 :: Int)+  where+    go 0 e = e+    go fuel e =+        let e' = step e+         in if e' == e then e else go (fuel - 1) e'++step :: SRExpr -> SRExpr+step (SUn op e) = simplifyUn op (step e)+step (SBin op a b) = simplifyBin op (step a) (step b)+step e = e++simplifyUn :: UnOp -> SRExpr -> SRExpr+simplifyUn SNeg (SUn SNeg e) = e+simplifyUn op (SConst c) = SConst (foldUn op c)+simplifyUn op e = SUn op e++simplifyBin :: BinOp -> SRExpr -> SRExpr -> SRExpr+simplifyBin op (SConst a) (SConst b) = SConst (foldBin op a b)+simplifyBin SAdd a (SConst 0) = a+simplifyBin SAdd (SConst 0) b = b+simplifyBin SSub a (SConst 0) = a+simplifyBin SSub a b | a == b = SConst 0+simplifyBin SMul _ (SConst 0) = SConst 0+simplifyBin SMul (SConst 0) _ = SConst 0+simplifyBin SMul a (SConst 1) = a+simplifyBin SMul (SConst 1) b = b+simplifyBin SDiv a (SConst 1) = a+simplifyBin SDiv a b | a == b = SConst 1+simplifyBin op a b+    | commutative op && a > b = SBin op b a+    | otherwise = SBin op a b++commutative :: BinOp -> Bool+commutative SAdd = True+commutative SMul = True+commutative _ = False++foldBin :: BinOp -> Double -> Double -> Double+foldBin SAdd a b = a + b+foldBin SSub a b = a - b+foldBin SMul a b = a * b+foldBin SDiv a b = if abs b < 1e-9 then 1 else a / b++foldUn :: UnOp -> Double -> Double+foldUn SNeg = negate+foldUn SSin = sin+foldUn SCos = cos+foldUn SExp = exp . min 50+foldUn SLog = \x -> log (abs x + 1e-9)+foldUn SSqrt = sqrt . abs

@@ -1,483 +1,371 @@ {-# LANGUAGE BangPatterns #-}-{-# LANGUAGE ExplicitNamespaces #-}-{-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE FlexibleInstances #-}-{-# LANGUAGE GADTs #-} {-# LANGUAGE MultiParamTypeClasses #-}-{-# LANGUAGE OverloadedStrings #-}-{-# LANGUAGE RankNTypes #-} {-# LANGUAGE ScopedTypeVariables #-}-{-# LANGUAGE TypeApplications #-}-{-# LANGUAGE UndecidableInstances #-} -module DataFrame.Synthesis where+{- | Feature synthesis by bottom-up enumerative search with observational+equivalence — the canonical enumerative method from Solar-Lezama's+/Introduction to Program Synthesis/, hardened for a numeric, examples-only+setting. +Given a frame and a numeric target column, it searches for a small, interpretable+arithmetic expression over the other columns whose values track the target. The+specification is purely the example rows; there is no SMT solver and no logical+spec. Deterministic and pure.++The engine:++  * enumerates programs by increasing AST size (so the first representative of any+    behaviour is the smallest — interpretability for free);+  * evaluates each candidate /incrementally/ by combining the cached result+    vectors of its subprograms (one vector op), never re-interpreting the whole+    tree;+  * keeps exactly one program per /observational-equivalence/ class — candidates+    producing the same column (up to a float tolerance) are interchangeable, so+    duplicates are dropped rather than re-explored;+  * breaks commutative symmetry (never both @a+b@ and @b+a@) and uses protected+    operators (@sqrt|x|@, @log(|x|+1)@) plus a denominator guard so domain errors+    never arise;+  * caps each size layer by fit score when it grows large (a cost-guided+    tractability bound over /distinct/ behaviours, not a lossy beam over raw+    syntax).++'fit' returns the best 'SynthesizedFeature'; 'predict' is its expression.+'synthesizeFeatures' returns the whole ranked, deduplicated feature bank — useful+as automated feature engineering feeding a downstream model.++Deferred (documented next steps, not yet implemented): skeleton enumeration with+closed-form least-squares coefficient fitting, hard-row counterexample sampling+for very large frames, and piecewise (condition-abduction) features.+-}+module DataFrame.Synthesis (+    LossFunction (..),+    SynthesisConfig (..),+    defaultSynthesisConfig,+    SynthesizedFeature (..),+    synthesizeFeatures,+) where++import Data.Bits (xor)+import Data.Either (fromRight)+import Data.List (sortBy)+import qualified Data.Map.Strict as M+import Data.Maybe (fromMaybe)+import Data.Ord (Down (..), comparing)+import qualified Data.Text as T+import qualified Data.Vector.Unboxed as VU+import Data.Word (Word64)+import GHC.Float (castDoubleToWord64)++import DataFrame.Featurize.Internal (featureNames) import qualified DataFrame.Functions as F-import DataFrame.Internal.Column-import DataFrame.Internal.DataFrame (-    DataFrame (..),-    columnNames,- )-import DataFrame.Internal.Expression (-    Expr (..),-    eSize,-    eqExpr,+import DataFrame.Internal.DataFrame (DataFrame)+import DataFrame.Internal.Expression (Expr (..))+import DataFrame.Internal.Statistics (+    meanSquaredError,+    mutualInformationBinned,+    percentile',+    variance',  )-import DataFrame.Internal.Interpreter (interpret)-import DataFrame.Internal.Statistics+import DataFrame.Model (Fit (..), Predict (..)) import DataFrame.Operations.Core (columnAsDoubleVector)-import qualified DataFrame.Operations.Statistics as Stats-import DataFrame.Operations.Subset (exclude) -import Control.Exception (throw)-import Data.Function-import qualified Data.List as L-import qualified Data.Map as M-import Data.Maybe (listToMaybe)-import qualified Data.Text as T-import Data.Type.Equality-import qualified Data.Vector.Unboxed as VU-import DataFrame.Operators-import Debug.Trace (trace)-import Type.Reflection (typeRep)+-- | How a candidate's output column is scored against the target (higher is better).+data LossFunction+    = -- | Pearson @r²@: scale-invariant, the default for derived features.+      PearsonCorrelation+    | -- | Binned mutual information: captures nonlinear association.+      MutualInformation+    | -- | Negative mean squared error: for reproducing a target exactly.+      MeanSquaredError+    deriving (Eq, Show) -generateConditions ::-    TypedColumn Double -> [Expr Bool] -> [Expr Double] -> DataFrame -> [Expr Bool]-generateConditions labels conds ps df =-    let-        newConds =-            [ p .<= q-            | p <- filter (not . isLiteral) ps-            , q <- ps-            , Prelude.not (eqExpr p q)-            ]-                ++ [ F.not p-                   | p <- conds-                   ]-        expandedConds =-            conds-                ++ newConds-                ++ [p .&& q | p <- newConds, q <- conds, Prelude.not (eqExpr p q)]-                ++ [p .|| q | p <- newConds, q <- conds, Prelude.not (eqExpr p q)]-     in-        pickTopNBool df labels (deduplicate df expandedConds)+-- | Search hyperparameters.+data SynthesisConfig = SynthesisConfig+    { synMaxSize :: !Int+    -- ^ Largest AST (node count) to enumerate.+    , synBankCap :: !Int+    -- ^ Max observationally-distinct programs kept per size layer.+    , synLoss :: !LossFunction+    , synTopK :: !Int+    -- ^ How many ranked features to return in the bank.+    }+    deriving (Eq, Show) -generatePrograms ::-    Bool ->-    [Expr Bool] ->-    [Expr Double] ->-    [Expr Double] ->-    [Expr Double] ->-    [Expr Double]-generatePrograms _ _ vars' constants [] = vars' ++ constants-generatePrograms includeConds conds vars constants ps =-    let-        existingPrograms = ps ++ vars ++ constants-     in-        existingPrograms-            ++ [ transform p-               | p <- ps ++ vars-               , Prelude.not (isConditional p)-               , transform <--                    [ sqrt-                    , abs-                    , log . (+ Lit 1)-                    , exp-                    , sin-                    , cos-                    , F.relu-                    , signum-                    ]-               ]-            ++ [ F.pow p i-               | p <- existingPrograms-               , Prelude.not (isConditional p)-               , i <- [2 .. 6]-               ]-            ++ [ p + q-               | (i, p) <- zip [(0 :: Int) ..] existingPrograms-               , (j, q) <- zip [(0 :: Int) ..] existingPrograms-               , Prelude.not (isLiteral p && isLiteral q)-               , Prelude.not (isConditional p || isConditional q)-               , i >= j-               ]-            ++ [ p - q-               | (i, p) <- zip [(0 :: Int) ..] existingPrograms-               , (j, q) <- zip [(0 :: Int) ..] existingPrograms-               , Prelude.not (isLiteral p && isLiteral q)-               , Prelude.not (isConditional p || isConditional q)-               , i /= j-               ]-            ++ ( if includeConds-                    then-                        [ F.min p q-                        | (i, p) <- zip [(0 :: Int) ..] existingPrograms-                        , (j, q) <- zip [(0 :: Int) ..] existingPrograms-                        , Prelude.not (isLiteral p && isLiteral q)-                        , Prelude.not (isConditional p || isConditional q)-                        , Prelude.not (eqExpr p q)-                        , i > j-                        ]-                            ++ [ F.max p q-                               | (i, p) <- zip [(0 :: Int) ..] existingPrograms-                               , (j, q) <- zip [(0 :: Int) ..] existingPrograms-                               , Prelude.not (isLiteral p && isLiteral q)-                               , Prelude.not (isConditional p || isConditional q)-                               , Prelude.not (eqExpr p q)-                               , i > j-                               ]-                            ++ [ F.ifThenElse cond r s-                               | cond <- conds-                               , r <- existingPrograms-                               , s <- existingPrograms-                               , Prelude.not (isConditional r || isConditional s)-                               , Prelude.not (eqExpr r s)-                               ]-                    else []-               )-            ++ [ p * q-               | (i, p) <- zip [(0 :: Int) ..] existingPrograms-               , (j, q) <- zip [(0 :: Int) ..] existingPrograms-               , Prelude.not (isLiteral p && isLiteral q)-               , Prelude.not (isConditional p || isConditional q)-               , i >= j-               ]-            ++ [ p / q-               | p <- existingPrograms-               , q <- existingPrograms-               , Prelude.not (isLiteral p && isLiteral q)-               , Prelude.not (isConditional p || isConditional q)-               , Prelude.not (eqExpr p q)-               ]+defaultSynthesisConfig :: SynthesisConfig+defaultSynthesisConfig =+    SynthesisConfig+        { synMaxSize = 6+        , synBankCap = 500+        , synLoss = PearsonCorrelation+        , synTopK = 16+        } -isLiteral :: Expr a -> Bool-isLiteral (Lit _) = True-isLiteral _ = False+{- | A synthesized feature. 'sfExpr' is the best-scoring expression and 'sfFeatures'+is the ranked, observationally-distinct bank (expression and its score).+-}+data SynthesizedFeature = SynthesizedFeature+    { sfExpr :: !(Expr Double)+    , sfScore :: !Double+    , sfFeatures :: ![(Expr Double, Double)]+    } -isConditional :: Expr a -> Bool-isConditional (If{}) = True-isConditional _ = False+instance Fit SynthesisConfig (Expr Double) SynthesizedFeature where+    fit = synthesizeFeatures -deduplicate ::-    forall a.-    (Columnable a) =>-    DataFrame ->-    [Expr a] ->-    [(Expr a, TypedColumn a)]-deduplicate df = go [] . L.nubBy eqExpr . L.sortBy (\e1 e2 -> compare (eSize e1) (eSize e2))+instance Predict SynthesizedFeature Double where+    predict = sfExpr++-- | A candidate's evaluated column over the example rows.+type Output = VU.Vector Double++data Prog = Prog+    { progExpr :: !(Expr Double)+    , progSize :: !Int+    , progOut :: !Output+    }++-- | Search for expressions over the non-target columns that track @target@.+synthesizeFeatures ::+    SynthesisConfig -> Expr Double -> DataFrame -> SynthesizedFeature+synthesizeFeatures cfg target df+    | null leaves || VU.null tgt = SynthesizedFeature (Lit 0) (negate (1 / 0)) []+    | otherwise = SynthesizedFeature best bestScore ranked   where-    go _ [] = []-    go seen (x : xs)-        | hasInvalid = go seen xs-        | res `elem` seen = go seen xs-        | otherwise = (x, res) : go (res : seen) xs-      where-        res = case interpret @a df x of-            Left e -> throw e-            Right v -> v-        hasInvalid = case res of-            (TColumn (UnboxedColumn _ (column :: VU.Vector b))) -> case testEquality (typeRep @Double) (typeRep @b) of-                Just Refl -> VU.any (\n -> isNaN n || isInfinite n) column-                Nothing -> False-            _ -> False+    feats = featureNames target df+    tgt = fromRight VU.empty (columnAsDoubleVector target df)+    n = VU.length tgt+    leaves = mkLeaves df feats n+    bank = grow cfg tgt leaves+    scored =+        [ (progExpr p, progSize p, s)+        | p <- bank+        , Just s <- [scoreOf (synLoss cfg) tgt (progOut p)]+        ]+    sorted = sortBy (comparing (\(_, sz, s) -> (Down s, sz))) scored+    ranked = [(e, s) | (e, _, s) <- take (synTopK cfg) sorted]+    (best, bestScore) = case ranked of+        ((e, s) : _) -> (e, s)+        [] -> (Lit 0, negate (1 / 0)) --- | Checks if two programs generate the same outputs given all the same inputs.-equivalent :: DataFrame -> Expr Double -> Expr Double -> Bool-equivalent df p1 p2 = case (==) <$> interpret df p1 <*> interpret df p2 of-    Left e -> throw e-    Right v -> v+-- | Size-1 programs: numeric feature columns and a pool of constants, OE-deduped.+mkLeaves :: DataFrame -> [T.Text] -> Int -> [Prog]+mkLeaves df feats n = fst (dedupProgs M.empty candidates)+  where+    candidates =+        [ Prog (Col name) 1 o+        | name <- feats+        , Right o <- [columnAsDoubleVector (Col name :: Expr Double) df]+        ]+            ++ [Prog (Lit v) 1 (VU.replicate n v) | v <- constantPool df feats] -synthesizeFeatureExpr ::-    -- | Target expression-    T.Text ->-    BeamConfig ->-    DataFrame ->-    Either String (Expr Double)-synthesizeFeatureExpr target cfg df =-    let-        df' = exclude [target] df-        t = case interpret df (Col target) of-            Left e -> throw e-            Right v -> v-     in-        case beamSearch-            df'-            cfg-            t-            (percentiles df')-            []-            [] of-            Nothing -> Left "No programs found"-            Just p -> Right p+{- | Domain-informed constants: per-column quartiles, variance, and std, plus a few+small integers. (Duplicates collapse under observational equivalence.)+-}+constantPool :: DataFrame -> [T.Text] -> [Double]+constantPool df feats =+    [0, 1, 2, -1]+        ++ [ roundSig 2 v+           | name <- feats+           , Right c <- [columnAsDoubleVector (Col name :: Expr Double) df]+           , v <-+                [percentile' p c | p <- [1, 25, 75, 99]] ++ [variance' c, sqrt (variance' c)]+           ] -f1FromBinary :: VU.Vector Double -> VU.Vector Double -> Maybe Double-f1FromBinary trues preds =-    let (!tp, !fp, !fn) =-            VU.foldl' step (0 :: Int, 0 :: Int, 0 :: Int) $-                VU.zip (VU.map (> 0) preds) (VU.map (> 0) trues)-     in f1FromCounts tp fp fn+-- | Grow the bank one size layer at a time, keeping one program per OE class.+grow :: SynthesisConfig -> Output -> [Prog] -> [Prog]+grow cfg tgt leaves = go 2 leaves (foldr (seenInsert . progOut) M.empty leaves)   where-    step (!tp, !fp, !fn) (!p, !t) =-        case (p, t) of-            (True, True) -> (tp + 1, fp, fn)-            (True, False) -> (tp, fp + 1, fn)-            (False, True) -> (tp, fp, fn + 1)-            (False, False) -> (tp, fp, fn)+    go size bank seen+        | size > synMaxSize cfg = bank+        | otherwise =+            let (kept, seen') = absorb cfg tgt seen (layer size bank)+             in go (size + 1) (bank ++ kept) seen' -f1FromCounts :: Int -> Int -> Int -> Maybe Double-f1FromCounts tp fp fn =-    let tp' = fromIntegral tp-        fp' = fromIntegral fp-        fn' = fromIntegral fn-        precision = if tp' + fp' == 0 then 0 else tp' / (tp' + fp')-        recall = if tp' + fn' == 0 then 0 else tp' / (tp' + fn')-     in if precision + recall == 0-            then Nothing-            else Just (2 * precision * recall / (precision + recall))+-- | All candidate programs of exactly @size@ nodes, built from smaller ones.+layer :: Int -> [Prog] -> [Prog]+layer size bank = unaries ++ pows ++ comms ++ subs ++ divs+  where+    atSize s = filter ((== s) . progSize) bank+    args1 = atSize (size - 1)+    unaries =+        [ Prog (mk e) size (VU.map f o)+        | (mk, f) <- unaryProds+        , Prog e _ o <- args1+        ]+    pows =+        [ Prog (F.pow e k) size (VU.map (^ k) o)+        | Prog e _ o <- args1+        , k <- [2 .. 6 :: Int]+        ]+    comms =+        [ Prog (mk ea eb) size (VU.zipWith f oa ob)+        | (mk, f) <- commutativeProds+        , (Prog ea _ oa, Prog eb _ ob) <- unorderedPairs size bank+        ]+    subs =+        [ Prog (ea - eb) size (VU.zipWith (-) oa ob)+        | (Prog ea _ oa, Prog eb _ ob) <- orderedPairs size bank+        ]+    divs =+        [ Prog (ea / eb) size (VU.zipWith (/) oa ob)+        | (Prog ea _ oa, Prog eb _ ob) <- orderedPairs size bank+        , VU.all ((> 1e-9) . abs) ob+        ] -fitClassifier ::-    -- | Target expression-    T.Text ->-    -- | Depth of search (Roughly, how many terms in the final expression)-    Int ->-    -- | Beam size - the number of candidate expressions to consider at a time.-    Int ->-    DataFrame ->-    Either String (Expr Int)-fitClassifier target d b df =-    let-        df' = exclude [target] df-        t = case interpret df (Col target) of-            Left e -> throw e-            Right v -> v-     in-        case beamSearch-            df'-            (BeamConfig d b F1 True)-            t-            (percentiles df' ++ [Lit 1, Lit 0, Lit (-1)])-            []-            [] of-            Nothing -> Left "No programs found"-            Just p -> Right (F.ifThenElse (p .> (0 :: Expr Double)) 1 0)+-- | Protected unary operators: total on all reals (no NaN/domain errors).+unaryProds :: [(Expr Double -> Expr Double, Double -> Double)]+unaryProds =+    [ (sqrt . abs, sqrt . abs)+    , (abs, abs)+    , (\e -> log (abs e + 1), \x -> log (abs x + 1))+    , (exp, exp)+    , (sin, sin)+    , (cos, cos)+    , (F.relu, max 0)+    , (signum, signum)+    ] -percentiles :: DataFrame -> [Expr Double]-percentiles df =-    let-        doubleColumns =-            map-                (either throw id . ((`columnAsDoubleVector` df) . Col @Double))-                (columnNames df)-     in-        concatMap-            (\c -> map (Lit . roundTo2SigDigits . (`percentile'` c)) [1, 25, 75, 99])-            doubleColumns-            ++ map (Lit . roundTo2SigDigits . variance') doubleColumns-            ++ map (Lit . roundTo2SigDigits . sqrt . variance') doubleColumns+-- | Commutative binary operators (enumerated over unordered operand pairs).+commutativeProds ::+    [(Expr Double -> Expr Double -> Expr Double, Double -> Double -> Double)]+commutativeProds =+    [ ((+), (+))+    , ((*), (*))+    , (F.min, min)+    , (F.max, max)+    ] -roundToSigDigits :: Int -> Double -> Double-roundToSigDigits n x-    | x == 0 = 0-    | otherwise =-        let magnitude = floor (logBase 10 (abs x))-            scale = 10 ** fromIntegral (n - 1 - magnitude)-         in fromIntegral (round (x * scale) :: Int) / scale+-- | Ordered operand pairs whose sizes sum to @size-1@ (for non-commutative ops).+orderedPairs :: Int -> [Prog] -> [(Prog, Prog)]+orderedPairs size bank =+    [ (a, b)+    | sa <- [1 .. size - 2]+    , let sb = size - 1 - sa+    , sb >= 1+    , a <- atSize sa+    , b <- atSize sb+    ]+  where+    atSize s = filter ((== s) . progSize) bank -roundTo2SigDigits :: Double -> Double-roundTo2SigDigits = roundToSigDigits 2+-- | Unordered operand pairs (for commutative ops): each pair once.+unorderedPairs :: Int -> [Prog] -> [(Prog, Prog)]+unorderedPairs size bank =+    [ (a, b)+    | sa <- [1 .. size - 2]+    , let sb = size - 1 - sa+    , sb >= 1+    , sa <= sb+    , (i, a) <- zip [0 :: Int ..] (atSize sa)+    , (j, b) <- zip [0 :: Int ..] (atSize sb)+    , sa < sb || i <= j+    ]+  where+    atSize s = filter ((== s) . progSize) bank -fitRegression ::-    -- | Target expression-    T.Text ->-    -- | Depth of search (Roughly, how many terms in the final expression)-    Int ->-    -- | Beam size - the number of candidate expressions to consider at a time.-    Int ->-    DataFrame ->-    Either String (Expr Double)-fitRegression target d b df =-    let-        df' = exclude [target] df-        targetMean = Stats.mean (Col @Double target) df-        t = case interpret df (Col target) of-            Left e -> throw e-            Right v -> v-        cfg = BeamConfig d b MeanSquaredError True-        constants =-            percentiles df'-                ++ [Lit targetMean]-                ++ [ F.pow p i-                   | i <- [1 .. 6]-                   , p <- [Lit 10, Lit 1, Lit 0.1]-                   ]-     in-        case beamSearch df' cfg t constants [] [] of-            Nothing -> Left "No programs found"-            Just p -> Right p+{- | Keep the valid, observationally-novel candidates of a layer, then cap by fit+score (cost-guided). Returns the kept programs and the updated OE-class set.+-}+absorb ::+    SynthesisConfig -> Output -> Seen -> [Prog] -> ([Prog], Seen)+absorb cfg tgt seen0 cands = (capLayer cfg tgt fresh, seen')+  where+    (fresh, seen') = dedupProgs seen0 cands -data LossFunction-    = PearsonCorrelation-    | MutualInformation-    | MeanSquaredError-    | F1+-- | When a layer has more distinct programs than the cap, keep the best-scoring.+capLayer :: SynthesisConfig -> Output -> [Prog] -> [Prog]+capLayer cfg tgt progs+    | length progs <= synBankCap cfg = progs+    | otherwise = take (synBankCap cfg) (sortBy (comparing (Down . rank)) progs)+  where+    rank p = fromMaybe (negate (1 / 0)) (scoreOf (synLoss cfg) tgt (progOut p)) -getLossFunction ::-    LossFunction -> (VU.Vector Double -> VU.Vector Double -> Maybe Double)-getLossFunction f = case f of-    MutualInformation ->-        ( \l r ->-            mutualInformationBinned-                (Prelude.max 10 (ceiling (sqrt (fromIntegral (VU.length l) :: Double))))-                l-                r-        )-    PearsonCorrelation -> (\l r -> (^ (2 :: Int)) <$> correlation' l r)-    MeanSquaredError -> (\l r -> fmap negate (meanSquaredError l r))-    F1 -> f1FromBinary+-- | Fit score of an output against the target (higher is better), or @Nothing@.+scoreOf :: LossFunction -> Output -> Output -> Maybe Double+scoreOf lf tgt out+    | VU.length out /= VU.length tgt = Nothing+    | otherwise = finite $ case lf of+        PearsonCorrelation -> pearsonR2 tgt out+        MutualInformation -> mutualInformationBinned bins tgt out+        MeanSquaredError -> negate <$> meanSquaredError tgt out+  where+    bins = max 10 (ceiling (sqrt (fromIntegral (VU.length tgt) :: Double)))+    -- Belt-and-suspenders: drop any non-finite score so it cannot win the ranking.+    finite (Just s) | isNaN s || isInfinite s = Nothing+    finite ms = ms -data BeamConfig = BeamConfig-    { searchDepth :: Int-    , beamLength :: Int-    , lossFunction :: LossFunction-    , includeConditionals :: Bool-    }+{- | Pearson @r²@ via the numerically stable centered two-pass formula. Returns+'Nothing' when the feature (or target) is constant, and is bounded by+Cauchy–Schwarz to @[0,1]@ — unlike the one-pass @n·Σxy − Σx·Σy@ form, which+cancels catastrophically for low-variance features and can report @r² > 1@.+-}+pearsonR2 :: Output -> Output -> Maybe Double+pearsonR2 ys xs+    | n < 2 = Nothing+    | sxx <= 0 || syy <= 0 = Nothing+    | otherwise = Just (min 1 (sxy * sxy / (sxx * syy)))+  where+    n = VU.length xs+    nf = fromIntegral n+    mx = VU.sum xs / nf+    my = VU.sum ys / nf+    sxy = VU.sum (VU.zipWith (\x y -> (x - mx) * (y - my)) xs ys)+    sxx = VU.sum (VU.map (\x -> (x - mx) * (x - mx)) xs)+    syy = VU.sum (VU.map (\y -> (y - my) * (y - my)) ys) -defaultBeamConfig :: BeamConfig-defaultBeamConfig = BeamConfig 2 100 PearsonCorrelation False+-- | An output is usable iff it is non-empty and free of NaN/±Inf.+valid :: Output -> Bool+valid o = not (VU.null o) && VU.all (\x -> not (isNaN x || isInfinite x)) o -beamSearch ::-    DataFrame ->-    -- | Parameters of the beam search.-    BeamConfig ->-    -- | Examples-    TypedColumn Double ->-    -- | Constants-    [Expr Double] ->-    -- | Conditions-    [Expr Bool] ->-    -- | Programs-    [Expr Double] ->-    Maybe (Expr Double)-beamSearch df cfg outputs constants conds programs-    | searchDepth cfg == 0 = case ps of-        [] -> Nothing-        (x : _) -> Just x-    | otherwise =-        beamSearch-            df-            (cfg{searchDepth = searchDepth cfg - 1})-            outputs-            constants-            conditions-            (generatePrograms (includeConditionals cfg) conditions vars constants ps)+{- | Quantize an output to nine significant digits, so float noise (@x*2@ vs+@x+x@) collapses while genuinely distinct features stay apart.+-}+quantize :: Output -> Output+quantize = VU.map (\x -> if x == 0 then 0 else signum x * roundSig 9 (abs x))++{- | The observational-equivalence class set: a map from a 64-bit FNV-1a+fingerprint of the quantized output to the (usually one) quantized outputs with+that fingerprint. Bucketing on the fingerprint keeps membership cheap, and+verifying exact equality within the bucket makes a hash collision harmless — two+genuinely different columns that happen to collide are kept apart, not merged.+-}+type Seen = M.Map Int [Output]++-- | FNV-1a fingerprint of an already-quantized output's bit patterns.+fpOf :: Output -> Int+fpOf = fromIntegral . VU.foldl' step (1469598103934665603 :: Word64)   where-    vars = map Col names-    conditions = generateConditions outputs conds (vars ++ constants) df-    ps = pickTopN df outputs cfg $ deduplicate df programs-    names = (map fst . L.sortBy (compare `on` snd) . M.toList . columnIndices) df+    step !h x = (h `xor` castDoubleToWord64 x) * 1099511628211 -pickTopN ::-    DataFrame ->-    TypedColumn Double ->-    BeamConfig ->-    [(Expr Double, TypedColumn a)] ->-    [Expr Double]-pickTopN _ _ _ [] = []-pickTopN df (TColumn column) cfg ps =-    let-        l = case toVector @Double @VU.Vector column of-            Left e -> throw e-            Right v -> v-        ordered =-            Prelude.take-                (beamLength cfg)-                ( map fst $-                    L.sortBy-                        ( \(_, c2) (_, c1) ->-                            if maybe False isInfinite c1-                                || maybe False isInfinite c2-                                || maybe False isNaN c1-                                || maybe False isNaN c2-                                then LT-                                else compare c1 c2-                        )-                        ( map-                            (\(e, res) -> (e, getLossFunction (lossFunction cfg) l (asDoubleVector res)))-                            ps-                        )-                )-        asDoubleVector c =-            let-                (TColumn col') = c-             in-                case toVector @Double @VU.Vector col' of-                    Left e -> throw e-                    Right v -> VU.convert v-        interpretDoubleVector e' =-            let-                (TColumn col') = case interpret df e' of-                    Left err -> throw err-                    Right v -> v-             in-                case toVector @Double @VU.Vector col' of-                    Left err -> throw err-                    Right v -> VU.convert v-     in-        trace-            ( "Best loss: "-                ++ show-                    ( getLossFunction (lossFunction cfg) l . interpretDoubleVector-                        <$> listToMaybe ordered-                    )-                ++ " "-                ++ (if null ordered then "empty" else show (listToMaybe ordered))-            )-            ordered+-- | Record an output's observational-equivalence class.+seenInsert :: Output -> Seen -> Seen+seenInsert o = M.insertWith (++) (fpOf q) [q]+  where+    q = quantize o -pickTopNBool ::-    DataFrame ->-    TypedColumn Double ->-    [(Expr Bool, TypedColumn Bool)] ->-    [Expr Bool]-pickTopNBool _ _ [] = []-pickTopNBool _df (TColumn column) ps =-    let-        l = case toVector @Double @VU.Vector column of-            Left e -> throw e-            Right v -> v-        ordered =-            Prelude.take-                10-                ( map fst $-                    L.sortBy-                        ( \(_, c2) (_, c1) ->-                            if maybe False isInfinite c1-                                || maybe False isInfinite c2-                                || maybe False isNaN c1-                                || maybe False isNaN c2-                                then LT-                                else compare c1 c2-                        )-                        ( map-                            (\(e, res) -> (e, getLossFunction MutualInformation l (asDoubleVector res)))-                            ps-                        )-                )-        asDoubleVector c =-            let-                (TColumn col') = c-             in-                case toVector @Bool @VU.Vector col' of-                    Left e -> throw e-                    Right v -> VU.map (fromIntegral @Int @Double . fromEnum) v-     in-        ordered+-- | Round a positive double to @n@ significant digits.+roundSig :: Int -> Double -> Double+roundSig n x+    | x == 0 = 0+    | otherwise =+        let magnitude = floor (logBase 10 (abs x)) :: Int+            scale = 10 ** fromIntegral (n - 1 - magnitude)+         in fromIntegral (round (x * scale) :: Integer) / scale -satisfiesExamples :: DataFrame -> TypedColumn Double -> Expr Double -> Bool-satisfiesExamples df column expr =-    let-        result = case interpret df expr of-            Left e -> throw e-            Right v -> v-     in-        result == column+{- | Keep the first valid program of each observational-equivalence class,+preserving order; returns the kept programs and the grown class set.+-}+dedupProgs :: Seen -> [Prog] -> ([Prog], Seen)+dedupProgs = go []+  where+    go acc s [] = (reverse acc, s)+    go acc s (p : ps)+        | not (valid o) = go acc s ps+        | member = go acc s ps+        | otherwise = go (p : acc) (M.insertWith (++) fp [q] s) ps+      where+        o = progOut p+        q = quantize o+        fp = fpOf q+        member = maybe False (q `elem`) (M.lookup fp s)