diff --git a/README.md b/README.md new file mode 100644 --- /dev/null +++ b/README.md @@ -0,0 +1,473 @@ + + +# 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)) +``` + +> +> 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)) +``` + +> +> 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) ++ ")") +``` + +> +> 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" ]) +``` + +> +> 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) +``` + +> +> ["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 ]) +``` + +> +> 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 +``` + +> +> 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 +``` + +> + +**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) ]) +``` + +> +> 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) ]) +``` + +> +> 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))) +``` + +> +> 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)) +``` + +> +> 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) +``` + +> +> ["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 +``` + +> +> [[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) + ]) +``` + +> +> 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.) diff --git a/dataframe-learn.cabal b/dataframe-learn.cabal --- a/dataframe-learn.cabal +++ b/dataframe-learn.cabal @@ -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 diff --git a/src/DataFrame/Boosting.hs b/src/DataFrame/Boosting.hs new file mode 100644 --- /dev/null +++ b/src/DataFrame/Boosting.hs @@ -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 diff --git a/src/DataFrame/Boosting/AdaBoost.hs b/src/DataFrame/Boosting/AdaBoost.hs new file mode 100644 --- /dev/null +++ b/src/DataFrame/Boosting/AdaBoost.hs @@ -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) diff --git a/src/DataFrame/Boosting/GBM.hs b/src/DataFrame/Boosting/GBM.hs new file mode 100644 --- /dev/null +++ b/src/DataFrame/Boosting/GBM.hs @@ -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 diff --git a/src/DataFrame/DBSCAN.hs b/src/DataFrame/DBSCAN.hs new file mode 100644 --- /dev/null +++ b/src/DataFrame/DBSCAN.hs @@ -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)))] diff --git a/src/DataFrame/DecisionTree/Cart.hs b/src/DataFrame/DecisionTree/Cart.hs --- a/src/DataFrame/DecisionTree/Cart.hs +++ b/src/DataFrame/DecisionTree/Cart.hs @@ -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 diff --git a/src/DataFrame/DecisionTree/Categorical.hs b/src/DataFrame/DecisionTree/Categorical.hs --- a/src/DataFrame/DecisionTree/Categorical.hs +++ b/src/DataFrame/DecisionTree/Categorical.hs @@ -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) diff --git a/src/DataFrame/DecisionTree/CondVec.hs b/src/DataFrame/DecisionTree/CondVec.hs --- a/src/DataFrame/DecisionTree/CondVec.hs +++ b/src/DataFrame/DecisionTree/CondVec.hs @@ -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)) diff --git a/src/DataFrame/DecisionTree/Fit.hs b/src/DataFrame/DecisionTree/Fit.hs --- a/src/DataFrame/DecisionTree/Fit.hs +++ b/src/DataFrame/DecisionTree/Fit.hs @@ -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] diff --git a/src/DataFrame/DecisionTree/Linear.hs b/src/DataFrame/DecisionTree/Linear.hs --- a/src/DataFrame/DecisionTree/Linear.hs +++ b/src/DataFrame/DecisionTree/Linear.hs @@ -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 [] -> "" --- | 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]) diff --git a/src/DataFrame/DecisionTree/Model.hs b/src/DataFrame/DecisionTree/Model.hs new file mode 100644 --- /dev/null +++ b/src/DataFrame/DecisionTree/Model.hs @@ -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 diff --git a/src/DataFrame/DecisionTree/Numeric.hs b/src/DataFrame/DecisionTree/Numeric.hs --- a/src/DataFrame/DecisionTree/Numeric.hs +++ b/src/DataFrame/DecisionTree/Numeric.hs @@ -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 = [] diff --git a/src/DataFrame/DecisionTree/Pool.hs b/src/DataFrame/DecisionTree/Pool.hs --- a/src/DataFrame/DecisionTree/Pool.hs +++ b/src/DataFrame/DecisionTree/Pool.hs @@ -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 [] -> "" @@ -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) diff --git a/src/DataFrame/DecisionTree/Predict.hs b/src/DataFrame/DecisionTree/Predict.hs --- a/src/DataFrame/DecisionTree/Predict.hs +++ b/src/DataFrame/DecisionTree/Predict.hs @@ -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) diff --git a/src/DataFrame/DecisionTree/Prune.hs b/src/DataFrame/DecisionTree/Prune.hs --- a/src/DataFrame/DecisionTree/Prune.hs +++ b/src/DataFrame/DecisionTree/Prune.hs @@ -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) diff --git a/src/DataFrame/DecisionTree/Regression.hs b/src/DataFrame/DecisionTree/Regression.hs new file mode 100644 --- /dev/null +++ b/src/DataFrame/DecisionTree/Regression.hs @@ -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) diff --git a/src/DataFrame/DecisionTree/Tao.hs b/src/DataFrame/DecisionTree/Tao.hs --- a/src/DataFrame/DecisionTree/Tao.hs +++ b/src/DataFrame/DecisionTree/Tao.hs @@ -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) diff --git a/src/DataFrame/DecisionTree/Types.hs b/src/DataFrame/DecisionTree/Types.hs --- a/src/DataFrame/DecisionTree/Types.hs +++ b/src/DataFrame/DecisionTree/Types.hs @@ -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 diff --git a/src/DataFrame/Featurize/Internal.hs b/src/DataFrame/Featurize/Internal.hs new file mode 100644 --- /dev/null +++ b/src/DataFrame/Featurize/Internal.hs @@ -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) diff --git a/src/DataFrame/GMM.hs b/src/DataFrame/GMM.hs new file mode 100644 --- /dev/null +++ b/src/DataFrame/GMM.hs @@ -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) diff --git a/src/DataFrame/KMeans.hs b/src/DataFrame/KMeans.hs new file mode 100644 --- /dev/null +++ b/src/DataFrame/KMeans.hs @@ -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] diff --git a/src/DataFrame/LinearAlgebra.hs b/src/DataFrame/LinearAlgebra.hs new file mode 100644 --- /dev/null +++ b/src/DataFrame/LinearAlgebra.hs @@ -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 diff --git a/src/DataFrame/LinearAlgebra/Eigen.hs b/src/DataFrame/LinearAlgebra/Eigen.hs new file mode 100644 --- /dev/null +++ b/src/DataFrame/LinearAlgebra/Eigen.hs @@ -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) diff --git a/src/DataFrame/LinearAlgebra/Solve.hs b/src/DataFrame/LinearAlgebra/Solve.hs new file mode 100644 --- /dev/null +++ b/src/DataFrame/LinearAlgebra/Solve.hs @@ -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) diff --git a/src/DataFrame/LinearModel.hs b/src/DataFrame/LinearModel.hs new file mode 100644 --- /dev/null +++ b/src/DataFrame/LinearModel.hs @@ -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 diff --git a/src/DataFrame/LinearModel/Logistic.hs b/src/DataFrame/LinearModel/Logistic.hs new file mode 100644 --- /dev/null +++ b/src/DataFrame/LinearModel/Logistic.hs @@ -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)) diff --git a/src/DataFrame/LinearModel/Regression.hs b/src/DataFrame/LinearModel/Regression.hs new file mode 100644 --- /dev/null +++ b/src/DataFrame/LinearModel/Regression.hs @@ -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) diff --git a/src/DataFrame/LinearSolver.hs b/src/DataFrame/LinearSolver.hs --- a/src/DataFrame/LinearSolver.hs +++ b/src/DataFrame/LinearSolver.hs @@ -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) diff --git a/src/DataFrame/LinearSolver/Loss.hs b/src/DataFrame/LinearSolver/Loss.hs new file mode 100644 --- /dev/null +++ b/src/DataFrame/LinearSolver/Loss.hs @@ -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 diff --git a/src/DataFrame/Metrics.hs b/src/DataFrame/Metrics.hs new file mode 100644 --- /dev/null +++ b/src/DataFrame/Metrics.hs @@ -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 diff --git a/src/DataFrame/Metrics/Report.hs b/src/DataFrame/Metrics/Report.hs new file mode 100644 --- /dev/null +++ b/src/DataFrame/Metrics/Report.hs @@ -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) diff --git a/src/DataFrame/Model.hs b/src/DataFrame/Model.hs new file mode 100644 --- /dev/null +++ b/src/DataFrame/Model.hs @@ -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 _ = [] diff --git a/src/DataFrame/ModelSelection.hs b/src/DataFrame/ModelSelection.hs new file mode 100644 --- /dev/null +++ b/src/DataFrame/ModelSelection.hs @@ -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) diff --git a/src/DataFrame/PCA.hs b/src/DataFrame/PCA.hs new file mode 100644 --- /dev/null +++ b/src/DataFrame/PCA.hs @@ -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 _ _ _ _ = [] diff --git a/src/DataFrame/PCA/Kernel.hs b/src/DataFrame/PCA/Kernel.hs new file mode 100644 --- /dev/null +++ b/src/DataFrame/PCA/Kernel.hs @@ -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] diff --git a/src/DataFrame/Random.hs b/src/DataFrame/Random.hs new file mode 100644 --- /dev/null +++ b/src/DataFrame/Random.hs @@ -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') diff --git a/src/DataFrame/SVM.hs b/src/DataFrame/SVM.hs new file mode 100644 --- /dev/null +++ b/src/DataFrame/SVM.hs @@ -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))) diff --git a/src/DataFrame/SVM/RFF.hs b/src/DataFrame/SVM/RFF.hs new file mode 100644 --- /dev/null +++ b/src/DataFrame/SVM/RFF.hs @@ -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) diff --git a/src/DataFrame/SymbolicRegression.hs b/src/DataFrame/SymbolicRegression.hs new file mode 100644 --- /dev/null +++ b/src/DataFrame/SymbolicRegression.hs @@ -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 diff --git a/src/DataFrame/SymbolicRegression/Expr.hs b/src/DataFrame/SymbolicRegression/Expr.hs new file mode 100644 --- /dev/null +++ b/src/DataFrame/SymbolicRegression/Expr.hs @@ -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'') diff --git a/src/DataFrame/SymbolicRegression/GP.hs b/src/DataFrame/SymbolicRegression/GP.hs new file mode 100644 --- /dev/null +++ b/src/DataFrame/SymbolicRegression/GP.hs @@ -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 diff --git a/src/DataFrame/SymbolicRegression/Optimize.hs b/src/DataFrame/SymbolicRegression/Optimize.hs new file mode 100644 --- /dev/null +++ b/src/DataFrame/SymbolicRegression/Optimize.hs @@ -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) diff --git a/src/DataFrame/SymbolicRegression/Simplify.hs b/src/DataFrame/SymbolicRegression/Simplify.hs new file mode 100644 --- /dev/null +++ b/src/DataFrame/SymbolicRegression/Simplify.hs @@ -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 diff --git a/src/DataFrame/Synthesis.hs b/src/DataFrame/Synthesis.hs --- a/src/DataFrame/Synthesis.hs +++ b/src/DataFrame/Synthesis.hs @@ -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) diff --git a/src/DataFrame/Transform.hs b/src/DataFrame/Transform.hs new file mode 100644 --- /dev/null +++ b/src/DataFrame/Transform.hs @@ -0,0 +1,103 @@ +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE TypeApplications #-} + +{- | Fitted column transforms as a composable monoid. A 'Transform' is a list of +named output expressions; @s <> t@ means \"apply @s@, then @t@\", fusing @t@'s +references to @s@'s outputs by simultaneous substitution. 'applyTransform' runs +one against a frame; 'compileThrough' folds a transform into a model's +prediction expression so the result is a single expression over the raw inputs. + +Every right-hand side must be row-wise (no aggregation/window), and within one +transform each expression reads the original frame. +-} +module DataFrame.Transform ( + Transform (..), + applyTransform, + compileThrough, + ScalerModel (..), + standardScaler, + scalerTransform, +) where + +import Control.Exception (throw) +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 qualified DataFrame.Functions as F +import DataFrame.Internal.Column (Columnable) +import DataFrame.Internal.DataFrame (DataFrame) +import DataFrame.Internal.Expression ( + Expr (..), + NamedExpr, + UExpr (..), + substituteColumns, + ) +import DataFrame.Operations.Core (columnAsDoubleVector) +import DataFrame.Operations.Transformations (deriveMany) +import DataFrame.Operators ((.-.), (./.)) + +-- | A fitted transform: named output columns derived from the input frame. +newtype Transform = Transform {transformOutputs :: [NamedExpr]} + +instance Semigroup Transform where + Transform s <> Transform t = + Transform (s ++ map (subst (M.fromList s)) t) + where + subst :: M.Map T.Text UExpr -> NamedExpr -> NamedExpr + subst m (nm, UExpr e) = (nm, UExpr (substituteColumns m e)) + +instance Monoid Transform where + mempty = Transform [] + +-- | Apply a transform to a frame (deriving its outputs in order). +applyTransform :: Transform -> DataFrame -> DataFrame +applyTransform (Transform os) = deriveMany os + +{- | Fold a preprocessing transform into a model's prediction expression, +yielding one expression over the transform's input columns. +-} +compileThrough :: (Columnable a) => Transform -> Expr a -> Expr a +compileThrough (Transform os) = substituteColumns (M.fromList os) + +-- | A fitted standardizer: per-column means and standard deviations. +data ScalerModel = ScalerModel + { smColumns :: !(V.Vector T.Text) + , smMeans :: !(VU.Vector Double) + , smStds :: !(VU.Vector Double) + } + deriving (Eq, Show) + +-- | Fit a standard scaler over the named columns. +standardScaler :: [T.Text] -> DataFrame -> ScalerModel +standardScaler names df = + ScalerModel (V.fromList names) (VU.fromList means) (VU.fromList stds) + where + cols = map column names + column n = case columnAsDoubleVector (F.col @Double n) df of + Right v -> v + Left e -> throw e + means = [VU.sum c / fromIntegral (max 1 (VU.length c)) | c <- cols] + stds = + [ let mu = mean + v = + VU.sum (VU.map (\x -> (x - mu) ^ (2 :: Int)) c) + / fromIntegral (max 1 (VU.length c)) + s = sqrt v + in if s < 1e-12 then 1 else s + | (mean, c) <- zip means cols + ] + +-- | The scaler as a 'Transform': @(col - μ) / σ@ per column. +scalerTransform :: ScalerModel -> Transform +scalerTransform m = + Transform + [ (n, UExpr ((F.col @Double n .-. F.lit mu) ./. F.lit sigma)) + | (n, mu, sigma) <- + zip3 + (V.toList (smColumns m)) + (VU.toList (smMeans m)) + (VU.toList (smStds m)) + ]