-- SPDX-License-Identifier: MPL-2.0

{- | Dense decomposition and manifold reference algorithms. PCA signs and
degenerate eigenspaces are not unique. t-SNE implements exact all-pairs
Student-t KL gradients (van der Maaten & Hinton, JMLR 2008); UMAP uses the
fuzzy union and full cross-entropy gradient (McInnes et al., arXiv:1802.03426)
with caller-supplied positive a,b curve parameters. Neither is a Barnes-Hut,
ANN, or stochastic negative-sampling implementation. Budgets are explicit.
-}
module HMath.MachineLearning.Decomposition (
    PCA,
    pcaMean,
    pcaComponents,
    pcaVariances,
    fitPCA,
    transformPCA,
    inversePCA,
    IncrementalPCA,
    incrementalPCA,
    updatePCA,
    mergePCA,
    finishPCA,
    KernelPCA,
    fitKernelPCA,
    transformKernelPCA,
    truncatedSVD,
    ICA,
    fitICA,
    transformICA,
    FactorAnalysis,
    factorLoadings,
    factorNoise,
    fitFactorAnalysis,
    CCA,
    fitCCA,
    transformCCA,
    canonicalCorrelations,
    PLS,
    fitPLS,
    predictPLS,
    randomProjection,
    classicalMDS,
    isomap,
    tsne,
    umap,
) where

import Control.Monad (foldM)
import Data.List (sortOn, transpose)
import HMath.Error
import HMath.MachineLearning.Internal
import HMath.MachineLearning.Regression (Kernel, kernel)
import qualified HMath.Matrix.LAPACK as L
import qualified HMath.Random as R

{- | Fitted feature center,principal-component directions,and their explained variances; use
transform functions for new observations.
-}
data PCA = PCA
    { PCA -> [Double]
pcaMean :: [Double]
    -- ^ Training feature means subtracted before projection.
    , PCA -> Dataset
pcaComponents :: Dataset
    -- ^ Principal component loading vectors in component order.
    , PCA -> [Double]
pcaVariances :: [Double]
    -- ^ Explained variances associated with the retained principal components.
    }
    deriving (PCA -> PCA -> Bool
(PCA -> PCA -> Bool) -> (PCA -> PCA -> Bool) -> Eq PCA
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: PCA -> PCA -> Bool
== :: PCA -> PCA -> Bool
$c/= :: PCA -> PCA -> Bool
/= :: PCA -> PCA -> Bool
Eq, Int -> PCA -> ShowS
[PCA] -> ShowS
PCA -> String
(Int -> PCA -> ShowS)
-> (PCA -> String) -> ([PCA] -> ShowS) -> Show PCA
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> PCA -> ShowS
showsPrec :: Int -> PCA -> ShowS
$cshow :: PCA -> String
show :: PCA -> String
$cshowList :: [PCA] -> ShowS
showList :: [PCA] -> ShowS
Show)

-- | Fit count principal axes from at least two observations; uses sample covariance (denominator n-1).
fitPCA :: Int -> Dataset -> Either MathError PCA
fitPCA :: Int -> Dataset -> Either MathError PCA
fitPCA Int
count Dataset
xs = do
    (n, d) <- Dataset -> Either MathError (Int, Int)
dataset Dataset
xs
    require (n >= 2 && count >= 1 && count <= d) "PCA requires n>=2 and 1<=components<=features"
    let mu = Dataset -> [Double]
meanV Dataset
xs
        centered = ([Double] -> [Double]) -> Dataset -> Dataset
forall a b. (a -> b) -> [a] -> [b]
map ([Double] -> [Double] -> [Double]
`sub` [Double]
mu) Dataset
xs
        cov = ([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
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1))) (Dataset -> Dataset -> Dataset
matmul (Dataset -> Dataset
forall a. [[a]] -> [[a]]
transpose Dataset
centered) Dataset
centered)
    pcaFromCov count mu cov

pcaFromCov :: Int -> [Double] -> Dataset -> Either MathError PCA
pcaFromCov :: Int -> [Double] -> Dataset -> Either MathError PCA
pcaFromCov Int
count [Double]
mu Dataset
cov = do
    (vals, vs) <- Dataset -> Either MathError ([Double], Dataset)
symmetricEigen Dataset
cov
    let components = Int -> Dataset -> Dataset
forall a. Int -> [a] -> [a]
take Int
count (Dataset -> Dataset
forall a. [[a]] -> [[a]]
transpose Dataset
vs); variance = Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
take Int
count [Double]
vals
    require (all (>= -1e-8) variance) "PCA covariance is not positive semidefinite"
    Right (PCA mu components (map (max 0) variance))

-- | Center with the training mean and project onto retained axes; no whitening is applied.
transformPCA :: PCA -> Dataset -> Either MathError Dataset
transformPCA :: PCA -> Dataset -> Either MathError Dataset
transformPCA (PCA [Double]
mu Dataset
components [Double]
_) Dataset
xs = do
    Int -> Dataset -> Either MathError ()
queries ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
mu) Dataset
xs
    let out :: Dataset
out = [([Double] -> Double) -> Dataset -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map ([Double] -> [Double] -> Double
`dot` [Double] -> [Double] -> [Double]
sub [Double]
r [Double]
mu) Dataset
components | [Double]
r <- Dataset
xs]
    [Double] -> Either MathError ()
checked (Dataset -> [Double]
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat Dataset
out)
    Dataset -> Either MathError Dataset
forall a b. b -> Either a b
Right Dataset
out

-- | Reconstruct original features from retained scores and restore the training mean.
inversePCA :: PCA -> Dataset -> Either MathError Dataset
inversePCA :: PCA -> Dataset -> Either MathError Dataset
inversePCA (PCA [Double]
mu Dataset
components [Double]
_) Dataset
xs = do
    Int -> Dataset -> Either MathError ()
queries (Dataset -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length Dataset
components) Dataset
xs
    let out :: Dataset
out = ([Double] -> [Double]) -> Dataset -> Dataset
forall a b. (a -> b) -> [a] -> [b]
map (\[Double]
r -> [Double] -> [Double] -> [Double]
add [Double]
mu (([Double] -> Double) -> Dataset -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map ([Double] -> [Double] -> Double
dot [Double]
r) (Dataset -> Dataset
forall a. [[a]] -> [[a]]
transpose Dataset
components))) Dataset
xs
    [Double] -> Either MathError ()
checked (Dataset -> [Double]
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat Dataset
out)
    Dataset -> Either MathError Dataset
forall a b. b -> Either a b
Right Dataset
out

{- | Exact sufficient statistics, including between-batch mean correction.
Sample count, running mean, and centered scatter sufficient statistics.
-}
data IncrementalPCA = IPCA Int [Double] Dataset deriving (IncrementalPCA -> IncrementalPCA -> Bool
(IncrementalPCA -> IncrementalPCA -> Bool)
-> (IncrementalPCA -> IncrementalPCA -> Bool) -> Eq IncrementalPCA
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: IncrementalPCA -> IncrementalPCA -> Bool
== :: IncrementalPCA -> IncrementalPCA -> Bool
$c/= :: IncrementalPCA -> IncrementalPCA -> Bool
/= :: IncrementalPCA -> IncrementalPCA -> Bool
Eq, Int -> IncrementalPCA -> ShowS
[IncrementalPCA] -> ShowS
IncrementalPCA -> String
(Int -> IncrementalPCA -> ShowS)
-> (IncrementalPCA -> String)
-> ([IncrementalPCA] -> ShowS)
-> Show IncrementalPCA
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> IncrementalPCA -> ShowS
showsPrec :: Int -> IncrementalPCA -> ShowS
$cshow :: IncrementalPCA -> String
show :: IncrementalPCA -> String
$cshowList :: [IncrementalPCA] -> ShowS
showList :: [IncrementalPCA] -> ShowS
Show)

-- | Create empty running covariance state for a positive feature dimension.
incrementalPCA :: Int -> Either MathError IncrementalPCA
incrementalPCA :: Int -> Either MathError IncrementalPCA
incrementalPCA Int
d = do
    Bool -> String -> Either MathError ()
require (Int
d Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
0 Bool -> Bool -> Bool
&& Int
d Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
100000) String
"invalid incremental PCA dimension"
    IncrementalPCA -> Either MathError IncrementalPCA
forall a b. b -> Either a b
Right (Int -> [Double] -> Dataset -> IncrementalPCA
IPCA Int
0 (Int -> Double -> [Double]
forall a. Int -> a -> [a]
replicate Int
d Double
0) (Int -> [Double] -> Dataset
forall a. Int -> a -> [a]
replicate Int
d (Int -> Double -> [Double]
forall a. Int -> a -> [a]
replicate Int
d Double
0)))

-- | Merge observations by Welford updates. Empty batches preserve the previous state.
updatePCA :: IncrementalPCA -> Dataset -> Either MathError IncrementalPCA
updatePCA :: IncrementalPCA -> Dataset -> Either MathError IncrementalPCA
updatePCA old :: IncrementalPCA
old@(IPCA Int
_ [Double]
mu Dataset
_) Dataset
xs = do
    Int -> Dataset -> Either MathError ()
queries ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
mu) Dataset
xs
    (IncrementalPCA -> [Double] -> Either MathError IncrementalPCA)
-> IncrementalPCA -> Dataset -> Either MathError IncrementalPCA
forall (t :: * -> *) (m :: * -> *) b a.
(Foldable t, Monad m) =>
(b -> a -> m b) -> b -> t a -> m b
foldM IncrementalPCA -> [Double] -> Either MathError IncrementalPCA
update IncrementalPCA
old Dataset
xs
  where
    update :: IncrementalPCA -> [Double] -> Either MathError IncrementalPCA
update (IPCA Int
n [Double]
m Dataset
s) [Double]
x = do
        Bool -> String -> Either MathError ()
require (Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
forall a. Bounded a => a
maxBound) String
"incremental PCA sample count overflow"
        let delta :: [Double]
delta = [Double] -> [Double] -> [Double]
sub [Double]
x [Double]
m
            m' :: [Double]
m' = [Double] -> [Double] -> [Double]
add [Double]
m (Double -> [Double] -> [Double]
scale (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 Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1)) [Double]
delta)
            delta' :: [Double]
delta' = [Double] -> [Double] -> [Double]
sub [Double]
x [Double]
m'
            s' :: Dataset
s' = [[Dataset
s 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
+ [Double]
delta [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i Double -> Double -> Double
forall a. Num a => a -> a -> a
* [Double]
delta' [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
j | Int
j <- [Int
0 .. [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
m Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]] | Int
i <- [Int
0 .. [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
m Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]]
        [Double] -> Either MathError ()
checked ([Double]
m' [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ Dataset -> [Double]
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat Dataset
s')
        IncrementalPCA -> Either MathError IncrementalPCA
forall a b. b -> Either a b
Right (Int -> [Double] -> Dataset -> IncrementalPCA
IPCA (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) [Double]
m' Dataset
s')

-- | Combine independent batch statistics with the between-mean scatter correction.
mergePCA :: IncrementalPCA -> IncrementalPCA -> Either MathError IncrementalPCA
mergePCA :: IncrementalPCA -> IncrementalPCA -> Either MathError IncrementalPCA
mergePCA a :: IncrementalPCA
a@(IPCA Int
n [Double]
ma Dataset
sa) b :: IncrementalPCA
b@(IPCA Int
m [Double]
mb Dataset
sb) = do
    Bool -> String -> Either MathError ()
require ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
ma Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
mb) String
"incremental PCA feature mismatch"
    Bool -> String -> Either MathError ()
require (Int -> Integer
forall a. Integral a => a -> Integer
toInteger Int
n Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Int -> Integer
forall a. Integral a => a -> Integer
toInteger Int
m Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<= Int -> Integer
forall a. Integral a => a -> Integer
toInteger (Int
forall a. Bounded a => a
maxBound :: Int)) String
"incremental PCA count overflow"
    if Int
n Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0
        then IncrementalPCA -> Either MathError IncrementalPCA
forall a b. b -> Either a b
Right IncrementalPCA
b
        else
            if Int
m Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0
                then IncrementalPCA -> Either MathError IncrementalPCA
forall a b. b -> Either a b
Right IncrementalPCA
a
                else do
                    let total :: Int
total = Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
m
                        delta :: [Double]
delta = [Double] -> [Double] -> [Double]
sub [Double]
mb [Double]
ma
                        mu :: [Double]
mu = [Double] -> [Double] -> [Double]
add [Double]
ma (Double -> [Double] -> [Double]
scale (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
m Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
total) [Double]
delta)
                        s :: Dataset
s =
                            [ [ Dataset
sa 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
sb 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
+ [Double]
delta [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i Double -> Double -> Double
forall a. Num a => a -> a -> a
* [Double]
delta [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
j Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
m Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
total)
                              | Int
j <- [Int
0 .. [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
ma Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]
                              ]
                            | Int
i <- [Int
0 .. [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
ma Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]
                            ]
                    [Double] -> Either MathError ()
checked ([Double]
mu [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ Dataset -> [Double]
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat Dataset
s)
                    IncrementalPCA -> Either MathError IncrementalPCA
forall a b. b -> Either a b
Right (Int -> [Double] -> Dataset -> IncrementalPCA
IPCA Int
total [Double]
mu Dataset
s)

-- | Extract principal axes from accumulated sample covariance; requires at least two observations.
finishPCA :: Int -> IncrementalPCA -> Either MathError PCA
finishPCA :: Int -> IncrementalPCA -> Either MathError PCA
finishPCA Int
count (IPCA Int
n [Double]
mu Dataset
scatter) = do
    Bool -> String -> Either MathError ()
require
        (Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
2 Bool -> Bool -> Bool
&& Int
count Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
0 Bool -> Bool -> Bool
&& Int
count Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
mu)
        String
"insufficient incremental PCA samples/components"
    Int -> [Double] -> Dataset -> Either MathError PCA
pcaFromCov Int
count [Double]
mu (([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
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1))) Dataset
scatter)

-- | Centered training kernel state and normalized positive eigenvectors.
data KernelPCA = KPCA Kernel Dataset [Double] Double Dataset deriving (KernelPCA -> KernelPCA -> Bool
(KernelPCA -> KernelPCA -> Bool)
-> (KernelPCA -> KernelPCA -> Bool) -> Eq KernelPCA
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: KernelPCA -> KernelPCA -> Bool
== :: KernelPCA -> KernelPCA -> Bool
$c/= :: KernelPCA -> KernelPCA -> Bool
/= :: KernelPCA -> KernelPCA -> Bool
Eq, Int -> KernelPCA -> ShowS
[KernelPCA] -> ShowS
KernelPCA -> String
(Int -> KernelPCA -> ShowS)
-> (KernelPCA -> String)
-> ([KernelPCA] -> ShowS)
-> Show KernelPCA
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> KernelPCA -> ShowS
showsPrec :: Int -> KernelPCA -> ShowS
$cshow :: KernelPCA -> String
show :: KernelPCA -> String
$cshowList :: [KernelPCA] -> ShowS
showList :: [KernelPCA] -> ShowS
Show)

-- | Fit kernel principal axes; count must be smaller than n and within positive numerical rank.
fitKernelPCA :: Int -> Kernel -> Dataset -> Either MathError KernelPCA
fitKernelPCA :: Int -> Kernel -> Dataset -> Either MathError KernelPCA
fitKernelPCA Int
count Kernel
kind Dataset
xs = do
    (n, _) <- Dataset -> Either MathError (Int, Int)
dataset Dataset
xs
    require (count > 0 && count < n) "kernel PCA needs fewer components than samples"
    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 means = ([Double] -> Double) -> Dataset -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (\[Double]
r -> [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Double]
r Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n) Dataset
gram
        overall = [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Double]
means Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n
        centered = [[Dataset
gram 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
- [Double]
means [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i Double -> Double -> Double
forall a. Num a => a -> a -> a
- [Double]
means [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
j Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
overall | 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]]
    (vals, vectors) <- symmetricEigen centered
    require (all (> 1e-12) (take count vals)) "kernel PCA requested null component"
    let coefficients' = [[Dataset
vectors 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. Fractional a => a -> a -> a
/ Double -> Double
forall a. Floating a => a -> a
sqrt ([Double]
vals [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
j) | Int
j <- [Int
0 .. Int
count 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]]
    Right (KPCA kind xs means overall coefficients')

-- | Center query kernels with training means before projection, without fitting query statistics.
transformKernelPCA :: KernelPCA -> Dataset -> Either MathError Dataset
transformKernelPCA :: KernelPCA -> Dataset -> Either MathError Dataset
transformKernelPCA (KPCA Kernel
kind Dataset
train [Double]
means Double
overall Dataset
weights) 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 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
        ( \[Double]
x -> do
            k <- ([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 average = [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Double]
k 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]
k)
                centered = (Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\Double
v Double
m -> Double
v Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
average Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
m Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
overall) [Double]
k [Double]
means
            let result = ([Double] -> Double) -> Dataset -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map ([Double] -> [Double] -> Double
dot [Double]
centered) (Dataset -> Dataset
forall a. [[a]] -> [[a]]
transpose Dataset
weights)
            checked result
            Right result
        )
        Dataset
xs

-- | Return reduced (U, singular values, V-transpose) from LAPACK, keeping the largest count values.
truncatedSVD :: Int -> Dataset -> Either MathError (Dataset, [Double], Dataset)
truncatedSVD :: Int -> Dataset -> Either MathError (Dataset, [Double], Dataset)
truncatedSVD Int
count Dataset
xs = do
    (n, d) <- Dataset -> Either MathError (Int, Int)
dataset Dataset
xs
    require (count > 0 && count <= min n d) "invalid truncated SVD count"
    m <- matrix xs
    s <- L.svd m
    Right
        (map (take count) (rows (L.svdU s)), take count (L.singularValues s), take count (rows (L.svdVT s)))

-- | Training center and whitened unmixing directions; component signs/permutations are arbitrary.
data ICA = ICA [Double] Dataset deriving (ICA -> ICA -> Bool
(ICA -> ICA -> Bool) -> (ICA -> ICA -> Bool) -> Eq ICA
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: ICA -> ICA -> Bool
== :: ICA -> ICA -> Bool
$c/= :: ICA -> ICA -> Bool
/= :: ICA -> ICA -> Bool
Eq, Int -> ICA -> ShowS
[ICA] -> ShowS
ICA -> String
(Int -> ICA -> ShowS)
-> (ICA -> String) -> ([ICA] -> ShowS) -> Show ICA
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> ICA -> ShowS
showsPrec :: Int -> ICA -> ShowS
$cshow :: ICA -> String
show :: ICA -> String
$cshowList :: [ICA] -> ShowS
showList :: [ICA] -> ShowS
Show)

{- | Symmetric FastICA with tanh nonlinearity (Hyvarinen 1999), whitening by
sample covariance. Convergence accounts for independent component signs.
-}
fitICA :: MLConfig -> Int -> Dataset -> Either MathError (Fit ICA)
fitICA :: MLConfig -> Int -> Dataset -> Either MathError (Fit ICA)
fitICA MLConfig
c Int
count Dataset
xs = do
    MLConfig -> Either MathError ()
config MLConfig
c
    p <- Int -> Dataset -> Either MathError PCA
fitPCA Int
count Dataset
xs
    require (all (> 1e-12) (pcaVariances p)) "ICA whitening requires full component rank"
    z0 <- transformPCA p xs
    let scales = (Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/) ([Double] -> [Double])
-> ([Double] -> [Double]) -> [Double] -> [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
sqrt ([Double] -> [Double]) -> [Double] -> [Double]
forall a b. (a -> b) -> a -> b
$ PCA -> [Double]
pcaVariances PCA
p
        z = ([Double] -> [Double]) -> Dataset -> Dataset
forall a b. (a -> b) -> [a] -> [b]
map ((Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Double -> Double -> Double
forall a. Num a => a -> a -> a
(*) [Double]
scales) Dataset
z0
        n = Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Dataset -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length Dataset
xs)
        decorrelate Dataset
w = do
            (vals, vs) <- Dataset -> Either MathError ([Double], Dataset)
symmetricEigen (Dataset -> Dataset -> Dataset
matmul Dataset
w (Dataset -> Dataset
forall a. [[a]] -> [[a]]
transpose Dataset
w))
            require (all (> 1e-14) vals) "FastICA unmixing lost rank"
            let invroot = Dataset -> Dataset -> Dataset
matmul (([Double] -> [Double]) -> Dataset -> Dataset
forall a b. (a -> b) -> [a] -> [b]
map (\[Double]
r -> (Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\Double
a Double
v -> Double
a Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double -> Double
forall a. Floating a => a -> a
sqrt Double
v) [Double]
r [Double]
vals) Dataset
vs) (Dataset -> Dataset
forall a. [[a]] -> [[a]]
transpose Dataset
vs)
            Right (matmul invroot w)
        step Dataset
w = do
            let one :: [Double] -> [Double]
one [Double]
row =
                    let projected :: [Double]
projected = ([Double] -> Double) -> Dataset -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map ([Double] -> [Double] -> Double
`dot` [Double]
row) Dataset
z
                        g :: [Double]
g = (Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map Double -> Double
forall a. Floating a => a -> a
tanh [Double]
projected
                        gp :: Double
gp = [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
v Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
v | Double
v <- [Double]
g] Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
n
                     in [Double] -> [Double] -> [Double]
sub (Double -> [Double] -> [Double]
scale (Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
n) (([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]
g Dataset
z))) (Double -> [Double] -> [Double]
scale Double
gp [Double]
row)
            Dataset -> Either MathError Dataset
decorrelate (([Double] -> [Double]) -> Dataset -> Dataset
forall a b. (a -> b) -> [a] -> [b]
map [Double] -> [Double]
one Dataset
w)
    startRaw <- normalMatrix (mlSeed c) count count
    start <- decorrelate startRaw
    fit <- iterateFit c step (\Dataset
a Dataset
b -> [Double] -> Double
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum [Double -> Double
forall a. Num a => a -> a
abs (Double -> Double
forall a. Num a => a -> a
abs ([Double] -> [Double] -> Double
dot [Double]
x [Double]
y) Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1) | ([Double]
x, [Double]
y) <- Dataset -> Dataset -> [([Double], [Double])]
forall a b. [a] -> [b] -> [(a, b)]
zip Dataset
a Dataset
b]) start
    let whitening = (Double -> [Double] -> [Double]) -> [Double] -> Dataset -> Dataset
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Double -> [Double] -> [Double]
scale [Double]
scales (PCA -> Dataset
pcaComponents PCA
p)
    Right fit{fitted = ICA (pcaMean p) (matmul (fitted fit) whitening)}

-- | Apply the fitted centering and unmixing directions to new observations.
transformICA :: ICA -> Dataset -> Either MathError Dataset
transformICA :: ICA -> Dataset -> Either MathError Dataset
transformICA (ICA [Double]
mu Dataset
unmixing) Dataset
xs = do
    Int -> Dataset -> Either MathError ()
queries ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
mu) Dataset
xs
    let answer :: Dataset
answer = [([Double] -> Double) -> Dataset -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map ([Double] -> [Double] -> Double
`dot` [Double] -> [Double] -> [Double]
sub [Double]
x [Double]
mu) Dataset
unmixing | [Double]
x <- Dataset
xs]
    [Double] -> Either MathError ()
checked (Dataset -> [Double]
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat Dataset
answer)
    Dataset -> Either MathError Dataset
forall a b. b -> Either a b
Right Dataset
answer

normalMatrix :: (Integral s) => s -> Int -> Int -> Either MathError Dataset
normalMatrix :: forall s. Integral s => s -> Int -> Int -> Either MathError Dataset
normalMatrix s
seed Int
n Int
d = do
    (values, _) <- 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
n Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
d) (Double -> Double -> RNG -> Either MathError (Double, RNG)
R.normal Double
0 Double
1) (Generator -> Word64 -> RNG
R.seedRNG Generator
R.PCG32 (s -> Word64
forall a b. (Integral a, Num b) => a -> b
fromIntegral s
seed))
    Right [take d (drop (i * d) values) | i <- [0 .. n - 1]]

-- | Fitted feature center,latent factor loadings,and diagonal feature-noise variances.
data FactorAnalysis = FA
    { FactorAnalysis -> [Double]
factorCenter :: [Double]
    , FactorAnalysis -> Dataset
factorLoadings :: Dataset
    -- ^ Fitted feature-by-factor loading matrix.
    , FactorAnalysis -> [Double]
factorNoise :: [Double]
    -- ^ Per-feature residual noise variances in the factor model.
    }
    deriving (FactorAnalysis -> FactorAnalysis -> Bool
(FactorAnalysis -> FactorAnalysis -> Bool)
-> (FactorAnalysis -> FactorAnalysis -> Bool) -> Eq FactorAnalysis
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: FactorAnalysis -> FactorAnalysis -> Bool
== :: FactorAnalysis -> FactorAnalysis -> Bool
$c/= :: FactorAnalysis -> FactorAnalysis -> Bool
/= :: FactorAnalysis -> FactorAnalysis -> Bool
Eq, Int -> FactorAnalysis -> ShowS
[FactorAnalysis] -> ShowS
FactorAnalysis -> String
(Int -> FactorAnalysis -> ShowS)
-> (FactorAnalysis -> String)
-> ([FactorAnalysis] -> ShowS)
-> Show FactorAnalysis
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> FactorAnalysis -> ShowS
showsPrec :: Int -> FactorAnalysis -> ShowS
$cshow :: FactorAnalysis -> String
show :: FactorAnalysis -> String
$cshowList :: [FactorAnalysis] -> ShowS
showList :: [FactorAnalysis] -> ShowS
Show)

{- | Gaussian factor-analysis EM with diagonal noise floor. Unlike PCA, each
observed feature has its own learned residual variance.
-}
fitFactorAnalysis :: MLConfig -> Int -> Double -> Dataset -> Either MathError (Fit FactorAnalysis)
fitFactorAnalysis :: MLConfig
-> Int
-> Double
-> Dataset
-> Either MathError (Fit FactorAnalysis)
fitFactorAnalysis MLConfig
c Int
count Double
floorNoise Dataset
xs = do
    MLConfig -> Either MathError ()
config MLConfig
c
    (n, d) <- Dataset -> Either MathError (Int, Int)
dataset Dataset
xs
    require
        (n >= 2 && count > 0 && count < d && finite floorNoise && floorNoise > 0)
        "invalid factor-analysis rank/noise floor"
    p <- fitPCA count xs
    let mu = PCA -> [Double]
pcaMean PCA
p
        centered = ([Double] -> [Double]) -> Dataset -> Dataset
forall a b. (a -> b) -> [a] -> [b]
map ([Double] -> [Double] -> [Double]
`sub` [Double]
mu) Dataset
xs
        covariance = ([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
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n)) (Dataset -> Dataset -> Dataset
matmul (Dataset -> Dataset
forall a. [[a]] -> [[a]]
transpose Dataset
centered) Dataset
centered)
        noise = Int -> Double -> [Double]
forall a. Int -> a -> [a]
replicate Int
d Double
floorNoise
        loading =
            Dataset -> Dataset
forall a. [[a]] -> [[a]]
transpose ((Double -> [Double] -> [Double]) -> [Double] -> Dataset -> Dataset
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\Double
v [Double]
r -> Double -> [Double] -> [Double]
scale (Double -> Double
forall a. Floating a => a -> a
sqrt (Double -> Double -> Double
forall a. Ord a => a -> a -> a
max Double
floorNoise Double
v)) [Double]
r) (PCA -> [Double]
pcaVariances PCA
p) (PCA -> Dataset
pcaComponents PCA
p))
        step (FA [Double]
_ Dataset
load [Double]
psi) = do
            let weighted :: Dataset
weighted = ([Double] -> Double -> [Double]) -> Dataset -> [Double] -> Dataset
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\[Double]
r Double
v -> Double -> [Double] -> [Double]
scale (Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
v) [Double]
r) Dataset
load [Double]
psi
                information :: Dataset
information = ([Double] -> [Double] -> [Double]) -> Dataset -> Dataset -> Dataset
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith [Double] -> [Double] -> [Double]
add (Int -> Dataset
ident Int
count) (Dataset -> Dataset -> Dataset
matmul (Dataset -> Dataset
forall a. [[a]] -> [[a]]
transpose Dataset
load) Dataset
weighted)
            invInfo <- Dataset -> Either MathError Dataset
inverse Dataset
information
            let beta = Dataset -> Dataset -> Dataset
matmul Dataset
invInfo (Dataset -> Dataset
forall a. [[a]] -> [[a]]
transpose Dataset
weighted)
                sxz = Dataset -> Dataset -> Dataset
matmul Dataset
covariance (Dataset -> Dataset
forall a. [[a]] -> [[a]]
transpose Dataset
beta)
                szz = ([Double] -> [Double] -> [Double]) -> Dataset -> Dataset -> Dataset
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith [Double] -> [Double] -> [Double]
add Dataset
invInfo (Dataset -> Dataset -> Dataset
matmul Dataset
beta Dataset
sxz)
            invZZ <- inverse szz
            let newLoad = Dataset -> Dataset -> Dataset
matmul Dataset
sxz Dataset
invZZ
                newNoise = [Double -> Double -> Double
forall a. Ord a => a -> a -> a
max Double
floorNoise (Dataset
covariance 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] -> [Double] -> Double
dot (Dataset
newLoad Dataset -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! Int
j) (Dataset
sxz Dataset -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! Int
j)) | Int
j <- [Int
0 .. Int
d Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]]
            checked (concat newLoad ++ newNoise)
            Right (FA mu newLoad newNoise)
    iterateFit
        c
        step
        ( \FactorAnalysis
a FactorAnalysis
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 (FactorAnalysis -> Dataset
factorLoadings FactorAnalysis
a) (FactorAnalysis -> Dataset
factorLoadings FactorAnalysis
b)))
                Double -> Double -> Double
forall a. Num a => a -> a -> a
+ [Double] -> Double
norm ([Double] -> [Double] -> [Double]
sub (FactorAnalysis -> [Double]
factorNoise FactorAnalysis
a) (FactorAnalysis -> [Double]
factorNoise FactorAnalysis
b))
        )
        (FA mu loading noise)

-- | Paired training centers, canonical directions, and regularized canonical correlations.
data CCA = CCA [Double] [Double] Dataset Dataset [Double] deriving (CCA -> CCA -> Bool
(CCA -> CCA -> Bool) -> (CCA -> CCA -> Bool) -> Eq CCA
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: CCA -> CCA -> Bool
== :: CCA -> CCA -> Bool
$c/= :: CCA -> CCA -> Bool
/= :: CCA -> CCA -> Bool
Eq, Int -> CCA -> ShowS
[CCA] -> ShowS
CCA -> String
(Int -> CCA -> ShowS)
-> (CCA -> String) -> ([CCA] -> ShowS) -> Show CCA
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> CCA -> ShowS
showsPrec :: Int -> CCA -> ShowS
$cshow :: CCA -> String
show :: CCA -> String
$cshowList :: [CCA] -> ShowS
showList :: [CCA] -> ShowS
Show)

-- | Singular values of the whitened cross-covariance, in descending order.
canonicalCorrelations :: CCA -> [Double]
canonicalCorrelations :: CCA -> [Double]
canonicalCorrelations (CCA [Double]
_ [Double]
_ Dataset
_ Dataset
_ [Double]
correlations) = [Double]
correlations

invsqrt :: Dataset -> Either MathError Dataset
invsqrt :: Dataset -> Either MathError Dataset
invsqrt Dataset
cov = do
    (vals, vs) <- Dataset -> Either MathError ([Double], Dataset)
symmetricEigen Dataset
cov
    require (all (> 1e-14) vals) "positive covariance required"
    Right (matmul (map (\[Double]
r -> (Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\Double
a Double
v -> Double
a Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double -> Double
forall a. Floating a => a -> a
sqrt Double
v) [Double]
r [Double]
vals) vs) (transpose vs))

-- | Fit paired canonical directions with a strictly positive diagonal covariance ridge.
fitCCA :: Int -> Double -> Dataset -> Dataset -> Either MathError CCA
fitCCA :: Int -> Double -> Dataset -> Dataset -> Either MathError CCA
fitCCA Int
count Double
ridge Dataset
xs Dataset
ys = do
    (n, dx) <- Dataset -> Either MathError (Int, Int)
dataset Dataset
xs
    (ny, dy) <- dataset ys
    require
        (n == ny && n >= 2 && count > 0 && count <= min dx dy && finite ridge && ridge > 0)
        "invalid CCA shape/rank/ridge"
    let mx = Dataset -> [Double]
meanV Dataset
xs
        my = Dataset -> [Double]
meanV Dataset
ys
        x = ([Double] -> [Double]) -> Dataset -> Dataset
forall a b. (a -> b) -> [a] -> [b]
map ([Double] -> [Double] -> [Double]
`sub` [Double]
mx) Dataset
xs
        y = ([Double] -> [Double]) -> Dataset -> Dataset
forall a b. (a -> b) -> [a] -> [b]
map ([Double] -> [Double] -> [Double]
`sub` [Double]
my) Dataset
ys
        covariance Dataset
a Dataset
b = ([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
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1))) (Dataset -> Dataset -> Dataset
matmul (Dataset -> Dataset
forall a. [[a]] -> [[a]]
transpose Dataset
a) Dataset
b)
        regularize Dataset
z =
            [[Dataset
z 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 .. Dataset -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length Dataset
z Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]] | Int
i <- [Int
0 .. Dataset -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length Dataset
z Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]]
    wx <- invsqrt (regularize (covariance x x))
    wy <- invsqrt (regularize (covariance y y))
    let whitened = Dataset -> Dataset -> Dataset
matmul (Dataset -> Dataset -> Dataset
matmul Dataset
wx (Dataset -> Dataset -> Dataset
covariance Dataset
x Dataset
y)) Dataset
wy
    (u, s, vt) <- truncatedSVD count whitened
    Right (CCA mx my (matmul wx u) (matmul wy (transpose vt)) s)

-- | Project paired query rows with fitted training centers and canonical directions.
transformCCA :: CCA -> Dataset -> Dataset -> Either MathError (Dataset, Dataset)
transformCCA :: CCA -> Dataset -> Dataset -> Either MathError (Dataset, Dataset)
transformCCA (CCA [Double]
mx [Double]
my Dataset
wx Dataset
wy [Double]
_) Dataset
xs Dataset
ys = do
    Int -> Dataset -> Either MathError ()
queries ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
mx) Dataset
xs
    Int -> Dataset -> Either MathError ()
queries ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
my) Dataset
ys
    Bool -> String -> Either MathError ()
require (Dataset -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length Dataset
xs Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Dataset -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length Dataset
ys) String
"CCA paired sample mismatch"
    let a :: Dataset
a = Dataset -> Dataset -> Dataset
matmul (([Double] -> [Double]) -> Dataset -> Dataset
forall a b. (a -> b) -> [a] -> [b]
map ([Double] -> [Double] -> [Double]
`sub` [Double]
mx) Dataset
xs) Dataset
wx; b :: Dataset
b = Dataset -> Dataset -> Dataset
matmul (([Double] -> [Double]) -> Dataset -> Dataset
forall a b. (a -> b) -> [a] -> [b]
map ([Double] -> [Double] -> [Double]
`sub` [Double]
my) Dataset
ys) Dataset
wy
    [Double] -> Either MathError ()
checked (Dataset -> [Double]
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat Dataset
a [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ Dataset -> [Double]
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat Dataset
b)
    (Dataset, Dataset) -> Either MathError (Dataset, Dataset)
forall a b. b -> Either a b
Right (Dataset
a, Dataset
b)

-- | PLS2 regression state with predictor/response centers and a feature-by-response coefficient matrix.
data PLS = PLS [Double] [Double] Dataset deriving (PLS -> PLS -> Bool
(PLS -> PLS -> Bool) -> (PLS -> PLS -> Bool) -> Eq PLS
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: PLS -> PLS -> Bool
== :: PLS -> PLS -> Bool
$c/= :: PLS -> PLS -> Bool
/= :: PLS -> PLS -> Bool
Eq, Int -> PLS -> ShowS
[PLS] -> ShowS
PLS -> String
(Int -> PLS -> ShowS)
-> (PLS -> String) -> ([PLS] -> ShowS) -> Show PLS
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> PLS -> ShowS
showsPrec :: Int -> PLS -> ShowS
$cshow :: PLS -> String
show :: PLS -> String
$cshowList :: [PLS] -> ShowS
showList :: [PLS] -> ShowS
Show)

{- | PLS2 NIPALS with regression deflation; count cannot exceed the predictor
rank. Inner fixed-point iteration obeys the supplied convergence budget.
-}
fitPLS :: MLConfig -> Int -> Dataset -> Dataset -> Either MathError (Fit PLS)
fitPLS :: MLConfig -> Int -> Dataset -> Dataset -> Either MathError (Fit PLS)
fitPLS MLConfig
c Int
count Dataset
xs Dataset
ys = do
    MLConfig -> Either MathError ()
config MLConfig
c
    (n, dx) <- Dataset -> Either MathError (Int, Int)
dataset Dataset
xs
    (ny, dy) <- dataset ys
    require (n == ny && n >= 2 && count > 0 && count <= min (n - 1) dx) "invalid PLS shape/rank"
    let mx = Dataset -> [Double]
meanV Dataset
xs
        my = Dataset -> [Double]
meanV Dataset
ys
        centerX = ([Double] -> [Double]) -> Dataset -> Dataset
forall a b. (a -> b) -> [a] -> [b]
map ([Double] -> [Double] -> [Double]
`sub` [Double]
mx) Dataset
xs
        centerY = ([Double] -> [Double]) -> Dataset -> Dataset
forall a b. (a -> b) -> [a] -> [b]
map ([Double] -> [Double] -> [Double]
`sub` [Double]
my) Dataset
ys
        component (Dataset
x, Dataset
y, Dataset
ws, Dataset
ps, Dataset
qs, Int
total, Double
res, Bool
converged) p
_ = do
            let start :: [Double]
start = Dataset -> [Double]
maximumByNorm (Dataset -> Dataset
forall a. [[a]] -> [[a]]
transpose Dataset
y)
                iteration :: [Double] -> Either MathError [Double]
iteration [Double]
u = do
                    let w0 :: [Double]
w0 = ([Double] -> Double) -> Dataset -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map ([Double] -> [Double] -> Double
`dot` [Double]
u) (Dataset -> Dataset
forall a. [[a]] -> [[a]]
transpose Dataset
x); wn :: Double
wn = [Double] -> Double
norm [Double]
w0
                    Bool -> String -> Either MathError ()
require (Double
wn Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
1e-14) String
"PLS predictor rank exhausted"
                    let w :: [Double]
w = Double -> [Double] -> [Double]
scale (Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
wn) [Double]
w0; t :: [Double]
t = ([Double] -> Double) -> Dataset -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map ([Double] -> [Double] -> Double
`dot` [Double]
w) Dataset
x; q :: [Double]
q = ([Double] -> Double) -> Dataset -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map ([Double] -> [Double] -> Double
`dot` [Double]
t) (Dataset -> Dataset
forall a. [[a]] -> [[a]]
transpose Dataset
y); u' :: [Double]
u' = ([Double] -> Double) -> Dataset -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map ([Double] -> [Double] -> Double
`dot` [Double]
q) Dataset
y
                    Bool -> String -> Either MathError ()
require ([Double] -> Double
norm [Double]
u' Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
1e-14) String
"PLS response rank exhausted"
                    [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] -> Double
norm [Double]
u') [Double]
u')
            inner <-
                MLConfig
-> ([Double] -> Either MathError [Double])
-> ([Double] -> [Double] -> Double)
-> [Double]
-> Either MathError (Fit [Double])
forall a.
MLConfig
-> (a -> Either MathError a)
-> (a -> a -> Double)
-> a
-> Either MathError (Fit a)
iterateFit
                    MLConfig
c
                    [Double] -> Either MathError [Double]
iteration
                    (\[Double]
a [Double]
b -> Double -> Double -> Double
forall a. Ord a => a -> a -> a
min ([Double] -> Double
norm ([Double] -> [Double] -> [Double]
sub [Double]
a [Double]
b)) ([Double] -> Double
norm ([Double] -> [Double] -> [Double]
add [Double]
a [Double]
b)))
                    (Double -> [Double] -> [Double]
scale (Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double -> Double -> Double
forall a. Ord a => a -> a -> a
max Double
1e-300 ([Double] -> Double
norm [Double]
start)) [Double]
start)
            let u = Fit [Double] -> [Double]
forall a. Fit a -> a
fitted Fit [Double]
inner
                w0 = ([Double] -> Double) -> Dataset -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map ([Double] -> [Double] -> Double
`dot` [Double]
u) (Dataset -> Dataset
forall a. [[a]] -> [[a]]
transpose Dataset
x)
                w = Double -> [Double] -> [Double]
scale (Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ [Double] -> Double
norm [Double]
w0) [Double]
w0
                t = ([Double] -> Double) -> Dataset -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map ([Double] -> [Double] -> Double
`dot` [Double]
w) Dataset
x
                tt = [Double] -> [Double] -> Double
dot [Double]
t [Double]
t
                p = Double -> [Double] -> [Double]
scale (Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
tt) (([Double] -> Double) -> Dataset -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map ([Double] -> [Double] -> Double
`dot` [Double]
t) (Dataset -> Dataset
forall a. [[a]] -> [[a]]
transpose Dataset
x))
                q = Double -> [Double] -> [Double]
scale (Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
tt) (([Double] -> Double) -> Dataset -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map ([Double] -> [Double] -> Double
`dot` [Double]
t) (Dataset -> Dataset
forall a. [[a]] -> [[a]]
transpose Dataset
y))
                x' = ([Double] -> Double -> [Double]) -> Dataset -> [Double] -> Dataset
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\[Double]
r Double
z -> [Double] -> [Double] -> [Double]
sub [Double]
r (Double -> [Double] -> [Double]
scale Double
z [Double]
p)) Dataset
x [Double]
t
                y' = ([Double] -> Double -> [Double]) -> Dataset -> [Double] -> Dataset
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\[Double]
r Double
z -> [Double] -> [Double] -> [Double]
sub [Double]
r (Double -> [Double] -> [Double]
scale Double
z [Double]
q)) Dataset
y [Double]
t
            checked (concat x' ++ concat y')
            Right
                ( x'
                , y'
                , ws ++ [w]
                , ps ++ [p]
                , qs ++ [q]
                , total + fitIterations inner
                , max res (fitResidual inner)
                , converged && fitConverged inner
                )
    (_, _, ws, ps, qs, it, res, conv) <-
        foldM component (centerX, centerY, [], [], [], 0, 0, True) [1 .. count]
    let w = Dataset -> Dataset
forall a. [[a]] -> [[a]]
transpose Dataset
ws; p = Dataset -> Dataset
forall a. [[a]] -> [[a]]
transpose Dataset
ps; q = Dataset -> Dataset
forall a. [[a]] -> [[a]]
transpose Dataset
qs
    rotationInv <- inverse (matmul (transpose p) w)
    let b = Dataset -> Dataset -> Dataset
matmul (Dataset -> Dataset -> Dataset
matmul Dataset
w Dataset
rotationInv) (Dataset -> Dataset
forall a. [[a]] -> [[a]]
transpose Dataset
q)
    require (length b == dx && all ((== dy) . length) b) "PLS coefficient shape"
    checked (concat b)
    Right (Fit (PLS mx my b) it res conv)
  where
    maximumByNorm :: Dataset -> [Double]
maximumByNorm Dataset
cols = Dataset
cols Dataset -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! [Double] -> Int
forall a. Ord a => [a] -> Int
argmax (([Double] -> Double) -> Dataset -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map [Double] -> Double
norm Dataset
cols)

-- | Predict response rows using fitted PLS coefficients and training centers.
predictPLS :: PLS -> Dataset -> Either MathError Dataset
predictPLS :: PLS -> Dataset -> Either MathError Dataset
predictPLS (PLS [Double]
mx [Double]
my Dataset
b) Dataset
xs = do
    Int -> Dataset -> Either MathError ()
queries ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
mx) Dataset
xs
    let answer :: Dataset
answer = ([Double] -> [Double]) -> Dataset -> Dataset
forall a b. (a -> b) -> [a] -> [b]
map ([Double] -> [Double] -> [Double]
add [Double]
my) (Dataset -> Dataset -> Dataset
matmul (([Double] -> [Double]) -> Dataset -> Dataset
forall a b. (a -> b) -> [a] -> [b]
map ([Double] -> [Double] -> [Double]
`sub` [Double]
mx) Dataset
xs) Dataset
b)
    [Double] -> Either MathError ()
checked (Dataset -> [Double]
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat Dataset
answer)
    Dataset -> Either MathError Dataset
forall a b. b -> Either a b
Right Dataset
answer

-- | Seeded Gaussian projection with variance 1/dimension per coefficient. Returns transformed rows.
randomProjection :: MLConfig -> Int -> Dataset -> Either MathError Dataset
randomProjection :: MLConfig -> Int -> Dataset -> Either MathError Dataset
randomProjection MLConfig
c Int
dimension Dataset
xs = do
    MLConfig -> Either MathError ()
config MLConfig
c
    (_, d) <- Dataset -> Either MathError (Int, Int)
dataset Dataset
xs
    require (dimension > 0 && dimension <= 10000) "invalid random-projection dimension"
    projection <- normalMatrix (mlSeed c) d dimension
    let out = Dataset -> Dataset -> Dataset
matmul Dataset
xs (([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 -> Double
forall a. Floating a => a -> a
sqrt (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
dimension))) Dataset
projection)
    checked (concat out)
    Right out

{- | Classical metric MDS. Negative eigenvalues are omitted; requesting a
dimension beyond positive rank fails instead of creating imaginary axes.
-}
classicalMDS :: Int -> Dataset -> Either MathError Dataset
classicalMDS :: Int -> Dataset -> Either MathError Dataset
classicalMDS Int
dimension Dataset
distances = do
    (n, d) <- Dataset -> Either MathError (Int, Int)
dataset Dataset
distances
    require
        (n == d && dimension > 0 && dimension < n && all (>= 0) (concat distances))
        "invalid MDS distance matrix"
    require
        ( and
            [ abs (distances !! i !! j - distances !! j !! i) < 1e-10 && abs (distances !! i !! i) < 1e-10
            | i <- [0 .. n - 1]
            , j <- [0 .. n - 1]
            ]
        )
        "MDS distances must be symmetric with zero diagonal"
    let squares = ([Double] -> [Double]) -> Dataset -> Dataset
forall a b. (a -> b) -> [a] -> [b]
map ((Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
2 :: Int))) Dataset
distances
        means = ([Double] -> Double) -> Dataset -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (\[Double]
r -> [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Double]
r Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n) Dataset
squares
        grand = [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Double]
means Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n
        gram =
            [ [-Double
0.5 Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Dataset
squares 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
- [Double]
means [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i Double -> Double -> Double
forall a. Num a => a -> a -> a
- [Double]
means [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
j Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
grand) | 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]
            ]
    (vals, vs) <- symmetricEigen gram
    require (all (> 1e-12) (take dimension vals)) "MDS positive rank insufficient"
    Right [[vs !! i !! j * sqrt (vals !! j) | j <- [0 .. dimension - 1]] | i <- [0 .. n - 1]]

-- | Build a symmetric k-neighbor graph, compute geodesics, then metric MDS. Disconnected graphs fail.
isomap :: Int -> Int -> Dataset -> Either MathError Dataset
isomap :: Int -> Int -> Dataset -> Either MathError Dataset
isomap Int
neighbors Int
dimension Dataset
xs = do
    (n, _) <- Dataset -> Either MathError (Int, Int)
dataset Dataset
xs
    require (neighbors > 0 && neighbors < n) "Isomap neighbor count invalid"
    let ds = [[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]
        knn Int
i = Int -> [Int] -> [Int]
forall a. Int -> [a] -> [a]
take Int
neighbors (((Double, Int) -> Int) -> [(Double, Int)] -> [Int]
forall a b. (a -> b) -> [a] -> [b]
map (Double, Int) -> Int
forall a b. (a, b) -> b
snd (((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
ds 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], Int
j Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
i]))
        graph =
            [ [ if Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
j then Double
0 else if 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 -> [Int]
knn Int
i Bool -> Bool -> Bool
|| 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 -> [Int]
knn Int
j then Dataset
ds Dataset -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! Int
i [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
j else Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
0
              | 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]
            ]
        shortest =
            (Dataset -> Int -> Dataset) -> Dataset -> [Int] -> 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
a Int
k -> [[Double -> Double -> Double
forall a. Ord a => a -> a -> a
min (Dataset
a Dataset -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! Int
i [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
j) (Dataset
a 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
+ Dataset
a Dataset -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! Int
k [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
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]])
                Dataset
graph
                [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]
    require (all finite (concat shortest)) "Isomap neighborhood graph disconnected"
    classicalMDS dimension shortest

embeddingStart :: MLConfig -> Int -> Int -> Either MathError Dataset
embeddingStart :: MLConfig -> Int -> Int -> Either MathError Dataset
embeddingStart MLConfig
c Int
n Int
dimension = ([Double] -> [Double]) -> Dataset -> Dataset
forall a b. (a -> b) -> [a] -> [b]
map ((Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (Double
1e-4 Double -> Double -> Double
forall a. Num a => a -> a -> a
*)) (Dataset -> Dataset)
-> Either MathError Dataset -> Either MathError Dataset
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Word64 -> Int -> Int -> Either MathError Dataset
forall s. Integral s => s -> Int -> Int -> Either MathError Dataset
normalMatrix (MLConfig -> Word64
mlSeed MLConfig
c) Int
n Int
dimension
center :: Dataset -> Dataset
center :: Dataset -> Dataset
center Dataset
x = let mu :: [Double]
mu = Dataset -> [Double]
meanV Dataset
x in ([Double] -> [Double]) -> Dataset -> Dataset
forall a b. (a -> b) -> [a] -> [b]
map ([Double] -> [Double] -> [Double]
`sub` [Double]
mu) Dataset
x

-- | Exact all-pairs t-SNE; perplexity must be strictly between one and n.
tsne :: MLConfig -> Int -> Double -> Dataset -> Either MathError (Fit Dataset)
tsne :: MLConfig
-> Int -> Double -> Dataset -> Either MathError (Fit Dataset)
tsne MLConfig
c Int
dimension Double
perplexity Dataset
xs = do
    MLConfig -> Either MathError ()
config MLConfig
c
    (n, _) <- Dataset -> Either MathError (Int, Int)
dataset Dataset
xs
    require
        ( dimension > 0
            && dimension <= 10
            && perplexity > 1
            && perplexity < fromIntegral n
            && finite perplexity
        )
        "invalid t-SNE dimension/perplexity"
    let squared = [[[Double] -> [Double] -> Double
dist2 [Double]
x [Double]
y | [Double]
y <- Dataset
xs] | [Double]
x <- Dataset
xs]
        conditional Int
i =
            let distances :: [Double]
distances = [Dataset
squared 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
i]
                minimumD :: Double
minimumD = [Double] -> Double
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
minimum [Double]
distances
                search :: t -> Double -> Maybe Double -> Maybe Double -> [Double]
search t
k Double
beta Maybe Double
lo Maybe Double
hi =
                    let ex :: [Double]
ex = (Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (\Double
d -> Double -> Double
forall a. Floating a => a -> a
exp (-Double
beta Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
d Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
minimumD))) [Double]
distances
                        total :: Double
total = [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Double]
ex
                        probs :: [Double]
probs = (Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
total) [Double]
ex
                        entropy :: Double
entropy = 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 [Double
prob Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
log Double
prob | Double
prob <- [Double]
probs, Double
prob Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0])
                        diff :: Double
diff = Double
entropy Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double -> Double
forall a. Floating a => a -> a
log Double
perplexity
                        next :: Double
next =
                            if Double
diff Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0
                                then Double -> (Double -> Double) -> Maybe Double -> Double
forall b a. b -> (a -> b) -> Maybe a -> b
maybe (Double
beta Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
2) (\Double
h -> (Double
beta Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
h) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2) Maybe Double
hi
                                else Double -> (Double -> Double) -> Maybe Double -> Double
forall b a. b -> (a -> b) -> Maybe a -> b
maybe (Double
beta Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2) (\Double
l -> (Double
beta Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
l) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2) Maybe Double
lo
                     in if t
k t -> t -> Bool
forall a. Ord a => a -> a -> Bool
> t
100 Bool -> Bool -> Bool
|| Double -> Double
forall a. Num a => a -> a
abs Double
diff Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
1e-8
                            then [Double]
probs
                            else if Double
diff Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0 then t -> Double -> Maybe Double -> Maybe Double -> [Double]
search (t
k t -> t -> t
forall a. Num a => a -> a -> a
+ t
1) Double
next (Double -> Maybe Double
forall a. a -> Maybe a
Just Double
beta) Maybe Double
hi else t -> Double -> Maybe Double -> Maybe Double -> [Double]
search (t
k t -> t -> t
forall a. Num a => a -> a -> a
+ t
1) Double
next Maybe Double
lo (Double -> Maybe Double
forall a. a -> Maybe a
Just Double
beta)
                probabilities :: [Double]
probabilities = Int -> Double -> Maybe Double -> Maybe Double -> [Double]
forall {t}.
(Ord t, Num t) =>
t -> Double -> Maybe Double -> Maybe Double -> [Double]
search (Int
0 :: Int) Double
1 Maybe Double
forall a. Maybe a
Nothing Maybe Double
forall a. Maybe a
Nothing
             in Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
take Int
i [Double]
probabilities [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ [Double
0] [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
drop Int
i [Double]
probabilities
        conditionals = (Int -> [Double]) -> [Int] -> Dataset
forall a b. (a -> b) -> [a] -> [b]
map Int -> [Double]
conditional [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]
        p =
            [ [(Dataset
conditionals 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
conditionals Dataset -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! Int
j [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (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
n) | 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]
            ]
        step Dataset
ys = do
            let numerator :: Dataset
numerator =
                    [ [if Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
j then Double
0 else Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ [Double] -> [Double] -> Double
dist2 (Dataset
ys Dataset -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! Int
i) (Dataset
ys 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]
                    ]
                total :: Double
total = [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum (Dataset -> [Double]
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat Dataset
numerator)
                grad :: Dataset
grad =
                    [ ([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]
scale
                            (Double
4 Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Dataset
p 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
numerator 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. Fractional a => a -> a -> a
/ Double
total) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Dataset
numerator 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]
sub (Dataset
ys Dataset -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! Int
i) (Dataset
ys 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]
                    ]
                next :: Dataset
next = Dataset -> Dataset
center (([Double] -> [Double] -> [Double]) -> Dataset -> Dataset -> Dataset
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\[Double]
y [Double]
g -> [Double] -> [Double] -> [Double]
sub [Double]
y (Double -> [Double] -> [Double]
scale (MLConfig -> Double
mlRate MLConfig
c) [Double]
g)) Dataset
ys Dataset
grad)
            [Double] -> Either MathError ()
checked (Dataset -> [Double]
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat Dataset
next)
            Dataset -> Either MathError Dataset
forall a b. b -> Either a b
Right Dataset
next
    start <- embeddingStart c n dimension
    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

{- | Exact full fuzzy-set cross entropy, caller supplies fitted curve a,b.
Local connectivity rho is the first positive neighbor distance. This dense
reference resolves every repulsive pair instead of sampling negatives.
-}
umap :: MLConfig -> Int -> Int -> Double -> Double -> Dataset -> Either MathError (Fit Dataset)
umap :: MLConfig
-> Int
-> Int
-> Double
-> Double
-> Dataset
-> Either MathError (Fit Dataset)
umap MLConfig
c Int
neighbors Int
dimension Double
a Double
b Dataset
xs = do
    MLConfig -> Either MathError ()
config MLConfig
c
    (n, _) <- Dataset -> Either MathError (Int, Int)
dataset Dataset
xs
    require
        ( neighbors >= 2
            && neighbors < n
            && dimension > 0
            && dimension <= 10
            && 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) [a, b]
        )
        "invalid UMAP configuration"
    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]
        nearest Int
i = Int -> [(Double, Int)] -> [(Double, Int)]
forall a. Int -> [a] -> [a]
take Int
neighbors (((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], Int
j Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
i])
        memberships Int
i =
            let near :: [(Double, Int)]
near = Int -> [(Double, Int)]
nearest Int
i
                positive :: [Double]
positive = (Double -> Bool) -> [Double] -> [Double]
forall a. (a -> Bool) -> [a] -> [a]
filter (Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0) (((Double, Int) -> Double) -> [(Double, Int)] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (Double, Int) -> Double
forall a b. (a, b) -> a
fst [(Double, Int)]
near)
                rho :: Double
rho = if [Double] -> Bool
forall a. [a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [Double]
positive then Double
0 else [Double] -> Double
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
minimum [Double]
positive
                target :: Double
target = 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
neighbors)
                value :: Double -> Double
value Double
bandwidth = [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Double -> Double
forall a. Floating a => a -> a
exp (-Double -> Double -> Double
forall a. Ord a => a -> a -> a
max Double
0 (Double
d Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
rho) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
bandwidth) | (Double
d, Int
_) <- [(Double, Int)]
near]
                search :: t -> Double -> Double -> Double
search t
k Double
lo Double
hi =
                    let mid :: Double
mid = (Double
lo Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
hi) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2
                     in if t
k t -> t -> Bool
forall a. Ord a => a -> a -> Bool
> t
80 then Double
mid else if Double -> Double
value Double
mid Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
target then t -> Double -> Double -> Double
search (t
k t -> t -> t
forall a. Num a => a -> a -> a
+ t
1) Double
lo Double
mid else t -> Double -> Double -> Double
search (t
k t -> t -> t
forall a. Num a => a -> a -> a
+ t
1) Double
mid Double
hi
                sigma :: Double
sigma = Int -> Double -> Double -> Double
forall {t}. (Ord t, Num t) => t -> Double -> Double -> Double
search (Int
0 :: Int) Double
1e-12 (Double -> Double -> Double
forall a. Ord a => a -> a -> a
max Double
1 ([Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum (((Double, Int) -> Double) -> [(Double, Int)] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (Double, Int) -> Double
forall a b. (a, b) -> a
fst [(Double, Int)]
near) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
100))
             in [ case Int -> [(Int, Double)] -> Maybe Double
forall a b. Eq a => a -> [(a, b)] -> Maybe b
lookup Int
j [(Int
index, Double -> Double
forall a. Floating a => a -> a
exp (-Double -> Double -> Double
forall a. Ord a => a -> a -> a
max Double
0 (Double
d Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
rho) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
sigma)) | (Double
d, Int
index) <- [(Double, Int)]
near] of
                    Just Double
v -> Double
v
                    Maybe Double
Nothing -> Double
0
                | Int
j <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]
                ]
        directed = (Int -> [Double]) -> [Int] -> Dataset
forall a b. (a -> b) -> [a] -> [b]
map Int -> [Double]
memberships [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]
        p =
            [ [let u :: Double
u = Dataset
directed Dataset -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! Int
i [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
j; v :: Double
v = Dataset
directed Dataset -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! Int
j [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i in Double
u Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
v Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
u Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
v | 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]
            ]
        step Dataset
ys = do
            let grad :: Dataset
grad =
                    [ ([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
                        [ if Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
j
                            then Int -> Double -> [Double]
forall a. Int -> a -> [a]
replicate Int
dimension Double
0
                            else
                                let delta :: [Double]
delta = [Double] -> [Double] -> [Double]
sub (Dataset
ys Dataset -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! Int
i) (Dataset
ys Dataset -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! Int
j); r2 :: Double
r2 = Double -> Double -> Double
forall a. Ord a => a -> a -> a
max Double
1e-8 ([Double] -> [Double] -> Double
dot [Double]
delta [Double]
delta); q :: Double
q = Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
r2 Double -> Double -> Double
forall a. Floating a => a -> a -> a
** Double
b)
                                 in (Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (\Double
v -> Double -> Double -> Double
forall a. Ord a => a -> a -> a
max (-Double
4) (Double -> Double -> Double
forall a. Ord a => a -> a -> a
min Double
4 Double
v)) (Double -> [Double] -> [Double]
scale (Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Dataset
p 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
- Double
q) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
r2) [Double]
delta)
                        | 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]
                    ]
                next :: Dataset
next = Dataset -> Dataset
center (([Double] -> [Double] -> [Double]) -> Dataset -> Dataset -> Dataset
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\[Double]
y [Double]
g -> [Double] -> [Double] -> [Double]
sub [Double]
y (Double -> [Double] -> [Double]
scale (MLConfig -> Double
mlRate MLConfig
c) [Double]
g)) Dataset
ys Dataset
grad)
            [Double] -> Either MathError ()
checked (Dataset -> [Double]
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat Dataset
next)
            Dataset -> Either MathError Dataset
forall a b. b -> Either a b
Right Dataset
next
    start <- embeddingStart c n dimension
    iterateFit c step (\Dataset
u Dataset
v -> [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
u Dataset
v))) start