module Data.Eigen.Matrix (
Matrix(..),
valid,
fromList,
toList,
generate,
empty,
null,
square,
zero,
ones,
identity,
constant,
random,
cols,
rows,
(!),
coeff,
unsafeCoeff,
col,
row,
block,
topRows,
bottomRows,
leftCols,
rightCols,
sum,
prod,
mean,
minCoeff,
maxCoeff,
trace,
norm,
squaredNorm,
blueNorm,
hypotNorm,
determinant,
all,
any,
count,
add,
sub,
mul,
diagonal,
transpose,
inverse,
adjoint,
conjugate,
normalize,
modify,
thaw,
freeze,
unsafeThaw,
unsafeFreeze,
unsafeWith,
) where
import Prelude hiding (null, sum, all, any)
import Data.List (intercalate)
import Data.Tuple
import Foreign.Ptr
import Foreign.C.Types
import Foreign.C.String
import Text.Printf
import Control.Monad
import Control.Monad.ST
import Control.Monad.Primitive
import Control.Applicative hiding (empty)
import qualified Data.Vector.Storable as VS
import qualified Data.Vector.Storable.Mutable as VSM
import qualified Data.Eigen.Internal as I
import qualified Data.Eigen.Matrix.Mutable as M
data Matrix = Matrix {
m_rows :: Int,
m_cols :: Int,
m_vals :: VS.Vector CDouble
};
instance Show Matrix where
show m@Matrix{..} = concat [
"Matrix ", show m_rows, "x", show m_cols,
"\n", intercalate "\n" $ map (intercalate "\t" . map show) $ toList m, "\n"]
instance Num Matrix where
(*) = mul
(+) = add
() = sub
fromInteger = constant 1 1 . fromInteger
signum m@Matrix{..} = m { m_vals = VS.map signum m_vals }
abs m@Matrix{..} = m { m_vals = VS.map abs m_vals }
empty :: Matrix
empty = Matrix 0 0 VS.empty
null :: Matrix -> Bool
null Matrix{..} = m_rows == 0 && m_cols == 0
square :: Matrix -> Bool
square Matrix{..} = m_rows == m_cols
constant :: Int -> Int -> Double -> Matrix
constant rows cols val = Matrix rows cols $ VS.replicate (rows * cols) (I.cast val)
zero :: Int -> Int -> Matrix
zero rows cols = constant rows cols 0
ones :: Int -> Int -> Matrix
ones rows cols = constant rows cols 1
identity :: Int -> Matrix
identity size = Matrix size size $ VS.create $ do
vm <- VSM.replicate (size * size) 0
forM_ [0..pred size] $ \n ->
VSM.write vm (n * size + n) 1
return vm
random :: Int -> Int -> IO Matrix
random rows cols = do
m <- M.new rows cols
I.call $ M.unsafeWith m I.c_random
unsafeFreeze m
rows :: Matrix -> Int
rows = m_rows
cols :: Matrix -> Int
cols = m_cols
dims :: Matrix -> (Int, Int)
dims Matrix{..} = (m_rows, m_cols)
(!) :: Matrix -> (Int,Int) -> Double
(!) m (row,col) = coeff row col m
coeff :: Int -> Int -> Matrix -> Double
coeff row col m@Matrix{..}
| not (valid m) = error "matrix is not valid"
| row < 0 || row >= m_rows = error $ printf "Matrix.coeff: row %d is out of bounds [0..%d)" row m_rows
| col < 0 || col >= m_cols = error $ printf "Matrix.coeff: col %d is out of bounds [0..%d)" col m_cols
| otherwise = unsafeCoeff row col m
unsafeCoeff :: Int -> Int -> Matrix -> Double
unsafeCoeff row col Matrix{..} = I.cast $ VS.unsafeIndex m_vals $ col * m_rows + row
col :: Int -> Matrix -> [Double]
col c m@Matrix{..} = [coeff r c m | r <- [0..pred m_rows]]
row :: Int -> Matrix -> [Double]
row r m@Matrix{..} = [coeff r c m | c <- [0..pred m_cols]]
block :: Int -> Int -> Int -> Int -> Matrix -> Matrix
block startRow startCol blockRows blockCols m =
generate blockRows blockCols $ \row col ->
coeff (startRow + row) (startCol + col) m
valid :: Matrix -> Bool
valid Matrix{..} = m_rows >= 0 && m_cols >= 0 && VS.length m_vals == m_rows * m_cols
maxCoeff :: Matrix -> Double
maxCoeff Matrix{..} = I.cast $ VS.maximum m_vals
minCoeff :: Matrix -> Double
minCoeff Matrix{..} = I.cast $ VS.minimum m_vals
topRows :: Int -> Matrix -> Matrix
topRows rows m@Matrix{..} = block 0 0 rows m_cols m
bottomRows :: Int -> Matrix -> Matrix
bottomRows rows m@Matrix{..} = block (m_rows rows) 0 rows m_cols m
leftCols :: Int -> Matrix -> Matrix
leftCols cols m@Matrix{..} = block 0 0 m_rows cols m
rightCols :: Int -> Matrix -> Matrix
rightCols cols m@Matrix{..} = block 0 (m_cols cols) m_rows cols m
fromList :: [[Double]] -> Matrix
fromList list = Matrix rows cols vals where
rows = length list
cols = maximum $ map length list
vals = VS.create $ do
vm <- VSM.replicate (rows * cols) 0
forM_ (zip [0..] list) $ \(row, vals) ->
forM_ (zip [0..] vals) $ \(col, val) ->
VSM.write vm (col * rows + row) (I.cast val)
return vm
toList :: Matrix -> [[Double]]
toList Matrix{..} = [[I.cast $ m_vals VS.! (col * m_rows + row) | col <- [0..pred m_cols]] | row <- [0..pred m_rows]]
generate :: Int -> Int -> (Int -> Int -> Double) -> Matrix
generate rows cols f = Matrix rows cols $ VS.create $ do
vals <- VSM.new (rows * cols)
forM_ [0..pred rows] $ \row ->
forM_ [0..pred cols] $ \col ->
VSM.write vals (col * rows + row) (I.cast $ f row col)
return vals
sum :: Matrix -> Double
sum = _prop I.c_sum
prod :: Matrix -> Double
prod = _prop I.c_prod
mean :: Matrix -> Double
mean = _prop I.c_prod
trace :: Matrix -> Double
trace = _prop I.c_trace
all :: (Double -> Bool) -> Matrix -> Bool
all f = VS.all (f . I.cast) . m_vals
any :: (Double -> Bool) -> Matrix -> Bool
any f = VS.any (f . I.cast) . m_vals
count :: (Double -> Bool) -> Matrix -> Int
count f = VS.foldl' (\n x -> if f (I.cast x) then succ n else n) 0 . m_vals
norm :: Matrix -> Double
norm = _prop I.c_norm
squaredNorm :: Matrix -> Double
squaredNorm = _prop I.c_squaredNorm
blueNorm :: Matrix -> Double
blueNorm = _prop I.c_blueNorm
hypotNorm :: Matrix -> Double
hypotNorm = _prop I.c_hypotNorm
determinant :: Matrix -> Double
determinant m
| square m = _prop I.c_determinant m
| otherwise = error "Matrix.determinant: non-square matrix"
add :: Matrix -> Matrix -> Matrix
add m1 m2
| dims m1 == dims m2 = _binop const I.c_add m1 m2
| otherwise = error "Matrix.add: matrix should have the same size"
sub :: Matrix -> Matrix -> Matrix
sub m1 m2
| dims m1 == dims m2 = _binop const I.c_sub m1 m2
| otherwise = error "Matrix.add: matrix should have the same size"
mul :: Matrix -> Matrix -> Matrix
mul m1 m2
| m_cols m1 == m_rows m2 = _binop (\(rows, _) (_, cols) -> (rows, cols)) I.c_mul m1 m2
| otherwise = error "Matrix.mul: number of columns for lhs matrix should be the same as number of rows for rhs matrix"
diagonal :: Matrix -> Matrix
diagonal = _unop swap I.c_diagonal
inverse :: Matrix -> Matrix
inverse m@Matrix{..}
| m_rows == m_cols = _unop id I.c_inverse m
| otherwise = error "Matrix.inverse: non-square matrix"
adjoint :: Matrix -> Matrix
adjoint = _unop swap I.c_adjoint
transpose :: Matrix -> Matrix
transpose = _unop swap I.c_transpose
conjugate :: Matrix -> Matrix
conjugate = _unop id I.c_conjugate
normalize :: Matrix -> Matrix
normalize Matrix{..} = I.performIO $ do
vals <- VS.thaw m_vals
VSM.unsafeWith vals $ \p ->
I.call $ I.c_normalize p (I.cast m_rows) (I.cast m_cols)
Matrix m_rows m_cols <$> VS.unsafeFreeze vals
modify :: (forall s. M.MMatrix s -> ST s ()) -> Matrix -> Matrix
modify f m@Matrix{..} = m { m_vals = VS.modify (f . M.MMatrix m_rows m_cols) m_vals } where
freeze :: PrimMonad m => M.MMatrix (PrimState m) -> m Matrix
freeze M.MMatrix{..} = VS.freeze mm_vals >>= return . Matrix mm_rows mm_cols
thaw :: PrimMonad m => Matrix -> m (M.MMatrix (PrimState m))
thaw Matrix{..} = VS.thaw m_vals >>= return . M.MMatrix m_rows m_cols
unsafeFreeze :: PrimMonad m => M.MMatrix (PrimState m) -> m Matrix
unsafeFreeze M.MMatrix{..} = VS.unsafeFreeze mm_vals >>= return . Matrix mm_rows mm_cols
unsafeThaw :: PrimMonad m => Matrix -> m (M.MMatrix (PrimState m))
unsafeThaw Matrix{..} = VS.unsafeThaw m_vals >>= return . M.MMatrix m_rows m_cols
unsafeWith :: Matrix -> (Ptr CDouble -> CInt -> CInt -> IO a) -> IO a
unsafeWith m@Matrix{..} f
| not (valid m) = fail "matrix layout is invalid"
| otherwise = VS.unsafeWith m_vals $ \p -> f p (I.cast m_rows) (I.cast m_cols)
_prop :: (Ptr CDouble -> CInt -> CInt -> IO CDouble) -> Matrix -> Double
_prop f m = I.cast $ I.performIO $ unsafeWith m f
_binop :: ((Int, Int) -> (Int, Int) -> (Int, Int)) -> (Ptr CDouble -> CInt -> CInt -> Ptr CDouble -> CInt -> CInt -> Ptr CDouble -> CInt -> CInt -> IO CString) -> Matrix -> Matrix -> Matrix
_binop f g m1 m2 = I.performIO $ do
m0 <- uncurry M.new $ f (dims m1) (dims m2)
M.unsafeWith m0 $ \vals0 rows0 cols0 ->
unsafeWith m1 $ \vals1 rows1 cols1 ->
unsafeWith m2 $ \vals2 rows2 cols2 ->
I.call $ g
vals0 rows0 cols0
vals1 rows1 cols1
vals2 rows2 cols2
unsafeFreeze m0
_unop :: ((Int,Int) -> (Int,Int)) -> (Ptr CDouble -> CInt -> CInt -> Ptr CDouble -> CInt -> CInt -> IO CString) -> Matrix -> Matrix
_unop f g m1 = I.performIO $ do
m0 <- uncurry M.new $ f (dims m1)
M.unsafeWith m0 $ \vals0 rows0 cols0 ->
unsafeWith m1 $ \vals1 rows1 cols1 ->
I.call $ g
vals0 rows0 cols0
vals1 rows1 cols1
unsafeFreeze m0