-- SPDX-License-Identifier: MPL-2.0

{- | Exact branch-and-bound KD/ball trees and budgeted approximate search.
IDs always refer to original training rows. Neighbors sort by distance then
ID. LSH is Euclidean p-stable Gaussian projection hashing (Datar et al. 2004).
Isolation forest follows Liu et al. ICDM 2008, with seeded subsampling and the
harmonic-number correction for unsplit leaves. Larger anomaly scores indicate
more isolated observations. One-class SVM returns positive inlier decisions.
-}
module HMath.MachineLearning.Neighbors (
    IndexKind (..),
    SpatialIndex,
    buildIndex,
    NeighborResult (..),
    nearestNeighbors,
    knnRegress,
    knnClassify,
    LSH,
    buildLSH,
    queryLSH,
    localOutlierFactor,
    IsolationForest,
    fitIsolationForest,
    scoreIsolationForest,
    OneClassSVM,
    fitOneClassSVM,
    decisionOneClassSVM,
) where

import Control.Monad (foldM)
import Data.List (nub, sort, sortOn)
import Data.Word (Word64)
import HMath.Error
import HMath.MachineLearning.Internal
import HMath.MachineLearning.Regression (Kernel, kernel)
import qualified HMath.Random as R

-- | Axis-aligned KD partitions or enclosing Euclidean ball partitions.
data IndexKind = KDTree | BallTree deriving (IndexKind -> IndexKind -> Bool
(IndexKind -> IndexKind -> Bool)
-> (IndexKind -> IndexKind -> Bool) -> Eq IndexKind
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: IndexKind -> IndexKind -> Bool
== :: IndexKind -> IndexKind -> Bool
$c/= :: IndexKind -> IndexKind -> Bool
/= :: IndexKind -> IndexKind -> Bool
Eq, Int -> IndexKind -> ShowS
[IndexKind] -> ShowS
IndexKind -> String
(Int -> IndexKind -> ShowS)
-> (IndexKind -> String)
-> ([IndexKind] -> ShowS)
-> Show IndexKind
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> IndexKind -> ShowS
showsPrec :: Int -> IndexKind -> ShowS
$cshow :: IndexKind -> String
show :: IndexKind -> String
$cshowList :: [IndexKind] -> ShowS
showList :: [IndexKind] -> ShowS
Show)

data Tree = Leaf [(Int, [Double])] | KDNode Int Double Tree Tree | BallNode [Double] Double Tree Tree
    deriving (Tree -> Tree -> Bool
(Tree -> Tree -> Bool) -> (Tree -> Tree -> Bool) -> Eq Tree
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: Tree -> Tree -> Bool
== :: Tree -> Tree -> Bool
$c/= :: Tree -> Tree -> Bool
/= :: Tree -> Tree -> Bool
Eq, Int -> Tree -> ShowS
[Tree] -> ShowS
Tree -> String
(Int -> Tree -> ShowS)
-> (Tree -> String) -> ([Tree] -> ShowS) -> Show Tree
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> Tree -> ShowS
showsPrec :: Int -> Tree -> ShowS
$cshow :: Tree -> String
show :: Tree -> String
$cshowList :: [Tree] -> ShowS
showList :: [Tree] -> ShowS
Show)

-- | Immutable tree retaining original observation IDs and its validated feature dimension.
data SpatialIndex = SpatialIndex Int Int Tree deriving (SpatialIndex -> SpatialIndex -> Bool
(SpatialIndex -> SpatialIndex -> Bool)
-> (SpatialIndex -> SpatialIndex -> Bool) -> Eq SpatialIndex
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: SpatialIndex -> SpatialIndex -> Bool
== :: SpatialIndex -> SpatialIndex -> Bool
$c/= :: SpatialIndex -> SpatialIndex -> Bool
/= :: SpatialIndex -> SpatialIndex -> Bool
Eq, Int -> SpatialIndex -> ShowS
[SpatialIndex] -> ShowS
SpatialIndex -> String
(Int -> SpatialIndex -> ShowS)
-> (SpatialIndex -> String)
-> ([SpatialIndex] -> ShowS)
-> Show SpatialIndex
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> SpatialIndex -> ShowS
showsPrec :: Int -> SpatialIndex -> ShowS
$cshow :: SpatialIndex -> String
show :: SpatialIndex -> String
$cshowList :: [SpatialIndex] -> ShowS
showList :: [SpatialIndex] -> ShowS
Show)

-- | Build an exact search tree with at most leafSize observations per leaf.
buildIndex :: IndexKind -> Int -> Dataset -> Either MathError SpatialIndex
buildIndex :: IndexKind -> Int -> Dataset -> Either MathError SpatialIndex
buildIndex IndexKind
kind Int
leafSize Dataset
xs = do
    (n, d) <- Dataset -> Either MathError (Int, Int)
dataset Dataset
xs
    require (leafSize > 0) "tree leaf size must be positive"
    let build [(Int, [Double])]
points
            | [(Int, [Double])] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [(Int, [Double])]
points Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
leafSize = [(Int, [Double])] -> Tree
Leaf [(Int, [Double])]
points
            | Bool
otherwise = case IndexKind
kind of
                IndexKind
KDTree ->
                    let spans :: [Double]
spans = [[Double] -> Double
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum [Double]
col Double -> Double -> Double
forall a. Num a => a -> a -> a
- [Double] -> Double
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
minimum [Double]
col | [Double]
col <- Dataset -> Dataset
forall {a}. [[a]] -> [[a]]
transposeLocal (((Int, [Double]) -> [Double]) -> [(Int, [Double])] -> Dataset
forall a b. (a -> b) -> [a] -> [b]
map (Int, [Double]) -> [Double]
forall a b. (a, b) -> b
snd [(Int, [Double])]
points)]
                        axis :: Int
axis = [Double] -> Int
forall a. Ord a => [a] -> Int
argmax [Double]
spans
                        ordered :: [(Int, [Double])]
ordered = ((Int, [Double]) -> Double)
-> [(Int, [Double])] -> [(Int, [Double])]
forall b a. Ord b => (a -> b) -> [a] -> [a]
sortOn (\(Int
_, [Double]
r) -> [Double]
r [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
axis) [(Int, [Double])]
points
                        middle :: Int
middle = [(Int, [Double])] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [(Int, [Double])]
points Int -> Int -> Int
forall a. Integral a => a -> a -> a
`div` Int
2
                        cutoff :: Double
cutoff = (Int, [Double]) -> [Double]
forall a b. (a, b) -> b
snd ([(Int, [Double])]
ordered [(Int, [Double])] -> Int -> (Int, [Double])
forall a. HasCallStack => [a] -> Int -> a
!! Int
middle) [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
axis
                     in Int -> Double -> Tree -> Tree -> Tree
KDNode Int
axis Double
cutoff ([(Int, [Double])] -> Tree
build (Int -> [(Int, [Double])] -> [(Int, [Double])]
forall a. Int -> [a] -> [a]
take Int
middle [(Int, [Double])]
ordered)) ([(Int, [Double])] -> Tree
build (Int -> [(Int, [Double])] -> [(Int, [Double])]
forall a. Int -> [a] -> [a]
drop Int
middle [(Int, [Double])]
ordered))
                IndexKind
BallTree ->
                    let center' :: [Double]
center' = Dataset -> [Double]
meanV (((Int, [Double]) -> [Double]) -> [(Int, [Double])] -> Dataset
forall a b. (a -> b) -> [a] -> [b]
map (Int, [Double]) -> [Double]
forall a b. (a, b) -> b
snd [(Int, [Double])]
points)
                        radius :: Double
radius = [Double] -> Double
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum [Double -> Double
forall a. Floating a => a -> a
sqrt ([Double] -> [Double] -> Double
dist2 [Double]
center' [Double]
p) | (Int
_, [Double]
p) <- [(Int, [Double])]
points]
                        pivot :: [Double]
pivot = (Int, [Double]) -> [Double]
forall a b. (a, b) -> b
snd ([(Int, [Double])]
points [(Int, [Double])] -> Int -> (Int, [Double])
forall a. HasCallStack => [a] -> Int -> a
!! [Double] -> Int
forall a. Ord a => [a] -> Int
argmax [[Double] -> [Double] -> Double
dist2 [Double]
center' [Double]
p | (Int
_, [Double]
p) <- [(Int, [Double])]
points])
                        ordered :: [(Int, [Double])]
ordered = ((Int, [Double]) -> Double)
-> [(Int, [Double])] -> [(Int, [Double])]
forall b a. Ord b => (a -> b) -> [a] -> [a]
sortOn (\(Int
_, [Double]
p) -> [Double] -> [Double] -> Double
dist2 [Double]
pivot [Double]
p) [(Int, [Double])]
points
                        middle :: Int
middle = [(Int, [Double])] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [(Int, [Double])]
points Int -> Int -> Int
forall a. Integral a => a -> a -> a
`div` Int
2
                     in [Double] -> Double -> Tree -> Tree -> Tree
BallNode [Double]
center' Double
radius ([(Int, [Double])] -> Tree
build (Int -> [(Int, [Double])] -> [(Int, [Double])]
forall a. Int -> [a] -> [a]
take Int
middle [(Int, [Double])]
ordered)) ([(Int, [Double])] -> Tree
build (Int -> [(Int, [Double])] -> [(Int, [Double])]
forall a. Int -> [a] -> [a]
drop Int
middle [(Int, [Double])]
ordered))
    let tree = [(Int, [Double])] -> Tree
build ([Int] -> Dataset -> [(Int, [Double])]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
0 ..] Dataset
xs)
        validate (Leaf [(Int, [Double])]
_) = () -> Either MathError ()
forall a b. b -> Either a b
Right ()
        validate (KDNode Int
_ Double
cutoff Tree
left Tree
right) = [Double] -> Either MathError ()
checked [Double
cutoff] Either MathError () -> Either MathError () -> Either MathError ()
forall a b.
Either MathError a -> Either MathError b -> Either MathError b
forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> Tree -> Either MathError ()
validate Tree
left Either MathError () -> Either MathError () -> Either MathError ()
forall a b.
Either MathError a -> Either MathError b -> Either MathError b
forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> Tree -> Either MathError ()
validate Tree
right
        validate (BallNode [Double]
center' Double
radius Tree
left Tree
right) = [Double] -> Either MathError ()
checked (Double
radius Double -> [Double] -> [Double]
forall a. a -> [a] -> [a]
: [Double]
center') Either MathError () -> Either MathError () -> Either MathError ()
forall a b.
Either MathError a -> Either MathError b -> Either MathError b
forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> Tree -> Either MathError ()
validate Tree
left Either MathError () -> Either MathError () -> Either MathError ()
forall a b.
Either MathError a -> Either MathError b -> Either MathError b
forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> Tree -> Either MathError ()
validate Tree
right
    validate tree
    Right (SpatialIndex d n tree)
  where
    transposeLocal :: [[a]] -> [[a]]
transposeLocal [] = []
    transposeLocal [[a]]
rows' = let d :: Int
d = [a] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length ([[a]] -> [a]
forall a. [[a]] -> [a]
firstRow [[a]]
rows') in [[[a]
r [a] -> Int -> a
forall a. HasCallStack => [a] -> Int -> a
!! Int
j | [a]
r <- [[a]]
rows'] | Int
j <- [Int
0 .. Int
d Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]]

{- | Neighbor indices/distances plus visited-leaf and exhaustive-search diagnostics.
Budget-limited candidates are explicitly marked approximate.
-}
data NeighborResult = NeighborResult
    { NeighborResult -> [(Int, Double)]
neighbors :: [(Int, Double)]
    -- ^ Training-row indices paired with returned Euclidean distances, ordered by distance and then index.
    , NeighborResult -> Int
visitedLeaves :: Int
    -- ^ Number of leaves visited against the explicit query budget.
    , NeighborResult -> Bool
exhaustive :: Bool
    -- ^ True if the neighbor search completed its exact stopping condition; False if unexplored nodes may improve candidates.
    }
    deriving (NeighborResult -> NeighborResult -> Bool
(NeighborResult -> NeighborResult -> Bool)
-> (NeighborResult -> NeighborResult -> Bool) -> Eq NeighborResult
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: NeighborResult -> NeighborResult -> Bool
== :: NeighborResult -> NeighborResult -> Bool
$c/= :: NeighborResult -> NeighborResult -> Bool
/= :: NeighborResult -> NeighborResult -> Bool
Eq, Int -> NeighborResult -> ShowS
[NeighborResult] -> ShowS
NeighborResult -> String
(Int -> NeighborResult -> ShowS)
-> (NeighborResult -> String)
-> ([NeighborResult] -> ShowS)
-> Show NeighborResult
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> NeighborResult -> ShowS
showsPrec :: Int -> NeighborResult -> ShowS
$cshow :: NeighborResult -> String
show :: NeighborResult -> String
$cshowList :: [NeighborResult] -> ShowS
showList :: [NeighborResult] -> ShowS
Show)

{- | maxLeaves is a query work budget. A False exhaustive flag means unexplored
nodes might improve the result. Even approximate results have exact distances.
-}
nearestNeighbors :: Int -> Int -> SpatialIndex -> [Double] -> Either MathError NeighborResult
nearestNeighbors :: Int
-> Int
-> SpatialIndex
-> [Double]
-> Either MathError NeighborResult
nearestNeighbors Int
k Int
maxLeaves (SpatialIndex Int
dimension Int
count Tree
tree) [Double]
point = do
    Int -> Dataset -> Either MathError ()
queries Int
dimension [[Double]
point]
    Bool -> String -> Either MathError ()
require (Int
k Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
0 Bool -> Bool -> Bool
&& Int
k Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
count Bool -> Bool -> Bool
&& Int
maxLeaves Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
0) String
"invalid neighbor count/leaf budget"
    let trim :: [(Int, Double)] -> [(Int, Double)]
trim = Int -> [(Int, Double)] -> [(Int, Double)]
forall a. Int -> [a] -> [a]
take Int
k ([(Int, Double)] -> [(Int, Double)])
-> ([(Int, Double)] -> [(Int, Double)])
-> [(Int, Double)]
-> [(Int, Double)]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. ((Int, Double) -> (Double, Int))
-> [(Int, Double)] -> [(Int, Double)]
forall b a. Ord b => (a -> b) -> [a] -> [a]
sortOn (\(Int
i, Double
d) -> (Double
d, Int
i))
        bound :: [(a, a)] -> a
bound [(a, a)]
best = if [(a, a)] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [(a, a)]
best Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
k then a
1 a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
0 else (a, a) -> a
forall a b. (a, b) -> b
snd ([(a, a)] -> (a, a)
forall a. HasCallStack => [a] -> a
last [(a, a)]
best)
        search :: Tree -> [(Int, Double)] -> Int -> ([(Int, Double)], Int, Bool)
search Tree
node [(Int, Double)]
best Int
used
            | Int
used Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
maxLeaves = ([(Int, Double)]
best, Int
used, Bool
False)
            | Bool
otherwise = case Tree
node of
                Leaf [(Int, [Double])]
pts -> ([(Int, Double)] -> [(Int, Double)]
trim ([(Int, Double)]
best [(Int, Double)] -> [(Int, Double)] -> [(Int, Double)]
forall a. [a] -> [a] -> [a]
++ [(Int
i, Double -> Double
forall a. Floating a => a -> a
sqrt ([Double] -> [Double] -> Double
dist2 [Double]
point [Double]
p)) | (Int
i, [Double]
p) <- [(Int, [Double])]
pts]), Int
used Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1, Bool
True)
                BallNode [Double]
center' Double
radius Tree
left Tree
right
                    | Double -> Double
forall a. Floating a => a -> a
sqrt ([Double] -> [Double] -> Double
dist2 [Double]
point [Double]
center') Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
radius Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> [(Int, Double)] -> Double
forall {a} {a}. Fractional a => [(a, a)] -> a
bound [(Int, Double)]
best -> ([(Int, Double)]
best, Int
used, Bool
True)
                    | Bool
otherwise ->
                        let ([(Int, Double)]
a, Int
u, Bool
e1) = Tree -> [(Int, Double)] -> Int -> ([(Int, Double)], Int, Bool)
search Tree
left [(Int, Double)]
best Int
used; ([(Int, Double)]
b, Int
v, Bool
e2) = Tree -> [(Int, Double)] -> Int -> ([(Int, Double)], Int, Bool)
search Tree
right [(Int, Double)]
a Int
u in ([(Int, Double)]
b, Int
v, Bool
e1 Bool -> Bool -> Bool
&& Bool
e2)
                KDNode Int
axis Double
cutoff Tree
left Tree
right ->
                    let (Tree
near, Tree
far) = if [Double]
point [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
axis Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
cutoff then (Tree
left, Tree
right) else (Tree
right, Tree
left)
                        ([(Int, Double)]
a, Int
u, Bool
e1) = Tree -> [(Int, Double)] -> Int -> ([(Int, Double)], Int, Bool)
search Tree
near [(Int, Double)]
best Int
used
                     in if Double -> Double
forall a. Num a => a -> a
abs ([Double]
point [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
axis Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
cutoff) Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> [(Int, Double)] -> Double
forall {a} {a}. Fractional a => [(a, a)] -> a
bound [(Int, Double)]
a
                            then ([(Int, Double)]
a, Int
u, Bool
e1)
                            else
                                let ([(Int, Double)]
b, Int
v, Bool
e2) = Tree -> [(Int, Double)] -> Int -> ([(Int, Double)], Int, Bool)
search Tree
far [(Int, Double)]
a Int
u in ([(Int, Double)]
b, Int
v, Bool
e1 Bool -> Bool -> Bool
&& Bool
e2)
        ([(Int, Double)]
answer, Int
leafCount, Bool
complete) = Tree -> [(Int, Double)] -> Int -> ([(Int, Double)], Int, Bool)
search Tree
tree [] Int
0
    [Double] -> Either MathError ()
checked (((Int, Double) -> Double) -> [(Int, Double)] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (Int, Double) -> Double
forall a b. (a, b) -> b
snd [(Int, Double)]
answer)
    NeighborResult -> Either MathError NeighborResult
forall a b. b -> Either a b
Right ([(Int, Double)] -> Int -> Bool -> NeighborResult
NeighborResult [(Int, Double)]
answer Int
leafCount Bool
complete)

-- | Average the k exact nearest training responses for each query row.
knnRegress :: Int -> SpatialIndex -> [Double] -> Dataset -> Either MathError [Double]
knnRegress :: Int
-> SpatialIndex -> [Double] -> Dataset -> Either MathError [Double]
knnRegress Int
k index :: SpatialIndex
index@(SpatialIndex Int
_ Int
n Tree
_) [Double]
targets Dataset
xs = do
    Bool -> String -> Either MathError ()
require ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
targets Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
n) String
"KNN response count mismatch"
    [Double] -> Either MathError ()
checked [Double]
targets
    ([Double] -> Either MathError Double)
-> Dataset -> Either MathError [Double]
forall (t :: * -> *) (f :: * -> *) a b.
(Traversable t, Applicative f) =>
(a -> f b) -> t a -> f (t b)
forall (f :: * -> *) a b.
Applicative f =>
(a -> f b) -> [a] -> f [b]
traverse
        ( \[Double]
x -> do
            r <- Int
-> Int
-> SpatialIndex
-> [Double]
-> Either MathError NeighborResult
nearestNeighbors Int
k Int
forall a. Bounded a => a
maxBound SpatialIndex
index [Double]
x
            Right (sum [targets !! i | (i, _) <- neighbors r] / fromIntegral k)
        )
        Dataset
xs

-- | Majority vote among k exact neighbors; targets follow original training row order.
knnClassify :: Int -> SpatialIndex -> [Int] -> Dataset -> Either MathError [Int]
knnClassify :: Int -> SpatialIndex -> [Int] -> Dataset -> Either MathError [Int]
knnClassify Int
k index :: SpatialIndex
index@(SpatialIndex Int
_ Int
n Tree
_) [Int]
targets Dataset
xs = do
    Bool -> String -> Either MathError ()
require ([Int] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Int]
targets Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
n) String
"KNN label count mismatch"
    ([Double] -> Either MathError Int)
-> Dataset -> Either MathError [Int]
forall (t :: * -> *) (f :: * -> *) a b.
(Traversable t, Applicative f) =>
(a -> f b) -> t a -> f (t b)
forall (f :: * -> *) a b.
Applicative f =>
(a -> f b) -> [a] -> f [b]
traverse
        ( \[Double]
x -> do
            r <- Int
-> Int
-> SpatialIndex
-> [Double]
-> Either MathError NeighborResult
nearestNeighbors Int
k Int
forall a. Bounded a => a
maxBound SpatialIndex
index [Double]
x
            let labels = [[Int]
targets [Int] -> Int -> Int
forall a. HasCallStack => [a] -> Int -> a
!! Int
i | (Int
i, Double
_) <- NeighborResult -> [(Int, Double)]
neighbors NeighborResult
r]
                classes = [Int] -> [Int]
forall a. Ord a => [a] -> [a]
sort ([Int] -> [Int]
forall a. Eq a => [a] -> [a]
nub [Int]
labels)
                counts = [[Int] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length ((Int -> Bool) -> [Int] -> [Int]
forall a. (a -> Bool) -> [a] -> [a]
filter (Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
c) [Int]
labels) | Int
c <- [Int]
classes]
            Right (classes !! argmax counts)
        )
        Dataset
xs

-- | Seeded p-stable Euclidean hash tables retaining Gaussian projections and offset buckets.
data LSH = LSH Dataset Double [[([Double], Double)]] [[([Integer], Int)]] deriving (LSH -> LSH -> Bool
(LSH -> LSH -> Bool) -> (LSH -> LSH -> Bool) -> Eq LSH
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: LSH -> LSH -> Bool
== :: LSH -> LSH -> Bool
$c/= :: LSH -> LSH -> Bool
/= :: LSH -> LSH -> Bool
Eq, Int -> LSH -> ShowS
[LSH] -> ShowS
LSH -> String
(Int -> LSH -> ShowS)
-> (LSH -> String) -> ([LSH] -> ShowS) -> Show LSH
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> LSH -> ShowS
showsPrec :: Int -> LSH -> ShowS
$cshow :: LSH -> String
show :: LSH -> String
$cshowList :: [LSH] -> ShowS
showList :: [LSH] -> ShowS
Show)

-- | Build tables with projections scalar hashes each and positive bucket width. Storage budget is explicit.
buildLSH :: Word64 -> Int -> Int -> Double -> Dataset -> Either MathError LSH
buildLSH :: Word64 -> Int -> Int -> Double -> Dataset -> Either MathError LSH
buildLSH Word64
seed Int
tables Int
projections Double
width Dataset
xs = do
    (_, d) <- Dataset -> Either MathError (Int, Int)
dataset Dataset
xs
    require
        ( tables > 0
            && tables <= 10000
            && projections > 0
            && projections <= 10000
            && finite width
            && width > 0
        )
        "invalid LSH dimensions/width"
    require
        (toInteger tables * toInteger projections * toInteger d <= 10000000)
        "LSH projection budget exceeded"
    let generate RNG
r p
_ = do
            (direction, r') <- Int
-> (RNG -> Either MathError (Double, RNG))
-> RNG
-> Either MathError ([Double], RNG)
forall a.
Int
-> (RNG -> Either MathError (a, RNG))
-> RNG
-> Either MathError ([a], RNG)
R.sampleN Int
d (Double -> Double -> RNG -> Either MathError (Double, RNG)
R.normal Double
0 Double
1) RNG
r
            let (u, r'') = R.uniform r'
            Right (r'', (direction, u * width))
    (_, flat) <-
        foldM
            (\(RNG
r, [([Double], Double)]
acc) Int
i -> do (next, value) <- RNG -> Int -> Either MathError (RNG, ([Double], Double))
forall {p}. RNG -> p -> Either MathError (RNG, ([Double], Double))
generate RNG
r Int
i; Right (next, acc ++ [value]))
            (R.seedRNG R.PCG32 seed, [])
            [1 .. tables * projections]
    let planes = [Int -> [([Double], Double)] -> [([Double], Double)]
forall a. Int -> [a] -> [a]
take Int
projections (Int -> [([Double], Double)] -> [([Double], Double)]
forall a. Int -> [a] -> [a]
drop (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
projections) [([Double], Double)]
flat) | Int
i <- [Int
0 .. Int
tables Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]]
    checked [(dot direction row + offset) / width | ps <- planes, (direction, offset) <- ps, row <- xs]
    let signature [([Double], Double)]
ps [Double]
row = [Double -> a
forall b. Integral b => Double -> b
forall a b. (RealFrac a, Integral b) => a -> b
floor (([Double] -> [Double] -> Double
dot [Double]
direction [Double]
row Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
offset) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
width) | ([Double]
direction, Double
offset) <- [([Double], Double)]
ps]
        buckets = [[([([Double], Double)] -> [Double] -> [Integer]
forall {a}. Integral a => [([Double], Double)] -> [Double] -> [a]
signature [([Double], Double)]
ps [Double]
row, Int
i) | (Int
i, [Double]
row) <- [Int] -> Dataset -> [(Int, [Double])]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
0 ..] Dataset
xs] | [([Double], Double)]
ps <- [[([Double], Double)]]
planes]
    Right (LSH xs width planes buckets)

-- | Rank matching-bucket candidates by exact distance; may return fewer than k results, including none.
queryLSH :: Int -> LSH -> [Double] -> Either MathError [(Int, Double)]
queryLSH :: Int -> LSH -> [Double] -> Either MathError [(Int, Double)]
queryLSH Int
k (LSH Dataset
xs Double
width [[([Double], Double)]]
planes [[([Integer], Int)]]
buckets) [Double]
point = do
    Int -> Dataset -> Either MathError ()
queries ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length (Dataset -> [Double]
forall a. [[a]] -> [a]
firstRow Dataset
xs)) [[Double]
point]
    Bool -> String -> Either MathError ()
require (Int
k Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
0 Bool -> Bool -> Bool
&& Int
k Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Dataset -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length Dataset
xs) String
"invalid LSH requested neighbor count"
    [Double] -> Either MathError ()
checked [([Double] -> [Double] -> Double
dot [Double]
direction [Double]
point Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
offset) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
width | [([Double], Double)]
ps <- [[([Double], Double)]]
planes, ([Double]
direction, Double
offset) <- [([Double], Double)]
ps]
    let signature :: [([Double], Double)] -> [a]
signature [([Double], Double)]
ps = [Double -> a
forall b. Integral b => Double -> b
forall a b. (RealFrac a, Integral b) => a -> b
floor (([Double] -> [Double] -> Double
dot [Double]
direction [Double]
point Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
offset) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
width) | ([Double]
direction, Double
offset) <- [([Double], Double)]
ps]
        candidates :: [Int]
candidates = [Int] -> [Int]
forall a. Eq a => [a] -> [a]
nub [Int
i | ([([Double], Double)]
ps, [([Integer], Int)]
bucket) <- [[([Double], Double)]]
-> [[([Integer], Int)]]
-> [([([Double], Double)], [([Integer], Int)])]
forall a b. [a] -> [b] -> [(a, b)]
zip [[([Double], Double)]]
planes [[([Integer], Int)]]
buckets, ([Integer]
key, Int
i) <- [([Integer], Int)]
bucket, [Integer]
key [Integer] -> [Integer] -> Bool
forall a. Eq a => a -> a -> Bool
== [([Double], Double)] -> [Integer]
forall {a}. Integral a => [([Double], Double)] -> [a]
signature [([Double], Double)]
ps]
        answer :: [(Int, Double)]
answer = Int -> [(Int, Double)] -> [(Int, Double)]
forall a. Int -> [a] -> [a]
take Int
k (((Int, Double) -> (Double, Int))
-> [(Int, Double)] -> [(Int, Double)]
forall b a. Ord b => (a -> b) -> [a] -> [a]
sortOn (\(Int
i, Double
d) -> (Double
d, Int
i)) [(Int
i, Double -> Double
forall a. Floating a => a -> a
sqrt ([Double] -> [Double] -> Double
dist2 [Double]
point (Dataset
xs Dataset -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! Int
i))) | Int
i <- [Int]
candidates])
    [Double] -> Either MathError ()
checked (((Int, Double) -> Double) -> [(Int, Double)] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (Int, Double) -> Double
forall a b. (a, b) -> b
snd [(Int, Double)]
answer)
    [(Int, Double)] -> Either MathError [(Int, Double)]
forall a b. b -> Either a b
Right [(Int, Double)]
answer

{- | Training-set LOF with ties included at the kth neighbor distance. Identical
neighborhoods have LOF 1, including the infinite-density 0/0 limiting case.
-}
localOutlierFactor :: Int -> Dataset -> Either MathError [Double]
localOutlierFactor :: Int -> Dataset -> Either MathError [Double]
localOutlierFactor Int
k Dataset
xs = do
    (n, _) <- Dataset -> Either MathError (Int, Int)
dataset Dataset
xs
    require (k > 0 && k < n) "LOF k must exclude the point itself"
    let distances = [[Double -> Double
forall a. Floating a => a -> a
sqrt ([Double] -> [Double] -> Double
dist2 [Double]
x [Double]
y) | [Double]
y <- Dataset
xs] | [Double]
x <- Dataset
xs]
        kth = [[Double] -> [Double]
forall a. Ord a => [a] -> [a]
sort [Dataset
distances Dataset -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! Int
i [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
j | Int
j <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1], Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
j] [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! (Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) | Int
i <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]]
        near Int
i = [Int
j | Int
j <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1], Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
j, Dataset
distances Dataset -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! Int
i [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
j Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= [Double]
kth [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i]
        reach Int
i = [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Double -> Double -> Double
forall a. Ord a => a -> a -> a
max ([Double]
kth [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
j) (Dataset
distances Dataset -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! Int
i [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
j) | Int
j <- Int -> [Int]
near Int
i] Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral ([Int] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length (Int -> [Int]
near Int
i))
        density = (Int -> Double) -> [Int] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map Int -> Double
reach [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]
        ratio a
a a
b | a
a a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0 Bool -> Bool -> Bool
&& a
b a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0 = a
1 | a
b a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0 = a
1 a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
0 | Bool
otherwise = a
a a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
b
        scores =
            [ [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Double -> Double -> Double
forall {a}. (Eq a, Fractional a) => a -> a -> a
ratio ([Double]
density [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i) ([Double]
density [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
j) | Int
j <- Int -> [Int]
near Int
i] Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral ([Int] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length (Int -> [Int]
near Int
i))
            | Int
i <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]
            ]
    checked scores
    Right scores

data IsolationTree = IsolationLeaf Int | IsolationNode Int Double IsolationTree IsolationTree
    deriving (IsolationTree -> IsolationTree -> Bool
(IsolationTree -> IsolationTree -> Bool)
-> (IsolationTree -> IsolationTree -> Bool) -> Eq IsolationTree
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: IsolationTree -> IsolationTree -> Bool
== :: IsolationTree -> IsolationTree -> Bool
$c/= :: IsolationTree -> IsolationTree -> Bool
/= :: IsolationTree -> IsolationTree -> Bool
Eq, Int -> IsolationTree -> ShowS
[IsolationTree] -> ShowS
IsolationTree -> String
(Int -> IsolationTree -> ShowS)
-> (IsolationTree -> String)
-> ([IsolationTree] -> ShowS)
-> Show IsolationTree
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> IsolationTree -> ShowS
showsPrec :: Int -> IsolationTree -> ShowS
$cshow :: IsolationTree -> String
show :: IsolationTree -> String
$cshowList :: [IsolationTree] -> ShowS
showList :: [IsolationTree] -> ShowS
Show)

-- | Owned random isolation trees with fixed training dimension and subsample size.
data IsolationForest = IsolationForest Int Int [IsolationTree] deriving (IsolationForest -> IsolationForest -> Bool
(IsolationForest -> IsolationForest -> Bool)
-> (IsolationForest -> IsolationForest -> Bool)
-> Eq IsolationForest
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: IsolationForest -> IsolationForest -> Bool
== :: IsolationForest -> IsolationForest -> Bool
$c/= :: IsolationForest -> IsolationForest -> Bool
/= :: IsolationForest -> IsolationForest -> Bool
Eq, Int -> IsolationForest -> ShowS
[IsolationForest] -> ShowS
IsolationForest -> String
(Int -> IsolationForest -> ShowS)
-> (IsolationForest -> String)
-> ([IsolationForest] -> ShowS)
-> Show IsolationForest
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> IsolationForest -> ShowS
showsPrec :: Int -> IsolationForest -> ShowS
$cshow :: IsolationForest -> String
show :: IsolationForest -> String
$cshowList :: [IsolationForest] -> ShowS
showList :: [IsolationForest] -> ShowS
Show)

-- | Seeded independent isolation trees using subsamples of 2..2048 observations.
fitIsolationForest :: Word64 -> Int -> Int -> Dataset -> Either MathError IsolationForest
fitIsolationForest :: Word64 -> Int -> Int -> Dataset -> Either MathError IsolationForest
fitIsolationForest Word64
seed Int
treeCount Int
sampleCount Dataset
xs = do
    (n, d) <- Dataset -> Either MathError (Int, Int)
dataset Dataset
xs
    require
        (treeCount > 0 && treeCount <= 10000 && sampleCount >= 2 && sampleCount <= min 2048 n)
        "invalid isolation forest tree/subsample budget"
    let maxDepth = Double -> Int
forall b. Integral b => Double -> b
forall a b. (RealFrac a, Integral b) => a -> b
ceiling (Double -> Double -> Double
forall a. Floating a => a -> a -> a
logBase Double
2 (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
sampleCount) :: Double)
        grow Int
depth Dataset
rows' RNG
rng
            | Dataset -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length Dataset
rows' Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
1 Bool -> Bool -> Bool
|| Int
depth Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
maxDepth = (IsolationTree, RNG) -> Either MathError (IsolationTree, RNG)
forall a b. b -> Either a b
Right (Int -> IsolationTree
IsolationLeaf (Dataset -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length Dataset
rows'), RNG
rng)
            | Bool
otherwise = do
                let dimensions :: [Int]
dimensions = [Int
j | Int
j <- [Int
0 .. Int
d Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1], [Double] -> Double
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum (([Double] -> Double) -> Dataset -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map ([Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
j) Dataset
rows') Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> [Double] -> Double
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
minimum (([Double] -> Double) -> Dataset -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map ([Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
j) Dataset
rows')]
                if [Int] -> Bool
forall a. [a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [Int]
dimensions
                    then (IsolationTree, RNG) -> Either MathError (IsolationTree, RNG)
forall a b. b -> Either a b
Right (Int -> IsolationTree
IsolationLeaf (Dataset -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length Dataset
rows'), RNG
rng)
                    else do
                        (choice, next) <- Int -> Int -> RNG -> Either MathError (Int, RNG)
R.uniformInt Int
0 ([Int] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Int]
dimensions Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) RNG
rng
                        let axis = [Int]
dimensions [Int] -> Int -> Int
forall a. HasCallStack => [a] -> Int -> a
!! Int
choice
                            lo = [Double] -> Double
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
minimum (([Double] -> Double) -> Dataset -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map ([Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
axis) Dataset
rows')
                            hi = [Double] -> Double
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum (([Double] -> Double) -> Dataset -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map ([Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
axis) Dataset
rows')
                            (u, next') = R.uniformOpen next
                            cutoff = Double
lo Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (Double
hi Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
lo) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
u
                        require (finite cutoff && cutoff > lo && cutoff < hi) "isolation split not representable"
                        (left, leftR) <- grow (depth + 1) (filter (\[Double]
r -> [Double]
r [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
axis Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
cutoff) rows') next'
                        (right, rightR) <- grow (depth + 1) (filter (\[Double]
r -> [Double]
r [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
axis Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
cutoff) rows') leftR
                        Right (IsolationNode axis cutoff left right, rightR)
        one (RNG
rng, [IsolationTree]
trees) p
_ = do
            (shuffled, next) <- Dataset -> RNG -> Either MathError (Dataset, RNG)
forall a. [a] -> RNG -> Either MathError ([a], RNG)
R.shuffle Dataset
xs RNG
rng
            (tree, next') <- grow (0 :: Int) (take sampleCount shuffled) next
            Right (next', tree : trees)
    (_, trees) <- foldM one (R.seedRNG R.PCG32 seed, []) [1 .. treeCount]
    Right (IsolationForest d sampleCount trees)

-- | Return 2^(-mean corrected path length / expected random-tree path length); larger means more anomalous.
scoreIsolationForest :: IsolationForest -> Dataset -> Either MathError [Double]
scoreIsolationForest :: IsolationForest -> Dataset -> Either MathError [Double]
scoreIsolationForest (IsolationForest Int
d Int
sampleCount [IsolationTree]
trees) Dataset
xs = do
    Int -> Dataset -> Either MathError ()
queries Int
d Dataset
xs
    let correction :: a -> a
correction a
n
            | a
n a -> a -> Bool
forall a. Ord a => a -> a -> Bool
<= a
1 = a
0
            | Bool
otherwise =
                a
2 a -> a -> a
forall a. Num a => a -> a -> a
* [a] -> a
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [a
1 a -> a -> a
forall a. Fractional a => a -> a -> a
/ a -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral a
j | a
j <- [a
1 .. a
n a -> a -> a
forall a. Num a => a -> a -> a
- a
1]] a -> a -> a
forall a. Num a => a -> a -> a
- a
2 a -> a -> a
forall a. Num a => a -> a -> a
* a -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral (a
n a -> a -> a
forall a. Num a => a -> a -> a
- a
1) a -> a -> a
forall a. Fractional a => a -> a -> a
/ a -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral a
n
        path :: t -> IsolationTree -> [Double] -> t
path t
depth (IsolationLeaf Int
size) [Double]
_ = t -> t
forall a b. (Integral a, Num b) => a -> b
fromIntegral t
depth t -> t -> t
forall a. Num a => a -> a -> a
+ Int -> t
forall {a} {a}. (Fractional a, Integral a) => a -> a
correction Int
size
        path t
depth (IsolationNode Int
axis Double
cutoff IsolationTree
left IsolationTree
right) [Double]
point = t -> IsolationTree -> [Double] -> t
path (t
depth t -> t -> t
forall a. Num a => a -> a -> a
+ t
1) (if [Double]
point [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
axis Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
cutoff then IsolationTree
left else IsolationTree
right) [Double]
point
        answer :: [Double]
answer =
            [ Double
2
                Double -> Double -> Double
forall a. Floating a => a -> a -> a
** ( Double -> Double
forall a. Num a => a -> a
negate ([Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Int -> IsolationTree -> [Double] -> Double
forall {t} {t}.
(Integral t, Fractional t) =>
t -> IsolationTree -> [Double] -> t
path (Int
0 :: Int) IsolationTree
t [Double]
x | IsolationTree
t <- [IsolationTree]
trees] Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral ([IsolationTree] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [IsolationTree]
trees))
                        Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall {a} {a}. (Fractional a, Integral a) => a -> a
correction Int
sampleCount
                   )
            | [Double]
x <- Dataset
xs
            ]
    [Double] -> Either MathError ()
checked [Double]
answer
    [Double] -> Either MathError [Double]
forall a b. b -> Either a b
Right [Double]
answer

-- | Kernel novelty detector with normalized nonnegative dual weights and learned offset.
data OneClassSVM = OneClassSVM Kernel Dataset [Double] Double deriving (OneClassSVM -> OneClassSVM -> Bool
(OneClassSVM -> OneClassSVM -> Bool)
-> (OneClassSVM -> OneClassSVM -> Bool) -> Eq OneClassSVM
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: OneClassSVM -> OneClassSVM -> Bool
== :: OneClassSVM -> OneClassSVM -> Bool
$c/= :: OneClassSVM -> OneClassSVM -> Bool
/= :: OneClassSVM -> OneClassSVM -> Bool
Eq, Int -> OneClassSVM -> ShowS
[OneClassSVM] -> ShowS
OneClassSVM -> String
(Int -> OneClassSVM -> ShowS)
-> (OneClassSVM -> String)
-> ([OneClassSVM] -> ShowS)
-> Show OneClassSVM
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> OneClassSVM -> ShowS
showsPrec :: Int -> OneClassSVM -> ShowS
$cshow :: OneClassSVM -> String
show :: OneClassSVM -> String
$cshowList :: [OneClassSVM] -> ShowS
showList :: [OneClassSVM] -> ShowS
Show)

{- | nu-SVM dual coordinate descent: sum alpha=1, 0<=alpha<=1/(nu*n).
The reported gap is the pairwise box-constrained KKT violation.
-}
fitOneClassSVM :: MLConfig -> Kernel -> Double -> Dataset -> Either MathError (Fit OneClassSVM)
fitOneClassSVM :: MLConfig
-> Kernel
-> Double
-> Dataset
-> Either MathError (Fit OneClassSVM)
fitOneClassSVM MLConfig
c Kernel
kind Double
nu Dataset
xs = do
    MLConfig -> Either MathError ()
config MLConfig
c
    (n, _) <- Dataset -> Either MathError (Int, Int)
dataset Dataset
xs
    require (nu > 0 && nu <= 1 && finite nu) "one-class SVM nu must lie in (0,1]"
    gram <- traverse (\[Double]
x -> ([Double] -> Either MathError Double)
-> Dataset -> Either MathError [Double]
forall (t :: * -> *) (f :: * -> *) a b.
(Traversable t, Applicative f) =>
(a -> f b) -> t a -> f (t b)
forall (f :: * -> *) a b.
Applicative f =>
(a -> f b) -> [a] -> f [b]
traverse (Kernel -> [Double] -> [Double] -> Either MathError Double
kernel Kernel
kind [Double]
x) Dataset
xs) xs
    let cap = Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
nu Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n)
        initial = Int -> Double -> [Double]
forall a. Int -> a -> [a]
replicate Int
n (Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n)
        go Int
iteration [Double]
alpha = do
            let gradient :: [Double]
gradient = ([Double] -> Double) -> Dataset -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map ([Double] -> [Double] -> Double
`dot` [Double]
alpha) Dataset
gram
                canIncrease :: [Int]
canIncrease = [Int
i | Int
i <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1], [Double]
alpha [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
cap Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1e-14]
                canDecrease :: [Int]
canDecrease = [Int
i | Int
i <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1], [Double]
alpha [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
1e-14]
                free :: [Double]
free = [[Double]
gradient [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i | Int
i <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1], [Double]
alpha [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
1e-10 Bool -> Bool -> Bool
&& [Double]
alpha [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
cap Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1e-10]
                rho :: Double
rho =
                    if Bool -> Bool
not ([Double] -> Bool
forall a. [a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [Double]
free)
                        then [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Double]
free Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
free)
                        else
                            let lower :: Double
lower = if [Int] -> Bool
forall a. [a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [Int]
canDecrease then [Double] -> Double
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
minimum [Double]
gradient else [Double] -> Double
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum [[Double]
gradient [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i | Int
i <- [Int]
canDecrease]
                                upper :: Double
upper = if [Int] -> Bool
forall a. [a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [Int]
canIncrease then [Double] -> Double
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum [Double]
gradient else [Double] -> Double
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
minimum [[Double]
gradient [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i | Int
i <- [Int]
canIncrease]
                             in (Double
lower Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
upper) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2
                result :: Double -> Bool -> Either a (Fit OneClassSVM)
result Double
residual Bool
conv = Fit OneClassSVM -> Either a (Fit OneClassSVM)
forall a b. b -> Either a b
Right (OneClassSVM -> Int -> Double -> Bool -> Fit OneClassSVM
forall a. a -> Int -> Double -> Bool -> Fit a
Fit (Kernel -> Dataset -> [Double] -> Double -> OneClassSVM
OneClassSVM Kernel
kind Dataset
xs [Double]
alpha Double
rho) Int
iteration Double
residual Bool
conv)
            if [Int] -> Bool
forall a. [a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [Int]
canIncrease Bool -> Bool -> Bool
|| [Int] -> Bool
forall a. [a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [Int]
canDecrease
                then Double -> Bool -> Either MathError (Fit OneClassSVM)
forall {a}. Double -> Bool -> Either a (Fit OneClassSVM)
result Double
0 Bool
True
                else do
                    let i :: Int
i = [Int]
canIncrease [Int] -> Int -> Int
forall a. HasCallStack => [a] -> Int -> a
!! [Double] -> Int
forall a. Ord a => [a] -> Int
argmin [[Double]
gradient [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
q | Int
q <- [Int]
canIncrease]
                        j :: Int
j = [Int]
canDecrease [Int] -> Int -> Int
forall a. HasCallStack => [a] -> Int -> a
!! [Double] -> Int
forall a. Ord a => [a] -> Int
argmax [[Double]
gradient [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
q | Int
q <- [Int]
canDecrease]
                        gap :: Double
gap = Double -> Double -> Double
forall a. Ord a => a -> a -> a
max Double
0 ([Double]
gradient [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
j Double -> Double -> Double
forall a. Num a => a -> a -> a
- [Double]
gradient [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i)
                    if Double
gap Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= MLConfig -> Double
mlTolerance MLConfig
c Bool -> Bool -> Bool
|| Int
iteration Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= MLConfig -> Int
mlIterations MLConfig
c
                        then Double -> Bool -> Either MathError (Fit OneClassSVM)
forall {a}. Double -> Bool -> Either a (Fit OneClassSVM)
result Double
gap (Double
gap Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= MLConfig -> Double
mlTolerance MLConfig
c)
                        else do
                            let curvature :: Double
curvature = Dataset
gram Dataset -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! Int
i [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Dataset
gram Dataset -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! Int
j [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
j Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Dataset
gram Dataset -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! Int
i [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
j
                                room :: Double
room = Double -> Double -> Double
forall a. Ord a => a -> a -> a
min (Double
cap Double -> Double -> Double
forall a. Num a => a -> a -> a
- [Double]
alpha [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i) ([Double]
alpha [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
j)
                                delta :: Double
delta = if Double
curvature Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
1e-14 then Double
room else Double -> Double -> Double
forall a. Ord a => a -> a -> a
min Double
room (Double
gap Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
curvature)
                                next :: [Double]
next = Int -> Double -> [Double] -> [Double]
forall a. Int -> a -> [a] -> [a]
replace Int
i ([Double]
alpha [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
delta) (Int -> Double -> [Double] -> [Double]
forall a. Int -> a -> [a] -> [a]
replace Int
j ([Double]
alpha [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
j Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
delta) [Double]
alpha)
                            [Double] -> Either MathError ()
checked [Double]
next
                            Int -> [Double] -> Either MathError (Fit OneClassSVM)
go (Int
iteration Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) [Double]
next
    go 0 initial

-- | Return signed novelty margins: positive indicates an inlier under the fitted kernel boundary.
decisionOneClassSVM :: OneClassSVM -> Dataset -> Either MathError [Double]
decisionOneClassSVM :: OneClassSVM -> Dataset -> Either MathError [Double]
decisionOneClassSVM (OneClassSVM Kernel
kind Dataset
train [Double]
alpha Double
rho) Dataset
xs = do
    Int -> Dataset -> Either MathError ()
queries ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length (Dataset -> [Double]
forall a. [[a]] -> [a]
firstRow Dataset
train)) Dataset
xs
    ([Double] -> Either MathError Double)
-> Dataset -> Either MathError [Double]
forall (t :: * -> *) (f :: * -> *) a b.
(Traversable t, Applicative f) =>
(a -> f b) -> t a -> f (t b)
forall (f :: * -> *) a b.
Applicative f =>
(a -> f b) -> [a] -> f [b]
traverse
        ( \[Double]
x -> do
            values <- ([Double] -> Either MathError Double)
-> Dataset -> Either MathError [Double]
forall (t :: * -> *) (f :: * -> *) a b.
(Traversable t, Applicative f) =>
(a -> f b) -> t a -> f (t b)
forall (f :: * -> *) a b.
Applicative f =>
(a -> f b) -> [a] -> f [b]
traverse (Kernel -> [Double] -> [Double] -> Either MathError Double
kernel Kernel
kind [Double]
x) Dataset
train
            let value = [Double] -> [Double] -> Double
dot [Double]
alpha [Double]
values Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
rho
            checked [value]
            Right value
        )
        Dataset
xs