-- |
-- Module      :  Statistics.Covariance.GraphicalLasso
-- Description :  Graphical lasso
-- Copyright   :  2021 Dominik Schrempf
-- License     :  GPL-3.0-or-later
--
-- Maintainer  :  dominik.schrempf@gmail.com
-- Stability   :  experimental
-- Portability :  portable
--
-- Creation date: Wed Sep 15 09:23:19 2021.
module Statistics.Covariance.GraphicalLasso
  ( graphicalLasso,
  )
where

import Algorithms.GLasso
import Data.Bifunctor
import qualified Numeric.LinearAlgebra as L

-- | Gaussian graphical model based estimator.
--
-- This function estimates both, the covariance and the precision matrices. It
-- is best suited for sparse covariance matrices.
--
-- For now, this is just a wrapper around 'glasso'.
--
-- See Friedman, J., Hastie, T., & Tibshirani, R., Sparse inverse covariance
-- estimation with the graphical lasso, Biostatistics, 9(3), 432–441 (2007).
-- http://dx.doi.org/10.1093/biostatistics/kxm045.
--
-- Return 'Left' if
--
-- - the regularization parameter is out of bounds \([0, \infty)\).
--
-- - only one sample is available.
--
-- - no parameters are available.
--
-- NOTE: This function may call 'error' due to partial library functions.
graphicalLasso ::
  -- | Regularization or lasso parameter; penalty for non-zero covariances. The
  -- higher the lasso parameter, the sparser the estimated inverse covariance
  -- matrix. Must be non-negative.
  Double ->
  -- | Sample data matrix of dimension \(n \times p\), where \(n\) is the number
  -- of samples (rows), and \(p\) is the number of parameters (columns).
  L.Matrix Double ->
  -- | @Either ErrorString (Covariance matrix, Precision matrix)@.
  Either String (L.Herm Double, L.Herm Double)
graphicalLasso :: Double -> Matrix Double -> Either String (Herm Double, Herm Double)
graphicalLasso Double
l Matrix Double
xs
  | Double
l Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
0 = String -> Either String (Herm Double, Herm Double)
forall a b. a -> Either a b
Left String
"graphicalLasso: Regularization parameter is negative."
  | Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
2 = String -> Either String (Herm Double, Herm Double)
forall a b. a -> Either a b
Left String
"graphicalLasso: Need more than one sample."
  | Int
p Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
1 = String -> Either String (Herm Double, Herm Double)
forall a b. a -> Either a b
Left String
"graphicalLasso: Need at least one parameter."
  | Bool
otherwise =
      (Herm Double, Herm Double)
-> Either String (Herm Double, Herm Double)
forall a b. b -> Either a b
Right ((Herm Double, Herm Double)
 -> Either String (Herm Double, Herm Double))
-> (Herm Double, Herm Double)
-> Either String (Herm Double, Herm Double)
forall a b. (a -> b) -> a -> b
$
        (Vector Double -> Herm Double)
-> (Vector Double -> Herm Double)
-> (Vector Double, Vector Double)
-> (Herm Double, Herm Double)
forall (p :: * -> * -> *) a b c d.
Bifunctor p =>
(a -> b) -> (c -> d) -> p a c -> p b d
bimap Vector Double -> Herm Double
convert Vector Double -> Herm Double
convert ((Vector Double, Vector Double) -> (Herm Double, Herm Double))
-> (Vector Double, Vector Double) -> (Herm Double, Herm Double)
forall a b. (a -> b) -> a -> b
$ Int -> Vector Double -> Double -> (Vector Double, Vector Double)
glasso Int
p (Matrix Double -> Vector Double
forall t. Element t => Matrix t -> Vector t
L.flatten (Matrix Double -> Vector Double) -> Matrix Double -> Vector Double
forall a b. (a -> b) -> a -> b
$ Herm Double -> Matrix Double
forall t. Herm t -> Matrix t
L.unSym Herm Double
sigma) Double
l
  where
    n :: Int
n = Matrix Double -> Int
forall t. Matrix t -> Int
L.rows Matrix Double
xs
    p :: Int
p = Matrix Double -> Int
forall t. Matrix t -> Int
L.cols Matrix Double
xs
    (Vector Double
_, Herm Double
sigma) = Matrix Double -> (Vector Double, Herm Double)
L.meanCov Matrix Double
xs
    convert :: Vector Double -> Herm Double
convert = Matrix Double -> Herm Double
forall t. Matrix t -> Herm t
L.trustSym (Matrix Double -> Herm Double)
-> (Vector Double -> Matrix Double) -> Vector Double -> Herm Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> Vector Double -> Matrix Double
forall t. Storable t => Int -> Vector t -> Matrix t
L.reshape Int
p