{-# LANGUAGE DataKinds #-}
{-# LANGUAGE GADTs #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE KindSignatures #-}
{-# LANGUAGE LambdaCase #-}
{-# LANGUAGE MagicHash #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE StandaloneDeriving #-}
module Eigen.Solver.LA
( Decomposition(..)
, solve
, relativeError
, rank
, kernel
, image
, linearRegression
) where
import Eigen.Internal (Elem, Cast(..))
import Eigen.Matrix
import Foreign.C.Types (CInt)
import Foreign.Marshal.Alloc (alloca)
import Foreign.Storable (Storable(..))
import GHC.TypeLits (KnownNat)
import GHC.Types
import Prelude
import qualified Eigen.Internal as Internal
import qualified Eigen.Matrix.Mutable as MM
import qualified Eigen.Matrix as M
import qualified Data.List as List
import qualified Data.Vector.Storable as VS
import qualified Foreign.Concurrent as FC
data Decomposition
= PartialPivLU
| FullPivLU
| HouseholderQR
| ColPivHouseholderQR
| FullPivHouseholderQR
| LLT
| LDLT
| JacobiSVD
deriving (Enum, Eq, Show, Read)
con2CTag :: Decomposition -> CInt
{-# INLINE con2CTag #-}
con2CTag = \case
PartialPivLU -> 0
FullPivLU -> 1
HouseholderQR -> 2
ColPivHouseholderQR -> 3
FullPivHouseholderQR -> 4
LLT -> 5
LDLT -> 6
JacobiSVD -> 7
solve :: (KnownNat n, KnownNat m, KnownNat n1, KnownNat m1, Elem a)
=> Decomposition
-> Matrix n m a
-> Matrix n1 m1 a
-> Matrix m 1 a
solve d a b = Internal.performIO $ do
x :: MM.IOMatrix m 1 a <- MM.new
MM.unsafeWith x $ \x_vals x_rows x_cols ->
unsafeWith a $ \a_vals a_rows a_cols ->
unsafeWith b $ \b_vals b_rows b_cols ->
Internal.call $ Internal.solve (con2CTag d)
x_vals x_rows x_cols
a_vals a_rows a_cols
b_vals b_rows b_cols
unsafeFreeze x
relativeError :: (KnownNat n, KnownNat m, KnownNat n1, KnownNat m1, KnownNat n2, KnownNat m2, Elem a)
=> Matrix n m a
-> Matrix n1 m1 a
-> Matrix n2 m2 a
-> a
relativeError x a b = Internal.performIO $
unsafeWith x $ \x_vals x_rows x_cols ->
unsafeWith a $ \a_vals a_rows a_cols ->
unsafeWith b $ \b_vals b_rows b_cols ->
alloca $ \pe -> do
Internal.call $ Internal.relativeError pe
x_vals x_rows x_cols
a_vals a_rows a_cols
b_vals b_rows b_cols
fromC <$> peek pe
rank :: (KnownNat n, KnownNat m, Elem a)
=> Decomposition
-> Matrix n m a
-> Int
rank d m = Internal.performIO $ alloca $ \pr -> do
Internal.call $ unsafeWith m $ Internal.rank (con2CTag d) pr
fromC <$> peek pr
kernel :: forall a n m. (Elem a, KnownNat n, KnownNat m)
=> Decomposition
-> Matrix n m a
-> Matrix n m a
kernel d m = Internal.performIO $
alloca $ \pvals ->
alloca $ \prows ->
alloca $ \pcols ->
unsafeWith m $ \vals1 rows1 cols1 -> do
Internal.call $ Internal.kernel (con2CTag d)
pvals prows pcols
vals1 rows1 cols1
vals <- peek pvals
rs <- fromC <$> peek prows
cs <- fromC <$> peek pcols
fp <- FC.newForeignPtr vals $ Internal.free vals
pure $ Matrix . Vec $ VS.unsafeFromForeignPtr0 fp (rs * cs)
image :: forall a n m. (Elem a, KnownNat n, KnownNat m)
=> Decomposition
-> Matrix n m a
-> Matrix n m a
image d m = Internal.performIO $
alloca $ \pvals ->
alloca $ \prows ->
alloca $ \pcols ->
unsafeWith m $ \vals1 rows1 cols1 -> do
Internal.call $ Internal.image (con2CTag d)
pvals prows pcols
vals1 rows1 cols1
vals <- peek pvals
rs <- fromC <$> peek prows
cs <- fromC <$> peek pcols
fp <- FC.newForeignPtr vals $ Internal.free vals
pure $ Matrix . Vec $ VS.unsafeFromForeignPtr0 fp (rs * cs)
linearRegression :: forall r. (KnownNat r)
=> Internal.Row r
-> [[Double]]
-> Maybe ([Double], Double)
linearRegression _ points = do
_a :: MatrixXd r 2 <- M.fromList $ List.map ((1:) . tail) points
_b :: MatrixXd r 1 <- M.fromList $ List.map ((:[]) . head) points
let _x = solve ColPivHouseholderQR _a _b
let e = relativeError _x _a _b
let coeffs = List.map head $ M.toList _x
return (coeffs, e)