-- SPDX-License-Identifier: MPL-2.0

{- | Dense, finite-budget clustering. Distances are Euclidean, labels zero
based; DBSCAN noise is -1. OPTICS preserves undefined reachability as Nothing.
Affinity propagation uses Frey & Dueck's responsibility/availability updates
(Science 2007, doi:10.1126/science.1136800), with explicit damping.
-}
module HMath.MachineLearning.Clustering (
    ClusterModel,
    clusterCenters,
    predictClusters,
    kMeans,
    kMeansPlusPlus,
    miniBatchKMeans,
    kMedoids,
    GaussianMixture,
    mixtureWeights,
    mixtureMeans,
    mixtureCovariances,
    fitGaussianMixture,
    mixtureResponsibilities,
    dbscan,
    optics,
    extractOPTICS,
    meanShift,
    Linkage (..),
    Merge (..),
    hierarchical,
    cutHierarchy,
    spectralClustering,
    affinityPropagation,
) where

import Control.Monad (foldM)
import Data.List (nub, sortOn, transpose)
import HMath.Error
import HMath.MachineLearning.Internal
import qualified HMath.Random as R

{- | Fitted cluster centroids as cluster-by-feature observations; prediction assigns new
queries to the learned centers.
-}
newtype ClusterModel = Clusters
    { ClusterModel -> Dataset
clusterCenters :: Dataset
    -- ^ Cluster centroids as cluster-by-feature rows.
    }
    deriving (ClusterModel -> ClusterModel -> Bool
(ClusterModel -> ClusterModel -> Bool)
-> (ClusterModel -> ClusterModel -> Bool) -> Eq ClusterModel
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: ClusterModel -> ClusterModel -> Bool
== :: ClusterModel -> ClusterModel -> Bool
$c/= :: ClusterModel -> ClusterModel -> Bool
/= :: ClusterModel -> ClusterModel -> Bool
Eq, Int -> ClusterModel -> ShowS
[ClusterModel] -> ShowS
ClusterModel -> String
(Int -> ClusterModel -> ShowS)
-> (ClusterModel -> String)
-> ([ClusterModel] -> ShowS)
-> Show ClusterModel
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> ClusterModel -> ShowS
showsPrec :: Int -> ClusterModel -> ShowS
$cshow :: ClusterModel -> String
show :: ClusterModel -> String
$cshowList :: [ClusterModel] -> ShowS
showList :: [ClusterModel] -> ShowS
Show)

-- | Assign each row to its closest squared-Euclidean center; ties use the first center.
predictClusters :: ClusterModel -> Dataset -> Either MathError [Int]
predictClusters :: ClusterModel -> Dataset -> Either MathError [Int]
predictClusters (Clusters Dataset
centers) 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
centers)) Dataset
xs
    let ds :: Dataset
ds = [[[Double] -> [Double] -> Double
dist2 [Double]
x [Double]
c | [Double]
c <- Dataset
centers] | [Double]
x <- Dataset
xs]
    [Double] -> Either MathError ()
checked (Dataset -> [Double]
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat Dataset
ds)
    [Int] -> Either MathError [Int]
forall a b. b -> Either a b
Right (([Double] -> Int) -> Dataset -> [Int]
forall a b. (a -> b) -> [a] -> [b]
map [Double] -> Int
forall a. Ord a => [a] -> Int
argmin Dataset
ds)

initialCenters :: Bool -> MLConfig -> Int -> Dataset -> Either MathError Dataset
initialCenters :: Bool -> MLConfig -> Int -> Dataset -> Either MathError Dataset
initialCenters Bool
plus MLConfig
c Int
k Dataset
xs = do
    (n, _) <- Dataset -> Either MathError (Int, Int)
dataset Dataset
xs
    require (k > 0 && k <= n) "cluster count must be between one and sample count"
    let rng = Generator -> Word64 -> RNG
R.seedRNG Generator
R.PCG32 (MLConfig -> Word64
mlSeed MLConfig
c)
    if not plus
        then do
            (shuffled, _) <- R.shuffle xs rng
            Right (take k shuffled)
        else do
            (first, rng') <- R.uniformInt 0 (n - 1) rng
            let go [Int]
chosen RNG
r
                    | [Int] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Int]
chosen Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
k = Dataset -> Either MathError Dataset
forall a b. b -> Either a b
Right ((Int -> [Double]) -> [Int] -> Dataset
forall a b. (a -> b) -> [a] -> [b]
map (Dataset
xs Dataset -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!!) [Int]
chosen)
                    | Bool
otherwise = do
                        let weights :: [Double]
weights =
                                [if Int
i Int -> [Int] -> Bool
forall a. Eq a => a -> [a] -> Bool
forall (t :: * -> *) a. (Foldable t, Eq a) => a -> t a -> Bool
`elem` [Int]
chosen then Double
0 else [Double] -> Double
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
minimum [[Double] -> [Double] -> Double
dist2 [Double]
x (Dataset
xs Dataset -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! Int
j) | Int
j <- [Int]
chosen] | (Int
i, [Double]
x) <- [Int] -> Dataset -> [(Int, [Double])]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
0 ..] Dataset
xs]
                            total :: Double
total = [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Double]
weights
                        [Double] -> Either MathError ()
checked [Double]
weights
                        let (Double
u, RNG
r') = RNG -> (Double, RNG)
R.uniform RNG
r
                            remaining :: [Int]
remaining = [Int
i | Int
i <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1], Int
i Int -> [Int] -> Bool
forall (t :: * -> *) a. (Foldable t, Eq a) => a -> t a -> Bool
`notElem` [Int]
chosen]
                            index :: Int
index =
                                if Double
total Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
0
                                    then [Int]
remaining [Int] -> Int -> Int
forall a. HasCallStack => [a] -> Int -> a
!! Int
0
                                    else Int -> Int -> Int
forall a. Ord a => a -> a -> a
min (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length ((Double -> Bool) -> [Double] -> [Double]
forall a. (a -> Bool) -> [a] -> [a]
takeWhile (Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
u Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
total) (Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
drop Int
1 ((Double -> Double -> Double) -> Double -> [Double] -> [Double]
forall b a. (b -> a -> b) -> b -> [a] -> [b]
scanl Double -> Double -> Double
forall a. Num a => a -> a -> a
(+) Double
0 [Double]
weights))))
                        [Int] -> RNG -> Either MathError Dataset
go ([Int]
chosen [Int] -> [Int] -> [Int]
forall a. [a] -> [a] -> [a]
++ [Int
index]) RNG
r'
            go [first] rng'

-- | Seeded random-center Lloyd iteration; residual is total center displacement.
kMeans :: MLConfig -> Int -> Dataset -> Either MathError (Fit ClusterModel)
kMeans :: MLConfig -> Int -> Dataset -> Either MathError (Fit ClusterModel)
kMeans = Bool
-> MLConfig
-> Int
-> Dataset
-> Either MathError (Fit ClusterModel)
lloyd Bool
False

-- | D-squared k-means++ initialization followed by Lloyd center updates.
kMeansPlusPlus :: MLConfig -> Int -> Dataset -> Either MathError (Fit ClusterModel)
kMeansPlusPlus :: MLConfig -> Int -> Dataset -> Either MathError (Fit ClusterModel)
kMeansPlusPlus = Bool
-> MLConfig
-> Int
-> Dataset
-> Either MathError (Fit ClusterModel)
lloyd Bool
True

lloyd :: Bool -> MLConfig -> Int -> Dataset -> Either MathError (Fit ClusterModel)
lloyd :: Bool
-> MLConfig
-> Int
-> Dataset
-> Either MathError (Fit ClusterModel)
lloyd Bool
plus MLConfig
c Int
k Dataset
xs = do
    MLConfig -> Either MathError ()
config MLConfig
c
    start <- Bool -> MLConfig -> Int -> Dataset -> Either MathError Dataset
initialCenters Bool
plus MLConfig
c Int
k Dataset
xs
    let step Dataset
centers = do
            labels <- ClusterModel -> Dataset -> Either MathError [Int]
predictClusters (Dataset -> ClusterModel
Clusters Dataset
centers) Dataset
xs
            let assigned Int
j = [[Double]
x | ([Double]
x, Int
l) <- Dataset -> [Int] -> [([Double], Int)]
forall a b. [a] -> [b] -> [(a, b)]
zip Dataset
xs [Int]
labels, Int
l Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
j]
                -- Keep an empty center at its prior location. This is stable on
                -- duplicate points and preserves the caller's requested k.
                next = [if Dataset -> Bool
forall a. [a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null (Int -> Dataset
assigned Int
j) then Dataset
centers Dataset -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! Int
j else Dataset -> [Double]
meanV (Int -> Dataset
assigned Int
j) | Int
j <- [Int
0 .. Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]]
            checked (concat next)
            Right next
    fit <- iterateFit c step (\Dataset
a Dataset
b -> [Double] -> Double
norm (Dataset -> [Double]
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat (([Double] -> [Double] -> [Double]) -> Dataset -> Dataset -> Dataset
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith [Double] -> [Double] -> [Double]
sub Dataset
a Dataset
b))) start
    Right fit{fitted = Clusters (fitted fit)}

{- | Online center updates with per-cluster counts and seeded shuffled batches.
One outer iteration consumes batchSize observations, cycling shuffled data.
-}
miniBatchKMeans :: MLConfig -> Int -> Int -> Dataset -> Either MathError (Fit ClusterModel)
miniBatchKMeans :: MLConfig
-> Int -> Int -> Dataset -> Either MathError (Fit ClusterModel)
miniBatchKMeans MLConfig
c Int
k Int
batchSize Dataset
xs = do
    MLConfig -> Either MathError ()
config MLConfig
c
    Bool -> String -> Either MathError ()
require (Int
batchSize Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
0 Bool -> Bool -> Bool
&& Int
batchSize 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 minibatch size"
    start <- Bool -> MLConfig -> Int -> Dataset -> Either MathError Dataset
initialCenters Bool
True MLConfig
c Int
k Dataset
xs
    let step (Dataset
centers, [a]
counts, RNG
rng) = do
            (batch, rng') <- Dataset -> RNG -> Either MathError (Dataset, RNG)
forall a. [a] -> RNG -> Either MathError ([a], RNG)
R.shuffle Dataset
xs RNG
rng
            (next, newCounts) <-
                foldM
                    ( \(Dataset
cs, [a]
ns) [Double]
x -> do
                        let j :: Int
j = [Double] -> Int
forall a. Ord a => [a] -> Int
argmin (([Double] -> Double) -> Dataset -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map ([Double] -> [Double] -> Double
dist2 [Double]
x) Dataset
cs)
                            count :: a
count = [a]
ns [a] -> Int -> a
forall a. HasCallStack => [a] -> Int -> a
!! Int
j a -> a -> a
forall a. Num a => a -> a -> a
+ a
1
                            value :: [Double]
value = [Double] -> [Double] -> [Double]
add (Dataset
cs Dataset -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! Int
j) (Double -> [Double] -> [Double]
scale (Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ a -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral a
count) ([Double] -> [Double] -> [Double]
sub [Double]
x (Dataset
cs Dataset -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! Int
j)))
                        [Double] -> Either MathError ()
checked [Double]
value
                        (Dataset, [a]) -> Either MathError (Dataset, [a])
forall a b. b -> Either a b
Right (Int -> [Double] -> Dataset -> Dataset
forall a. Int -> a -> [a] -> [a]
replace Int
j [Double]
value Dataset
cs, Int -> a -> [a] -> [a]
forall a. Int -> a -> [a] -> [a]
replace Int
j a
count [a]
ns)
                    )
                    (centers, counts)
                    (take batchSize batch)
            Right (next, newCounts, rng')
    fit <-
        iterateFit
            c
            step
            (\(Dataset
a, [Int]
_, RNG
_) (Dataset
b, [Int]
_, RNG
_) -> [Double] -> Double
norm (Dataset -> [Double]
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat (([Double] -> [Double] -> [Double]) -> Dataset -> Dataset -> Dataset
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith [Double] -> [Double] -> [Double]
sub Dataset
a Dataset
b)))
            (start, replicate k (0 :: Int), R.seedRNG R.PCG32 (mlSeed c))
    let (centers, _, _) = fitted fit
    Right fit{fitted = Clusters centers}

-- | PAM exhaustive medoid swaps. Return original row indices; residual is the objective decrease.
kMedoids :: MLConfig -> Int -> Dataset -> Either MathError (Fit [Int])
kMedoids :: MLConfig -> Int -> Dataset -> Either MathError (Fit [Int])
kMedoids MLConfig
c Int
k Dataset
xs = do
    MLConfig -> Either MathError ()
config MLConfig
c
    (n, _) <- Dataset -> Either MathError (Int, Int)
dataset Dataset
xs
    require (k > 0 && k <= n) "invalid medoid count"
    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]
        cost [Int]
ids = [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [[Double] -> Double
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
minimum [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]
ids] | Int
i <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]]
        step [Int]
ids = do
            let candidates :: [[Int]]
candidates = [Int]
ids [Int] -> [[Int]] -> [[Int]]
forall a. a -> [a] -> [a]
: [Int -> Int -> [Int] -> [Int]
forall a. Int -> a -> [a] -> [a]
replace Int
p Int
j [Int]
ids | Int
p <- [Int
0 .. Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1], Int
j <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1], Int
j Int -> [Int] -> Bool
forall (t :: * -> *) a. (Foldable t, Eq a) => a -> t a -> Bool
`notElem` [Int]
ids]
                costs :: [Double]
costs = ([Int] -> Double) -> [[Int]] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map [Int] -> Double
cost [[Int]]
candidates
            [Double] -> Either MathError ()
checked [Double]
costs
            [Int] -> Either MathError [Int]
forall a b. b -> Either a b
Right ([[Int]]
candidates [[Int]] -> Int -> [Int]
forall a. HasCallStack => [a] -> Int -> a
!! [Double] -> Int
forall a. Ord a => [a] -> Int
argmin [Double]
costs)
    iterateFit c step (\[Int]
a [Int]
b -> Double -> Double
forall a. Num a => a -> a
abs ([Int] -> Double
cost [Int]
a Double -> Double -> Double
forall a. Num a => a -> a -> a
- [Int] -> Double
cost [Int]
b)) [0 .. k - 1]

{- | Fitted component probabilities,mean vectors,and covariance matrices in shared component
order.
-}
data GaussianMixture = GMM
    { GaussianMixture -> [Double]
mixtureWeights :: [Double]
    -- ^ Gaussian-mixture component probabilities in component order.
    , GaussianMixture -> Dataset
mixtureMeans :: Dataset
    -- ^ Component mean vectors, matching mixtureWeights order.
    , GaussianMixture -> [Dataset]
mixtureCovariances :: [Dataset]
    -- ^ One feature-by-feature covariance matrix per Gaussian component.
    }
    deriving (GaussianMixture -> GaussianMixture -> Bool
(GaussianMixture -> GaussianMixture -> Bool)
-> (GaussianMixture -> GaussianMixture -> Bool)
-> Eq GaussianMixture
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: GaussianMixture -> GaussianMixture -> Bool
== :: GaussianMixture -> GaussianMixture -> Bool
$c/= :: GaussianMixture -> GaussianMixture -> Bool
/= :: GaussianMixture -> GaussianMixture -> Bool
Eq, Int -> GaussianMixture -> ShowS
[GaussianMixture] -> ShowS
GaussianMixture -> String
(Int -> GaussianMixture -> ShowS)
-> (GaussianMixture -> String)
-> ([GaussianMixture] -> ShowS)
-> Show GaussianMixture
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> GaussianMixture -> ShowS
showsPrec :: Int -> GaussianMixture -> ShowS
$cshow :: GaussianMixture -> String
show :: GaussianMixture -> String
$cshowList :: [GaussianMixture] -> ShowS
showList :: [GaussianMixture] -> ShowS
Show)

logsumexp :: [Double] -> Double
logsumexp :: [Double] -> Double
logsumexp [Double]
xs = let m :: Double
m = [Double] -> Double
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum [Double]
xs in Double
m Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double -> Double
forall a. Floating a => a -> a
log ([Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ((Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (Double -> Double
forall a. Floating a => a -> a
exp (Double -> Double) -> (Double -> Double) -> Double -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Double -> Double -> Double
forall a. Num a => a -> a -> a
subtract Double
m) [Double]
xs))

-- | Compute normalized posterior component probabilities with log-sum-exp stabilization.
mixtureResponsibilities :: GaussianMixture -> Dataset -> Either MathError Dataset
mixtureResponsibilities :: GaussianMixture -> Dataset -> Either MathError Dataset
mixtureResponsibilities (GMM [Double]
weights Dataset
means [Dataset]
covariances) Dataset
xs = do
    let d :: Int
d = [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length (Dataset -> [Double]
forall a. [[a]] -> [a]
firstRow Dataset
means)
    Int -> Dataset -> Either MathError ()
queries Int
d Dataset
xs
    inverses <- (Dataset -> Either MathError Dataset)
-> [Dataset] -> Either MathError [Dataset]
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 Dataset -> Either MathError Dataset
inverse [Dataset]
covariances
    eigenvalues <- traverse symmetricEigen covariances
    require (all (> 0) (concatMap fst eigenvalues)) "mixture covariance must be positive definite"
    let logdet = (([Double], Dataset) -> Double)
-> [([Double], Dataset)] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map ([Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ([Double] -> Double)
-> (([Double], Dataset) -> [Double])
-> ([Double], Dataset)
-> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map Double -> Double
forall a. Floating a => a -> a
log ([Double] -> [Double])
-> (([Double], Dataset) -> [Double])
-> ([Double], Dataset)
-> [Double]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. ([Double], Dataset) -> [Double]
forall a b. (a, b) -> a
fst) [([Double], Dataset)]
eigenvalues
        logDensity [Double]
row Int
j =
            let delta :: [Double]
delta = [Double] -> [Double] -> [Double]
sub [Double]
row (Dataset
means Dataset -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! Int
j)
             in Double -> Double
forall a. Floating a => a -> a
log ([Double]
weights [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
j)
                    Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
0.5 Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
d Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
log (Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
forall a. Floating a => a
pi) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ [Double]
logdet [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
j Double -> Double -> Double
forall a. Num a => a -> a -> a
+ [Double] -> [Double] -> Double
dot [Double]
delta (([Double] -> Double) -> Dataset -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map ([Double] -> [Double] -> Double
`dot` [Double]
delta) ([Dataset]
inverses [Dataset] -> Int -> Dataset
forall a. HasCallStack => [a] -> Int -> a
!! Int
j)))
        answer =
            [ let logs :: [Double]
logs = [[Double] -> Int -> Double
logDensity [Double]
row Int
j | Int
j <- [Int
0 .. [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
weights Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]]; z :: Double
z = [Double] -> Double
logsumexp [Double]
logs
               in (Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (\Double
v -> Double -> Double
forall a. Floating a => a -> a
exp (Double
v Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
z)) [Double]
logs
            | [Double]
row <- Dataset
xs
            ]
    checked (concat answer)
    Right answer

-- | Full-covariance EM initialized by k-means++; positive ridge is added at each covariance update. Collapsed component mass fails.
fitGaussianMixture :: MLConfig -> Int -> Double -> Dataset -> Either MathError (Fit GaussianMixture)
fitGaussianMixture :: MLConfig
-> Int
-> Double
-> Dataset
-> Either MathError (Fit GaussianMixture)
fitGaussianMixture MLConfig
c Int
k Double
ridge Dataset
xs = do
    MLConfig -> Either MathError ()
config MLConfig
c
    (n, d) <- Dataset -> Either MathError (Int, Int)
dataset Dataset
xs
    require (finite ridge && ridge > 0) "mixture covariance ridge must be positive"
    initial <- kMeansPlusPlus c k xs
    let centers = ClusterModel -> Dataset
clusterCenters (Fit ClusterModel -> ClusterModel
forall a. Fit a -> a
fitted Fit ClusterModel
initial)
        n' = Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n
        initialCov = ([Double] -> [Double]) -> Dataset -> Dataset
forall a b. (a -> b) -> [a] -> [b]
map (Double -> [Double] -> [Double]
scale (Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
n')) (Dataset -> Dataset -> Dataset
matmul (Dataset -> Dataset
forall a. [[a]] -> [[a]]
transpose (([Double] -> [Double]) -> Dataset -> Dataset
forall a b. (a -> b) -> [a] -> [b]
map ([Double] -> [Double] -> [Double]
`sub` Dataset -> [Double]
meanV Dataset
xs) Dataset
xs)) (([Double] -> [Double]) -> Dataset -> Dataset
forall a b. (a -> b) -> [a] -> [b]
map ([Double] -> [Double] -> [Double]
`sub` Dataset -> [Double]
meanV Dataset
xs) Dataset
xs))
        regularize Dataset
cov = [[Dataset
cov Dataset -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! Int
i [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
j Double -> Double -> Double
forall a. Num a => a -> a -> a
+ if Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
j then Double
ridge else Double
0 | Int
j <- [Int
0 .. Int
d Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]] | Int
i <- [Int
0 .. Int
d Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]]
        step GaussianMixture
model = do
            responsibilities <- GaussianMixture -> Dataset -> Either MathError Dataset
mixtureResponsibilities GaussianMixture
model Dataset
xs
            let totals = ([Double] -> Double) -> Dataset -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum (Dataset -> Dataset
forall a. [[a]] -> [[a]]
transpose Dataset
responsibilities)
            require (all (> 1e-12) totals) "Gaussian mixture component collapsed"
            let means =
                    [ Double -> [Double] -> [Double]
scale (Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
total) (([Double] -> [Double] -> [Double]) -> Dataset -> [Double]
forall a. (a -> a -> a) -> [a] -> a
forall (t :: * -> *) a. Foldable t => (a -> a -> a) -> t a -> a
foldl1 [Double] -> [Double] -> [Double]
add ((Double -> [Double] -> [Double]) -> [Double] -> Dataset -> Dataset
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Double -> [Double] -> [Double]
scale [Double]
weights Dataset
xs))
                    | ([Double]
weights, Double
total) <- Dataset -> [Double] -> [([Double], Double)]
forall a b. [a] -> [b] -> [(a, b)]
zip (Dataset -> Dataset
forall a. [[a]] -> [[a]]
transpose Dataset
responsibilities) [Double]
totals
                    ]
                covariances =
                    [ Dataset -> Dataset
regularize
                        [ [ [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum
                                [ (Dataset
responsibilities Dataset -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! Int
i [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
j) Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Dataset
xs Dataset -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! Int
i [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
a Double -> Double -> Double
forall a. Num a => a -> a -> a
- Dataset
means Dataset -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! Int
j [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
a) Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Dataset
xs Dataset -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! Int
i [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
b Double -> Double -> Double
forall a. Num a => a -> a -> a
- Dataset
means Dataset -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! Int
j [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
b)
                                | Int
i <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]
                                ]
                                Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ ([Double]
totals [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
j)
                          | Int
b <- [Int
0 .. Int
d Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]
                          ]
                        | Int
a <- [Int
0 .. Int
d Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]
                        ]
                    | Int
j <- [Int
0 .. Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]
                    ]
            checked (concat means ++ concatMap concat covariances)
            Right (GMM (map (/ n') totals) means covariances)
        difference GaussianMixture
a GaussianMixture
b =
            [Double] -> Double
norm ([Double] -> [Double] -> [Double]
sub (GaussianMixture -> [Double]
mixtureWeights GaussianMixture
a) (GaussianMixture -> [Double]
mixtureWeights GaussianMixture
b))
                Double -> Double -> Double
forall a. Num a => a -> a -> a
+ [Double] -> Double
norm (Dataset -> [Double]
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat (([Double] -> [Double] -> [Double]) -> Dataset -> Dataset -> Dataset
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith [Double] -> [Double] -> [Double]
sub (GaussianMixture -> Dataset
mixtureMeans GaussianMixture
a) (GaussianMixture -> Dataset
mixtureMeans GaussianMixture
b)))
                Double -> Double -> Double
forall a. Num a => a -> a -> a
+ [Double] -> Double
norm ((Dataset -> [Double]) -> [Dataset] -> [Double]
forall (t :: * -> *) a b. Foldable t => (a -> [b]) -> t a -> [b]
concatMap Dataset -> [Double]
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat ((Dataset -> Dataset -> Dataset)
-> [Dataset] -> [Dataset] -> [Dataset]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (([Double] -> [Double] -> [Double]) -> Dataset -> Dataset -> Dataset
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith [Double] -> [Double] -> [Double]
sub) (GaussianMixture -> [Dataset]
mixtureCovariances GaussianMixture
a) (GaussianMixture -> [Dataset]
mixtureCovariances GaussianMixture
b)))
    iterateFit
        c
        step
        difference
        (GMM (replicate k (1 / fromIntegral k)) centers (replicate k (regularize initialCov)))

-- | Density expansion using a positive Euclidean radius; minimumPoints includes the observation itself.
dbscan :: Double -> Int -> Dataset -> Either MathError [Int]
dbscan :: Double -> Int -> Dataset -> Either MathError [Int]
dbscan Double
epsilon Int
minimumPoints Dataset
xs = do
    (n, _) <- Dataset -> Either MathError (Int, Int)
dataset Dataset
xs
    require
        (finite epsilon && epsilon > 0 && minimumPoints > 0 && minimumPoints <= n)
        "invalid DBSCAN parameters"
    let distances = [[[Double] -> Double
norm ([Double] -> [Double] -> [Double]
sub [Double]
x [Double]
y) | [Double]
y <- Dataset
xs] | [Double]
x <- Dataset
xs]
    checked (concat distances)
    let neighbors Int
i = [Int
j | Int
j <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1], 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
epsilon]
        expand t
_ [] [Int]
visited [t]
labels = ([Int]
visited, [t]
labels)
        expand t
label (Int
j : [Int]
todo) [Int]
visited [t]
labels =
            let unseen :: Bool
unseen = Int
j Int -> [Int] -> Bool
forall (t :: * -> *) a. (Foldable t, Eq a) => a -> t a -> Bool
`notElem` [Int]
visited
                visited' :: [Int]
visited' = if Bool
unseen then Int
j Int -> [Int] -> [Int]
forall a. a -> [a] -> [a]
: [Int]
visited else [Int]
visited
                todo' :: [Int]
todo' = if Bool
unseen Bool -> Bool -> Bool
&& [Int] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length (Int -> [Int]
neighbors Int
j) Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
minimumPoints then [Int] -> [Int]
forall a. Eq a => [a] -> [a]
nub ([Int]
todo [Int] -> [Int] -> [Int]
forall a. [a] -> [a] -> [a]
++ Int -> [Int]
neighbors Int
j) else [Int]
todo
                labels' :: [t]
labels' = if [t]
labels [t] -> Int -> t
forall a. HasCallStack => [a] -> Int -> a
!! Int
j t -> t -> Bool
forall a. Ord a => a -> a -> Bool
< t
0 then Int -> t -> [t] -> [t]
forall a. Int -> a -> [a] -> [a]
replace Int
j t
label [t]
labels else [t]
labels
             in t -> [Int] -> [Int] -> [t] -> ([Int], [t])
expand t
label [Int]
todo' [Int]
visited' [t]
labels'
        go Int
i t
label [Int]
visited [t]
labels
            | Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
n = [t]
labels
            | Int
i Int -> [Int] -> Bool
forall a. Eq a => a -> [a] -> Bool
forall (t :: * -> *) a. (Foldable t, Eq a) => a -> t a -> Bool
`elem` [Int]
visited = Int -> t -> [Int] -> [t] -> [t]
go (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) t
label [Int]
visited [t]
labels
            | [Int] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length (Int -> [Int]
neighbors Int
i) Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
minimumPoints = Int -> t -> [Int] -> [t] -> [t]
go (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) t
label (Int
i Int -> [Int] -> [Int]
forall a. a -> [a] -> [a]
: [Int]
visited) [t]
labels
            | Bool
otherwise =
                let ([Int]
v, [t]
l) = t -> [Int] -> [Int] -> [t] -> ([Int], [t])
forall {t}.
(Ord t, Num t) =>
t -> [Int] -> [Int] -> [t] -> ([Int], [t])
expand t
label (Int -> [Int]
neighbors Int
i) (Int
i Int -> [Int] -> [Int]
forall a. a -> [a] -> [a]
: [Int]
visited) (Int -> t -> [t] -> [t]
forall a. Int -> a -> [a] -> [a]
replace Int
i t
label [t]
labels)
                 in Int -> t -> [Int] -> [t] -> [t]
go (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) (t
label t -> t -> t
forall a. Num a => a -> a -> a
+ t
1) [Int]
v [t]
l
    Right (go 0 0 [] (replicate n (-1)))

{- | OPTICS (Ankerst et al., SIGMOD 1999): ordered point index, core distance,
reachability distance. Undefined reachability starts each connected expansion.
-}
optics :: Double -> Int -> Dataset -> Either MathError [(Int, Maybe Double, Maybe Double)]
optics :: Double
-> Int
-> Dataset
-> Either MathError [(Int, Maybe Double, Maybe Double)]
optics Double
epsilon Int
minimumPoints Dataset
xs = do
    (n, _) <- Dataset -> Either MathError (Int, Int)
dataset Dataset
xs
    require
        (finite epsilon && epsilon > 0 && minimumPoints > 0 && minimumPoints <= n)
        "invalid OPTICS parameters"
    let distances = [[[Double] -> Double
norm ([Double] -> [Double] -> [Double]
sub [Double]
x [Double]
y) | [Double]
y <- Dataset
xs] | [Double]
x <- Dataset
xs]
    checked (concat distances)
    let near Int
i = ((Double, Int) -> Double) -> [(Double, Int)] -> [(Double, Int)]
forall b a. Ord b => (a -> b) -> [a] -> [a]
sortOn (Double, Int) -> Double
forall a b. (a, b) -> a
fst [(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
j <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1], 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
epsilon]
        core Int
i = if [(Double, Int)] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length (Int -> [(Double, Int)]
near Int
i) Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
minimumPoints then Double -> Maybe Double
forall a. a -> Maybe a
Just ((Double, Int) -> Double
forall a b. (a, b) -> a
fst (Int -> [(Double, Int)]
near Int
i [(Double, Int)] -> Int -> (Double, Int)
forall a. HasCallStack => [a] -> Int -> a
!! (Int
minimumPoints Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1))) else Maybe Double
forall a. Maybe a
Nothing
        visit Int
i Maybe Double
reached [Int]
processed [(Int, Double)]
queue [(Int, Maybe Double, Maybe Double)]
result =
            let processed' :: [Int]
processed' = Int
i Int -> [Int] -> [Int]
forall a. a -> [a] -> [a]
: [Int]
processed
                update :: b -> [(Int, b)] -> (b, Int) -> [(Int, b)]
update b
cd [(Int, b)]
q (b
distance, Int
j)
                    | Int
j Int -> [Int] -> Bool
forall a. Eq a => a -> [a] -> Bool
forall (t :: * -> *) a. (Foldable t, Eq a) => a -> t a -> Bool
`elem` [Int]
processed' = [(Int, b)]
q
                    | Bool
otherwise =
                        let r :: b
r = b -> b -> b
forall a. Ord a => a -> a -> a
max b
cd b
distance
                         in case Int -> [(Int, b)] -> Maybe b
forall a b. Eq a => a -> [(a, b)] -> Maybe b
lookup Int
j [(Int, b)]
q of
                                Maybe b
Nothing -> (Int
j, b
r) (Int, b) -> [(Int, b)] -> [(Int, b)]
forall a. a -> [a] -> [a]
: [(Int, b)]
q
                                Just b
old -> (Int
j, b -> b -> b
forall a. Ord a => a -> a -> a
min b
old b
r) (Int, b) -> [(Int, b)] -> [(Int, b)]
forall a. a -> [a] -> [a]
: ((Int, b) -> Bool) -> [(Int, b)] -> [(Int, b)]
forall a. (a -> Bool) -> [a] -> [a]
filter ((Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
j) (Int -> Bool) -> ((Int, b) -> Int) -> (Int, b) -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Int, b) -> Int
forall a b. (a, b) -> a
fst) [(Int, b)]
q
                queue' :: [(Int, Double)]
queue' = case Int -> Maybe Double
core Int
i of
                    Maybe Double
Nothing -> [(Int, Double)]
queue
                    Just Double
cd -> ([(Int, Double)] -> (Double, Int) -> [(Int, Double)])
-> [(Int, Double)] -> [(Double, Int)] -> [(Int, Double)]
forall b a. (b -> a -> b) -> b -> [a] -> b
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl (Double -> [(Int, Double)] -> (Double, Int) -> [(Int, Double)]
forall {b}. Ord b => b -> [(Int, b)] -> (b, Int) -> [(Int, b)]
update Double
cd) [(Int, Double)]
queue (Int -> [(Double, Int)]
near Int
i)
                result' :: [(Int, Maybe Double, Maybe Double)]
result' = [(Int, Maybe Double, Maybe Double)]
result [(Int, Maybe Double, Maybe Double)]
-> [(Int, Maybe Double, Maybe Double)]
-> [(Int, Maybe Double, Maybe Double)]
forall a. [a] -> [a] -> [a]
++ [(Int
i, Int -> Maybe Double
core Int
i, Maybe Double
reached)]
             in case ((Int, Double) -> Double) -> [(Int, Double)] -> [(Int, Double)]
forall b a. Ord b => (a -> b) -> [a] -> [a]
sortOn (Int, Double) -> Double
forall a b. (a, b) -> b
snd [(Int, Double)]
queue' of
                    [] -> ([Int]
processed', [(Int, Maybe Double, Maybe Double)]
result')
                    (Int
j, Double
r) : [(Int, Double)]
rest -> Int
-> Maybe Double
-> [Int]
-> [(Int, Double)]
-> [(Int, Maybe Double, Maybe Double)]
-> ([Int], [(Int, Maybe Double, Maybe Double)])
visit Int
j (Double -> Maybe Double
forall a. a -> Maybe a
Just Double
r) [Int]
processed' [(Int, Double)]
rest [(Int, Maybe Double, Maybe Double)]
result'
        go Int
i [Int]
processed [(Int, Maybe Double, Maybe Double)]
result
            | Int
i Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
n = [(Int, Maybe Double, Maybe Double)]
result
            | Int
i Int -> [Int] -> Bool
forall a. Eq a => a -> [a] -> Bool
forall (t :: * -> *) a. (Foldable t, Eq a) => a -> t a -> Bool
`elem` [Int]
processed = Int
-> [Int]
-> [(Int, Maybe Double, Maybe Double)]
-> [(Int, Maybe Double, Maybe Double)]
go (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) [Int]
processed [(Int, Maybe Double, Maybe Double)]
result
            | Bool
otherwise = let ([Int]
p, [(Int, Maybe Double, Maybe Double)]
r) = Int
-> Maybe Double
-> [Int]
-> [(Int, Double)]
-> [(Int, Maybe Double, Maybe Double)]
-> ([Int], [(Int, Maybe Double, Maybe Double)])
visit Int
i Maybe Double
forall a. Maybe a
Nothing [Int]
processed [] [(Int, Maybe Double, Maybe Double)]
result in Int
-> [Int]
-> [(Int, Maybe Double, Maybe Double)]
-> [(Int, Maybe Double, Maybe Double)]
go (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) [Int]
p [(Int, Maybe Double, Maybe Double)]
r
    Right (go 0 [] [])

-- | Extract DBSCAN-style labels at epsilon from a checked OPTICS ordering; noise is -1.
extractOPTICS :: Double -> [(Int, Maybe Double, Maybe Double)] -> Either MathError [Int]
extractOPTICS :: Double
-> [(Int, Maybe Double, Maybe Double)] -> Either MathError [Int]
extractOPTICS Double
epsilon [(Int, Maybe Double, Maybe Double)]
ordered = do
    let n :: Int
n = [(Int, Maybe Double, Maybe Double)] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [(Int, Maybe Double, Maybe Double)]
ordered
    Bool -> String -> Either MathError ()
require
        (Double -> Bool
finite Double
epsilon Bool -> Bool -> Bool
&& Double
epsilon Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0 Bool -> Bool -> Bool
&& (Int -> Int) -> [Int] -> [Int]
forall b a. Ord b => (a -> b) -> [a] -> [a]
sortOn Int -> Int
forall a. a -> a
id (((Int, Maybe Double, Maybe Double) -> Int)
-> [(Int, Maybe Double, Maybe Double)] -> [Int]
forall a b. (a -> b) -> [a] -> [b]
map (\(Int
i, Maybe Double
_, Maybe Double
_) -> Int
i) [(Int, Maybe Double, Maybe Double)]
ordered) [Int] -> [Int] -> Bool
forall a. Eq a => a -> a -> Bool
== [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1])
        String
"invalid OPTICS ordering/threshold"
    Bool -> String -> Either MathError ()
require
        (((Int, Maybe Double, Maybe Double) -> Bool)
-> [(Int, Maybe Double, Maybe Double)] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (\(Int
_, Maybe Double
c, Maybe Double
r) -> (Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (\Double
v -> Double -> Bool
finite Double
v Bool -> Bool -> Bool
&& Double
v Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
0) ([Double] -> (Double -> [Double]) -> Maybe Double -> [Double]
forall b a. b -> (a -> b) -> Maybe a -> b
maybe [] (Double -> [Double] -> [Double]
forall a. a -> [a] -> [a]
: []) Maybe Double
c [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ [Double] -> (Double -> [Double]) -> Maybe Double -> [Double]
forall b a. b -> (a -> b) -> Maybe a -> b
maybe [] (Double -> [Double] -> [Double]
forall a. a -> [a] -> [a]
: []) Maybe Double
r)) [(Int, Maybe Double, Maybe Double)]
ordered)
        String
"invalid OPTICS distances"
    let (Int
_, [(Int, Int)]
assignments) =
            ((Int, [(Int, Int)])
 -> (Int, Maybe Double, Maybe Double) -> (Int, [(Int, Int)]))
-> (Int, [(Int, Int)])
-> [(Int, Maybe Double, Maybe Double)]
-> (Int, [(Int, Int)])
forall b a. (b -> a -> b) -> b -> [a] -> b
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl
                ( \(Int
cluster, [(Int, Int)]
labels) (Int
i, Maybe Double
cd, Maybe Double
rd) ->
                    if Bool -> (Double -> Bool) -> Maybe Double -> Bool
forall b a. b -> (a -> b) -> Maybe a -> b
maybe Bool
True (Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
epsilon) Maybe Double
rd
                        then
                            if Bool -> (Double -> Bool) -> Maybe Double -> Bool
forall b a. b -> (a -> b) -> Maybe a -> b
maybe Bool
False (Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
epsilon) Maybe Double
cd
                                then (Int
cluster Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1, (Int
i, Int
cluster Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) (Int, Int) -> [(Int, Int)] -> [(Int, Int)]
forall a. a -> [a] -> [a]
: [(Int, Int)]
labels)
                                else (Int
cluster, (Int
i, -Int
1) (Int, Int) -> [(Int, Int)] -> [(Int, Int)]
forall a. a -> [a] -> [a]
: [(Int, Int)]
labels)
                        else (Int
cluster, (Int
i, Int
cluster) (Int, Int) -> [(Int, Int)] -> [(Int, Int)]
forall a. a -> [a] -> [a]
: [(Int, Int)]
labels)
                )
                (-Int
1, [])
                [(Int, Maybe Double, Maybe Double)]
ordered
    [Int] -> Either MathError [Int]
forall a b. b -> Either a b
Right [Int -> (Int -> Int) -> Maybe Int -> Int
forall b a. b -> (a -> b) -> Maybe a -> b
maybe (-Int
1) Int -> Int
forall a. a -> a
id (Int -> [(Int, Int)] -> Maybe Int
forall a b. Eq a => a -> [(a, b)] -> Maybe b
lookup Int
i [(Int, Int)]
assignments) | Int
i <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]]

-- | Gaussian-kernel mean shift from every training point; modes within bandwidth/2 are merged.
meanShift :: MLConfig -> Double -> Dataset -> Either MathError (Fit ClusterModel)
meanShift :: MLConfig
-> Double -> Dataset -> Either MathError (Fit ClusterModel)
meanShift MLConfig
c Double
bandwidth Dataset
xs = do
    MLConfig -> Either MathError ()
config MLConfig
c
    _ <- Dataset -> Either MathError (Int, Int)
dataset Dataset
xs
    require (finite bandwidth && bandwidth > 0) "mean-shift bandwidth must be positive"
    let step t [Double]
centers =
            ([Double] -> Either MathError [Double])
-> t [Double] -> Either MathError (t [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) -> t a -> f (t b)
traverse
                ( \[Double]
mu -> do
                    let weights :: [Double]
weights = ([Double] -> Double) -> Dataset -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (\[Double]
x -> Double -> Double
forall a. Floating a => a -> a
exp (-[Double] -> [Double] -> Double
dist2 [Double]
x [Double]
mu Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
bandwidth Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
bandwidth))) Dataset
xs; total :: Double
total = [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Double]
weights
                    Bool -> String -> Either MathError ()
require (Double
total Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0 Bool -> Bool -> Bool
&& Double -> Bool
finite Double
total) String
"mean-shift kernel mass vanished"
                    [Double] -> Either MathError [Double]
forall a b. b -> Either a b
Right (Double -> [Double] -> [Double]
scale (Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
total) (([Double] -> [Double] -> [Double]) -> Dataset -> [Double]
forall a. (a -> a -> a) -> [a] -> a
forall (t :: * -> *) a. Foldable t => (a -> a -> a) -> t a -> a
foldl1 [Double] -> [Double] -> [Double]
add ((Double -> [Double] -> [Double]) -> [Double] -> Dataset -> Dataset
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Double -> [Double] -> [Double]
scale [Double]
weights Dataset
xs)))
                )
                t [Double]
centers
    fit <- iterateFit c step (\Dataset
a Dataset
b -> [Double] -> Double
norm (Dataset -> [Double]
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat (([Double] -> [Double] -> [Double]) -> Dataset -> Dataset -> Dataset
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith [Double] -> [Double] -> [Double]
sub Dataset
a Dataset
b))) xs
    let uniqueModes =
            (Dataset -> [Double] -> Dataset) -> Dataset -> Dataset -> Dataset
forall b a. (b -> a -> b) -> b -> [a] -> b
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl
                (\Dataset
m [Double]
x -> if ([Double] -> Bool) -> Dataset -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
any (\[Double]
y -> [Double] -> [Double] -> Double
dist2 [Double]
x [Double]
y Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
bandwidth Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
bandwidth Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
4) Dataset
m then Dataset
m else Dataset
m Dataset -> Dataset -> Dataset
forall a. [a] -> [a] -> [a]
++ [[Double]
x])
                []
                (Fit Dataset -> Dataset
forall a. Fit a -> a
fitted Fit Dataset
fit)
    Right fit{fitted = Clusters uniqueModes}

-- | Single, complete, average Euclidean distance or Ward variance-increase linkage.
data Linkage = Single | Complete | Average | Ward deriving (Linkage -> Linkage -> Bool
(Linkage -> Linkage -> Bool)
-> (Linkage -> Linkage -> Bool) -> Eq Linkage
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: Linkage -> Linkage -> Bool
== :: Linkage -> Linkage -> Bool
$c/= :: Linkage -> Linkage -> Bool
/= :: Linkage -> Linkage -> Bool
Eq, Int -> Linkage -> ShowS
[Linkage] -> ShowS
Linkage -> String
(Int -> Linkage -> ShowS)
-> (Linkage -> String) -> ([Linkage] -> ShowS) -> Show Linkage
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> Linkage -> ShowS
showsPrec :: Int -> Linkage -> ShowS
$cshow :: Linkage -> String
show :: Linkage -> String
$cshowList :: [Linkage] -> ShowS
showList :: [Linkage] -> ShowS
Show)

{- | One hierarchical-clustering merge, retaining child IDs,linkage distance,and merged
observation count.
-}
data Merge = Merge
    { Merge -> Int
leftCluster :: Int
    -- ^ Left child cluster ID in a hierarchical merge.
    , Merge -> Int
rightCluster :: Int
    -- ^ Right child cluster ID in a hierarchical merge.
    , Merge -> Double
mergeDistance :: Double
    -- ^ Linkage distance at which the two child clusters merge.
    , Merge -> Int
mergedSize :: Int
    -- ^ Number of original observations in the merged cluster.
    }
    deriving (Merge -> Merge -> Bool
(Merge -> Merge -> Bool) -> (Merge -> Merge -> Bool) -> Eq Merge
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: Merge -> Merge -> Bool
== :: Merge -> Merge -> Bool
$c/= :: Merge -> Merge -> Bool
/= :: Merge -> Merge -> Bool
Eq, Int -> Merge -> ShowS
[Merge] -> ShowS
Merge -> String
(Int -> Merge -> ShowS)
-> (Merge -> String) -> ([Merge] -> ShowS) -> Show Merge
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> Merge -> ShowS
showsPrec :: Int -> Merge -> ShowS
$cshow :: Merge -> String
show :: Merge -> String
$cshowList :: [Merge] -> ShowS
showList :: [Merge] -> ShowS
Show)

-- | Agglomerate the closest active clusters to produce all n-1 linkage rows.
hierarchical :: Linkage -> Dataset -> Either MathError [Merge]
hierarchical :: Linkage -> Dataset -> Either MathError [Merge]
hierarchical Linkage
method Dataset
xs = do
    (n, _) <- Dataset -> Either MathError (Int, Int)
dataset Dataset
xs
    let distance [Int]
as [Int]
bs = case Linkage
method of
            Linkage
Single -> [Double] -> Double
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
minimum [Double]
ds
            Linkage
Complete -> [Double] -> Double
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum [Double]
ds
            Linkage
Average -> [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Double]
ds 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]
ds)
            Linkage
Ward ->
                Double -> Double
forall a. Floating a => a -> a
sqrt
                    ( Double
2
                        Double -> Double -> Double
forall a. Num 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]
as Int -> Int -> Int
forall a. Num a => a -> a -> a
* [Int] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Int]
bs)
                        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]
as Int -> Int -> Int
forall a. Num a => a -> a -> a
+ [Int] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Int]
bs)
                        Double -> Double -> Double
forall a. Num a => a -> a -> a
* [Double] -> [Double] -> Double
dist2 (Dataset -> [Double]
meanV ((Int -> [Double]) -> [Int] -> Dataset
forall a b. (a -> b) -> [a] -> [b]
map (Dataset
xs Dataset -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!!) [Int]
as)) (Dataset -> [Double]
meanV ((Int -> [Double]) -> [Int] -> Dataset
forall a b. (a -> b) -> [a] -> [b]
map (Dataset
xs Dataset -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!!) [Int]
bs))
                    )
          where
            ds :: [Double]
ds = [Double -> Double
forall a. Floating a => a -> a
sqrt ([Double] -> [Double] -> Double
dist2 (Dataset
xs Dataset -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! Int
i) (Dataset
xs Dataset -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! Int
j)) | Int
i <- [Int]
as, Int
j <- [Int]
bs]
        go Int
_ [(Int, [Int])
_] [Merge]
result = [Merge] -> Either MathError [Merge]
forall a b. b -> Either a b
Right [Merge]
result
        go Int
next [(Int, [Int])]
clusters [Merge]
result = do
            let pairs :: [(Double, Int, Int, [Int])]
pairs = [([Int] -> [Int] -> Double
distance [Int]
as [Int]
bs, Int
ca, Int
cb, [Int]
as [Int] -> [Int] -> [Int]
forall a. [a] -> [a] -> [a]
++ [Int]
bs) | (Int
ca, [Int]
as) <- [(Int, [Int])]
clusters, (Int
cb, [Int]
bs) <- [(Int, [Int])]
clusters, Int
ca Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
cb]
                (Double
height, Int
a, Int
b, [Int]
members) = ((Double, Int, Int, [Int]) -> Double)
-> [(Double, Int, Int, [Int])] -> [(Double, Int, Int, [Int])]
forall b a. Ord b => (a -> b) -> [a] -> [a]
sortOn (\(Double
v, Int
_, Int
_, [Int]
_) -> Double
v) [(Double, Int, Int, [Int])]
pairs [(Double, Int, Int, [Int])] -> Int -> (Double, Int, Int, [Int])
forall a. HasCallStack => [a] -> Int -> a
!! Int
0
            [Double] -> Either MathError ()
checked [Double
height]
            Int -> [(Int, [Int])] -> [Merge] -> Either MathError [Merge]
go
                (Int
next Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1)
                ((Int
next, [Int]
members) (Int, [Int]) -> [(Int, [Int])] -> [(Int, [Int])]
forall a. a -> [a] -> [a]
: ((Int, [Int]) -> Bool) -> [(Int, [Int])] -> [(Int, [Int])]
forall a. (a -> Bool) -> [a] -> [a]
filter (\(Int
j, [Int]
_) -> Int
j Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
a Bool -> Bool -> Bool
&& Int
j Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
b) [(Int, [Int])]
clusters)
                ([Merge]
result [Merge] -> [Merge] -> [Merge]
forall a. [a] -> [a] -> [a]
++ [Int -> Int -> Double -> Int -> Merge
Merge Int
a Int
b Double
height ([Int] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Int]
members)])
    go n [(i, [i]) | i <- [0 .. n - 1]] []

-- | Validate the entire hierarchy and return k cluster labels in original observation order.
cutHierarchy :: Int -> Int -> [Merge] -> Either MathError [Int]
cutHierarchy :: Int -> Int -> [Merge] -> Either MathError [Int]
cutHierarchy Int
n Int
k [Merge]
merges = do
    Bool -> String -> Either MathError ()
require (Int
n 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
0 Bool -> Bool -> Bool
&& Int
k Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
n Bool -> Bool -> Bool
&& [Merge] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Merge]
merges Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) String
"invalid hierarchy cut"
    let step :: [(Int, [a])] -> (Int, Merge) -> Either MathError [(Int, [a])]
step [(Int, [a])]
clusters (Int
j, Merge Int
a Int
b Double
height Int
size) = do
            Bool -> String -> Either MathError ()
require (Int
a Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
b Bool -> Bool -> Bool
&& Double -> Bool
finite Double
height Bool -> Bool -> Bool
&& Double
height Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
0) String
"invalid hierarchy merge"
            as <- Either MathError [a]
-> ([a] -> Either MathError [a])
-> Maybe [a]
-> Either MathError [a]
forall b a. b -> (a -> b) -> Maybe a -> b
maybe (MathError -> Either MathError [a]
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"ML: unknown hierarchy cluster")) [a] -> Either MathError [a]
forall a b. b -> Either a b
Right (Int -> [(Int, [a])] -> Maybe [a]
forall a b. Eq a => a -> [(a, b)] -> Maybe b
lookup Int
a [(Int, [a])]
clusters)
            bs <- maybe (Left (InvalidArgument "ML: unknown hierarchy cluster")) Right (lookup b clusters)
            require (size == length as + length bs) "invalid hierarchy size"
            Right ((n + j, as ++ bs) : filter (\(Int
i, [a]
_) -> Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
a Bool -> Bool -> Bool
&& Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
b) clusters)
    -- Validate the complete hierarchy, including merges after the requested cut.
    _ <- ([(Int, [Int])] -> (Int, Merge) -> Either MathError [(Int, [Int])])
-> [(Int, [Int])]
-> [(Int, Merge)]
-> Either MathError [(Int, [Int])]
forall (t :: * -> *) (m :: * -> *) b a.
(Foldable t, Monad m) =>
(b -> a -> m b) -> b -> t a -> m b
foldM [(Int, [Int])] -> (Int, Merge) -> Either MathError [(Int, [Int])]
forall {a}.
[(Int, [a])] -> (Int, Merge) -> Either MathError [(Int, [a])]
step [(Int
i, [Int
i]) | Int
i <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]] ([Int] -> [Merge] -> [(Int, Merge)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
0 ..] [Merge]
merges)
    clusters <- foldM step [(i, [i]) | i <- [0 .. n - 1]] (take (n - k) (zip [0 ..] merges))
    Right
        [ maybe
            (-1)
            id
            ( lookup
                i
                [(point, label) | (label, (_, members)) <- zip [0 ..] (sortOn fst clusters), point <- members]
            )
        | i <- [0 .. n - 1]
        ]

-- | Cluster row-normalized leading eigenvectors of the normalized Gaussian affinity matrix.
spectralClustering :: MLConfig -> Int -> Double -> Dataset -> Either MathError (Fit [Int])
spectralClustering :: MLConfig
-> Int -> Double -> Dataset -> Either MathError (Fit [Int])
spectralClustering MLConfig
c Int
k Double
gamma Dataset
xs = do
    MLConfig -> Either MathError ()
config MLConfig
c
    (n, _) <- Dataset -> Either MathError (Int, Int)
dataset Dataset
xs
    require (k > 0 && k < n && finite gamma && gamma > 0) "invalid spectral clustering parameters"
    let affinity =
            [ [if Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
j then Double
0 else Double -> Double
forall a. Floating a => a -> a
exp (-Double
gamma Double -> Double -> Double
forall a. Num a => a -> a -> a
* [Double] -> [Double] -> Double
dist2 (Dataset
xs Dataset -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! Int
i) (Dataset
xs Dataset -> 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
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]
            ]
        degrees = ([Double] -> Double) -> Dataset -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum Dataset
affinity
    require (all (> 0) degrees) "spectral affinity graph has isolated vertices"
    (_, vectors) <-
        symmetricEigen
            [[affinity !! i !! j / sqrt (degrees !! i * degrees !! j) | j <- [0 .. n - 1]] | i <- [0 .. n - 1]]
    let embedding = [let r :: [Double]
r = Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
take Int
k [Double]
v; s :: Double
s = [Double] -> Double
norm [Double]
r in if Double
s Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0 then [Double]
r else Double -> [Double] -> [Double]
scale (Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
s) [Double]
r | [Double]
v <- Dataset
vectors]
    fit <- kMeansPlusPlus c k embedding
    labels <- predictClusters (fitted fit) embedding
    Right fit{fitted = labels}

{- | The preference is the common diagonal similarity; damping is in [0.5,1).
Stops only after exemplar identities remain unchanged for ten iterations.
-}
affinityPropagation :: MLConfig -> Double -> Double -> Dataset -> Either MathError (Fit [Int])
affinityPropagation :: MLConfig
-> Double -> Double -> Dataset -> Either MathError (Fit [Int])
affinityPropagation MLConfig
c Double
preference Double
damping Dataset
xs = do
    MLConfig -> Either MathError ()
config MLConfig
c
    (n, _) <- Dataset -> Either MathError (Int, Int)
dataset Dataset
xs
    require
        (n >= 2 && finite preference && finite damping && damping >= 0.5 && damping < 1)
        "invalid affinity propagation parameters"
    let similarity =
            [ [if Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
j then Double
preference else Double -> Double
forall a. Num a => a -> a
negate ([Double] -> [Double] -> Double
dist2 (Dataset
xs Dataset -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! Int
i) (Dataset
xs Dataset -> 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
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]
            ]
        mix Double
old Double
new = Double
damping Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
old Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
damping) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
new
        go Int
iteration Dataset
responsibilities Dataset
availabilities [Int]
previous t
stable = do
            let rs :: Dataset
rs =
                    [ [ Double -> Double -> Double
mix
                            (Dataset
responsibilities Dataset -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! Int
i [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
k)
                            ( Dataset
similarity Dataset -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! Int
i [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
k
                                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
maximum [Dataset
availabilities Dataset -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! Int
i [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
j Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Dataset
similarity 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
j Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
k]
                            )
                      | Int
k <- [Int
0 .. Int
n 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]
                    ]
                av :: Dataset
av =
                    [ [ Double -> Double -> Double
mix
                            (Dataset
availabilities Dataset -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! Int
i [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
k)
                            ( if Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
k
                                then [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
0 (Dataset
rs Dataset -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! Int
j [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
k) | Int
j <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1], Int
j Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
k]
                                else Double -> Double -> Double
forall a. Ord a => a -> a -> a
min Double
0 (Dataset
rs Dataset -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! Int
k [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
k Double -> Double -> Double
forall a. Num a => a -> a -> a
+ [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
0 (Dataset
rs Dataset -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! Int
j [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
k) | Int
j <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1], Int
j Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
i, Int
j Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
k])
                            )
                      | Int
k <- [Int
0 .. Int
n 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]
                    ]
                exemplars :: [Int]
exemplars = [Int
k | Int
k <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1], Dataset
rs Dataset -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! Int
k [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
k Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Dataset
av Dataset -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! Int
k [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
k Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0]
                stable' :: t
stable' = if [Int]
exemplars [Int] -> [Int] -> Bool
forall a. Eq a => a -> a -> Bool
== [Int]
previous Bool -> Bool -> Bool
&& Bool -> Bool
not ([Int] -> Bool
forall a. [a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [Int]
exemplars) then t
stable t -> t -> t
forall a. Num a => a -> a -> a
+ t
1 else t
0
            [Double] -> Either MathError ()
checked (Dataset -> [Double]
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat Dataset
rs [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ Dataset -> [Double]
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat Dataset
av)
            if t
stable' t -> t -> Bool
forall a. Ord a => a -> a -> Bool
>= t
10 Bool -> Bool -> Bool
|| Int
iteration Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= MLConfig -> Int
mlIterations MLConfig
c
                then do
                    Bool -> String -> Either MathError ()
require (Bool -> Bool
not ([Int] -> Bool
forall a. [a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [Int]
exemplars)) String
"affinity propagation found no exemplar within budget"
                    let labels :: [Int]
labels = [[Double] -> Int
forall a. Ord a => [a] -> Int
argmax [Dataset
similarity 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]
exemplars] | Int
i <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]]
                    Fit [Int] -> Either MathError (Fit [Int])
forall a b. b -> Either a b
Right ([Int] -> Int -> Double -> Bool -> Fit [Int]
forall a. a -> Int -> Double -> Bool -> Fit a
Fit [Int]
labels Int
iteration (t -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (t -> t -> t
forall a. Ord a => a -> a -> a
max t
0 (t
10 t -> t -> t
forall a. Num a => a -> a -> a
- t
stable'))) (t
stable' t -> t -> Bool
forall a. Ord a => a -> a -> Bool
>= t
10))
                else Int
-> Dataset -> Dataset -> [Int] -> t -> Either MathError (Fit [Int])
go (Int
iteration Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Dataset
rs Dataset
av [Int]
exemplars t
stable'
    go 1 (replicate n (replicate n 0)) (replicate n (replicate n 0)) [] (0 :: Int)