module Data.Eigen.Matrix (
Matrix(..),
fromList,
toList,
empty,
zero,
ones,
identity,
constant,
cols,
rows,
coeff,
minCoeff,
maxCoeff,
col,
row,
block,
topRows,
bottomRows,
leftCols,
rightCols,
norm,
squaredNorm,
determinant,
inverse,
adjoint,
transpose,
normalize,
freeze,
thaw,
modify,
with
) where
import Data.List (intercalate)
import Text.Printf
import Foreign.Ptr
import Foreign.Storable
import Foreign.ForeignPtr
import Foreign.C.Types
import Foreign.Marshal.Array
import Control.Monad
import Control.Monad.ST
import System.IO.Unsafe
import qualified Data.Vector.Storable as VS
import qualified Data.Vector.Storable.Mutable as VSM
import qualified Data.Eigen.Matrix.Mutable as MM
import Data.Eigen.Internal
data Matrix = Matrix {
m_rows :: Int,
m_cols :: Int,
m_vals :: VS.Vector CDouble
};
instance Show Matrix where
show m@Matrix{..} = printf "Matrix %dx%d\n%s\n" m_rows m_cols $
intercalate "\n" $ map (intercalate "\t" . map show) $ toList m
instance Num Matrix where
(*) = binop MM.mul
(+) = binop MM.add
() = binop MM.sub
fromInteger = undefined
signum = undefined
abs = undefined
empty :: Matrix
empty = Matrix 0 0 VS.empty
constant :: Int -> Int -> Double -> Matrix
constant rows cols val = Matrix rows cols $ VS.replicate (rows * cols) (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 $ runST $ do
vm <- VSM.replicate (size * size) 0
forM_ [0..pred size] $ \n ->
VSM.write vm (n * size + n) 1
VS.unsafeFreeze vm
rows :: Matrix -> Int
rows = m_rows
cols :: Matrix -> Int
cols = m_cols
coeff :: Int -> Int -> Matrix -> Double
coeff row col Matrix{..} = cast $ m_vals VS.! (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 = fromList $
[[coeff row col m | col <- take blockCols [startCol..]] | row <- take blockRows [startRow..]]
maxCoeff :: Matrix -> Double
maxCoeff Matrix{..} = cast $ VS.maximum m_vals
minCoeff :: Matrix -> Double
minCoeff Matrix{..} = 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
binop :: (MM.MMatrix -> MM.MMatrix -> MM.MMatrix -> IO a) -> Matrix -> Matrix -> Matrix
binop f lhs rhs = unsafePerformIO $ do
ret <- MM.new 0 0
lhs <- thaw lhs
rhs <- thaw rhs
_ <- f ret lhs rhs
freeze ret
fromList :: [[Double]] -> Matrix
fromList list = Matrix rows cols vals where
rows = length list
cols = maximum $ map length list
vals = runST $ do
vm <- VSM.replicate (rows * cols) 0
zipWithM_ (\row vals ->
zipWithM_ (\col val ->
VSM.write vm (col * rows + row) (cast val)) [0..] vals) [0..] list
VS.unsafeFreeze vm
toList :: Matrix -> [[Double]]
toList Matrix{..} = [[cast $ m_vals VS.! (col * m_rows + row) | col <- [0..pred m_cols]] | row <- [0..pred m_rows]]
norm :: Matrix -> Double
norm m = unsafePerformIO $ thaw m >>= MM.norm
squaredNorm :: Matrix -> Double
squaredNorm m = unsafePerformIO $ thaw m >>= MM.squaredNorm
determinant :: Matrix -> Double
determinant m = unsafePerformIO $ thaw m >>= MM.determinant
inverse :: Matrix -> Matrix
inverse = modify MM.inverse
adjoint :: Matrix -> Matrix
adjoint = modify MM.adjoint
transpose :: Matrix -> Matrix
transpose = modify MM.transpose
normalize :: Matrix -> Matrix
normalize = modify MM.normalize
freeze :: MM.MMatrix -> IO Matrix
freeze mm = MM.with mm $ \pm -> do
rows <- fmap cast $ c_rows pm
cols <- fmap cast $ c_cols pm
let len = rows * cols
src <- c_data pm
fp <- mallocForeignPtrArray len
withForeignPtr fp $ \dst -> copyArray dst src len
return Matrix {
m_rows = rows,
m_cols = cols,
m_vals = VS.unsafeFromForeignPtr fp 0 len
}
thaw :: Matrix -> IO MM.MMatrix
thaw Matrix{..} = withForeignPtr fp $ \src -> do
let len = m_rows * m_cols
pm <- c_create (cast m_rows) (cast m_cols)
dst <- c_data pm
copyArray dst (plusPtr src $ off * sizeOf (undefined :: CDouble)) len
fmap MM.MMatrix $ newForeignPtr c_destroy pm
where (fp, off, _) = VS.unsafeToForeignPtr m_vals
modify :: (MM.MMatrix -> IO ()) -> Matrix -> Matrix
modify f m = unsafePerformIO $ thaw m >>= \mm -> f mm >> freeze mm
with :: Matrix -> (Ptr C_MatrixXd -> IO a) -> IO a
with m f = thaw m >>= \mm -> MM.with mm f